Cross-Comparison and Methodological Improvement in GPS Tomography

: GPS tomography has been investigated since 2000 as an attractive tool for retrieving the 3D ﬁeld of water vapour and wet refractivity. However, this observational technique still remains a challenging task that requires improvement of its methodology. This was the purpose of this study, and for this, GPS data from the Australian Continuously Operating Research Station (CORS) network during a severe weather event were used. Sensitivity tests and statistical cross-comparisons of tomography retrievals with independent observations from radiosonde and radio-occultation proﬁles showed improved results using the presented methodology. The initial conditions, which were associated with di ﬀ erent time-convergence of tomography inversion, play a critical role in GPS tomography. The best strategy can reduce the normalised root mean square (RMS) of the tomography solution by more than 3 with respect to radiosonde estimates. Data stacking and pseudo-slant observations can also signiﬁcantly improve tomography retrievals with respect to non-stacked solutions. A normalised RMS improvement up to 17% in the 0–8 km layer was found by using 30 min data stacking, and RMS values were divided by 5 for all the layers by using pseudo-observations. This result was due to a better geometrical distribution of mid- and low-tropospheric parts (a 30% coverage improvement). Our study of the impact of the uncertainty of GPS observations shows that there is an interest in evaluating tomography retrievals in comparison to independent external measurements and in estimating simultaneously the quality of weather forecasts. Finally, a comparison of multi-model tomography with numerical weather prediction shows the relevant use of tomography retrievals to improving the understanding of such severe weather conditions. alternative is to find m′ , which minimises the L 2 norm using the minimisation factor ξ , which is expressed Different standard techniques exist for solving such linear inverse problems, e.g., singular value decomposition (SVD), truncated singular value decomposition (TSVD), weighted and damped least-square adjustment (LS), Kalman-filter, simultaneous algebraic reconstruction technique (SART), conjugate gradient method, Tikhonov-regularisation, stacking, and double-difference–spline interpolation. This study uses five models with different inversion techniques, i.e., SVD, TSVD, LS, and SART. Retrievals were obtained mainly in 3D and for one model in 2D (the Technical University of Ostrava (TUO) model), and these were proceeded with for 𝑁 (cid:3050) and 𝜌 (cid:3050)(cid:3049) , except for the Technische Universität Wien (TUW) model ( 𝑁 (cid:3050) only) and the University of Beira Interior (UBI) model ( 𝜌 (cid:3050)(cid:3049) only). GNSS for Severe Weather and Climate monitoring) which aims to study the use of GNSS tropospheric products for high resolution NWP and severe weather forecasting (http: // gnss4swec.knmi.nl / wg2). This investigation was also a contribution to the International Association of Geodesy (http: // www.iag-aig.org) and the Solar-Terrestrial Centre of Excellence. We thank Roeland Van Malderen and Alex Deckmyn from the Royal Meteorological Institute of Belgium for providing radiosonde data and ALARO NWP outputs.

can take place seasonally, i.e., in the northern equatorial region, which has a tropical climate, and in the southern sub-tropical region, which has a with temperate climate. This study focusses on torrential rainfall, which happens occasionally during the autumn season in the southeastern Australian Victoria region, and more specifically on the meteorological situation of the March 2010 Melbourne storm.
Heavy rainfall affected the state of Victoria from 6 to 8 March, 2010, and caused consequent flash floods that disrupted the life of the population of the greater Melbourne region, with several injuries and damages resulting [27]. This event was probably the most severe of the past decade for this region [28]. Such an event, acting as a mesoscale convective system (MCS) with a horizontal extension of generally over 10 km, a long life-time (sometimes of a few hours), and a strong vertical vorticity, largely caused by wind shear, is called a supercell storm. The mean monthly precipitation of Melbourne (50 mm) fell in only one day (61 mm on 6 March, 2010), including hailstones with a size of up to 5 cm, as reported by the Australian Bureau of Meteorology [29].
Considering estimations of cloud top altitudes from Global Ozone Monitoring Experiment-2 (GOME-2) and Ozone Monitoring Instrument (OMI) sensors, a vertical extension which was at least higher than 10 km was induced close to Melbourne by a supercell storm and can be observed in Figure 1. Cloud top altitudes were evaluated using cloud top pressure retrievals [30,31] combined with the assumption of an atmosphere in a hydrostatic equilibrium. This application of hypsometric equations required inputs of ground pressure and mean temperature of the column of air for each footprint considered. These two parameters were obtained from the outputs of the analysis of the Australian Community Climate and Earth-System Simulator (ACCESS-A), with a focus on Australian area, and the operational Numerical Weather Prediction (NWP) of the Australian Bureau of Meteorology [32]. Figure 1a,b show the cloud top altitudes obtained by GOME-2 and OMI during the overpass over Melbourne at 23:56 UTC (on 5 March, 2010) and 04:07 UTC (on 6 March, 2010), respectively. In Figure  1a, the MCS (with a vertical extension of over 9 km) covers an area of 80 × 80 km² (the size of a pixel from the GOME-2 instrument is 80 × 40 km²), and the area with stratiform precipitation (with a vertical extension of over 6 km) is about 160 × 160 km². In Figure 1b, convective precipitation with a high vertical extension can be seen to concern an area of about 24 × 84 km² (a pixel from OMI is 24 × 13 km²) and stratiform precipitation in an area of 168 × 195 km². The background of Figure 1b is an image (visible channel) from the Moderate Resolution Imaging Spectroradiometer (MODIS). This imager, on board the Aqua satellite, simultaneously looks at the same area of the Earth as OMI/Aura. The dimension of the MCS detected by GOME-2 and OMI are respectively confirmed by heavy and moderate precipitation recorded by In Figure 1b, convective precipitation with a high vertical extension can be seen to concern an area of about 24 × 84 km 2 (a pixel from OMI is 24 × 13 km 2 ) and stratiform precipitation in an area of 168 × 195 km 2  by GOME-2 and OMI are respectively confirmed by heavy and moderate precipitation recorded by the weather radar of Melbourne; see Choy et al. [27] and Manning et al. [33] for more details regarding this meteorological situation.
A typical configuration of the atmosphere with instability in the troposphere (e.g., that triggered by the topography or due to low level convergence) associated with the opposition of dry and moist air was required to initiate deep convection during the 6-8 March superstorm of greater Melbourne [27]. At the end of the day on 5 March, 2010 (UTC time), the lifting of moist air coming from the south and southeast sides of Melbourne and warm dry air coming from westwards of Melbourne triggered the convection. The water vapour content of the troposphere is a key parameter which initiates and keeps running deep convection. Its monitoring is critical to evaluating the severity of severe weather events (nowcasting) and to forecasting meteorological situations. Figure 2 shows the 2D field of integrated water vapour (IWV) from ACCESS-A NWP during the supercell storm of Melbourne in March 2010. Even a horizontal resolution of ACCESS-A outputs of 12 km, the weather prediction used in this study, was not sufficient to properly forecast this storm and the creation of such a supercell. Strong differences with the IWV fields from GPS are observed in Figure 2 (see Section 3.1 for more details regarding estimates of water vapour content from GPS).
Remote Sens. 2019, 11, x FOR PEER REVIEW 4 of 36 the weather radar of Melbourne; see Choy et al. [27] and Manning et al. [33] for more details regarding this meteorological situation. A typical configuration of the atmosphere with instability in the troposphere (e.g., that triggered by the topography or due to low level convergence) associated with the opposition of dry and moist air was required to initiate deep convection during the 6-8 March superstorm of greater Melbourne [27]. At the end of the day on 5 March, 2010 (UTC time), the lifting of moist air coming from the south and southeast sides of Melbourne and warm dry air coming from westwards of Melbourne triggered the convection. The water vapour content of the troposphere is a key parameter which initiates and keeps running deep convection. Its monitoring is critical to evaluating the severity of severe weather events (nowcasting) and to forecasting meteorological situations. Figure 2 shows the 2D field of integrated water vapour (IWV) from ACCESS-A NWP during the supercell storm of Melbourne in March 2010. Even a horizontal resolution of ACCESS-A outputs of 12 km, the weather prediction used in this study, was not sufficient to properly forecast this storm and the creation of such a supercell. Strong differences with the IWV fields from GPS are observed in Figure 2 (see Section 3.1 for more details regarding estimates of water vapour content from GPS).  The flux of water vapour from the south of Melbourne is indicated by horizontal IWV gradients (Figure 2c), aligned along south north in a southeast-northwest direction in the central section of Victoria and north south in the west part of Victoria. These opposite gradients illustrate the critical condition that took place in the region surrounding Melbourne which maintained a convective process inside the supercell. Note that IWV gradients were obtained using a κ conversion factor [34] (see Equation (6)) applied on wet gradients (see Equation (4)). Monitoring by GPS techniques using a dense network can be essential for improving the capability to better understand and forecast such events. For this purpose, this study focusses on testing five different types of tomography software for estimating reliable tomography retrievals.

Data Inputs for GPS Tomography Using the Continuously Operating Reference Station (CORS) Network
During the Melbourne storm in March 2010, 69 stations of the Victorian GPS CORS network were recording signals emitted from 31 GPS satellites. The locations of these stations are presented in Figure 3a. The GPS data were processed by the Royal Melbourne Institute of Technology University using the double difference technique with Bernese 5.0 software [35]. The full details of the processing strategy are available in a paper by Rohm et al. [36]. The data used in this study were derived using the "shortest baselines network solution" with the following parameters: precise final International GNSS Service (IGS) orbits in an IGS05 reference frame with an L3/L5 strategy, minimum constraints applied on the translation parameters for Helmert transformation based on IGS fiducial stations, elevation dependent weighting, and a 5 • -cutoff angle. Zenith total delays of the neutral atmosphere (ZTDs) had relative constraints applied and were estimated in 30 min steps with gradients modelled by tilting every 6 h. The GMF mapping function [37] is employed for geodetic software analysis and retrievals of tropospheric parameters. Hence, ZTDs and horizontal delay gradients → G = (G NS , G EW ) were retrieved every 30 min (a piece-wise linear function was evaluated every 30 min) covering the period from 3 March, 2010 to 10 March, 2010.
The understanding of meteorological situations by using GPS tomography techniques can rely on retrieving the 3D fields of wet refractivity and water vapour density. Initial GPS tropospheric parameters are converted into slant wet delays (SWDs) and slant-integrated Water Vapour contents (SIWVs) to retrieve, respectively, these fields. Estimates of ground pressure and mean temperature (of the column) are required for each station of the GPS CORS network. These parameters were obtained using outputs from the ACCESS-A weather model. The ground pressures and temperatures of the four surrounding ACCESS-A grid points of a GPS station were converted to the altitude of the GPS site (using a hypsometric equation and linear vertical interpolation, respectively). Then, pressures and temperatures at the location of the GPS site were obtained using bi-linear horizontal interpolation. This strategy is best in the absence of a collocated meteorological station [38]. A time interpolation was also applied from 6 h of a weather model to a 30 min time resolution of GPS tropospheric parameters. The formula of Saastamoinen [39] was used to convert pressure data into zenith hydrostatic delay (ZHD). The zenith wet delay (ZWD) was obtained by simply subtracting ZHD from ZTD estimates (i.e., ZWD = ZTD − ZHD). Using the field of ground pressure from ACCESS-A, the hydrostatic part of the azimuthal contribution of GPS gradient to delays to slant delays was removed to allow only the wet contribution to remain [40,41].
The elevation (ie) and azimuth (α) dependency of the SWD is described by an isotropic (L wet sym ) and anisotropic (L wet az ) component in Equation (1).
Remote Sens. 2020, 12, 30 The isotropic SWD (L wet sym ) is obtained by mapping of the ZWD using a wet mapping function mf wet sym (in our case GMF from [37]) in the direction of the GPS satellite in a view above 5 • elevation (Equation (2)).
This equation uses a gradient mapping function mf wet az = 1/(sin ie tan ie + C) which depends on the satellite's elevation and on a constant (C = 0.0032) ( [42,43]), and the gradient wet components (G wet NS , G wet EW ) in the north-south and east-west direction, which are multiplied with the cosine or sine of the azimuth (α), respectively. Note that the one-way residual has not been considered in the L wet az and SWD retrievals here.
The gradient components (G NS , G EW ) retrieved by GPS technique describe the total asymmetric effect. This means there is no distinction between wet and hydrostatic gradients [1,43]. However, the wet gradient component can be expressed by the difference of the hydrostatic to the total component, i.e., To obtain the hydrostatic gradient components (G hydrostatic NS and G hydrostatic EW ), a characterisation of the surface pressure field around each GPS station is required. In that case, the hydrostatic gradient can be established by fitting a plane through the pressure measurements [40,41]. The spatial variations of the hydrostatic delay per unit of distance (km) in the north-south (Z hydrostatic NS ) and east-west (Z hydrostatic EW ) directions can be calculated from the pressure field for each GPS site. Surface pressure measurements around all GPS stations of the CORS network were not available during the supercell storm of 6-8 March 2010. For this reason, the pressure field of ACCESS-A NWP was considered. Assuming an exponential law in the hydrostatic refractivity and considering the scale height of the gradients in the hydrostatic delays as being set to H = 13 km, as suggested by [43], the spatial variations of the hydrostatic delay can be converted into hydrostatic gradients [1,44,45], i.e., The resulting slant wet delays SWD = N w ds are used as inputs for GPS tomography to retrieve 3D fields of wet refractivity (N w ). ds is the differential distance along the slant path travel of the GPS signal. The second type of slant measurements used is water vapour slant total column density, also called slant-integrated water vapour contents, which can be simply obtained from SWDs by multiplication with factor κ, as shown in Equation (6). Then, GPS tomography retrieves the 3D field of water vapour density (ρ wv ).
The expression of the κ factor [34] requires an estimation of the mean temperature of the column of air (T m ) above the GPS station. Hence, temperature profiles from ACCESS-A were interpolated in space and time. The obtained κ factors were applied in Equation (6)  atmosphere (ZTDs) had relative constraints applied and were estimated in 30 min steps with gradients modelled by tilting every 6 h. The GMF mapping function [37] is employed for geodetic software analysis and retrievals of tropospheric parameters. Hence, ZTDs and horizontal delay gradients ⃗ = ( , ) were retrieved every 30 min (a piece-wise linear function was evaluated every 30 min) covering the period from 3 March, 2010 to 10 March, 2010.

Profiles from Radiosonde
The most common and accurate technique used to measure the vertical profile of temperature, pressure, dew point temperature, and wind speed and direction is a sensor attached to an automatic radio-sounding balloon released either every 6, 12, or at a minimum 24 h. Usually, profiles reach up to 30 km and have a very high vertical resolution; a major drawback of such measurements is the horizontal movement of the balloon as it ascents vertically due to winds. The ascent time is around 1 h. In 2010 in Australia there were 32 radiosonde stations with only one located in Victoria (Melbourne airport-YMML) and another two within a short distance from the Victorian border (Wagga Wagga-YSWG and Mount Gambier-YMMG), as shown in Figure 3. Only the YMML radiosonde station, located at an altitude of 141 m and position of (37.67 • S, 144.83 • E), has been considered. This radiosonde (YMML-RS) data was available twice a day (at noon and midnight) for all days used in the study. The radiosonde technique is highly accurate for observations of the troposphere with a 1-2 hPa accuracy for pressure, 0.5 • C accuracy for temperature, and 5% accuracy for relative humidity [46]. Based on the radiosonde data, the following calculations were made to retrieve wet refractivity and water vapour content: (1) the formula used by Sargent [47] was applied to obtain the relative humidity from the dew point temperature [48]; (2) the formula used by Sonntag [49] was used to calculate the partial pressure of water vapour p w from RH and T; and (3) the atmospheric state parameters p w , p, and T were used to calculate wet refractivity N w and water density profiles ρ wv using equations proposed by Davis et al., [50], i.e., Remote Sens. 2020, 12, 30 where k 1 , k 2 , and k 3 denote empirically-derived refractivity coefficients, R d is a gas constant for dry air, and ρ is the atmospheric density.

Profiles from the Radio-Occultation Technique
The radio-occultation (RO) profiles used in this study are wet profiles which were acquired via the Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC) and re-processed in 2013 at the COSMIC Data Analysis Archive Center (CDAAC). The wet profile (wetPrf) is an interpolated product with a vertical resolution of 100 m which was obtained using a one-dimensional variational (1DVar) technique together with European Centre Medium Weather Forecast (ECMWF) low resolution analysis data, and it contains latitude and longitude of the perigee point, pressure, temperature (T), water vapour pressure (p w ), refractivity, and a mean sea level altitude of the perigee point. The atmospheric parameters were computed using a 1DVar approach where refractivity is weighted in such a way that the temperature is basically the same as the dry temperature in regions where the water vapour is insignificant [51]. Analyses of specific humidity estimated by GPS RO [52] have shown consistency of being within 0.1-0.3 g/kg in the median with Global Climate Observing System (GCOS) Reference Upper-Air Network (GRUAN) radiosonde data and with theoretical studies of accuracy [53,54]. The RO globally performs better than the radiosondes between 5 and 25 km of altitude when compared to the ECMWF global analysis [54].
For a comparison with tomography results, we converted the p w in ρ wv by using the gas equation for the four COSMIC RO reported in Table 1 and shown in Figure 3.

A Selection of Five Tomography Models
GPS tomography uses slant-integrated measurements (SWDs and SIWVs) inside a defined volume to retrieve N w or ρ wv estimates for each voxel N pqr , as illustrated in Figure 4. For highest consistency, the same GPS data set and the same 3D grid were used for each of the five tomography models considered in this study.
For each GPS station within the voxel model, the integrated slant observations SLANT GPS (in our case the SWDs or SIWVs) were converted into 3D fields of N w or ρ wv using Equation (9), i.e., where P is the position of a GPS station (latitude, longitude, and altitude), єis the elevation and α the azimuth angle of a GPS satellite, NB vx is the number of voxels crossed by each observation, and ∆s i is the ray length in each voxel i (x i ∈ [1, p], y i ∈ [1, q], z i ∈ [1, r]). Following this, axes x, y, and z are associated, respectively, with longitude, latitude, and altitude. For the tomographic grid, the horizontal size of the voxels in the inner grid is 0.5 •~5 0 km (thick black lines in Figure 3a), which is slightly denser than the mean baseline (~70 km) of the GPS network; see preconditions from Bender and Raabe [55]. Twelve voxels are considered longitudinally (p = 12) and six voxels latitudinally (q = 6). All around this inner grid, voxels of the outer grid are considered inside a band with a width of about 400 km (thin black lines of Figure 3a). The outer grid is used to compensate for water vapour contributions outside the tomographic model, and, therefore, keep stations in the solution which are located close to the edge of the inner grid. To deal with the topography under the tomographic grid, note that the lower layer of the grid is at the altitude above the sea level (altitude = 0 m a.s.l.). The voxels for which there is no GPS station inside or under are simply not retrieved by tomography models. Section 6.1 and Table 5 discuss the percentage of voxels retrieved in our calculations according to different types of calculations and settings.
Remote Sens. 2019, 11, x FOR PEER REVIEW 9 of 36 models. Section 6.1 and Table 5 discuss the percentage of voxels retrieved in our calculations according to different types of calculations and settings. Assuming the hypothesis of straight rays (the bending effect which can affect observations under 15° elevations ( [56]) has not been considered in our study), the inversion problem becomes linear. The formulation of the linear inverse problem reads where d represents the SWD/SIWV observations. This is connected to the model m of wet refractivity/water vapour density ( / ) (which we want to retrieve) by the geometrical matrix G and the observation errors defined by the covariance operator C . Positions of stations and satellites used to characterise the geometrical matrix are expressed in Cartesian coordinates. SWD/SIWV observations are associated with errors (∆SWD/∆SIWV) which are expressed using C . Positions and observations are the inputs of this inverse problem. The solution m does not exist because GPS tomography is an ill-posed problem. However, the alternative is to find m′, which minimises the L2norm using the minimisation factor ξ, which is expressed as Different standard techniques exist for solving such linear inverse problems, e.g., singular value decomposition (SVD), truncated singular value decomposition (TSVD), weighted and damped leastsquare adjustment (LS), Kalman-filter, simultaneous algebraic reconstruction technique (SART), conjugate gradient method, Tikhonov-regularisation, stacking, and double-difference-spline interpolation.
This study uses five models with different inversion techniques, i.e., SVD, TSVD, LS, and SART. Retrievals were obtained mainly in 3D and for one model in 2D (the Technical University of Ostrava (TUO) model), and these were proceeded with for and , except for the Technische Universität Wien (TUW) model ( only) and the University of Beira Interior (UBI) model ( only). Table 2 summarises the characterisation of the covariance operator of data and the a priori model, as well as a the quality check for tomographic retrievals. Assuming the hypothesis of straight rays (the bending effect which can affect observations under 15 • elevations ( [56]) has not been considered in our study), the inversion problem becomes linear. The formulation of the linear inverse problem reads d = Gm + C D (10) where d represents the SWD/SIWV observations. This is connected to the model m of wet refractivity/water vapour density (N w / ρ wv ) (which we want to retrieve) by the geometrical matrix G and the observation errors defined by the covariance operator C D . Positions of stations and satellites used to characterise the geometrical matrix are expressed in Cartesian coordinates. SWD/SIWV observations are associated with errors (∆SWD/∆SIWV) which are expressed using C D . Positions and observations are the inputs of this inverse problem. The solution m does not exist because GPS tomography is an ill-posed problem. However, the alternative is to find m , which minimises the L 2 -norm using the minimisation factor ξ, which is expressed as Different standard techniques exist for solving such linear inverse problems, e.g., singular value decomposition (SVD), truncated singular value decomposition (TSVD), weighted and damped least-square adjustment (LS), Kalman-filter, simultaneous algebraic reconstruction technique (SART), conjugate gradient method, Tikhonov-regularisation, stacking, and double-difference-spline interpolation.
This study uses five models with different inversion techniques, i.e., SVD, TSVD, LS, and SART. Retrievals were obtained mainly in 3D and for one model in 2D (the Technical University of Ostrava (TUO) model), and these were proceeded with for N w and ρ wv , except for the Technische Universität Wien (TUW) model (N w only) and the University of Beira Interior (UBI) model (ρ wv only). Table 2 summarises the characterisation of the covariance operator of data and the a priori model, as well as a the quality check for tomographic retrievals.
More information about the five tomography models can be found in the Appendices A-E. Note that for all the tomography models and all the tests performed in this study, the same a priori condition m 0 was considered. This condition was based on 6 h of NWP outputs from ACCESS-A interpolated to the tomography grid of Figure 3. The reader can take a look at Puri et al. [32] for more details regarding the observations considered in the assimilation chain of ACCESS-A, which is a four-dimensional observation variational assimilation method (4DVAR). In 2010, the radio-occultation observations were not assimilated in the ACCESS-A system.

Methodological Improvement in GNSS Tomography
The inverse problem treated by GNSS tomography was ill-posed due to a non-homogenous distribution of slant observations through the 3D grid (lack and redundancy), but it was also ill-conditioned due to the high number of parameters physically embedded. This explains why the stabilisation of the tomographic solution remains a challenging task, in spite of methodological improvements [25]. Two main limitations are currently identified in the methodology of the GNSS tomography technique: (1) the deficiency of geometric representativeness and (2) the problems of convergence related to a reliable solution in the inverse process. The present work looks at the limitation induced by a non-sufficient geometrical distribution of GPS observations (e.g., those due to a network of stations with a mean baseline >30 km) and investigates methods for improving it. This study also investigates the issue in the convergence process of the tomographic inversion, which is linked to geometrical representativeness, and the impact of the observation uncertainty on the quality of tomographic results.

Methodology for Improving of the Geometrical Representativeness
Two types of improvement of the geometrical distribution of tomography retrievals are suggested and tested, namely, stacking of slant GPS data and additional pseudo-slant observations.

Data Stacking of Slant Observations
The time resolution of the tropopsheric parameters (ZTD and gradients) used to assess the slant observations triggered the increment chosen for the stacking periods, i.e., 30 min. The resolution of our tomography grid (50 km), which was related to the resolution of the GPS network considered, was a limitation to improving the geometrical distribution. In our case, there was no interest in choosing a shorter stacking period (e.g., 15 min). The observation stacking procedure decreases, in theory, the condition number of the design matrix (see Equation (A3) and Equation (A4) in Appendix A) by improving the observation geometry; however, it can also decrease the accuracy of retrieval if the troposphere is active and changes rapidly within the stacking window. Stacking of GPS data (SWD/SIWV) was tested for three stacking periods (30, 60, and 120 min before and after the time of tomography reconstruction).

Use of Pseudo-Slant Observations
This section presents a methodology which is based on adding pseudo-slant observations to the systems of equations to be solved. The pseudo-slants were implemented according to the orientation of the delay gradient [57]. To obtain homogeneous repartition of sites (based on the CORS network), a distance of 20 km on both sides of a station in the direction defined by the GPS gradient was considered. This means that for two pseudo-sites, the wet delay and the water vapour content were propagated at a 20 km distance using the amplitude of the gradient, i.e., the additional wet delay equal to the multiplication of the gradient by the distance. The pseudo-SLANT GPS in the direction of the GPS satellites was estimated using an isotropic mapping function (GMF in this case; see [37]) and considered in tomography calculations, as illustrated in Figure 5. Note also that two other pseudo-stations located at 10 km from the station were used, considering gradient orientation and its opposite direction. This makes for a total of four pseudo-sites with four sets of additional pseudo-slants (only two additional sets are illustrated in Figure 5). The time resolution of the tropopsheric parameters (ZTD and gradients) used to assess the slant observations triggered the increment chosen for the stacking periods, i.e., 30 min. The resolution of our tomography grid (50 km), which was related to the resolution of the GPS network considered, was a limitation to improving the geometrical distribution. In our case, there was no interest in choosing a shorter stacking period (e.g., 15 min). The observation stacking procedure decreases, in theory, the condition number of the design matrix (see Equation (16) and Equation (17) in Appendix A) by improving the observation geometry; however, it can also decrease the accuracy of retrieval if the troposphere is active and changes rapidly within the stacking window. Stacking of GPS data (SWD/SIWV) was tested for three stacking periods (30, 60, and 120 min before and after the time of tomography reconstruction).

Use of Pseudo-Slant Observations
This section presents a methodology which is based on adding pseudo-slant observations to the systems of equations to be solved. The pseudo-slants were implemented according to the orientation of the delay gradient [57]. To obtain homogeneous repartition of sites (based on the CORS network), a distance of 20 km on both sides of a station in the direction defined by the GPS gradient was considered. This means that for two pseudo-sites, the wet delay and the water vapour content were propagated at a 20 km distance using the amplitude of the gradient, i.e., the additional wet delay equal to the multiplication of the gradient by the distance. The pseudo-SLANTGPS in the direction of the GPS satellites was estimated using an isotropic mapping function (GMF in this case; see [37]) and considered in tomography calculations, as illustrated in Figure 5. Note also that two other pseudostations located at 10 km from the station were used, considering gradient orientation and its opposite direction. This makes for a total of four pseudo-sites with four sets of additional pseudo-slants (only two additional sets are illustrated in Figure 5). The blue and dashed-blue arrows represent, respectively, the direction of the horizontal gradients and their opposites (e.g., for a defined distance of, e.g., 20 km), which were used to generate pseudo-slants (blue lines).

Figure 5.
Illustration of additional pseudo-slant observations in GPS tomography. The dashed black lines show slant observations. The blue and dashed-blue arrows represent, respectively, the direction of the horizontal gradients and their opposites (e.g., for a defined distance of, e.g., 20 km), which were used to generate pseudo-slants (blue lines).

A Priori Condition and Improvement of the Convergence in the Inversion Process
This section investigates the ability to obtain an optimal stabilisation of solutions for five types of tomography software based on different techniques (singular value decomposition, damped least-squares adjustment, algebraic reconstruction technique, and Kalman filter). Four tests related to the a priori conditions used in calculations were studied to determine the best strategy.
These four configurations were designed to test convergence in tomography calculations according to the a priori condition and time resolution of the tomography models (capital letters when the interval calculation was 30 min, i.e., tests A or B; lowercase letters when it was 6 h, i.e., tests a or b). An illustration of the configuration of these four tests is presented in Figure 6, in which arrows in red indicate when the a priori condition from ACCESS-A was considered in tomography calculations, and arrows in green indicate when this was the previous tomography retrieval (TR).

A Priori Condition and Improvement of the Convergence in the Inversion Process
This section investigates the ability to obtain an optimal stabilisation of solutions for five types of tomography software based on different techniques (singular value decomposition, damped leastsquares adjustment, algebraic reconstruction technique, and Kalman filter). Four tests related to the a priori conditions used in calculations were studied to determine the best strategy. These four configurations were designed to test convergence in tomography calculations according to the a priori condition and time resolution of the tomography models (capital letters when the interval calculation was 30 min, i.e., tests A or B; lowercase letters when it was 6 h, i.e., tests a or b). An illustration of the configuration of these four tests is presented in Figure 6, in which arrows in red indicate when the a priori condition from ACCESS-A was considered in tomography calculations, and arrows in green indicate when this was the previous tomography retrieval (TR).

Sensitivity Tests Based on the Uncertainty of Slant Observations
The quality of slant observations is critical for obtaining optimal GNSS tomography retrievals. The strategy of this section was to modify slant observations using minimal/maximal uncertainties and to look at the impact on the quality of retrievals. The uncertainties of parameters acting for obtaining SWDs and SIWVs are presented in Table 3. These parameters are the specific molar gas constant for dry air (Rd), the specific molar gas constant for water vapour (Rw), the refractivity coefficients (k1, k2, k3, and k2' = k2 -k1.Rd/Rw), the surface pressure (PS), the mean temperature of the column of air (Tm), the gravity in the centre of the atmospheric column pressure (gm), the GPS observation of ZTD, and estimations of ZHD, ZWD, factor κ, and IWV. More details regarding the links of these physical parameters with ZWD and IWV (based on preliminary results from [34,50,58]) and the impact of their uncertainties on simulations (using a high-resolution non-hydrostatic model during severe weather situations) are presented in Brenot et al. [40]. Finally, two kinds of empirical mapping functions (NMF of Niell [59] and GMF of Böhm et al. [37]) have been considered to estimate the impact on SLANTGPS and retrievals from GPS tomography. The impact on the precision of slant retrievals is presented by an absolute uncertainty (in the unit of the tested parameter) and by a

Sensitivity Tests Based on the Uncertainty of Slant Observations
The quality of slant observations is critical for obtaining optimal GNSS tomography retrievals. The strategy of this section was to modify slant observations using minimal/maximal uncertainties and to look at the impact on the quality of retrievals. The uncertainties of parameters acting for obtaining SWDs and SIWVs are presented in Table 3. These parameters are the specific molar gas constant for dry air (R d ), the specific molar gas constant for water vapour (R w ), the refractivity coefficients (k 1 , k 2 , k 3 , and k 2 ' = k 2 − k 1 .R d /R w ), the surface pressure (P S ), the mean temperature of the column of air (T m ), the gravity in the centre of the atmospheric column pressure (g m ), the GPS observation of ZTD, and estimations of ZHD, ZWD, factor κ, and IWV. More details regarding the links of these physical parameters with ZWD and IWV (based on preliminary results from [34,50,58]) and the impact of their uncertainties on simulations (using a high-resolution non-hydrostatic model during severe weather situations) are presented in Brenot et al. [40]. Finally, two kinds of empirical mapping functions (NMF of Niell [59] and GMF of Böhm et al. [37]) have been considered to estimate the impact on SLANT GPS and retrievals from GPS tomography. The impact on the precision of slant retrievals is presented by an absolute uncertainty (in the unit of the tested parameter) and by a relative error (in %). In the upper six lines of Table 3 it can be seen that the parameter with the highest relative error is the refractivity coefficient k' 2 . The uncertainties of refractivity coefficients (k 1 , k 2 , k 3 ) provided by Bevis et al. [60] were used, and the combination of higher and lower values for these coefficients show an absolute error of about 10% for k' 2 . Table 3. Absolute uncertainty and relative error of parameters related to GPS tomography. Legend: R d, specific molar gas constant for dry air; R w , specific molar gas constant for water vapour; k 1 , k 2 , k 3 , k' 2 , refractivity coefficients; P S , surface pressure; T m , mean temperature of the atmospheric column; g m , gravity in the centre of the atmospheric column; ZTD, zenith total delays; ZHD, zenith hydrostatic delay; ZWD, zenith wet delay; κ, conversion factor; IWV, integrated water vapour content; SWD, slant wet delay; SIWV, slant-IWV; L wet sym , isotropic wet component of SWD; L wet az , anisotropic wet component of SWD; G EW , East-West component of delay gradients; G NS , North-South component of delay gradients. The bottom part of Table 3 presents absolute uncertainty and relative error for typical values of the estimated parameters. In these columns, two types of uncertainty were considered: one which was more conservative (the higher values) and one which was more realistic (the lower values). The results related to the more realistic absolute uncertainty and relative error are presented in brackets (e.g., a conservative uncertainty of P S of 2 hPa against a more realistic uncertainty of 1 hPa). To estimate the uncertainty of g m , the bias between the formulation of Saastamoinen [39] and the formulation of Vedel et al. [61] was considered (a relative error of about 0.2%). The relative error in ZTD measurements and in P S , T m , and ZHD estimates are low (<0.5%). On the other hand, the relative errors of ZWD and IWV are high, being about 8.5% for both for the conservative error and 4.5% and 6.2%, respectively, for the more realistic error. The analytic expression of the conversion factor κ takes into account k 3 , k' 2 , and R d coefficients and T m parameters [34]. T m was estimated using output from the ACCESS-A model with an absolute uncertainty of 1 K (or 0.5 K). A moderate relative error was found for κ (<1%). Note that this parameter was able to obtain a higher relative error in the case of severe weather conditions and the occurrence of hydrometeors in the path travel of the GNSS signal [40]. The last parameter that considered was the uncertainty of the mapping function. The uncertainty of SLANT GPS was estimated using a satellite visible from one station during one epoch of measurements (10 satellites).

Parameter
To test the uncertainty of SWDs, the uncertainty of L wet sym of Equation (2) was considered using the absolute error of the NMF and GMF mapping function [37,59], as shown in Equation (12) L wet,MIN where GMF = mf wet,GMF sym (ie) and NMF = mf wet,NMF sym (ie). This shows a relative error of 9%. To obtain the uncertainty of L wet az of Equation (3), the formal error of gradients provided by geodetic software and the error estimated by [17] were considered. This brought about a relative error of 10% (or 5%, to be more realistic). Finally, the uncertainties of SWD and SIWV were found to be, respectively, 9% and 10%, with more realistic errors of 5% occurring for both.

Results of Methodological Improvement in GPS Tomography
Three GPS tomography methodologies were investigated and combined. The improvement of geometrical representativeness is linked to convergence in tomography calculations. Independently, a sensitivity test of the perturbation of the quality of slant observations, inputs of GNSS tomography, has been presented and was combined with the two others improvements to investigate the interest in this approach in meteorological applications. The first methodological improvement targets a better geometrical distribution. The second concerns a better convergence behaviour in tomography solutions with respect to the a priori condition applied and the time resolution of calculations (i.e., tests a, b, B, A in Table 4, as illustrated in Figure 6). The configurations of tests A and B were used as references and the improvement of the geometrical distribution was investigated using stacking data (e.g., tests A 30 , A 60 , and A 120 ) and/or pseudo-slant observations (e.g., tests A π and A π 30 ). In addition, via a third way of testing and again using the configurations of tests A and B as references, the impact of the observational uncertainty on the quality of tomographic results was studied (e.g., tests A + and A − ). Table 4 shows an overview of all tested model settings.
For all the tests presented in Table 4, spatial and temporal distributions of the N w or ρ wv parameters were determined for three types of data ("initial observations (IO)", "IO with additional pseudo-observations", and "observations modified considering uncertainties"). The interval calculation for all the tests in Table 4 was found to be 30 min (except for tests a and b). Using the configuration of test A or B, other tests with 30, 60, and 120 in the superscript (i.e., A 30 , A 60 , and A 120 ) were meant to test improvement of stacked data on solutions using 30, 60, and 120 min -stacked data, respectively. Tests with π in the subscript (i.e., tests A π and B π ) were set up to test the improvement of using pseudo-slants solutions and were combined with the test of stacked data (i.e., A π 30 , A π 60 , and A π 120 ). The impact of the uncertainty on solutions (columns "positive" and "negative" of Table 4) was investigated for two types of uncertainty: regular, with labels "+" and "-", and more realistic, with labels "(+)" and "(−)". The "positive" column indicates uncertainty increasing slant observations (superscript + or (+) used in the name of the test, e.g., test A + or A (+) ) and the "negative" column indicates uncertainty decreasing slant observations (superscript −or (−) , e.g., test A − or A (−) ). The impact of uncertainty on solutions was combined with the test of stacked data (i.e., A 30(+) , A 60(+) , and A 120(+) ). To visualise the results and validate our methodological improvements, tomography retrievals were compared with external observations (profiles of N w and ρ wv from RS and RO). Profiles of N w and ρ wv from ACCESS-A were also considered (as an example of forecasts) and were the a priori conditions used in our tomography calculations. ERA-Interim [62] profiles were also prepared and compared.
Three various comparisons are shown within this study: GPS tomography versus RS profiles, GPS tomography versus RO profiles, and individual GPS tomography solutions versus each other or ACCESS-A and ERA-Interim fields. Since all mentioned techniques sense tropospheric water vapour in a different way, the presented comparisons were based on different approaches. To compare RS profiles with GPS tomography results, or with ACCESS-A and ERA-Interim fields, RS profiles were linearly interpolated to individual heights of the tomography grid. Tomography retrievals and ACCESS-A/ERA-Interim estimates were then bilinearly interpolated to positions of RS stations; RS profiles were assumed to be vertical (wind field was not considered, given the assumption that tomography outputs and RS profiles are simultaneous). Slant RO profiles were derived at heights of the tomography grid. Then, tomography retrievals and ACCESS-A/ERA-Interim estimates were bilinearly interpolated to latitudes/longitudes of RO profiles. For comparisons between two individual GPS tomography solutions or between a GPS tomography solution and ACCESS-A/ERA-Interim fields, no interpolation had to be applied and a direct confrontation in 3D tomographic voxels was realised (for ACCESS-A/ERA-Interim the closest grid points to the tomography grid were considered).

Results Regarding the Improvement of Geometrical Distribution
An increased geometrical distribution can be a way to obtain a better stabilisation of tomography solutions. More voxels being crossed by slant GPS observations (SLANT GPS ) can improve the understanding of this meteorological situation. However, this means that more inputs are considered in tomography processing. Table 5 shows the status of the geometrical distribution (% of voxels crossed by SLANT GPS for the different tests considered) inside the inner tomography grid. The mean number of SLANT GPS is provided for each epoch of calculation. Note that the mean time of processing is also indicated (a 2.66 GHz processor with 24 GB RAM was used) as being a key parameter for nowcasting (non-linear quasi-quadratic dependency with the number of inputs). According to the distribution of rays from GPS stations inside the tomographic grid, for the 288 times of measurements from 3 to 10 March, 2010, the mean percentage of voxels retrieved by our classic tomography models was 67.8% for the whole tomographic grid (geometrical distribution obtained without stacking and pseudo-slant observations). Using stacked slant data (30, 60, or 120 min stacking), an increase in the geometrical distribution (+1.5%, +1.7% or +4.2%, respectively) and a consequent increase in the time of processing (×2, ×4, or ×16, respectively) was observed. Note that in these cases of stacked data, slants were only considered every 30 min, not continuously (this was not relevant for our tomography grid, which had a horizontal resolution of 50 km). The layers between 2 and 9 km show the main results with a high increase of geometrical distribution (about +10%). Using pseudo-slant observations, an increase of +25% can be observed for all layers on average, with an increase of +30% for the layers between 0 and 4 km and an increase of +20% for those between 4 and 9 km. The geometrical distribution only increased by +3% for the top layer of between 9 and 15 km (for more details regarding pseudo-slants, see Section 5.1.2). Three sets of stacked slant data were tested, namely, 30, 60, and 120 min stacking. Hence, the two processing strategies illustrated in Figure 6 by tests A and B were applied for tests (A 30 , A 60 , A 120 ) and (B 30 , B 60 , B 120 ), respectively.
In Figure 7a, the results for tests (A 30 , A 60 , A 120 ) are compared to test A, as this solution was also initiated with a priori data every 6 h. The three sets of solutions for the three models show different patterns. We can see that the 8-13 km layer was not well solved by the tomography models, which consequently impacted the overall result (respectively, lines/layers 4 and 7). Looking at the layers under 10 km, the Wroclaw University of Environmental and Life Sciences (WUELS) and University of Beira Interior (UBI) models show a similar performance for tests A and A 30 (shortest 30 min stacking) and an improved accuracy for the longer stacking period. WUELS reached a better performance in test A 60 , i.e., for 60 min stacking, with a normalised root mean square (RMS) of 0.52 g/m 3 compared with 0.96 g/m 3 for test A. Both UBI and the Royal Belgian Institute for Space Aeronomy (BIRA) model obtained comparable accuracy retrievals in the first two layers to ACCESS-A values and retrievals which were slightly better than ERA-Interim (normalised RMS of 0.16 g/m 3 versus 0.21 g/m 3 ). Note that normalised RMS values were computed from relative differences of two corresponding values from two solutions. The relative difference means that an absolute difference of two values is divided by one of the values. This method of statistics computation allows a reasonable comparison of studied N w and ρ wv parameters, whose absolute values change significantly with height.
The results shown in Figure 7b for tests (B 30 , B 60 , B 120 ) were compared to test B, as this solution was also initiated with a priori data from NWP (i.e., ACCESS-A) only in the first epoch, the so-called "stand-alone" mode. The WUELS and UBI models in tests (B 30 , B 60 , B 120 ) performed very similarly to test B as long as the troposphere below 4 km is considered; in the 4-8 km layer, both delivered improved solutions: WUELS showed in test B 60 a 30% improvement over test B and UBI performed in test B 30 about 20% better than in test B. However, both solutions were worse than ACCESS-A and ERA-Interim, especially above 8 km. The BIRA model produced somewhat different results: in the 0-1.5 km layer the solution for test B 30 was better than that in test B and in ERA-Interim (with a normalised RMS of 0.19 g/m 3 versus 0.20 g/m 3 and 0.21 g/m 3 , respectively). In all higher layers until 13 km test B 120 was much better than test B, while for the top-most layer it attained a better performance than ACCESS-A and ERA-Interim (0.90 g/m 3 versus 0.98 g/m 3 and 1.06 g/m 3 , respectively). Whatever the tomography model used, the performance in the 0-8 km layer was better in test (B 30 , B 60 , B 120 ) than in test B; however, above an 8 km height, the best performance was found (in the BIRA and WUELS models) for the longest stacking period of 120 min (test B 120 ). Note that for BIRA, test B 120 showed even better results than ACCESS-A (1.06 g/m 3 versus 0.90 g/m 3 ).
Remote Sens. 2019, 11, x FOR PEER REVIEW 17 of 36 respectively). Whatever the tomography model used, the performance in the 0-8 km layer was better in test (B 30 , B 60 , B 120 ) than in test B; however, above an 8 km height, the best performance was found (in the BIRA and WUELS models) for the longest stacking period of 120 min (test B 120 ). Note that for BIRA, test B 120 showed even better results than ACCESS-A (1.06 g/m³ versus 0.90 g/m³).

Interest in Pseudo-Slant Observations for Low-and Mid-Troposphere Retrievals
This section investigates the interest in adding pseudo-slants observations (SLANT GPS ) in tomography retrievals. As shown in Table 3, the use of pseudo-slants increased the geometrical distribution by +30% and +20% for the 0-4 km and 4-9 km layers, respectively. To test the implication of improving the geometry with pseudo-slant observations, test A (with a regular initialisation from NWP outputs; see Figure 6) and test B (only previous tomography retrievals used as a priori, making a Remote Sens. 2020, 12, 30 18 of 35 'stand-alone solution') were compared with tests A π and B π , respectively (tests that considered pseudo-observations); see Figure 8.
Tests A π 30 and B π 30 (with 30 min data stacking) are also shown. Tests with 60 and 120 min stacked pseudo-slants, as originally planned, were not processed because they are not relevant for nowcasting (they have time windows which are too long and involve a significant increase of time processing of about ×250 and ×4000, respectively). The comportments of the three models considered were different. While WUELS improved in tests A π and A π 30 over test A (an improvement of about 15%), a degradation was able to be observed in tests B π and B π 30 over test B (the normalised RMS was multiplied by 2.4 g/m 3 and 2.2 g/m 3 , respectively). This conclusion was mainly driven by analysis of the 8-13 km layer. Even with an improvement, the results were still far from the RS data, which had a normalised RMS of about 73 g/m 3 for WUELS and about 1 g/m 3 for ACCESS-A and ERA-Interim. For the UBI model, the same tropospheric layer drove the result, but with a degradation according to both tests A and B (the normalised RMS was multiplied by 6 and 12, respectively). For the BIRA model, tomography retrievals using pseudo-SLANT GPS with a regular initialisation from NWP showed a small degradation from a normalised RMS of 0.43 g/m 3 to 0.87 g/m 3 and 1.27 g/m 3 for tests A π and A π 30 , respectively (this was mainly driven by the high troposphere layer). However, the use of pseudo-SLANT GPS with 'stand-alone initialisation' (giving an overall normalised RMS of 2.84 g/m 3 for test B) obtained a consequent improvement with an overall RMS of 0.56 g/m 3 and 1.05 g/m 3 for tests B π and B π 30 , respectively. The improvement took place in five of the seven layers (altitude ranges) studied; no improvement was observed in the 0-1.5 km and 1.5-4 km layers. In the 8-13 km later, the BIRA model even showed a slightly better performance than ACCESS-A (which had a normalised RMS of 0.97 g/m 3 for BIRA compared with 0.98 g/m 3 for ACCESS-A). The recommendation of improving the geometrical matrix of tomography retrievals with pseudo-slant observations was applied using data from the Belgian dense network (baselines from 5 to 30 km) during the pouring rains of 15-17 August, 2010 (a weather depression called Yvette by German meteorologists). GAMIT software was used to calculate ZTDs and gradients with a respective time resolution of 15 and 30 min (see Brenot et al. [40] for more details regarding the settings of these measurements). SIWVs were obtained using the GMF mapping function and ground pressures and temperatures from synoptic stations (no post-fit residuals were considered). The horizontal resolution of the tomography grid was 10 × 10 km and a vertical resolution of 500 m (for up to 10 km) was set. Figure 9 presents a comparison of profiles from tomography retrievals and ALARO NWP outputs with RS data from the Uccle station (Brussels) at 12:00 UTC on 16 August, 2010. The a priori condition from the ALARO NWP model (4 km resolution) was considered and applied to retrieve water vapour density from GPS tomography (employing the classical method and using pseudo-SIWVs).
The normalised RMS between RS and a priori ALARO, RS and tomography (classical), and RS and tomography (pseudo-slants) were 1.00 g/m 3 , 0.97 g/m 3 , and 0.89 g/m 3 , respectively, for all layers. These normalised RMS values became 0.47 g/m 3 , 0.44 g/m 3 , and 0.39 g/m 3 in the 0-6 km layer and 0.11 g/m 3 , 0.09 g/m 3 , and 0.02 g/m 3 in the 0-3 km layer. This confirms interest in using pseudo-slant observations in GPS tomography.
showed a small degradation from a normalised RMS of 0.43 g/m³ to 0.87 g/m³ and 1.27 g/m³ for tests Aπ and Aπ 30 , respectively (this was mainly driven by the high troposphere layer). However, the use of pseudo-SLANTGPS with 'stand-alone initialisation' (giving an overall normalised RMS of 2.84 g/m³ for test B) obtained a consequent improvement with an overall RMS of 0.56 g/m³ and 1.05 g/m³ for tests Bπ and Bπ 30 , respectively. The improvement took place in five of the seven layers (altitude ranges) studied; no improvement was observed in the 0-1.5 km and 1.5-4 km layers. In the 8-13 km later, the BIRA model even showed a slightly better performance than ACCESS-A (which had a normalised RMS of 0.97 g/m³ for BIRA compared with 0.98 g/m³ for ACCESS-A). The recommendation of improving the geometrical matrix of tomography retrievals with pseudo-slant observations was applied using data from the Belgian dense network (baselines from 5 to 30 km) during the pouring rains of 15-17 August, 2010 (a weather depression called Yvette by German meteorologists). GAMIT software was used to calculate ZTDs and gradients with a respective time resolution of 15 and 30 min (see Brenot et al. [40] for more details regarding the settings of these measurements). SIWVs were obtained using the GMF mapping function and ground pressures and temperatures from synoptic stations (no post-fit residuals were considered). The horizontal resolution of the tomography grid was 10 × 10 km and a vertical resolution of 500 m (for up to 10 km) was set. Figure 9 presents a comparison of profiles from tomography retrievals and ALARO NWP outputs with RS data from the Uccle station (Brussels) at 12:00 UTC on 16 August, 2010. The a priori condition from the ALARO NWP model (4 km resolution) was considered and applied to retrieve water vapour density from GPS tomography (employing the classical method and using pseudo-SIWVs).  The recommendation of improving the geometrical matrix of tomography retrievals with pseudo-slant observations was applied using data from the Belgian dense network (baselines from 5 to 30 km) during the pouring rains of 15-17 August, 2010 (a weather depression called Yvette by German meteorologists). GAMIT software was used to calculate ZTDs and gradients with a respective time resolution of 15 and 30 min (see Brenot et al. [40] for more details regarding the settings of these measurements). SIWVs were obtained using the GMF mapping function and ground pressures and temperatures from synoptic stations (no post-fit residuals were considered). The horizontal resolution of the tomography grid was 10 × 10 km and a vertical resolution of 500 m (for up to 10 km) was set. Figure 9 presents a comparison of profiles from tomography retrievals and ALARO NWP outputs with RS data from the Uccle station (Brussels) at 12:00 UTC on 16 August, 2010. The a priori condition from the ALARO NWP model (4 km resolution) was considered and applied to retrieve water vapour density from GPS tomography (employing the classical method and using pseudo-SIWVs).

Convergence of Tomography Solutions with Respect to A Priori Conditions and Time Resolution of Processes
To look at the impact of a priori conditions on tomography retrievals (for 3-8 March 2010), this section focuses on tests A, b, and B, as presented in Figure 6. RS were launched at 00:00 and 12:00 UTC, so the results of tests A and a were the same. Figure 10 shows that in test A, wet refractivity estimates from GPS tomography models were usually the best in the bottom part of the troposphere (below 1.5 km), with two out of three models displaying a comparable performance (less than 10% degradation) to ACCESS-A, namely, BIRA and TUW. Tomography retrievals (test A for BIRA, WUELS, and TUW) in the bottom part of the model showed better agreement with RS than ERA-Interim results, and TUO displayed a similar performance (not shown in Figure 10). However, this shows that this positive impact of GPS observations on the estimated troposphere state was diminished once the a priori data introduced to the tomography model were no longer fed in every epoch (test b) or the time resolution of estimates became higher (i.e., 30 min for test B). In the second investigated layer (1.5-4 km) the tomography-retrieved refractivity displayed usually the same or slightly worse accuracy than the NWP-based solutions: BIRA and TUW displayed 0.02 ppm (test A). Once the a priori data were not fed into the solution, the quality of retrievals dropped substantially, by almost 40%. The retrievals in the higher layers, namely, 4-8 km and 8-13 km, for all models in this study, were substantially worse than the ACCESS-A retrieval in all three tests. It also should be noted that the TUO retrieval, even though less accurate than ACCESS-A, provided better solutions than ERA-Interim for all 10 km, producing a normalised RMS of 0. 36  The normalised RMS between RS and a priori ALARO, RS and tomography (classical), and RS and tomography (pseudo-slants) were 1.00 g/m 3 , 0.97 g/m 3 , and 0.89 g/m 3 , respectively, for all layers. These normalised RMS values became 0.47 g/m 3 , 0.44 g/m 3 , and 0.39 g/m 3 in the 0-6 km layer and 0.11 g/m 3 , 0.09 g/m 3 , and 0.02 g/m 3 in the 0-3 km layer. This confirms interest in using pseudo-slant observations in GPS tomography.

Convergence of Tomography Solutions with Respect to A Priori Conditions and Time Resolution of Processes
To look at the impact of a priori conditions on tomography retrievals (for 3-8 March 2010), this section focuses on tests A, b, and B, as presented in Figure 6. RS were launched at 00:00 and 12:00 UTC, so the results of tests A and a were the same. Figure 10 shows that in test A, wet refractivity estimates from GPS tomography models were usually the best in the bottom part of the troposphere (below 1.5 km), with two out of three models displaying a comparable performance (less than 10% degradation) to ACCESS-A, namely, BIRA and TUW. Tomography retrievals (test A for BIRA, WUELS, and TUW) in the bottom part of the model showed better agreement with RS than ERA-Interim results, and TUO displayed a similar performance (not shown in Figure 10). However, this shows that this positive impact of GPS observations on the estimated troposphere state was diminished once the a priori data introduced to the tomography model were no longer fed in every epoch (test b) or the time resolution of estimates became higher (i.e., 30 min for test B). In the second investigated layer (1.5-4 km) the tomographyretrieved refractivity displayed usually the same or slightly worse accuracy than the NWP-based solutions: BIRA and TUW displayed 0.02 ppm (test A). Once the a priori data were not fed into the solution, the quality of retrievals dropped substantially, by almost 40%. The retrievals in the higher layers, namely, 4-8 km and 8-13 km, for all models in this study, were substantially worse than the ACCESS-A retrieval in all three tests. It also should be noted that the TUO retrieval, even though less accurate than ACCESS-A, provided better solutions than ERA-Interim for all 10 km, producing a normalised RMS of 0.36 and 0.38.  Figure 7, with normalised RMS of the comparisons of wet refractivity ( ) profiles from RS with 11 datasets; two sets from NWP and nine sets from three tomography models (BIRA, WUELS, and TUW) for tests A, b, and B.
A similar pattern was obtained for models producing water vapour density (not presented). Compared to radiosonde in Melbourne, the tomography solution (test A) was of similar quality to the ACCESS-A model in the first two layers of the atmosphere (below 4 km), especially for the models BIRA and UBI, while WUELS and TUO showed a slightly lower performance. Between 4-8 km, only BIRA and TUO retrievals were of useful accuracy, while WUELS and UBI showed five and four times Figure 10. Same as in Figure 7, with normalised RMS of the comparisons of wet refractivity (N w ) profiles from RS with 11 datasets; two sets from NWP and nine sets from three tomography models (BIRA, WUELS, and TUW) for tests A, b, and B.
A similar pattern was obtained for models producing water vapour density (not presented). Compared to radiosonde in Melbourne, the tomography solution (test A) was of similar quality to the ACCESS-A model in the first two layers of the atmosphere (below 4 km), especially for the models BIRA and UBI, while WUELS and TUO showed a slightly lower performance. Between 4-8 km, only BIRA and TUO retrievals were of useful accuracy, while WUELS and UBI showed five and four times lower performance, respectively. Interestingly, the 8-13 km retrieval from TUO showed a better performance than ACCESS-A data (displaying a normalised RMS of 0.94 g/m 3 and 0.98 g/m 3 ). Test b and B showed a much worse solution, especially in the higher levels of the model. The BIRA model in test A reached a better agreement with RS than ERA-Interim results across the whole 0-10 km profile (producing a normalised RMS of 0.28 g/m 3 and 0.38 g/m 3 ).

Results Regarding the Impact on the Precision of Slant Retrievals
To investigate the impact on the precision of slant retrievals, the uncertainties of parameters acting for obtaining SWDs and SIWVs were considered; see Section 5.3. The higher and lower estimates of slants (as shown in Table 3) were used to look at the repercussions on tomography solutions. Figure 11 highlights differences of normalised RMS (related to ρ wv observations) from 13 datasets (compared to RS, e.g., RS-BIRA [A 30+ ]) with RS-ACCESS-A. A highlight of eight tests regarding the impact of uncertainty is presented (tests A + and A − , for higher and lower slants, respectively, without stacking; tests A 30+ and A 30with 30 min stacking; tests A 60+ and A 60with 60 min stacking; and tests A 120+ and A 120with 120 min stacking). Superscripted labels + and -and (+) and (−) refer to regular and more realistic uncertainties (see values in brackets in Table 3), respectively, shown in Figure 11a,b. These eight tests were used a priori from ACCESS-A every 6 h and can be compared to tests A, A 30 , A 60 , and A 120 . The difference in normalised RMS from RS-ERA-Interim with RS-ACCESS-A is also presented. To avoid differences in normalised RMS which were less or equal to 0.01 g/m 3 in the logarithm scale of the y-axis, an offset of +0.03 g/m 3 was considered.
Remote Sens. 2019, 11, x FOR PEER REVIEW 21 of 36 lower performance, respectively. Interestingly, the 8-13 km retrieval from TUO showed a better performance than ACCESS-A data (displaying a normalised RMS of 0.94 g/m³ and 0.98 g/m³). Test b and B showed a much worse solution, especially in the higher levels of the model. The BIRA model in test A reached a better agreement with RS than ERA-Interim results across the whole 0-10 km profile (producing a normalised RMS of 0.28 g/m³ and 0.38 g/m³).

Results Regarding the Impact on the Precision of Slant Retrievals
To investigate the impact on the precision of slant retrievals, the uncertainties of parameters acting for obtaining SWDs and SIWVs were considered; see Section 5.3. The higher and lower estimates of slants (as shown in Table 3) were used to look at the repercussions on tomography solutions. Figure 11 highlights differences of normalised RMS (related to observations) from 13 datasets (compared to RS, e.g., RS-BIRA [A 30+ ]) with RS-ACCESS-A. A highlight of eight tests regarding the impact of uncertainty is presented (tests A + and A -, for higher and lower slants, respectively, without stacking; tests A 30+ and A 30-with 30 min stacking; tests A 60+ and A 60-with 60 min stacking; and tests A 120+ and A 120-with 120 min stacking). Superscripted labels + and -and (+) and (-) refer to regular and more realistic uncertainties (see values in brackets in Table 3), respectively, shown in Figure 11a,b. These eight tests were used a priori from ACCESS-A every 6 h and can be compared to tests A, A 30 , A 60 , and A 120 . The difference in normalised RMS from RS-ERA-Interim with RS-ACCESS-A is also presented. To avoid differences in normalised RMS which were less or equal to 0.01 g/m³ in the logarithm scale of the y-axis, an offset of +0.03 g/m³ was considered. Figure 11. (left) (a) Difference in normalised RMS from 13 datasets of water vapour density ( ) profiles (one set from NWP, i.e., RS-ERA-Interim, and 12 sets from tomography, i.e., RS-BIRA A, A + , A -, A 30 , A 30+ , A 30-, A 60 , A 60+ , A 60-, A 120 , A 120+ , and A 120-); to avoid values of less than 0.01 g/m³ an offset of +0.03 g/m³ was used. (b) Same as (a) but using more realistic uncertainties. For each dataset, differences over seven altitude ranges have been presented. (right) Illustration of the altitude ranges covered by the seven layers.
We found consistent comportment of tomography retrievals according to the type of uncertainty (regular with strong impact and more realistic with weak impact). We also observed that stacked data amplified the impact of uncertainty in tomography retrievals. Considering RS-ACCESS-A as a Figure 11. (left) (a) Difference in normalised RMS from 13 datasets of water vapour density (ρ wv ) profiles (one set from NWP, i.e., RS-ERA-Interim, and 12 sets from tomography, i.e., RS-BIRA A, A + , A − , A 30 , A 30+ , A 30-, A 60 , A 60+ , A 60-, A 120 , A 120+ , and A 120-); to avoid values of less than 0.01 g/m 3 an offset of +0.03 g/m 3 was used. (b) Same as (a) but using more realistic uncertainties. For each dataset, differences over seven altitude ranges have been presented. (right) Illustration of the altitude ranges covered by the seven layers.
We found consistent comportment of tomography retrievals according to the type of uncertainty (regular with strong impact and more realistic with weak impact). We also observed that stacked data amplified the impact of uncertainty in tomography retrievals. Considering RS-ACCESS-A as a reference, the sensitivity of the bias (i.e., the difference of normalised RMS) between RS and tomography water vapour density for different kinds of slants (maximum/minimum uncertainty with regular/more realistic errors), in comparison to the RS-ERA-Interim, showed interesting results regarding the quality and possible improvements in SLANT GPS as well as the validation of the NWP model using tomography retrievals. Adding an uncertainty measure to the observation matrix increased the accuracy of stochastic modelling for tomography retrievals and using a priori every 6h (in Figure 11) and stacked solutions (e.g., A 30(−) , A 60(−) , and A 120(−) ) led to a better performance in the bottom 1.5 km.
In the upper layer (1.5-4 km), these solutions were of lower quality than in ACCESS-A, but still better than those of ERA-Interim. The 4-8 km layer was again better resolved using the BIRA model in tests A 30(−) and A 60(−) . The topmost layer was of lower quality in all tests. The same kinds of tests were performed using a priori from ACCESS-A only for the first epoch (to be compared with tests B, B 30 , B 60 , and B 120 ). In comparison to the results of Figure 11, the impact of the uncertainty was amplified due to the stand-alone situation of the tomography technique which was less constrained by the a priori condition and ACCESS-A model. To illustrate this, we looked at the 4-8 km layer and the comparison of ρ wv from RS with estimates from ACCESS-A, ERA-Interim, test A 30 , test A 30+ , and test A 30-. We found respective normalised RMS values of 0.50 g/m 3 , 0.67 g/m 3 , 0.61 g/m 3 , 0.72 g/m 3 (0.67 g/m 3 for test A 30(+) ), and 0.48 g/m 3 (0.55 g/m 3 for test A 30(−) ), and also looked at a comparison with test B 30 , test B 30+ , and test B 30-, which showed respective normalised RMS values of 3.49 g/m 3 , 3.50 g/m 3 (3.17 g/m 3 for test B 30(+) ), and 2.11 g/m 3 (2.49 g/m 3 for test B 30(−) ). Figure 12 illustrates the comportments of ρ wv retrievals from three tomography models (BIRA, UBI, and WUELS) with respect to Melbourne RS estimates during the whole five-day period; comportments of NWP (ACCESS-A or ERA-Interim) are also presented. Even ACCESS-A and ERA-Interim show ρ wv estimates closer to RS measurements than the tomography models, and it is possible to notice that the three tomography models show the same ρ wv pattern in the 0-4 km layer. A trend of moistening was observed in the low troposphere by GPS tomography and was not seen by RS and ACCESS-A NWP. Such features could be further studied by ensemble tomographic processing, which could be a way to identify anomalies in the low troposphere with respect to NWP or other measurement techniques. Note that tomography outputs and RS were assumed to be simultaneous (without considering the wind field that affected the trajectory of RS profiles). These assumptions could be the reason for the difference between tomography models and RS observations. Remote Sens. 2019, 11, x FOR PEER REVIEW 22 of 36 reference, the sensitivity of the bias (i.e., the difference of normalised RMS) between RS and tomography water vapour density for different kinds of slants (maximum/minimum uncertainty with regular/more realistic errors), in comparison to the RS-ERA-Interim, showed interesting results regarding the quality and possible improvements in SLANTGPS as well as the validation of the NWP model using tomography retrievals. Adding an uncertainty measure to the observation matrix increased the accuracy of stochastic modelling for tomography retrievals and using a priori every 6h (in Figure 11) and stacked solutions (e.g., A 30(-) , A 60(-) , and A 120(-) ) led to a better performance in the bottom 1.5 km. In the upper layer (1.5-4 km), these solutions were of lower quality than in ACCESS-A, but still better than those of ERA-Interim. The 4-8 km layer was again better resolved using the BIRA model in tests A 30(-) and A 60(-) . The topmost layer was of lower quality in all tests. The same kinds of tests were performed using a priori from ACCESS-A only for the first epoch (to be compared with tests B, B 30 , B 60 , and B 120 ). In comparison to the results of Figure 11, the impact of the uncertainty was amplified due to the stand-alone situation of the tomography technique which was less constrained by the a priori condition and ACCESS-A model. To illustrate this, we looked at the 4-8 km layer and the comparison of from RS with estimates from ACCESS-A, ERA-Interim, test A 30 , test A 30+ , and test A 30-. We found respective normalised RMS values of 0.50 g/m³, 0.67 g/m³, 0.61 g/m³, 0.72 g/m³ (0.67 g/m³ for test A 30(+) ), and 0.48 g/m³ (0.55 g/m³ for test A 30(-) ), and also looked at a comparison with test B 30 , test B 30+ , and test B 30-, which showed respective normalised RMS values of 3.49 g/m³, 3.50 g/m³ (3.17 g/m³ for test B 30(+) ), and 2.11 g/m³ (2.49 g/m³ for test B 30(-) ). Figure 12 illustrates the comportments of ρ retrievals from three tomography models (BIRA, UBI, and WUELS) with respect to Melbourne RS estimates during the whole five-day period; comportments of NWP (ACCESS-A or ERA-Interim) are also presented. Even ACCESS-A and ERA-Interim show ρ estimates closer to RS measurements than the tomography models, and it is possible to notice that the three tomography models show the same ρ pattern in the 0-4 km layer. A trend of moistening was observed in the low troposphere by GPS tomography and was not seen by RS and ACCESS-A NWP. Such features could be further studied by ensemble tomographic processing, which could be a way to identify anomalies in the low troposphere with respect to NWP or other measurement techniques. Note that tomography outputs and RS were assumed to be simultaneous (without considering the wind field that affected the trajectory of RS profiles). These assumptions could be the reason for the difference between tomography models and RS observations.    Figure 1) showed no significant difference between RO (mean normalised bias from 10 to 20%), tomography retrievals, and ACCESS-A and ERA-Interim results, and followed a typical exponential pattern. Hence, these profiles will not be discussed further. Figure 13 presents profiles which run from the south-west towards the north-east over the Australian Alps (see Figure 3) on 4 March, 2010 at 16:39 UTC (08:39 a.m. local time). Over 6 km the radio-occultation (RO) profiles are outside the inner tomography grid. For this reason, only RO and Numerical Weather Prediction (NWP) models show estimates of wet refractivity (N w ) over 6 km in Figure 13. The tomography solution for N w is the one referred to in Table 4 as "no stacking", with a priori data fed into the model only in the first epoch. The results of the WUELS, BIRA, and TUW tomography models were similar, distinguishing one inversion layer at heights between 4 and 6 km (we can define this as a N w anomaly). For the RO profile, there were also local minima at these heights; however, the largest inversion occurred below 4 km. The RO-and tomography-based values of N w and ρ wv significantly differed from the ACCESS-A results in the lowest layer, as the latter did not show any inversion. The ERA-Interim result for this particular location was much closer to the tomography solution (BIRA and WUELS) and RO retrievals, showing a slight inversion between 5 and 6 km. The N w anomaly retrieved by both tomography models is a relevant piece of information for nowcasting.  Figure 1) showed no significant difference between RO (mean normalised bias from 10 to 20%), tomography retrievals, and ACCESS-A and ERA-Interim results, and followed a typical exponential pattern. Hence, these profiles will not be discussed further. Figure 13 presents profiles which run from the south-west towards the north-east over the Australian Alps (see Figure 3) on 4 March, 2010 at 16:39 UTC (08:39 a.m. local time). Over 6 km the radio-occultation (RO) profiles are outside the inner tomography grid. For this reason, only RO and Numerical Weather Prediction (NWP) models show estimates of wet refractivity ( ) over 6 km in Figure 13. The tomography solution for is the one referred to in Table 4 as "no stacking", with a priori data fed into the model only in the first epoch. The results of the WUELS, BIRA, and TUW tomography models were similar, distinguishing one inversion layer at heights between 4 and 6 km (we can define this as a anomaly). For the RO profile, there were also local minima at these heights; however, the largest inversion occurred below 4 km. The RO-and tomography-based values of and significantly differed from the ACCESS-A results in the lowest layer, as the latter did not show any inversion. The ERA-Interim result for this particular location was much closer to the tomography solution (BIRA and WUELS) and RO retrievals, showing a slight inversion between 5 and 6 km. The anomaly retrieved by both tomography models is a relevant piece of information for nowcasting.

Selected Results -Focus on Tomography Network in the East and during Selected Epochs
Based on the day of the previous investigated RO profile, we decided to verify the performance of the tomography models with respect to NWP forecasts in the eastern part of this model (orange/red pixels with high top cloud detected by GOME-2 in Figure 1a). The investigation was focused on the most accurate processing test A 30 and times 12:00 and 18:00 UTC on 4 March 2010 (see Figure 14).

Selected Results-Focus on Tomography Network in the East and during Selected Epochs
Based on the day of the previous investigated RO profile, we decided to verify the performance of the tomography models with respect to NWP forecasts in the eastern part of this model (orange/red pixels with high top cloud detected by GOME-2 in Figure 1a). The investigation was focused on the most accurate processing test A 30 and times 12:00 and 18:00 UTC on 4 March 2010 (see Figure 14).
First of all, we can see in Figure 14 that the BIRA and TUW models were in better agreement with ACCESS-A in the nine voxels (the south-eastern part of the network) than for all voxels, with the lowest discrepancy in the lower 10 km of the atmosphere being 0.15 ppm and 0.14 ppm, respectively, (for all the voxels, as shown in Figure 14a), and 0.10 ppm and 0.08 ppm, respectively, (for the nine voxels in the south-eastern part of the network, as shown in Figure 14b). At the same time, there was a much larger discrepancy between ACCESS-A and ERA-Interim in the same height range (0.23 for all the voxels and 0.29 for the nine voxels). Taking into account this disagreement, it should also be mentioned that the discrepancy between ERA-Interim and BIRA as well as ERA-Interim and TUW was smaller (about 0.18 ppm for all the voxels and 0.20 ppm for the nine voxels) in the height range 0-1.5 km than the respective discrepancies between ACCESS-A and ERA-Interim (0.27 ppm for all the voxels and 0.28 ppm for the nine voxels). For the 4-8 km layer, similar behaviour can be observed for model TUW (0.16 ppm for all voxels and 0.11 ppm for the nine voxels) versus ERA-Interim (0.20 ppm for both tests with nine or all voxels). The 1.5-4 km layer in the tomography models (BIRA and TUW) did not align well with ERA-Interim but aligned much better with the ACCESS-A solution (about 0.24 ppm versus 0.11 ppm for all voxels and 0.30 ppm versus 0.07 for nine ppm voxels). Note that ACCESS-A was better at modelling the 1.5-4 km layer in the region around Melbourne than ERA-Interim, and that a substantial improvement in the tomography model solution was most visible in this layer. First of all, we can see in Figure 14 that the BIRA and TUW models were in better agreement with ACCESS-A in the nine voxels (the south-eastern part of the network) than for all voxels, with the lowest discrepancy in the lower 10 km of the atmosphere being 0.15 ppm and 0.14 ppm, respectively, (for all the voxels, as shown in Figure 14a), and 0.10 ppm and 0.08 ppm, respectively, (for the nine voxels in the south-eastern part of the network, as shown in Figure 14b). At the same time, there was a much larger discrepancy between ACCESS-A and ERA-Interim in the same height range (0.23 for all the voxels and 0.29 for the nine voxels). Taking into account this disagreement, it should also be mentioned that the discrepancy between ERA-Interim and BIRA as well as ERA-Interim and TUW was smaller (about 0.18 ppm for all the voxels and 0.20 ppm for the nine voxels) in the height range 0-1.5 km than the respective discrepancies between ACCESS-A and ERA-Interim (0.27 ppm for all the voxels and 0.28 ppm for the nine voxels). For the 4-8 km layer, similar behaviour

Summary, Conclusions, and Perspectives for Future Works
This study aimed to conduct sensitivity tests and methodological improvements by cross-comparing five independent tomography models. GPS data from the Australian CORS network were considered during a severe weather event, i.e., the 6-8 March superstorm of greater Melbourne. The same dataset of GPS slant observations was fed into tomography models and statistical results (bias, standard deviation, and root mean square error) have been shown in reference to independent observations from radiosonde and radio-occultation profiles. The inverse problem treated by GPS tomography was ill-posed due to non-homogenous distribution of slant observations through the 3D grid, but this was also ill-conditioned due to the high number of parameters physically embedded. This explains why the stabilisation of the tomographic solution remains a challenging task. In this study, software based on different techniques were considered to test improvements of the stabilisation of solutions. Three means of testing were investigated and combined: (1) the improvement of the geometrical representativeness of retrievals, (2) the influence of the a priori condition in the convergence process, and (3) the use of the uncertainty of inputs on tomographic solutions. This study has tested a number of variants for which the impact on the quality of tomographic results has been assessed. These variants are, notably, the interest in stacking data and the use of pseudo-slant observations, the a priori condition considered in the tomography process, and the quality of initial observations and the impact of observation uncertainty. We used five tomographic models, namely, BIRA, WUELS, TUO, UBI, and TUW (see the Appendices for more details about these models); the processed results were compared to data obtained from ACCESS-A and ERA-Interim outputs, RS profiles, and RO products.
Tomography models produce parameters that are usually of similar or worse quality than ACCESS-A retrievals (see Figures 7,8,10,11 and 14) with few exceptions, i.e., in the bottom part of the troposphere for the BIRA and TUW tomography models (0-1.5 km), wet refractivity retrievals were of 0.5% better quality compared to ACCESS-A for tomography models fed with 6 h a priori data. Time stacking improved the solution with respect to non-stacking in all investigated cases without a continuous feed of a priori data (i.e., in tests B 30 , B 60 , and B 120 ). In addition, a short stacking period of 30 min (test B 30 ) gave a better response in the BIRA, WUELS, and UBI models (with 5, 3, and 17% respective improvements in the 0-8 km layer in comparison to test B), while for all other layers longer stacking worked better (120 min, test B 120 ): the lowest discrepancy between the tomography model and RS was found in the 8-13 km layer and the stacking test B 120 , which showed an 8% improvement in the BIRA solution in comparison to ACCESS-A. The use of pseudo-slant observations showed variable results. Sometimes it did not bring about any visible improvement in terms of RMS, but in some cases (e.g., in the 0-10 km layer), the use of a pseudo-SLANT GPS with 'stand-alone initialisation' obtained a consequent improvement (see tests B π and B π 30 in comparison to test B in Figure 8b; 0.56 g/m 3 and 1.05 g/m 3 were obtained, respectively, in comparison to 2.84 g/m 3 , for all the layers).
The BIRA solution of test B π 30 was indeed 1% better than that of ACCESS-A for the 8-13 km layer.
Our sensitivity tests showed that when using 30 min stacking data, an improvement was obtained relative to RS profiles (a 5% improvement in the 0-1.5 km layer was observed for the BIRA model and an 18% improvement in the 4-8 km layer was observed for the UBI model) without bringing too old data into the tomography process, allowing for a proper understanding of the meteorological situation. However, the increase in the geometrical distribution (Table 3) was weak using only 30 min stacking data (+1.5 %). We found that the use of pseudo-slants can significantly improve the geometrical matrix of tomography retrievals (i.e., by +25 %, showing a better special representativeness) and can bring about, in some cases, an improvement (a five times lower RMS value was obtained for all the layers when comparing test B π with test B for the BIRA model). The a priori condition tests on tomography show that the best results (with respect to RS) can be obtained using the best a priori. This is why, assuming that the ACCESS-A a priori was reasonably good (and close to the RS profiles), tomography retrievals, using regular a priori from NWP, obtain the best results (a three times lower RMS value was obtained for the 0-8 km layer, considering solutions from the BIRA, WUELS, and UBI models), illustrating that the "stand-alone" strategy (using the previous tomography retrievals as the a priori of the next calculation) is not always successful. A divergence in results is obtained, especially when the a priori is far from the real state of the atmosphere. Figure 15 summarises the results of water vapour density retrievals obtained for three types of a priori conditions (tests b, B, and A; see Figure 6), the use of 30 min stacking (B 30 and A 30 ), the use of pseudo-observations (B π and A π ), and the combined models (B π 30 and A π 30 ). The impact of the uncertainty of tomography inputs (error for model A) is one order of magnitude lower than the use of stacking or pseudo-observations. The results for the 0-8 km, 8-13 km, and 0-13 km layers are shown. The mean and best results highlight the recommendations of using 30 min stacking or pseudo-observations for the 0-8 km layer. The use of combined stacking and pseudo-observations shows consequent instability in the solutions of the three types of tomography software considered. This is due to the high number of data inputs ingested by tomography models and an unsatisfied adjustment for the top layers of the tomography grid. For this reason, the overall mean result of the 8-13 km layer does not show an improvement when using stacking and pseudo-observations. using the best a priori. This is why, assuming that the ACCESS-A a priori was reasonably good (and close to the RS profiles), tomography retrievals, using regular a priori from NWP, obtain the best results (a three times lower RMS value was obtained for the 0-8 km layer, considering solutions from the BIRA, WUELS, and UBI models), illustrating that the "stand-alone" strategy (using the previous tomography retrievals as the a priori of the next calculation) is not always successful. A divergence in results is obtained, especially when the a priori is far from the real state of the atmosphere. Figure 15 summarises the results of water vapour density retrievals obtained for three types of a priori conditions (tests b, B, and A; see Figure 6), the use of 30 min stacking (B 30 and A 30 ), the use of pseudo-observations (Bπ and Aπ), and the combined models (Bπ 30 and Aπ 30 ). The impact of the uncertainty of tomography inputs (error for model A) is one order of magnitude lower than the use of stacking or pseudo-observations. The results for the 0-8 km, 8-13 km, and 0-13 km layers are shown. The mean and best results highlight the recommendations of using 30 min stacking or pseudo-observations for the 0-8 km layer. The use of combined stacking and pseudo-observations shows consequent instability in the solutions of the three types of tomography software considered. This is due to the high number of data inputs ingested by tomography models and an unsatisfied adjustment for the top layers of the tomography grid. For this reason, the overall mean result of the 8-13 km layer does not show an improvement when using stacking and pseudo-observations. The results regarding the quality of initial observations and the impact of observation uncertainty are the following. The relative error of parameters related to GPS tomography (e.g., refractivity coefficient k2', SWD, and SIWV) had a moderate impact on the GPS tomography solution (going from 4% without stacking to 24% with 2 h stacking), as shown in Figures 11 and 15. However, an interesting result of this sensitivity test is that, for some cases with lower (or higher) values of SWDs or SIWVs, closer tomography estimates were obtained in comparison to RS profiles. With the use of strong and more realistic uncertainties applied positively and negatively, respectively, this was, finally, the four sets of SLANTGPS that we tested to assess the sensitivity of tomography The results regarding the quality of initial observations and the impact of observation uncertainty are the following. The relative error of parameters related to GPS tomography (e.g., refractivity coefficient k 2 ', SWD, and SIWV) had a moderate impact on the GPS tomography solution (going from 4% without stacking to 24% with 2 h stacking), as shown in Figures 11 and 15. However, an interesting result of this sensitivity test is that, for some cases with lower (or higher) values of SWDs or SIWVs, closer tomography estimates were obtained in comparison to RS profiles. With the use of strong and more realistic uncertainties applied positively and negatively, respectively, this was, finally, the four sets of SLANT GPS that we tested to assess the sensitivity of tomography retrievals. The order of magnitude in observations was transferred to tomography retrievals (e.g., in Figure 11, an increase in observation amplitudes shows an increase in normalised RMS from tomography with respect to RS, meaning that lower SLANT GPS applied in tomography showed closer results to RS). A key result obtained in this study is that a moistening anomaly was found in tomography retrievals (and RO profiles) compared to the NWP forecast. This is all the more interesting if such an anomaly is identified by several tomography models, as shown in Figure 13. We recommend the use of ensemble tomography (from several models) in nowcasting systems.
Interest in using stacking slant data and optimal configuration of the tomography grid has been investigated by [17] and involved consideration of synthetic data from a BASCOE (Belgian Assimilation System of Chemical Observations from Environmental satellite) NWP model and external validation from satellite sensors and RS retrievals. Building on conclusions from previous studies, this paper focused on comparison with independent external observations from RS and RO techniques (measurements in all weather conditions) across five days (involving quiet and severe conditions) to validate methodological improvements. It is important to notice that this study did not contain an evaluation of SWD and SIWV quality; SWDs and SIWVs were computed from ZTDs and horizontal gradients estimated from GPS observations and obtained using outputs of ground pressure and mean temperature from ACCESS-A NWP. Note that an interpolation in time (from 6 h to 30 min) of ACCESS-A outputs was considered to obtain SWDs and SIWVs, rather than considering observations of ground pressures and temperatures from synoptic stations. Fortunately, when the a priori condition from ACCESS-A was considered every 6 h in GPS tomography, the time interpolation was not a problem in the comparison with RS data. Hence, the contribution of hydrometeors to the delay was not considered. Errors in ground pressure and hydrometeor contributions can impact ZTDs measurements by up to a few centimetres in extreme weather [40], which can lead to an error in IWV of up to more than 5 kg/m 2 .
Finally, to conclude, we used independent external observations to assess the impact of several tests. Even better results (in comparison to the a priori condition from ACCESS-A) were obtained only for a few configurations and layers (i.e., an up to 0.5% improvement of normalised RMS in the 0-1.5 km layer for tomography models fed with 6h a priori data was observed, as well as an up to 1% improvement in the 8-13 km layer when using pseudo-observations), and interest in improving the geometrical matrix (with a priori every 6 h or "stand-alone" solutions used) has been clearly highlighted by comparing improved solutions with classical retrievals. We recommend using stacking data (using a maximum time window of 30 min) and pseudo-slant observations in case studies and nowcasting applications.
A hypothesis of straight ray propagating was considered in this study for modelling the path delay from GPS stations to ground-based receivers in our tomography models. Because the bending effect is negligible over 10 • [3,63,64], the inversion problem becomes linear and can be formulated using the discrete theory. Post-fit residuals cleaned from systematic effects, which were not used in this study, could be beneficial for GPS slant observations under severe weather conditions. These two aspects need to be considered in future work.