Parameterisation of the GNSS troposphere tomography domain with optimisation of the nodes’ distribution

Water vapour is a highly variable constituent of the troposphere; thus, its high-resolution measurements are of great importance to weather prediction systems. The Global Navigation Satellite Systems (GNSS) are operationally used in the estimation of the tropospheric state and assimilation of the results into the weather models. One of the GNSS techniques of troposphere sensing is tomography which provides 3-D fields of wet refractivity. The tomographic results have been successfully assimilated into the numerical weather models, showing the great potential of this technique. The GNSS tomography can be based on two different approaches to the parameterisation of the model’s domain, i.e. block (voxel-based) or grid (node-based) approach. Regardless of the parameterisation approach, the tomographic domain should be discretised, which is usually performed in a regular manner, with a grid resolution depending on the mean distance between the GNSS receivers. In this work, we propose a new parameterisation approach based on the optimisation of the tomographic nodes’ location, taking into account the non-uniform distribution of the GNSS information in the troposphere. The experiment was performed using a dense network of 16 low-cost multi-GNSS receivers located in Wrocław and its suburbs, with a mean distance of 3 km. Cross-validation of four different parameterisation approaches is presented. The validation is performed based on the Weather Research and Forecasting model as well as radiosonde observations. The new approach improves the results of wet refractivity estimation by 0.5–2 ppm in terms of RMSE, especially for altitudes of 0.5–2.0 km.


Introduction
Water vapour strongly contributes to the basic physical and chemical processes in the atmosphere; thus, its spatial and temporal variability has a crucial impact on weather patterns (Turco 1992;Jacob 2001). In the troposphere, water vapour is a highly variable constituent; its concentration can vary by several orders of magnitude between locations in a distance less than 10 km (Couvreux et al. 2005 (Stensrud 2009). Since the conventional measurements are based on radiosonde stations, their spatial and temporal resolution is low (12 h, 250 km in Europe; Guerova et al. 2016), compared to the dynamics of the humidity changes. Recently, another observing system has been established, based on Global Navigation Satellite System (GNSS) observations (Guerova et al. 2016).
The main purpose of the GNSS is to provide information for navigation and positioning (Hofmann-Wellenhof et al. 2007). However, since the GNSS signals propagate through the atmosphere, they are affected by atmospheric components and thus can be utilised in atmospheric monitoring (Tralli and Lichten 1990;Böhm and Schuh 2013). The GNSS-based troposphere parameters are provided in all weather conditions, with a temporal resolution of 1 h or less (Tondaś et al. 2020). This makes the GNSS products a valuable source of information on the tropospheric state. The typical density of the GNSS network is on the order of tens of kilometres, which far exceeds the spatial resolution of the conventional radiosonde observations (Pacione et al. 2017). Recently, the low-cost GNSS receivers, based on mass-market chipsets capable of logging carrier phase measurements, were successfully applied in the troposphere profiling (Barindelli et al. 2018;Bramanto et al. 2018;Krietemeyer et al. 2018). Since this approach is much less expensive than the geodetic-grade receivers, it provides an opportunity for troposphere monitoring on the kilometre scale by further densification of the GNSS networks (Marut et al. 2022).
In the GNSS data processing, a tropospheric delay of the signal is estimated in the zenith direction as a zenith total delay (ZTD, Dach et al. 2015). Since the ZTD is a function of temperature, pressure, and humidity, it provides relevant information for the weather systems (Mendes 1999;Bevis et al. 1992;Wilgan et al. 2015). ZTD can be converted into integrated water vapour (IWV) which contains information about water vapour in the column above the GNSS site (Kleijer 2004;Benevides et al. 2019). Currently, GNSS meteorology is a well-established field in both research and operation (Bennitt and Jupp 2012;Dousa and Bennitt 2013;Li et al. 2015;Hadas and Hobiger 2020). The near-real-time (NRT) tropospheric parameters are operationally provided by European Meteorological Services Network (EUMETNET) under the GNSS Water Vapour Programme (E-GVAP; http:// egvap.dmi.dk) as a result of cooperation between meteorological services and GNSS analysis centres. The studies on the assimilation of the GNSS ZTD and PWV observations into the NWP models show a positive impact on the forecasts, especially during severe weather events, e.g. heavy precipitation or storm (Cucurull et al. 2004;Yan et al. 2009;Boniface et al. 2012;Rohm et al. 2019;Giannaros et al. 2020). It was found by Mahfouf et al. (2015) that although the ZTD observations represent a small fraction of the data assimilated into the NWP model (lower than 2%), their impact on the analysis is comparable to other humidity observing systems. The ZTD observations were assimilated into the numerical models with a horizontal grid resolution of the order of 10 km (Poli et al. 2007;Benjamin et al. 2010;Shoji et al. 2011) as well as a single kilometre scale (De Haan 2013; Arriola et al. 2016). Since the typical resolution of the GNSS network is on the order of tens of kilometres, the ability to detect small-scale weather events is limited (Stensrud 2009;Mateus et al. 2018). However, the studies show a link between the GNSS water vapour data and the development of convective systems that result in thunderstorms and lightning (Mazany et al. 2002).
Apart from the observations integrated in the zenith direction, also slant tropospheric delay (STD) and slant integrated water vapour (SIWV) can be estimated in the GNSS data processing using zenith observations and their gradients (Kačmařík et al. 2017). Since the observations integrated along the signal's path are asymmetrical, they have the potential to better describe the state of the troposphere. The first studies indicate their capability of obtaining information about storms (Kawabata and Shoji 2018). Although the assimilation of the slant observations into the NWP models is a challenging task, the first results show a positive impact on the weather forecasts (Zus et al. 2011(Zus et al. , 2015. GNSS tomography is another approach to the utilisation of the GNSS signal for troposphere monitoring. Using this technique, a 3-dimensional distribution of humidity in the troposphere is reconstructed based on the slant observations of the GNSS signal delay (Flores et al. 2000). Over the last 20 years, research has been conducted on the troposphere tomography in terms of different inversion strategies, constraints imposed on the model, discretisation of the model's domain, and parameterisation (Guerova et al. 2016). Also, an effort has been made to assimilate the tomographic results into the NWP models. Several studies show a positive impact on the models in terms of humidity forecasting, especially during extreme weather events (heavy precipitation; Hanna et al. 2019;Trzcina and Rohm 2019). Usually, a drying effect in the model is observed after assimilation (Möller et al. 2015). It was noticed that the tomographic data have a greater impact on the NWP model than the ZTD observations (Möller et al. 2015;Trzcina et al. 2020). Recently, an observation operator for the 3-D wet refractivity fields was introduced by Trzcina et al. (2020), so the tomography outputs can be directly assimilated into the Weather Research and Forecasting (WRF) model. The tomography-based vertical profiles of humidity were also successfully used for the prediction of storms using a machine learning approach (Łoś et al. 2020). These examples of the positive impact of the GNSS troposphere tomography outputs on weather prediction are encouraging for further development of this technique.
In order to reconstruct the 3-D humidity field in the troposphere, the tomographic model's domain has to be discretised. Usually, the discretisation is performed in a regular manner and the resolution of the horizontal grid depends on the mean distance between the GNSS stations (Bender and Raabe 2007;Troller et al. 2006). Also, an outer domain with a lower resolution can be appended in order to easily include low-elevation rays in the solution (Bosy et al. 2012;Brenot et al. 2019;Trzcina and Rohm 2019). Regarding the vertical layers, it was found by Manning (2013) that irregular distribution outperforms equidistant spacing.
There are two main strategies to parameterise the tomographic domain, i.e. block (voxel-based) parameterisation and grid (node-based) parameterisation. In the block approach, the estimated values (i.e. wet refractivity or water vapour humidity) are assumed to be constant within each small part of the 3-D tomographic domain (Flores et al. 2000;Möller 2017). On the contrary, in the grid parameterization, the values are estimated for particular nodes of the tomographic grid. In this approach, the values are interpolated between the grid and the GNSS rays' paths (Manning 2013;Perler et al. 2011). In the horizontal direction, linear interpolation is usually applied; in the vertical direction, it is recommended to use methods that better reflect humidity changes in the troposphere (e.g. natural splines or scale height of water vapour; Perler 2011; Ding et al. 2018). The two parameterisation approaches use different strategies for the modelling of the GNSS delay along the signal's path. In the block approach, the delay is modelled based on the multiplication of wet refractivity and distance that the ray travels in each voxel (Hirahara 2000;Rohm 2013). On the contrary, in the grid approach, the signal's delay is modelled using numerical integration techniques (Perler 2011;Zhang et al. 2020a).
The node-based parameterisation gives the opportunity to introduce new methods of tomographic nodes' distribution. Recently, it was proposed by Ding et al. (2018) to adjust the model's boundary to the distribution of the GNSS rays and to define tomographic nodes using meshing techniques; the results were encouraging since the RMSE of water vapour density decreased by 30-40% when compared to the standard method. Also, Zhang et al. (2020b) used a similar approach to define model boundaries in order to stabilise the tomographic solution. Although the new methods adjust the model's boundary to the geometry of the GNSS observations, the tomographic grid remains regular. This requirement of a regular grid is needed for the proper modelling of the GNSS delay. However, it is also a limitation; since the information about the tropospheric state comes from the GNSS signals, the quality of the wet refractivity estimation is not homogeneous in space. In this paper, we propose another approach that aims to identify the optimal location of the tomographic nodes, taking into account the spatial distribution of the GNSS rays inside the tomographic model's domain. In this new method, node-based parameterisation is applied, including both regularly and irregularly distributed points. In our study, the cross-validation of the different GNSS tomography parameterisation approaches is performed; the validation is based on the radiosonde measurements and NWP model.
The structure of this paper is as follows: Sect. 2 presents principles of the GNSS troposphere tomography technique. Section 3 describes the data and the tomographic model used in the performed experiment. Section 4 presents the results of the tomographic solutions and their validation with external data sources. The discussion of the results is presented in Sect. 5. The main conclusions of this work can be found in Sect. 6.

Tropospheric wet delay and wet refractivity
The GNSS troposphere tomography provides a 3-D distribution of humidity in the troposphere, based on observations of the GNSS signals' delays along the slant paths between the satellites and the receivers (Flores et al. 2000;Hirahara 2000). As the GNSS signal passes through the troposphere, its propagation is affected by the tropospheric components, which results in the signal's bending and delay. This tropospheric delay is related to the temperature, pressure, and humidity on the signal's path and can be estimated in the GNSS data analysis (Kleijer 2004). Usually, the tropospheric delay is modelled in the zenith direction above the GNSS receiver as a zenith total delay (ZTD; Dach et al. 2015). Davis et al. (1985) proposed a separation of two components of the ZTD, representing, respectively, the hydrostatic and non-hydrostatic effects of the troposphere on the GNSS signal. The two components are often referred to as ZHD (zenith hydrostatic delay) and ZWD (zenith wet delay). The relation between the ZTD and its two components is as follows: where ZTD refers to the total delay values, ZHD is the hydrostatic component, and ZWD is a non-hydrostatic component.
As the GNSS troposphere tomography estimates a distribution of humidity in the troposphere, it takes advantage of the wet delay. The ZWD values can be calculated by subtraction the ZHD from the total delay estimated in the GNSS data processing. For this purpose, Saastamoinen's formula is used (Kleijer 2004): where p is the total air pressure (hPa) at the GNSS station's location, φ denotes the location's latitude, and H is the altitude (m). The values of the wet delay in the zenith direction can be mapped into the directions of the GNSS satellites to get the values of slant wet delays (SWD) for each pair of the GNSS receiver and transmitter. The following mapping formula is applied: where e is an elevation angle of the tracked satellite, A is an azimuth, m w stands for a wet mapping function, G N w and G E w are tropospheric gradients in the north and east directions, respectively. In this work, Vienna Mapping Function (Boehm et al. 2009) was used for m w . The hydrostatic gradients were estimated based on the horizontal distribution of surface pressure using the methodology described in Shoji (2013) and subtracted from the GNSS-derived G N and G E ; thus, only the wet parts were used in Eq. 3. In the literature, there is a method of improving the resulting SWD values by adding post-fit residuals of the GNSS analysis in order to include the unmodelled parts of the signals' delays (Manning et al. 2014;Alber et al. 2000;Troller et al. 2006). However, other studies show that the post-fit residuals play a minor role in the estimation procedure (Hordyniec et al. 2018); also, they could contain errors of the GNSS measurements, e.g. multi-path (Böhm and Schuh 2013). In our work, the post-fit residuals are not considered.
The resulting slant wet delay (SWD) can be also expressed as an integral of wet refractivity values N w along the path s of the GNSS signal (Mendes 1999): In the GNSS troposphere tomography model, Eq. 4 links the SWD observations with the estimated parameters N w (Flores et al. 2000). The estimated wet refractivity is a part of the total tropospheric refractivity N . N w corresponds to the non-hydrostatic components of the troposphere, whereas the hydrostatic refractivity N h refers to the hydrostatic components (Mendes 1999): The two parts of N can be expressed based on the tropospheric parameters as follows: where R d is a mean specific gas constant for dry air (J kg −1 K −1 ), ρ is a total air density (kg m −3 ), e stands for water vapour partial pressure (hPa), T is temperature (K), Z −1 w is a compressibility factor of the wet air (-)]; k 1 , k 2 , and k 3 are empirically determined refractivity constants (Thayer 1974).

Tomographic solution
In the GNSS troposphere tomography model, a network of GNSS receivers is used to determine the wet refractivity distribution within the model's domain. The following system of equations is applied: where SWD G is an observation vector that contains SWD values for each pair of the GNSS transmitter and receiver (SWD 1 , SWD 2 ,…, SWD l ); N w is a vector of estimated N w parameters (N w 1 , N w 2 ,…, N w m ); A G is a projection matrix that contains derivatives of observations with respect to the estimated parameters: In this work, a priori wet refractivity field is used to stabilise the solution. This information is included in the observation vector and the projection matrix as follows: where elements of the N w apr vector are the a priori values of wet refractivity; elements of A apr matrix equal 0 (no apriori information) or 1 (a priori information is available). The final system of equations reads as follows: The solution is obtained by inversion of Eq. 11. In this paper, the least squares method with the More-Penrose pseudoinverse ( + ) is used: where P stands for the weighting matrix, i.e. inverse of the covariance matrix of the SWD observations and a priori data. Solving the final system of equations leads to the estimation of the 3-D wet refractivity field in the tomographic model's domain.

Parameterisation of the tomographic domain
The tomographic technique is based on a discretisation of the model's domain, i.e. the estimated wet refractivity field is not continuous in space. Apart from the wet refractivity field, also the tropospheric wet delay SWD along the GNSS ray's propagation path is discretised. Depending on the parameterisation approach, it can be divided into segments based on voxels' distribution (one segment within each voxel; Flores et al. 2000) or based on horizontal layers of the domain (one segment between two tomographic layers; Ding et al. 2017). Taking into account a relation between SWD and N w presented in Eq. 4, the SWD in a tomographic model can be expressed as follows: where k is a number of segments, swd i is a value of wet tropospheric delay for segment i, s i stands for a signal's path in a segment i. A single SWD observation is therefore modelled as a sum of the integrals of N w within all segments along the signal's path. As these integrals are not calculated directly, the method of modelling the swd value within each segment depends on the domain's parameterisation approach. In this work, two parameterisation approaches were applied, i.e. voxel-based and node-based. The former is based on a simplification stating that the N w within each voxel is constant (see Sect. 2.3.1). The latter uses interpolation methods to link wet refractivity at the domain's nodes with SWD observations (see Sect. 2.3.2).

Voxel-based parameterisation
In the voxel-based approach, an assumption is made that wet refractivity within each voxel is constant (Flores et al. 2000;Rohm and Bosy 2011;Möller 2017). Therefore, an integral from Eq. 13 for a voxel i can be represented by multiplication of wet refractivity N w i within the voxel by a distance d i that the GNSS ray travels in this volumetric element. As a result, the swd i inside of the voxel i is modelled as follows: with the following units: swd i (mm), d i (km), N wi (ppm).
In this approach, the number of estimated parameters equals the number of voxels in the model's domain. Each element in the N w vector (Eq. 11) represents wet refractivity within one voxel. The projection matrix A G (Eq. 8) consists of distances that the GNSS rays travel in particular voxels.

Node-based parameterisation
In the node-based parameterisation approach, wet refractivity at any point in the tomographic domain (e.g. at the GNSS signal's path) can be defined by interpolation of N w from tomographic nodes to the point of interest. The integral from Eq. 13 is resolved by numerical integration. In this work, Newton-Cotes formula is applied, based on evaluating an integrand (wet refractivity) at equally spaced points along the signal's path (Perler 2011). For each of the points, the weight of the integrand is calculated using Lagrange polynomials. Finally, the integral is estimated as a sum of the integrands multiplied by the corresponding weights. In this study, Lagrange polynomials of range 4 are used (Boole's rule), i.e. the path of the ray in segment i is divided into 4 equal parts (sections). Moreover, a composite rule is applied, which means that the Newton-Cotes formula is also used inside each section (4 subsections per each section). As a result, for each segment of the GNSS signal's path, the integrand (wet refractivity) is evaluated at 17 equidistant points. Based on these points, the slant wet delay swd i is expressed as follows: where K equals 4 (number of sections); N w i ,0...4K +1 are wet refractivity values at the signal's path in segment i, corresponding to the 17 evaluated points; ∂s is a distance that a ray travels in one subsection. Wet refractivity values N w i ,0...4K +1 are linked to the N w values at the tomographic nodes using interpolation methods. As a result, elements of the projection matrix A G (Eq. 8) correspond to the weights of the integrands, as well as the interpolation parameters. In this work, two methods of wet refractivity interpolation are tested, i.e. trilinear and spline/bilinear.

Trilinear interpolation
In a trilinear interpolation method, it is assumed that a distribution of wet refractivity in the troposphere is characterised by linear changes in all three directions. As a result, values of wet refractivity at the GNSS signals' paths are modelled by linear interpolation in both horizontal and vertical directions. An advantage of this approach is its simplicity in application. However, as wet refractivity in the troposphere usually does not change linearly with height (Böhm and Schuh 2013), this method can introduce a systematic error, especially in a case of sparse distribution of the horizontal tomographic layers (i.e. large vertical distances between the tomographic nodes).

Spline/bilinear interpolation
Another approach is to use different interpolation methods for horizontal and vertical directions to better reflect a real distribution of wet refractivity in the troposphere. As a result, wet refractivity values at the signal's path (N w i ,0...4K +1 in Eq. 15) are calculated in two steps. First, a vertical interpolation is performed, based on tomographic nodes below (h k ) and above (h k+1 ) the point of interest. In this work, Fig. 1 A scheme of the 4-step algorithm (a-d) to find the optimal locations of the tomographic nodes in the model's domain, based on the 3D distribution of the GNSS rays wet refractivity value at height h is expressed using natural splines, based on the following cubical polynomial given by Perler (2011): with: where N w (h) stands for wet refractivity value at a given height h, h k is the height of the tomographic node located below the point of interest, dh is a vertical distance between two tomographic points located below and above the point of interest. N w k and N w k +1 are wet refractivity values in the tomographic nodes at heights h k and h k+1 , respectively, whereas N w k and N w k +1 are their second derivatives. As the wet refractivity at height h is expressed by the natural splines for the neighbouring tomographic nodes, the second step is performed, i.e. a bilinear interpolation in a horizontal direction to the location of the point on the signal's path.

Optimisation of the tomographic nodes' distribution
In the GNSS troposphere tomography, wet refractivity in the model's domain is estimated based on the tropospheric wet delays along the GNSS rays. In a standard tomographic approach, wet refractivity is evaluated for regularly dis-tributed volume elements or nodes. The spatial distribution of the points of reference (voxels, nodes) is defined based on the location of the GNSS stations (Bender and Raabe 2007) or the GNSS rays (Ding et al. 2018) in order to optimise the geometry of the solution. However, even though the distribution of the voxels or nodes is adjusted to the geometry of the GNSS observations, the requirement of regularity is a limitation that does not allow for a flexible selection of the optimal locations for the estimation of wet refractivity. In this study, another approach is presented that aims to identify the best locations of the tomographic nodes. The algorithm is based on the GNSS rays' distribution in the model's domain (Fig. 1). The algorithm is based on 4 steps, presented on the panels a-d. First (Fig. 1a), the most useful GNSS rays are selected. In the figure, the GNSS rays are presented as straight lines. As the performance of the tomographic model is the most efficient at locations where multiple GNSS rays are utilised, the algorithm rejects the signals that do not cross with the others (panel a, red lines). In consequence, the tomographic nodes will not be located in spots where very little information can be derived from the GNSS data.
In the second step ( Fig. 1b), the locations of intersections of the remaining rays are detected (red points in panel b).
Since the GNSS signals contain information about the tropospheric state, the locations with the signals' crossing are supposed to have the greatest potential to reflect the actual values of wet refractivity in the tomographic solution. To find the intersection points, the path of each signal is sampled with a distance d st . At each of the sample points along the ray, the neighbouring points from the other rays are searched within a radius d r . If the neighbouring points are detected, the intersection point is defined as their geometric centre. To choose the sampling distance d st , we have tested values in a range of 2 km-10 m. The tests indicated that the sampling resolution higher than 100 m introduced unnatural distribu- Fig. 2 Flow chart of the experiment tion of the intersections, i.e. the resulting intersecting points were located mainly at the heights corresponding to multiple of d st ). The sampling resolution of 100 m and denser did not introduce this unnatural distribution; as a result, a 100 m distance was recognised as optimal. The radius of the neighbourhood for the sampling points from adjacent rays d r was set based on the mean distance between the GNSS stations.
In the next step ( Fig. 1c), the intersection points defined in step b are clustered. As the intersecting points from step b could occur close to each other, the clustering helps to avoid aggregation of the tomographic nodes in a small area, with low inter-distances between them. In Fig. 1c there is an example distribution of the intersection points resulting from step b; as the result of step c, the points are assigned to particular clusters (one colour denotes points assigned to one cluster). In our approach, a density-based algorithm for spatial applications DBSCAN (Density-Based Spatial Clustering of Applications with Noise) is used (Ester et al. 1996). The algorithm is based on two parameters, i.e. a radius of the neighbourhood R and a minimum number of points required to form a dense region (N ). In our approach, the radius R equals the distance d r used for the identification of the intersections of the GNSS rays (from step b). The number N is set to 1 in order to take into account all intersections, even those not grouped with any others. An advantage of the DBSCAN algorithm is that the number of clusters does not have to be defined initially; thus, the resulting number of the tomographic nodes is adjusted to the geometry of the GNSS observations. The geometric centres of the resultant clusters are considered the best locations for the GNSS tomography nodes.
However, as the SWD values cannot be properly modelled only based on these sparse and irregularly distributed points, also the regularly distributed part is involved (Fig. 1d). The regular nodes are distributed as in the standard node-based tomographic approach; however, only the points in close vicinity of the GNSS rays are used (those that are employed for the SWD modelling of at least one GNSS ray). As result, the GNSS tomography nodes consist of two subsets, i.e. the irregularly distributed part, with an optimal location for the estimation of N w (red colour in Fig. 1d), and the regularly dis-tributed part, used for proper modelling of the GNSS SWD (blue colour in Fig. 1d).
After the 4 steps ( Fig. 1a-d) are performed (i.e. the location of the tomographic nodes is defined), then the tomographic solution is run using node-based parameterisation with spline/bilinear interpolation (Eq. 12).

GNSS troposphere tomography model
The GNSS troposphere tomography was performed using the TOMO2 model developed at the Wrocław University and Environmental and Life Sciences (Rohm andBosy 2009, 2011;Rohm 2013). The model is based on the tomographic principles and allows for the estimation of the 3-dimensional distribution of wet refractivity or water vapour density in the troposphere. In this research, we reconstructed wet refractivity fields in the area of Wrocław for the period of March 15-28, 2021. The main scheme of the conducted research is presented in Fig. 2. The GNSS tomography was performed based on two types of tropospheric estimates: real (GNSS-based; see Sect. 3.2) and synthetic (ray-traced using the WRF model; see Sect. 3.3). In both cases, the a priori wet refractivity fields were derived from the Global Pressure and Temperature 2 (GPT2) model (see Sect. 3.4). The TOMO2 model was designed for voxel-based parameterisation. Within this study, we adjusted the model to use node-based parameterisation with irregular distribution of tomographic nodes. The tomography was performed with four different parameterisation approaches, i.e. (1) regular grid with voxel-based parameterisation (yellow line); (2) regular grid with node-based parameterisation and trilinear interpolation (green line); (3) regular grid with node-based parameterisation and spline/bilinear interpolation (purple line); (4) irregular nodes with node-based parameterisation and spline/bilinear interpolation (red line). The four tomographic solutions were validated with the radiosonde data (see Sect. 3.5) and the WRF model (see Sect. 3.3). Figure 3 presents the spatial domain of this study with the location of the GNSS receivers, nodes of the WRF model grid, and the RS station. Also, the tomographic grid is presented. In the regular approach, the tomography was based on the 4 × 4 voxel structure with a horizontal resolution of 5 km, in 11 irregularly distributed vertical layers ranging from 150 m to 11,750 m. In the node-based regular approach, the same grid was used (5 × 5 nodes with 12 vertical layers). In the node-based irregular approach, the methodology presented in Sect. 2.4 was applied for optimisation of the tomographic nodes' distribution.

GNSS tropospheric delays
The tropospheric tomography performed in this study was based on the ZTD observations from the network of 16 low-cost multi-GNSS receivers located in the urban area of Wrocław city, Poland, and its suburbs (Fig. 3). A detailed description of the receivers and the GNSS data processing strategy can be found in Marut et al. (2022). The receivers are based on the u-blox high-precision GNSS module (ZED-F9P), microcomputer Raspberry Pi 3B+, and the GSM modem. The signals from GPS, Galileo, and GLONASS systems were tracked. The data were processed in two strategies, i.e. real-time and final (post-processing). In our study, we used the product based on the CSRS-PPP software and IGS Final orbits. The GNSS network covers an area of 243 km 2 (14 km North-South, 24 km East-West); the mean distance between the receivers is about 3 km. In the GNSS data processing, the ZTD values and their gradients were estimated. As the GNSS tomography deals with the wet delays, the hydrostatic part was calculated based on the pressure values from the WRF model (see Sect. 3.3) using the Saastamoinen's formula (Eq. 2) and subtracted from the ZTD estimates. Next, the SWDs were calculated based on the ZWDs and tropospheric gradients using Vienna Mapping Functions (Eq. 3).
The SWD values were provided with a temporal resolution of 1 h.

WRF model
Apart from the GNSS estimates of the tropospheric delays, also the synthetic SWDs were derived from the Weather Research and Forecasting (WRF) model using a raytracing method (Lasota et al. 2020). The WRF model is a mesoscale numerical weather prediction system designed for atmospheric research (Hsiao et al. 2012;Givati et al. 2011;Nasrollahi et al. 2012) and operational weather forecasting (Benjamin et al. 2016). The WRF model consists of two dynamical solvers, i.e. the Advanced Research WRF (ARW) core (Skamarock et al. 2008) and the Non-hydrostatic Mesoscale Model (NMM) core (Janjic et al. 2005). The WRF-ARW is developed and maintained by U.S. NCAR's (National Center for Atmospheric Research) Mesoscale and Microscale Meteorology Laboratory. The code is open to public use. In this study, we use a WRF-ARW model in version 3.7 run by the Department of Climatology and Atmosphere Protection at the University of Wrocław (Kryza et al. 2013). The model is configured with two one-way nested domains with a spatial resolution of 12 km and 4 km, respec-tively. In the outer domain, the convection is parameterized, while in the inner domain it is explicitly resolved. In this model, 48-hour forecasts are run every 6 h (at 0, 6, 12, and 18 UTC) with a resolution of 1 h. Initial and boundary conditions are provided by the U.S. National Oceanic and Atmospheric Administration (NOAA) National Centers for Environmental Prediction (NCEP) Global Forecast System (GFS), available every 3 h with a spatial resolution of 0.5 • . The model consists of 35 terrain-following hydrostatic pressure levels, with the top fixed at 10 hPa. In this research, the ray-traced SWD values were used as input into the tomographic model in the synthetic approach. Also, wet refractivity fields were calculated based on the WRF forecasts using Eq. 7 for validation of the tomographic outputs.

GPT2 model
GPT2 (Global Pressure and Temperature 2) is an empirical model for the determination of slant tropospheric delays, dedicated to radio space geodetic techniques (Lagler et al. 2013). The model is based on monthly mean profiles from ERA-Interim reanalysis (ECMWF 2011) for the years 2001-2010. GPT2 provides pressure, temperature, temperature lapse rate, water vapour pressure, and mapping functions coefficients; all parameters contain mean, annual, and semi-annual terms. The model is based on the underlying horizontal resolution of 5 degrees and the mean ETOPO5 (NOAA 1993) heights. In this study, the GPT2 parameters were used to calculate values of wet refractivity (Eq. 7) for initialisation and stabilisation of the tomographic solution (a priori data in Eq. 10). As the tomographic model needs a rough approximation of the tropospheric state to initialise the solution, the GPT2 model was used in order to provide a mean climatology of the troposphere.

Meteorological radiosonde station
Radiosonde (RS) data were derived from the NOAA/ESLR Radiosonde Database managed by the Earth System Research Laboratory (ESLR) of the National Oceanic and Atmospheric Administration (NOAA: http://www.esrl.noaa.gov/). In the RS technique, a sensor attached to an automatic radiosounding balloon is released and measures meteorological parameters in the vertical profile. In this research, measurements from the station located in Wrocław (ID 12374) were used. The profiles of temperature, dew point temperature, and wind speed components are provided twice a day (at noon and midnight) for the pressure levels in a range of 1000-10 hPa. There are 14 mandatory levels; however, depending on the meteorological conditions, more data can be provided. Usually, measurements for 30-40 levels are available. The radiosonde measurement is the most common and accurate technique of troposphere profiling. According to the E-GVAP-II Product Requirements Document (Offiler et al. 2010), radiosonde stations are one of the standard data sources used as a reference in the measurement of the state of the troposphere. The accuracy of these measurements is 1-2 hPa for pressure, 0.5 • C for temperature, and 5% for relative humidity (Manning 2013).

Meteorological conditions
The experiment was conducted for March 15-28, 2021. Figure 4 presents time series of mean values and standard deviations of the ZTDs estimated for the network of the lowcost receivers in the given period. Also, ZTD values derived from the NWP model (WRF) are presented. In general, the ZTD values in this period are in the range of 2.3-2.4 m. The standard deviations for particular epochs are usually at the level of 10 mm. For a considerable part of the period, there is a good agreement between the GNSS estimates and the WRF model; the discrepancies are in the range of the standard deviations. The only exception is March 25-27 when the two data sources differ more significantly (by more than 10 mm). In the first part, the tropospheric delays do not exceed the value of 2.34 m, whereas in the second part of the experiment (after March 23) the values are noticeably higher (2.34-2.38 m). In the middle of the period (March 20-23), sharp changes in the ZTD values are noticed for both data sources (by more than 40 mm during one day).
According to the Interdisciplinary Centre for Mathematical and Computational Modelling (ICM) at the University of Warsaw (UW) service (www.meteo.pl), the weather conditions for March 15-20 were as follows. The North circulation prevailed, bringing cold air masses from Scandinavia, and temperature oscillated within 5 to 7 • C; low-level clouds, mild rain, and snow were predicted (which corresponds to the lower ZTD values in Fig. 4). Starting from March 21, more varying weather occurred with a warm front moving through the country from northwest to southeast (a spike in ZTD in Fig. 4), bringing more diverse air masses, including warmer and moist originating from the Norwegian Sea. After

Results
This section is organised as follows. In the beginning, a short description of the irregularly distributed nodes is introduced (4.1). Further, a comparison between the four tomographic solutions with different parameterisation approaches is presented. We refer to those approaches in the following way: (1) constant (voxel-based, regular grid); (2) trilinear (node-based, regular grid, trilinear interpolation); (3) spline (node-based, regular grid, bilinear/spline interpolation); (4) irregular (node-based, irregular grid, bilinear/spline interpolation).

Distribution of the irregular nodes
The locations of the regularly and irregularly distributed tomographic nodes are presented in Fig. 5. The plot presents information about the nodes from all tomographic epochs. In the horizontal direction, the irregular nodes are distributed rather uniformly, with the highest density in the centre of the tomographic area (panel a). A lower density of the irregu-larly distributed nodes can be noticed close to the boundary of the area, especially in the longitude direction. The side view (panel b) shows the inhomogeneous distribution of the nodes. Near ground (below 1 km) there are almost no irregular nodes since the GNSS rays do not cross in this altitude range. The highest density of the irregular nodes is noticed for the 2-4 km; in this range, about 40-50% of all irregular nodes occur. Above 4 km, the number of irregular nodes decreases with height.
The irregular distribution of the nodes was also investigated in terms of changes with time. Figure 6 shows a number of irregular nodes in the whole tomographic domain, for the particular epochs; the number of nodes is presented as a function of the hour of the day. In general, the number of irregularly distributed nodes is within the range of 0-50 points per epoch, with a mean of about 20 points. In the regular solution, there are 176 voxels in the tomographic domain; thus, the number of irregular nodes corresponds to 0-28% of the regular voxels, with a mean of 11%. Some diurnal patterns can be noticed in the number of irregular nodes. At the midnight, there are usually more than 20 points per epoch, up to 50. For the hours 2-6 AM this number does not exceed  25 points; later, around 8 AM there is another pick, up to 40 points per epoch. In the second part of the day, the number of the tomographic nodes is in the range of 5-30 nodes, with a noticeable decrease around 2 PM (not higher than 10 nodes in each day). Since the GPS constellation repeats its geometry each sidereal day (close to 24 h; Agnew and Larson 2007), the distribution of the GPS rays also follows this sidereal day pattern. As a consequence, the distribution of the irregular nodes, designated by the methodology given in Sect. 2.4, is to some extent regular in time. However, since the repetition time of the GPS satellites' positions is slightly shorter than 24 h, the particular graphs in Fig. 6 are shifted in time (compare, e.g. March 15 and March 28).

Validation of the reconstructed SWDs
As the parameterisation of the tomographic domain affects the way how the SWDs are modelled, in the first step we compare the reconstructed SWD values from the four different parameterisation approaches. The SWDs were calculated based on the wet refractivity fields from the WRF model and compared to the values provided by the ray-tracing technique. Figure 7 presents mean bias (μ) and standard deviation (σ ) of SWD in the consecutive epochs. The most common parameterisation approach (constant) is indicated by the yellow line. In this case, the mean bias is in general on the level of 5 mm; the lowest values are noticed at the beginning of the period (March 15-19; lower than 5 mm) and slightly higher values are observed at the end of the period (March 23-18; 5-7 mm).
In the period of March 20-23 (grey rectangle on the plot) the largest fluctuations of bias are noticed (between − 10 mm and 7 mm). The pattern of the bias fluctuations in the trilinear approach is similar to the constant one; however, the The statistics are presented for all observations from the analysed period (March 15-28, 2021) values are systematically larger by 2-5 mm. The spline and the irregular solutions are characterised by the lowest bias from all of the 4 approaches; the values are on the level of 0-2 mm at the beginning of the period and − 1 to 3 mm at the end of the period. These two approaches based on the spline vertical interpolation (spline and irregular) result in almost the same values of the reconstructed SWD observations. In terms of the standard deviation, the constant and trilinear solutions are similar, on the level of 0-5 mm. The spline and irregular approaches outperform the other two solutions; their standard deviation is in general in the range of 0-2 mm with some peaks up to 2.5-5 mm, especially during March 20-23. The overall statistics of the reconstructed SWD values validated based on the ray-traced data are presented in Table 1. The trilinear approach is characterised by the highest values of errors (mean bias, standard deviation, root-mean-square error) from all solutions. It is especially noticed in terms of the systematic error which is 2.7 mm for the constant solution, 5.9 mm for the trilinear one, and almost zero for the spline and irregular (− 0.2 mm and − 0.3 mm, respectively). The lowest standard deviation is reached for the spline and irregular solutions, i.e. 3.1 mm (decrease by 1.8 mm with respect to the constant solution). As a consequence, in terms of the rootmean-square error (RMSE) the quality of the reconstructed SWDs is the best for the spline and irregular parameterisation approaches (3.1 mm) and the worst for the trilinear approach (7.5 mm), while the RMSE of the constant solution equals 5.6 mm.

Validation of the wet refractivity fields: synthetic solution
The wet refractivity fields from the tomographic solutions based on the synthetic (ray-traced) SWD data were compared to the wet refractivity fields from the WRF model. The bias and standard deviation were calculated for each epoch and are presented in Fig. 8. Also, the statistics of the a priori fields derived from the GPT2 model were calculated. At the beginning of the analysed period (March 16 and 17), the quality of the a priori data was on the level of − 4 to − 6 ppm in terms of bias and 6-8 ppm in standard deviation. The GNSS observations improved the a priori fields and the errors Fig. 8 The statistics of the wet refractivity fields estimated using the four tomographic parameterisation approaches (synthetic solution) and derived from the GPT2 model (a priori data). The validation was performed based on the WRF model. The statistics are presented for the consecutive epochs of the analysed period (March 15-28, 2021). The shaded region corresponds to the period where the SWD values were poorly reconstructed for all tomographic approaches decreased to about − 2 ppm (bias) and about 4-6 ppm (standard deviation). In the period of March 18-19, the quality of the GPT2 model increased; the bias and standard deviation were on the level of − 2 to 1 ppm and 2-4 ppm, respectively. In this period, the tomographic technique did not improve the a priori wet refractivity fields for any of the parameterisation approaches. During the last week of the analysed period (March 23-28) the quality of the a priori data is lower, i.e. up to − 10 ppm in bias and 10 ppm in standard deviation; the GNSS tomography solutions were in general better than the GPT2 by 4-8 ppm in bias and 1-2 ppm in standard deviation. The largest discrepancies between the four tomographic solutions are noticed in the case of systematic error. The trilinear approach introduces the largest bias for the whole analysed period; it is usually 1-2 ppm larger than in the case of the constant parameterisation. The bias of the spline method is usually the same or worse than constant by 0-1 ppm. In the case of an irregular approach, the bias is the same or better than in the constant parameterisation (improvement by 0-1 ppm). In terms of the standard deviation, all tomographic solutions have similar results; the differences are usually not larger than 1 ppm. The trilinear and spline solutions are in general characterised by a slightly larger standard deviation than the constant and irregular approaches (differences not larger than 1 ppm). The cumulative statistics of the 3-D wet refractivity fields derived from the tomographic synthetic solutions and the GPT2 model, validated based on the WRF model, are presented in Table 2. The results show that all tomographic runs improved the a priori (GPT2) fields of the wet refractivity in terms of bias, standard deviation, and RMSE. In terms of bias, the largest improvement is noticed for the irregular (− 1.7 ppm) and constant (− 1.9 ppm) approaches. The largest bias was found for the trilinear solution (−2.4 ppm); however, even in this case, there is an improvement with respect to the a priori data (− 3.6 ppm, GPT2). In terms of the standard deviation, the smallest value was found for the irregular case (4.8 ppm). The regular approaches (constant, trilinear, spline) do not differ much from each other and are characterised by the standard deviation on the level of 5.1− 5.2 ppm. In consequence, the 3-D wet refractivity fields from the irregular approach have the lowest value of the RMSE (5.0 ppm) from all 4 tomographic runs.
The tomographic data were validated based on the WRF model within the particular height ranges (Fig. 9). In the case of the regular solutions (constant, trilinear, spline), the statistics were calculated for the consecutive tomographic layers. In the case of the irregular solution, the statistics were calculated for all observations between the two consecutive tomographic levels. The largest improvement of the a priori fields by the tomographic technique is noticed in the lower parts of the troposphere, i.e. between 0 and 3 km height. In this range, the GPT2 N w fields are characterised by the bias on the level from − 10 to − 1 ppm. The largest improvement is noticed in the case of the spline and irregular solutions; the bias is reduced by 5 ppm in the height range 0.5-1.5 km and by 2 ppm in the height range 1.5-2.5 km. From all tomographic parameterisation approaches, the largest bias is noticed in the trilinear approach (larger than spline and irregular by 2-3 ppm). In terms of RMSE, the irregular solution outperforms the other tomography parameterisation approaches in the height ranges 0.2-2 km. For the higher altitudes, the irregular and spline solutions are on the same RMSE level (2-4 ppm). The constant and trilinear approaches have RMSE error larger by 1-3 ppm than the two other tomographic solutions. In general, the RMSE of the a The validation was performed based on the WRF model. The statistics are presented for all observations from the analysed period (March 15-28, 2021)

Fig. 9
The statistics of the wet refractivity fields estimated using the four tomographic parameterisation approaches (synthetic solution) and derived from the GPT2 model (a priori data). The validation was performed based on the WRF model. The statistics are presented for the particular height ranges from 0 to 6 km priori GPT2 field is improved by the tomography solutions on the level of 1-5 ppm. The statistics calculated for the particular height ranges are also presented in Table 3. In the bottom level (150 m for trilinear, spline, irregular; 275 m for constant) the improvement of bias with respect to the a priori N w (GPT2) is the largest in the case of the constant solution (bias reduced by 2.6 ppm). Since in the voxel-based solution (constant) the bottom layer is higher than in the node-based solutions (trilinear, spline, irregular), more GNSS signals are affecting the results in this part of the troposphere. However, the improvement of standard deviation in the bottom layer is the largest in the node-based solutions (5.7-5.8 ppm for trilinear and irregular while 6.3 for constant). In the lower and middle troposphere (0.3-4.0 km) the largest impact of the tomographic technique on the a priori N w fields is noticed. In this height range, both bias and standard deviation are the lowest for the irregular parameterisation approach (− 2.3 ppm and 4.8 ppm, respectively). In the middle and upper troposphere (4.0-12.0 km) the impact of the tomographic technique on the a priori data is negligible (differences of bias and standard deviation on the level of 0.0-0.2 ppm).

Validation of the wet refractivity fields: real solution
In this section, the validation of the wet refractivity fields from the tomographic solution based on the real GNSS SWD estimates is presented. The validation is performed for the four tomographic parameterisation approaches and the a priori data from the GPT2 model. The comparison is based on the WRF model (Sect. 4.4.1) and the radiosonde observations (Sect. 4.4.2).

Validation based on the WRF model
The systematic error and the standard deviation of the tomographic runs and the GPT2 model, calculated based on the comparison with the WRF model, are presented in Fig. 10. The pattern of the tomographic statistics changes is similar to the synthetic solution (Fig. 8); however, the discrepancies between the a priori data (GPT2) and the tomographic estimates are smaller than in the synthetic case. The largest improvement is noticed for the first two days of the analysed period (March 15-16) and the last week (March 23-28). During these periods, the tomographic technique improved the a priori N w fields by 1-4 ppm in terms of systematic error The validation was performed based on the WRF model. The statistics are presented for the particular height ranges Fig. 10 The statistics of the wet refractivity fields estimated using the four tomographic parameterisation approaches (real solution) and derived from the GPT2 model (a priori data). The validation was performed based on the WRF model. The statistics are presented for the consecutive epochs of the analysed period (March 15-28, 2021). The shaded region corresponds to the period where the SWD values were poorly reconstructed for all tomographic approaches and 1-3 ppm in terms of standard deviation. The differences in bias between the particular parameterisation approaches are on the level of 1-2 ppm; the best results are found for the constant and irregular approaches, while the worst for the trilinear solution. The standard deviation is similar for all tomographic runs (differences not larger than 1 ppm). During March 17-23, no improvement of the a priori data is found; the tomographic estimates are worse than GPT2 by 1-5 ppm in bias and 1-6 ppm in standard deviation. The worst results are noticed for the period where the SWDs were poorly reconstructed by the tomographic model (March 20-23, indicated by the grey rectangle; see Fig. 7). The overall statistics indicate that all four tomographic solutions improved the a priori data from the GPT2 model (Table 4). The largest improvement is noticed in terms of bias in the irregular case (an improvement from the − 3.6 ppm GPT2 to the − 2.6 ppm irregular). The other tomographic parameterisation approaches result in a bias in the range from − 2.9 ppm (constant) to − 3.4 ppm (trilinear). The discrepancies in standard deviation are smaller; the best results are found for the constant and irregular solutions (5.8 ppm). The RMSE is the lowest for the irregular parameterisation approach (6.3 ppm).

Fig. 11
The statistics of the wet refractivity fields estimated using the four tomographic parameterisation approaches (real solution) and derived from the GPT2 model (a priori data). The validation was performed based on the WRF model. The statistics are presented for the particular height ranges from 0 to 6 km The validation of the tomographic results on the particular height ranges shows that the discrepancies between the a priori data and the estimated N w fields are in the range of 0-4 ppm (Fig. 11). The pattern of the statistics changes with height is similar to the synthetic solution (see Fig. 9); however, the impact of the GNSS observations is smaller than in the case of the ray-traced observations. The largest improvement of the bias with respect to the a priori data is found for the irregular solution (by 2-3 ppm in the range of 0.5-1.5 km) and the worst for the trilinear one (by 0-1 ppm in the same range). The RMSE is the lowest for the irregular and spline solution (improvement by 2-3 ppm with respect to the a priori data). Above the height of 4 km, the tomographic results are in general the same as the a priori fields.
In the bottom layer of the tomographic solutions (150 m for the trilinear, spline, and irregular; 275 m for the constant), the improvement with respect to the GPT2 data is the most visible in the case of irregular parameterisation approach (bias decreased by 0.9 ppm, standard deviation decreased by 1.1 ppm; Table 5). For the other tomographic solutions, changes in the statistics in the bottom layer are not larger Table 4 The statistics of the wet refractivity fields estimated using the four tomographic parameterisation approaches (real solution) and derived from the GPT2 model (a priori data) σ N w (mm) 6.6 5 .8 6 .1 6 .0 5 .8 RMSE N w (mm) 7.5 6 .5 7 .0 6 .7 6 .3 The validation was performed based on the WRF model. The statistics are presented for all observations from the analysed period (March 15-28, 2021) σ N w (ppm) 6.8 5.8 6.2 6.1 6.1 4000-12,000 m μN w (ppm) 0.9 1.0 0.8 0.9 0.9 σ N w (ppm) 1.1 1.0 0.8 0.8 0.8 The validation was performed based on the WRF model. The statistics are presented for the particular height ranges than 0.6 ppm. In the lower and middle troposphere (0.3-4.0 km) the lowest value of bias is noticed in the case of an irregular grid (improvement by 1.2 ppm). The best standard deviation is found for the constant approach (improvement by 1.0 ppm); however, differences in this parameter are very small between the four tomographic solutions (not larger than 0.3 ppm). In the middle and top troposphere (4.0-12.0 km) the discrepancies in the statistics between the GPT2 and the tomographic solutions are small, not larger than 0.3 ppm.

Validation based on the radiosonde observations
Another validation of the real (GNSS-based) tomographic solution was performed, using the radiosonde observations. The tomographic estimates of the wet refractivity from the whole domain were compared to the measurements from the one radiosonde station (see Fig. 3). Since the location of the radio-sounding balloon can change during the measurement by tens of kilometres (Wang and Zhang 2008;Seidel et al. 2011) and the dimensions of the tomographic area are relatively small (20 km), the observations from the whole tomographic domain were taken into account. The statistics of the comparison are presented in Fig. 12. The systematic error was improved with respect to the a priori data (GPT2) by 3-6 ppm at the beginning of the analysed period (March 15-17) and during the last week (March 23-28).
The largest improvement of bias is noticed for the irregular approach (1-2 ppm better than the other tomographic solutions). For March 18-20 the quality of the GPT2 model was much higher (bias in the range of − 2 ppm to 2 ppm) and the GNSS observations did not improve the a priori fields in this period. However, also during this period, the irregular tomographic solution was better than the others (bias − 2 to − 1 ppm). In terms of standard deviation, the tomographic technique improved the a priori fields by 1-2 ppm at the Fig. 12 The statistics of the wet refractivity fields estimated using the four tomographic parameterisation approaches (real solution) and derived from the GPT2 model (a priori data). The validation was performed based on the radiosonde measurements. The statistics are presented for the consecutive epochs of the analysed period (March 15-28, 2021). The shaded region corresponds to the period where the SWD values were poorly reconstructed for all tomographic approaches beginning and the end of the analysed period. In the epochs where the standard deviation of the GPT2 model was lower than 4 ppm (March 17-23), no improvement by the tomography was noticed. Almost through the whole period, the irregular solution had the lowest standard deviation from all tomographic solutions (lower by 1-2 ppm than constant and trilinear; lower by 0.5 ppm than spline in selected epochs). Only on March 26 the standard deviation of the constant solution was better than irregular (by 1 ppm). The quality of the wet refractivity fields from the tomographic solutions and the GPT2 model, compared to the RS data, is also presented in Table 6. The largest impact of the GNSS observations on the a priori fields is noticed in terms of RMSE N w (mm) 6.9 6 .4 6 .3 6 .0 5 .8 The validation was performed based on the radiosonde measurements. The statistics are presented for all observations from the analysed period (March 15-28, 2021) Fig. 13 The statistics of the wet refractivity fields estimated using the four tomographic parameterisation approaches (real solution) and derived from the GPT2 model (a priori data). The validation was performed based on the radiosonde measurements. The statistics are presented for the particular height ranges from 0 to 6 km bias in the irregular solution (improvement by 1.3 ppm, while for other tomographic solutions it is 0.3-0.9 ppm). The standard deviation was only slightly improved (by 0.8-0.5 ppm); the lowest value is found for the irregular solution (5.6 ppm); however, the differences between tomographic solutions are small (on the level of 0.1-0.3 ppm). In terms of RMSE, the best tomographic solution is the irregular one (5.8 ppm). In terms of bias, the irregular solution is the best from all tomographic runs, for the whole vertical profile (Fig. 13). In this parameterisation approach, the improvement of the bias with respect to the GPT2 model is on the level of 3.0-0.5 ppm in the height ranges of 0-3 km. The largest bias is noticed for the trilinear solution, in the whole vertical profile (improvement only by 0.0-1.0 ppm with respect to GPT2). However, in terms of RMSE, in the lowest parts of the troposphere (0-1.5 km) the constant solution is slightly better than the irregular (by 0.5 ppm), while in the higher parts (1.5-4.0 km) the irregular and spline solutions have the lowest RMSE.
In the bottom layer of the tomographic domain, the constant solution improved the GPT2 model only slightly in terms of bias (0.8 ppm) and did not change the standard deviation (Table 7). In the case of node-based approaches, the first tomographic layer is below the lowest radiosonde measurement; thus, there is no comparison for these cases. In the lower and middle troposphere (0.3-4.0 km) the bias was improved by 1.5 ppm for the irregular approach and by 0.4-1.3 ppm for the other tomographic runs. In terms of standard deviation, the constant solution was the best in this height range (6.1 ppm), while the node-based solutions did not differ much from each other (6.6-6.7 ppm). Above the 4 km altitude, no significant improvement of the GPT2 model was noticed (changes by 0.0-0.1 ppm).
In order to check whether the distance from the radiosonde station might impact the resulting statistics, another comparison was performed-based on the tomographic nodes not further than 4 km from the RS station. The results are similar to those presented in this chapter; the standard deviations and bias of the tomographic solutions are slightly lower than in the current comparison, however, with the same patterns of changes with time and height. The results can be found in Appendix A (Figures 15,16 and Tables 8,9).

Profiles of wet refractivity
The tomographic results were also examined for individual wet refractivity profiles in the vertical direction (Fig. 14). Three various cases were selected based on the characteristics of the RS wet refractivity profiles: a calm profile (panel a), a profile with several inversions (panel b), and a profile with one inversion (panel c). The first case presents wet refractivity data that reach 35 ppm near the ground and decrease regularly with height. In this case, the GPT2 model (a priori for the tomographic solution) indicates values of refractivity lower than 20 ppm near the ground, with a gradual decrease with height. The regular tomographic solutions outperform the GPT2 model by 10-15 ppm at 1 km altitude. In the case of the irregular solution, the irregular nodes occur at altitudes of 2.5-10 km; the wet refractivity values at these points do not differ significantly from the regular solutions (discrepancies by 1-2 ppm). In the second case (panel b.), the general shape of the wet refractivity profile is similar to the previous one, but with several inversions at altitudes 2-4 km. The GNSS observations improved the a priori fields; however, the inversions were not detected by any of the tomographic solutions. The last example (panel c.) presents a wet refractivity profile starting from 25 ppm near the ground, reaching 2-4 ppm at altitudes 2.5-3.5 km and 10 ppm at altitude 4 km. The validation was performed based on the radiosonde measurements. The statistics are presented for the particular height ranges In the lower troposphere, the best results are noticed for the spline tomographic approach, for which the wet refractivity reaches 25-27 ppm. The trilinear solution was slightly worse (22-23 ppm near the ground), whereas the constant approach did not change the a priori fields significantly (discrepancies between the constant solution and GPT2 model are not larger than 1-2 ppm). In the middle part of the profile, all regular tomographic solutions do not differ significantly from the a priori GPT2 data. In the case of the irregular solution, the wet refractivity estimated in the tomographic nodes is closer to the RS profile than the GPT2 model for the altitudes 2-3 km; although the values do not reach those from the RS data, they are the most successful from the all tomographic solutions.

Discussion
The conducted experiment was a base to analyse the tomographic results in the context of (1) the impact of the voxel-based or node-based parameterisation approach; (2) the impact of the interpolation method in the node-based approach; (3) the optimisation of the tomographic nodes' distribution; (4) use of the tropospheric estimates from a network of low-cost GNSS receivers. The impact of the voxel-based or the node-based parameterisation approach is clearly visible in the results of both SWD reconstruction (Sect. 4.2) and wet refractivity estimation (Sects. 4.3,4.4). In the case of the SWD reconstruction, a positive bias of about 2.7 mm is introduced by the voxelbased approach, while for the node-based it is − 0.2 mm (spline interpolation). The positive bias might be a result of the choice of the representative point within the voxel in order to model the SWD value (here, a central point of the voxel). As the wet refractivity in general (in calm weather conditions) decreases exponentially with height (Boehm et al. 2009), the choice of the central point within a voxel as a representative one might introduce a systematic error. This impact is also noticed in the final tomographic estimations of N w where a systematic error is about 1-2 ppm larger in the voxel-based than in the node-based (bilinear) approach. Similar results can be found in the work of Perler et al. (2011) where the voxel-based and the node-based estimations were compared. In that experiment, a systematic error of the voxelbased approach was more than two times larger than in the node-based (0.010 ppm and 0.004 ppm, respectively) for the synthetic solution. In the same work, a comparison based on a real GNSS network gave mixed results, i.e. lower bias in the voxel-based solution for the radiosonde measurements at midnight (0 UTC) but higher at noon (12 UTC). However, the experiment was conducted based on two validation epochs only. It is worth noticing that in our study not only the systematic errors but also the standard deviations of the reconstructed SWDs are smaller in the case of the node-based tomographic approach.
In the node-based approach, the wet refractivity values are estimated in the nodes of the tomographic grid. On the contrary, in the voxel-based approach, the wet refractivity values are estimated for the whole 3-D blocks in the tomographic domain. As a result, in the voxel-based approach, the N w fields have to be transformed from the block form into the grid form in order to validate them with external data sources. At the validation step, it introduces representation error of the tomographic data which might impact the assessment of their quality. Besides the validation issues, the same problem occurs in the assimilation step of the tomographic data. Usually, when the voxel-based tomographic results are assimilated into the NWP models, the N w within the voxel is assigned to the central point and then assimilated into the grid-based NWP model (Hanna et al. 2019;Trzcina and Rohm 2019;Trzcina et al. 2020). Since the representation error impacts the results of the voxel-based approach more than the node-based one, it would be reasonable to use the latter in the data assimilation applications.
Although the errors of the node-based parameterisation approach are lower than in the voxel-based solution, it is only in the case of the spline vertical interpolation. The linear vertical interpolation (i.e. trilinear solution in this study) introduces positive systematic error in the SWD reconstruction (5.9 mm) and negative in the N w estimation (− 2.4 mm synthetic solution, − 3.4 mm real solution). Also, the standard deviation of the trilinear approach (SWDs reconstruction, N w estimation) was the highest from the all solutions. To some extent, the results correspond well with the conclusions of Perler (2011) who found that the spline vertical interpolation gives better results than the trilinear. On the other hand, in the same work, the trilinear approach overperformed the voxel-based one, which is contrary to our results. However, in that research, the wet refractivity values were validated based on 1441 points distributed along one vertical profile in the centre of the tomographic domain. As a result, in the work of Perler (2011), the wet refractivity value estimated for one voxel is taken into account several times at the validation step, although it is constant for the whole voxel. In our case, we compared the estimated values at the heights of the tomographic layers (11 for the voxel-based and 12 for the node-based approach; see Fig. 3), for the whole domain, in all epochs. Thus, the differences between those results might be caused by the validation method. The results that we obtained correspond well with the physics of the troposphere where the humidity in general changes exponentially with height. If the linear interpolation is applied in the vertical direction, a positive bias should occur; since the tomography involves an inversion process, the final estimations are negatively biased.
In this work, we introduced a method of optimising the distribution of the tomographic nodes based on the distribution of the GNSS rays within the tomographic domain. In the literature, there are no works considering the irregular distribution of the tomographic nodes. However, there are works of Ding et al. (2017Ding et al. ( , 2018 and Zhang et al. (2020b) who propose to adjust the domain's boundaries to the distribution of the GNSS rays' location and then apply meshing techniques to define tomographic nodes. In our case, introducing the additional irregularly distributed nodes resulted in a decrease in the systematic error and standard deviation of the tomographic solution, for both synthetic and real data. The improvement was not high (up to 0.5 ppm when compared to the regular spline solution); however, it was consistent for all comparisons (synthetic and real data validated by the WRF model, and real data validated by the RS observations). Also, the detection of the inverse in the wet refractivity vertical profile was the most successful in the case of the irregular approach. In the other works concerning moving away from the regular cuboid grid in the GNSS troposphere tomography (Ding et al. 2017(Ding et al. , 2018Zhang et al. 2020b) the improvement was similar to our results, on the level of 30-40% (0.1-0.5 g m −3 ), mainly due to modifications in the lower troposphere (where regular cuboid grid introduces problems in the proper wet refractivity estimation). The mentioned works were focused on the adjustment of the domain's boundaries into the GNSS rays distribution, while our approach introduces additional irregularly distributed nodes, in the optimal location (close to the intersecting GNSS rays); thus, these two approaches are complementary and could be combined in future research. Also, the new method would be even more effective if observations from a number of different receivers (with different geometry) were added, e.g. mobile GNSS receivers, signals from the LEO satellites crossing the troposphere, and observations from High Altitude Platform Station (HAPS).
The GNSS troposphere tomography experiment was performed based on observations from a dense network of lowcost multi-GNSS receivers located in the area of Wrocław and its suburbs, with a mean distance between the receivers equal to 3 km. The presented work shows that using low-cost GNSS receivers for the 3-D reconstruction of wet refractivity gives positive results. The initial wet refractivity fields derived from the GPT2 model were improved by the GNSS observations in all presented cases (validation based on the WRF model and the RS observations). The best improvement (up to 4 ppm) was noticed in the altitude range of 0.2-3.0 km, contrary to the works based on the geodetic-grade GNSS stations, where the largest improvement is noticed for a range of 2.0-6.0 km (Möller 2017;Hanna et al. 2019) or even up to 8.0 km (Brenot et al. 2019). This discrepancy in altitude range comes from different characteristics of the GNSS rays' geometry within the tomographic domains based on the dense network (3 km) and the standard GNSS networks (tens of kilometres). In our case, the network of 16 receivers on an area of 14 × 24 km provides dense observations in the lower parts of the troposphere (signals' intersections are present mostly at the 0.2-3.0 km altitude) that enable to reconstruct wet refractivity on the height levels that are problematic in the sparse GNSS networks, i.e. below 2 km. The small interstation distance results in low variability of the water vapour content between the stations at altitudes higher than 3 km which might be the reason for a limited impact above this threshold. Also, the topography of Wrocław city is rather flat which results in similar altitudes of all GNSS receivers within the network (maximum height difference between the stations equal to 25 m); it is a limiting factor in the tomographic reconstruction of the wet refractivity fields (Rohm andBosy 2009, 2011). In the paper by Marut et al. (2022) it was proven that the low-cost GNSS receivers are able to provide 2-D information about water vapour distribution in the troposphere on the level of 1.0-1.5 kg m −2 IWV. It shows the ability of the low-cost receivers to provide water vapour distribution on the city scale, which was also shown in our research. The low-cost receivers could be used in future tomographic studies to densify the geodetic-grade GNSS networks in order to improve the solution in the lower parts of the troposphere (below 2.0 km altitude).

Conclusion
The GNSS tropospheric tomography provides 3-D fields of wet refractivity or water vapour density which are of great importance for weather prediction systems. The first results of the assimilation of the tomographic outputs into the NWP models are encouraging for the further development of this technique. Recently, several studies were conducted on the parameterisation of the tomographic domain, introducing two main approaches, i.e. voxel-based and node-based. Also, experiments with adjustment of the domain's boundaries to the distribution of the GNSS rays in space have been presented. However, all of these approaches are based on the regular distribution of the tomographic nodes within the domain. Since the distribution of the GNSS rays in the troposphere is not regular, we propose a new parameterisation approach with optimisation of the tomographic nodes' loca-tion. In this work, a cross-comparison of four tomographic solutions is performed in order to analyse them in the context of parameterisation. The experiment is conducted using a dense network of low-cost multi-GNSS receivers in Wrocław and its suburbs, for the 2-week period including calm weather and heavy precipitation event.
The GNSS troposphere tomography performed based on the dense network of the low-cost receivers improved the initial GPT2 model wet refractivity fields by up to 4.0 ppm. The best results were achieved for the altitudes of 0.2-3.0 km. The GNSS ray's geometry within the dense low-cost network (most intersections of the rays on the altitudes 0.2-3.0 km) enables the reconstruction of wet refractivity on the height levels below 2 km, which is usually not possible using sparse geodetic-grade networks. Future studies could be focused on the densification of the geodetic-grade GNSS networks by the low-cost receivers in order to better resolve lower parts of the troposphere.
The new approach of the GNSS troposphere tomography parameterisation with an optimisation of the tomographic nodes' distribution gave the best results from the all presented approaches. The improvement was noticed with respect to both voxel-based and regular node-based approaches by 0.5-2.0 ppm in terms of wet refractivity RMSE, especially for the altitudes 0.5-2.0 km. The new proposed approach shows potential for further use with a larger number of observations from various sensors, characterised by different geometries (e.g. signals from the LEO satellites) which could make an even larger impact on the optimal distribution of the tomographic nodes. Also, the new method could be combined with the approach based on the optimisation of the tomographic domain's boundaries. The experiment shows the potential of the node-based tomographic parameterisation which should be used in the context of assimilation of the tomographic results into the NWP models. National Science Centre project no. UMO-2020/37/B/ST10/03703 and the Wroclaw Center of Networking and Supercomputing (http://www. wcss.wroc.pl/) computational grant using MATLAB Software License no. 101979.

Data availability
The datasets generated or analysed during the current study are available from the corresponding author upon reasonable request.

Conflict of interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Ethics approval Not applicable.

Consent for publication Not applicable.
Code availability The code of the GNSS troposphere tomography model TOMO2 is available from the corresponding author upon reasonable request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Appendix A: Tomographic data validation based on the radiosonde measurements: nodes from the distance of 4 km from the RS station See Figures 15, 16 and Tables 8, 9. Fig. 15 The statistics of the wet refractivity fields estimated using the four tomographic parameterisation approaches (real solution) and derived from the GPT2 model (a priori data). The validation was performed based on the radiosonde measurements. The statistics are presented for the consecutive epochs of the analysed period (March 15-28, 2021), for the nodes located not more than 4 km from the radiosonde station. The shaded region corresponds to the period where the SWD values were poorly reconstructed for all tomographic approaches Fig. 16 The statistics of the wet refractivity fields estimated using the four tomographic parameterisation approaches (real solution) and derived from the GPT2 model (a priori data). The validation was performed based on the radiosonde measurements. The statistics are presented for the particular height ranges from 0 to 6 km, for the nodes located not more than 4 km from the radiosonde station RMSE N w (mm) 7.0 6 .3 7 .0 6 .5 5 .5 The validation was performed based on the radiosonde measurements. The statistics are presented for all observations from the analysed period (March 15-28, 2021), for the nodes located not more than 4 km from the radiosonde station Table 9 The statistics of the wet refractivity fields estimated using the four tomographic parameterisation approaches (real solution) and derived from the GPT2 model (a priori data) The validation was performed based on the radiosonde measurements. The statistics are presented for the particular height ranges, for the nodes located not more than 4 km from the radiosonde station