Atmospheric Water Vapor Variability over Houston: Continuous GNSS Tomography in the Year of Hurricane Harvey (2017)

Abstract

This study evaluates the capability of an unconstrained tomographic algorithm to capture 3D water vapor density variability throughout 2017 in Houston, U.S. The algorithm relies solely on Global Navigation Satellite System (GNSS) observations and does not require an initial guess or other specific constraints regarding water vapor density variability within the tomographic domain. The test domain, featuring 9 km horizontal, 500 m vertical, and 30 min temporal resolutions, yielded remarkable results when compared to data retrieved from the ECMWF Reanalysis v5 (ERA5), regional Weather Research and Forecasting Model (WRF) data, and GNSS-Radio Occultation (RO). For the first time, a time series of Precipitable Water Vapor maps derived from the Interferometric Synthetic Aperture Radar (InSAR) technique was used to validate the spatially integrated water vapor computed by GNSS tomography. Tomographic results clearly indicate the passage of Hurricane Harvey, with integrated water vapor peaking at 60 kg/m² and increased humidity at altitudes up to 7.5 km. Our findings suggest that GNSS tomography holds promise as a reliable source of atmospheric water vapor data for various applications. Future enhancements may arise from denser and multi-constellation networks.

1. Introduction

Water vapor is one of the most important greenhouse gases present in the troposphere. It plays a fundamental role in the hydrological cycle, Earth’s radiative budget, the atmosphere’s general circulation, and, more specifically, in the development of convective systems associated with extreme precipitation [1], which has been increasing with global warming. Its distribution in space and time is controlled by the atmospheric circulation on different scales (micro- to synoptic scale), making its observation by routine meteorological networks particularly challenging.

GNSS tomography is an innovative technique that attempts to use the Global Navigation Satellite System (GNSS) signals, from permanent and temporary GNSS networks, to reconstruct the 3D water vapor density distribution. It uses the interaction between the atmospheric water vapor and the electromagnetic waves that travel between the satellites and the GNSS antennas on Earth’s surface to infer the amount of water vapor along the signal path. This method utilizes numerous distinct views through the atmosphere within a specified time frame. The typical tomography system consists of using an integrated water vapor (or a related quantity) amount along the satellite–antenna path, retrieved from the GNSS processing, a 3D grid model that divides the atmosphere into voxels, and the discretization of signals (rays), which consists of finding the length of each ray inside the voxels crossed by it. In general, the solution consists of inverting an ill-posed problem, with incomplete data due to the heterogenous distribution of observations through the 3D tomography grid. One way to mitigate this limitation is to densify the GNSS network; however, that solution is costly and may not solve the problem depending on the grid horizontal resolution and the GNSS elevation cut-off mask used, typically 7–15° to avoid the delay error induced by incorporating low elevation observations. Another way consists of adding some inter-voxel physical constraints to add information to voxels not crossed by any ray (e.g., [2,3,4,5]) or adding a priori information from ground-based, air-borne, or space-borne instruments (e.g., [6,7,8,9]) from Numerical Weather Prediction models (NWP) data (e.g., [8,10,11,12,13]) or even creating virtual data (e.g., [14,15,16,17,18,19]).

Constraints and external information become necessary to overcome the rank deficiency problem introduced by empty voxels. However, the uncertainty introduced by that information can affect the tomography solution, restricting the mathematical model’s ability to fully simulate the dynamics of water vapor in the atmosphere and compromising the model’s reliability.

Another tomography challenge is the inversion algorithm. Some studies apply the generalized inverse method or a derivative [2,17,20,21,22] but these approaches can reveal numerical problems in the case of huge matrices and require large computational resources [23]. An alternative that avoids matrix inversion and even the creation of large matrices is the Algebraic Reconstruction Technique (ART) family. ART is an iterative algorithm that solves a system of linear equations by iteratively refining an initial guess [23]. A drawback of this technique is the requirement of solid initialization data, a general relaxation factor, and a criterion for stopping the iteration, which can depend on a specific configuration, leading to a degrading or oscillating solution [24]. Another technique applied in several studies is the Kalman Filter (KF), which combines the last tomographic solution and the observations using the corresponding error covariance matrices, providing an update of the error covariance for the tomographic solution, but realistic error estimates required in the KF are challenging to obtain [24,25,26,27]. An intercomparison between algorithms can be seen in [10,23,28,29,30].

To completely avoid the dependence on constraints and external data, in this study, we applied the unconstrained tomographic model proposed by [17], totally based on GNSS observations, to a one-year case study in Houston, U.S. We intend to validate the model’s performance in real-world scenarios and ensure that it meets reliability standards across various atmospheric conditions. The expected reliability standards involve metrics such as accuracy in capturing atmospheric variations, consistency across different temporal and spatial scales, and stability when applied in diverse weather scenarios.

In this approach, the rank of the design matrix is increased by applying observational angular interpolation and extrapolation. However, some changes were applied. The original inversion method uses the Moore–Penrose pseudoinverse algorithm [31], changed here to the Non-Negative Least Squares (NNLS) algorithm [32] to avoid unphysical negative values in the upper levels when the amount of water vapor is minimal. Next, we applied the approach implemented by [19] to increase the rank of the design matrix (contains information on the length of each ray per voxel), using the voxel geometric centers as a target to create a new angular domain (in both azimuth and elevation) instead of the regularly spaced angular domain initially proposed by [17].

Four different datasets were used to validate the tomography results. The European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis v5 (ERA5) data provides hourly estimates of a large number of atmospheric variables with a horizontal grid of 31 km [33]. The Weather Research and Forecasting Model (WRF), initialized by ERA5 data, was used here to improve the horizontal and vertical resolution of ERA5 [34]. Vertical profiles were retrieved from the Global Navigation Satellite System Radio Occultation (GNSS-RO) [35], and for the first time, we used precipitable water vapor (PWV) maps estimated by the Interferometric Synthetic-Aperture Radar (InSAR) technique with a horizontal resolution close to 300 m to evaluate the 2D integrated water vapor estimated by tomography, enabling the validation in areas lacking GNSS stations [36,37].

Using 1 year of data and the same tomography configuration, we evaluate the flexibility of the proposed algorithm to compute the 3D water vapor density variability for a full annual cycle in a humid subtropical climate with tropical influences. Additionally, we assess the algorithm’s capability to capture the signature of Hurricane Harvey. This devastating category 4 hurricane made landfall in Texas and Louisiana in August 2017, causing catastrophic flooding and resulting in over 100 deaths. While satellites primarily observed this event, in situ observations can fail due to harsh conditions. GNSS tomography, however, presents a breakthrough in evaluating the dynamics of severe weather, such as major hurricanes, providing detailed insights into atmospheric moisture conditions and offering a valuable tool for monitoring and understanding extreme weather events.

This study is organized as follows. Section 2 briefly describes the ERA5 and WRF model configuration, the GNSS-RO dataset, the corresponding InSAR dataset, and some details of the tomographic analysis methods. In Section 3, we present the results and discussion. A summary and concluding remarks are given in Section 4.

2. Materials and Methods

2.1. ERA5 Dataset

ERA5 data are produced by the Copernicus Climate Change Service (C3S) at ECMWF and is the fifth-generation ECMWF atmospheric reanalysis of the global climate covering the period from January 1940 to the present [33]. Compared to the discontinued ERA-Interim reanalysis model, ERA5 assimilates a larger number of observations, rising from an average of 0.75 million per day in 1979 to about 24 million in 2018, boosted mainly due to the increase in satellite radiances over time, as well as recent additions like GNSS-RO, scatterometer ocean vector wind and altimeter wave height data, ozone products, and ground-based radar observations [38]. It is freely available online by the C3S at https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-complete?tab=overview (accessed on 1 May 2024) at a horizontal resolution of 0.25° × 0.25° on 37 pressure levels from 1000 hPa to 1 hPa. Given the available ERA5 variables, the vertical water vapor density (${\rho }_v$) profiles, can be calculated by

(1)${\rho }_v={\displaystyle \frac{P·q}{287.055·T_v}}$

where $P$ is the atmospheric pressure (Pa), $T_v$ is the virtual temperature (K), and $q$ is the specific humidity (kg/kg) [39]. Hundreds of studies involving the ERA5 model have established it as one of the best reanalyses ever produced, providing an unparalleled dataset for validation in the absence of in situ observations.

2.2. WRF Dataset

The WRF model version 4.3 [34], available at https://github.com/wrf-model/WRF (accessed on 24 January 2024), was used to improve the grid resolution of ERA5. It was set up with two domains. The parent domain, centered on Houston, with a horizontal grid spacing of 12 × 12 km, covers 2,872,800 km² (1800 km in the east–west direction and 1596 km in the north–south direction). It is large enough to include the southeast part of the U.S., making it possible to simulate large hurricanes (like Hurricane Harvey) without creating discontinuities in the boundaries. The inner domain has a grid spacing of 4 × 4 km covering 1,044,975 km² (1084 km in the east–west direction and 964 km in the north–south direction). The model top reaches 10 hPa (about 30 km) and has 100 vertical levels. At 4 km horizontal resolution, the convection is explicitly solved. The remaining setup and the main physics options selected in the WRF model can be seen in [36]. As previously stated, the initial and lateral boundary conditions used in the parent domain are derived from the ERA5 data. Using the WRF output, ${\rho }_v$ is given by

(2)${\rho }_v={\displaystyle \frac{e}{461.51·T}}$

where $T$ is the air temperature (K), and $e$ is the water vapor pressure (Pa), given by

(3)$e={\displaystyle \frac{r·P}{0.622+r}}$

where $r$ is the water vapor mixing ratio (kg/kg) [39].

2.3. GNSS-RO Dataset

The GNSS-RO dataset used in this study includes data from several R.O. missions (COSMIC-1, KOMPSAT-5, MetOp-A, -B GRACE, TanDEM-X, and TerraSAR-X), all freely available from the COSMIC Data Analysis and Archive Center (CDAAC) (https://www.cosmic.ucar.edu/what-we-do/data-processing-center/data) (accessed on 12 April 2024). This study used the atmospheric profiles with moisture (wetPrf format) from the post-processed datasets (level 2). The CDAAC uses the one-dimensional variational approach (1D var) and the ERA-Interim model reanalysis as a first guest to inverse the observed refractivity to vertical profiles of temperature, atmospheric pressure, and water vapor pressure [40,41]. The vertical resolution is 100 m from the Earth’s surface to 40 km altitude. ${\rho }_v$ is calculated using Equation (2), and after checking the quality control flag, 30 profiles were selected to assess the tomography solution.

2.4. InSAR Dataset

The interferometric phase ($∆\varphi$) is a combination of several components given by

(4)$\Delta \varphi ={\Delta \varphi }_{suf}+{\Delta \varphi }_{iono}+{\Delta \varphi }_{hyd}+{\Delta \varphi }_{wet}+{\Delta \varphi }_{liq}+{\Delta \varphi }_{orb}+{\Delta \varphi }_{topo}+{\Delta \varphi }_{noise}$

where ${\Delta \varphi }_{suf}$ is the surface displacement contribution (due to tectonic activity or even anthropogenic effects), ${\Delta \varphi }_{iono}$ is the phase component related to the ionosphere, ${\Delta \varphi }_{hyd}$, ${\Delta \varphi }_{wet}$, and ${\Delta \varphi }_{liq}$ are the hydrostatic, non-hydrostatic (wet), and liquid contributions, respectively, ${\Delta \varphi }_{orb}$ and ${\Delta \varphi }_{topo}$ are the phase contributions related to inaccuracies in orbits and topographic distortions, respectively, and ${\Delta \varphi }_{noise}$ is the phase noise (decorrelation effects). The $\Delta$ symbol refers to a temporal variation [42]. In the absence of surface displacements, orbital and topographic effects (compensated by using precise orbits and a digital elevation model), and ionospheric contribution (can be neglected in the C-band SAR missions) and after removing the hydrostatic temporal variations (more stable in space and time, can be removed using NWP data), the remaining phase is related to the wet and liquid contribution. The liquid phase contribution is due to the liquid water mass in clouds and falling droplets (hydrometeors). Separating the wet (vapor) from the liquid part is challenging, requiring external measurements at the time of SAR acquisitions, e.g., by a meteorological weather radar. Under normal atmospheric conditions, this contribution can be disregarded; however, this may not hold when acquiring SAR images during a convective cell associated with extreme rainfall events [43,44].

In this study, 85 interferograms were generated using 112 Single-Look Complex (SLC) images acquired in 2017 by the C-band Sentinel-1 A&B Synthetic Aperture Radar (SAR) mission. Fifty SLC images were taken along the ascending orbit (path 34) and 62 along the descending orbit (path 143) between 12:00 and 12:30 UTC. Figure 1 shows both footprints. These images were processed using the SAR interferometry technique applied using the Sentinel Application Platform (SNAP) toolbox, available at https://step.esa.int/main/download/snap-download/ (accessed on 20 January 2024). The interferograms were processed based on the short temporal baseline method. This method consists of processing a time series of N SLC images to generate (N − 1) interferograms by selecting couples of subsequent SLC images with the shortest temporal baseline. Consequently, the phase contribution due to temporal decorrelation, especially in vegetated areas, and the Earth’s surface deformation, is diminished. Precise orbit ephemeris and a high-resolution digital elevation model (~30 m) were used to mitigate orbital and topographic phase contributions, respectively. A visual analysis of the interferometric fringes revealed no significant Earth surface deformation, indicating that the resulting interferometric phase primarily comprises contributions from ionospheric, tropospheric (hydrostatic, wet, and liquid phase components), and noise from other sources. The complex interferogram values are known only modulo-2π and must be unwrapped to achieve a quantitative interpretation. Phase unwrapping involves restoring the correct multiple of 2π to each point in the interferometric phase image [42]. In fact, successfully applying an unwrapping algorithm is one of the most challenging aspects of InSAR. In this way, to improve the success rate of the phase unwrapping, a Goldstein filter was applied to enhance the signal-to-noise ratio of the image [45]. We also used a spatial filter (multilook) to reduce the intensity of residues due to phase noise. The image will have less noise but less spatial resolution, changing from about 5 × 20 m² to 300 × 300 m² (in our case). Both filters are available in the SNAP toolbox. The unwrapping interferometry phase was computed using the algorithm proposed by [46].

On the hypothesis that all the phase contributions, except for the wet contribution, can be neglected or have been mitigated (see Equation (4)), we can estimate differential PWV maps from the unwrapping phase as

(5)$\Delta PWV=\mathsf{\Pi }·{\displaystyle \frac{\lambda }{4\pi }}·{\Delta \varphi }_{wet}·\mathrm{c}\mathrm{o}\mathrm{s}({\vartheta }_{look})$

where $\lambda$ is the radar wavelength corresponding to the central frequency of the Sentinel-1 signal, ${\vartheta }_{look}$ is the look angle along the swath, and $\mathsf{\Pi }$ is a proportionality constant that transforms signal delay in a corresponding PWV amount [47], given by

(6)$\mathsf{\Pi }={\displaystyle \frac{{10}^6}{\rho ·R_w·\left[\frac{k_3}{T_m}+k_2^{\prime }\right]}}$

where $\rho$ is the density of liquid water, $R_w$ is the specific gas constant for water vapor, $k_2^{\prime }$ and $k_3$ are refractivity constants [48], and $T_m$ is the weighted mean temperature of the atmosphere, given by

(7)$T_m=\int \left({\displaystyle \frac{e}{T}}\right)\mathrm{dz}/\int \left({\displaystyle \frac{e}{T^2}}\right)\mathrm{dz}$

with $e$ in hPa and $T$ in K. We used the ERA5 data to calculate $T_m$; details about this calculation can be seen in [49]. Note that while “mm” is commonly used for PWV, “kg/m²” is the correct physical unit. It aligns with the water budget equation and avoids confusion with ZWD, which is a delay computed as a distance. Our choice is consistent with [50].

Two essential aspects of InSAR data processing should be noted. Firstly, any geophysical quantity derived from the interferometric phase is inherently differential. Secondly, the unwrapping algorithm introduces an arbitrary constant to the phase. In order to obtain an absolute and calibrated quantity, the methodology proposed by [51] and improved by [52] was adopted here. This method is based on the least-squares approach to estimate single epochs and uses data from an NWP (in our case, from the ERA5 model) and GNSS observations (highlighted in red in Figure 1). It is important to note that there is no perfect method for deriving an absolute map from the InSAR time series; each technique has its own set of limitations. A comprehensive comparison of different methods, including their drawbacks, is detailed in [53].

It is now well established from a variety of studies that the InSAR technique can provide accurate PWV maps, achieving an accuracy of 1–2 mm at a spatial resolution of about 300 × 300 m² [42,43,49,54,55,56,57]. Furthermore, this technique uses microwave signals (most common are L-, C-, and X-band), showing advantages over other techniques (like visible light sensors) since it can operate day and night and is not affected by weather conditions such as clouds, making it a reliable data source under various atmospheric conditions [42]. It is sensitive enough to detect small-scale atmospheric phenomena, such as localized water vapor variations, which other observation methods often miss. It is used here to validate and complement other atmospheric observations, enhancing this study’s overall reliability and robustness. Details of the entire procedure, such as all the formulations used, are described in [58].

2.5. GNSS Dataset

A total of 27 Global Positioning System (GPS) permanent stations belonging to the Continuously Operating Reference Stations (CORS) network, managed by the National Oceanic and Atmospheric Administration (NOAA)/National Geodetic Survey (NGS), were used in this study. Figure 1 displays the station’s location (blue triangles for the tomographic grid). GPS observations in RINEX format are available online at https://geodesy.noaa.gov/CORS/ (accessed on 13 January 2024). We used the GAMIT/GLOBK (GG) package to process and analyze GPS data [59], available at https://geoweb.mit.edu/gg/ (accessed on 15 December 2023). The specific processing strategies are summarized in Table 1.

The GG output includes total, dry, and wet delays, as well as equivalent PWV values in the zenith direction and north–south ($G_{NS}$) and east–west ($G_{EW}$) gradients at a 10-degree elevation. The next step involves mapping the PWV to Slant Integrated Water Vapor (SIWV), considering the azimuth and elevation of each satellite and the effects of azimuthal asymmetry in the atmospheric delay. The contribution of the azimuthal asymmetry delay ($L_{az}$) can be calculated by

(8)$L_{az}(\epsilon ,\alpha )=G_{NS}·m_{az}\left(\epsilon \right)·\cos \left(\alpha \right)+G_{EW}·m_{az}\left(\epsilon \right)·\mathrm{s}\mathrm{i}\mathrm{n}\left(\alpha \right)$

where $\epsilon$ and $\alpha$ are the satellite elevation and azimuth angle, respectively, and $m_{az}$ is the gradient mapping function [60] that depends only on $\epsilon$, given by

(9)$m_{az}(\epsilon )={\displaystyle \frac{1}{\sin \left(\epsilon \right)·\tan \left(\epsilon \right)+0.0007}}$

Finally, the $\mathrm{S}\mathrm{I}\mathrm{W}\mathrm{V}$ is computed as

(10)$\mathrm{S}\mathrm{I}\mathrm{W}\mathrm{V}=PWV·m_w+L_{az}$

where $m_w$ is the wet mapping function (in our case, the VMF3).

To avoid day-boundary discontinuities, we adopt a multi-day approach. It consists of starting the run 6 h before (18 h of the previous day) and stopping 36 h later. Lastly, 12 h of processing are discarded (6 at the beginning and 6 h at the end), resulting in a seamless time series without discontinuities.

2.6. GNSS Tomographic Algorithm

The current study adopts the tomographic algorithm proposed by [17] (from now on, refer to MM21), with improvements reported in [18,19]. The advantages of MM21 include its simplicity of implementation, the absence of additional equations for imposing vertical and horizontal constraints, and the fact that there is no need for external observations for initialization. The tomography grid was configured with (7 × 10 × 20) 1400 voxels, featuring a horizontal resolution of 9 × 9 km² and a vertical resolution of 500 m (up to 10 km elevation). Figure 1 shows the tomographic grid footprint centered around Houston City. The tomographic domain was rotated by 30° to include more GPS stations and avoid a large number of empty voxels. The tomographic grid was configured with a vertical terrain-following coordinate [61]. This method transforms the vertical coordinate system to follow the contours of the Earth’s surface, ensuring that the vertical levels remain closer to the actual topography.

To overcome the matrix rank problem, new slant observations at different azimuths and elevations are computed by least-squares fitting a hypersurface to the SIWV observations [17], given by

(11)${\mathrm{S}\mathrm{I}\mathrm{W}\mathrm{V}}_{i,k}=a_i\cos \left({\alpha }_{i,k}\right)+b_i\sin \left({\alpha }_{i,k}\right)+{\displaystyle \frac{c_i}{\mathrm{s}\mathrm{i}\mathrm{n}\left({\epsilon }_{i,k}\right)}}+d_i (i=1,\dots ,n;k=1,\dots ,N)$

where ${\mathrm{S}\mathrm{I}\mathrm{W}\mathrm{V}}_{i,k}$ is the $k(=1,\dots ,N)$ slant observation at station $i$ within a specific time window, ${\alpha }_{i,k}$ and ${\epsilon }_{i,k}$ are the corresponding azimuth and elevation, and $N$ and $n$ are the number of slants (per time window) and stations, respectively. The coefficients $a_i, b_i, c_i, \mathrm{a}\mathrm{n}\mathrm{d} d_i (i=1,\dots ,n)$ are then used to create new slants at specific azimuths and elevations, as

(12)${\mathrm{S}\mathrm{I}\mathrm{W}\mathrm{V}}_j=a_i\cos \left({\alpha }_j\right)+b_i\sin \left({\alpha }_j\right)+{\displaystyle \frac{c_i}{\mathrm{s}\mathrm{i}\mathrm{n}\left({\epsilon }_j\right)}}+d_i (i=1,\dots ,n;j=1,\dots ,J)$

where $J$ is the number of targets (number of voxels); in our case, targeting the voxel centers from each station location is enough to obtain 100% rank. To exclude the slant contribution of the rays that enter the tomographic domain through its lateral boundaries, the ${\mathrm{S}\mathrm{I}\mathrm{W}\mathrm{V}}_j$ observation must be corrected. We used the equation proposed by [19] (an improvement of the original equation proposed by [17]), given by

(13)${\mathrm{S}\mathrm{I}\mathrm{W}\mathrm{V}}_j^{\prime }={\mathrm{S}\mathrm{I}\mathrm{W}\mathrm{V}}_j·e^{-z/H_v}$

where ${\mathrm{S}\mathrm{I}\mathrm{W}\mathrm{V}}_j^{\prime }$ is the corrected ${\mathrm{S}\mathrm{I}\mathrm{W}\mathrm{V}}_j$, $z$ is the height at which the rays enter the domain, and $H_v$ is the scale height of water vapor density. This formulation assumes that the background atmosphere is characterized by exponential decaying of the water vapor density. It was found in [19] that keeping $H_v$ constant (normally $H_v=2000$ [62]) does not produce the best tomographic results. Indeed, $H_v$ exhibits both seasonal and diurnal cycles, which influence the tomographic solution. In this study, we utilized ERA5 data to estimate an $H_v$ time series with a 1 h resolution; the downside is that our system is not entirely unconstrained. Finally, the solution is obtained by solving a determined (or underdetermined) system of linear equations as

(14)${\mathrm{S}\mathrm{I}\mathrm{W}\mathrm{V}}_{j\times 1}=D_{r\times j}·{\rho }_{v_{j\times 1}}$

where $D$ is the design matrix (contains the discretization of rays into voxels) and $r$ is the total number of rays. To solve Equation (14), we used the Non-Negative Least Squares (NNLS) algorithm [32] instead of the pseudoinverse. Although the NNLS algorithm is more computationally intensive, it ensures that the estimated parameters adhere to the natural or physical constraints of the problem domain.

2.7. Validation Tools

Four statistical metrics, including the Root Mean Square Error (RMSE), bias (BIAS), Pearson correlation coefficient (r), and Skill Score (SS), were used to evaluate the tomographic inversion performance. The first three are widely used and readily available. The SS is given by

(15)$\mathrm{SS}=100·\left(1-{\displaystyle \frac{{\mathrm{RMSE}}_{tomo,obs}}{{\mathrm{RMSE}}_{model,obs}}}\right)$

where ${\mathrm{RMSE}}_{tomo,obs}$ is the RMSE between the tomographic solution and the observations, and ${\mathrm{RMSE}}_{model,obs}$ is the RMSE between the model’s data and the observations. A positive SS value indicates an enhancement of the tomographic solution compared to the weather model data, with SS = 1 representing optimal performance. Conversely, a negative SS value signifies a deterioration of the tomographic accuracy, while SS = 0 indicates no change.

3. Results and Discussion

In this section, we evaluate the tomography solution for the year 2017, with a focused analysis on the month of August during Hurricane Harvey’s passage over Houston. The analysis will take three independent observational datasets as the reference: 30 min GNSS retrieved PWV at each station in the Houston network and, when available, high-resolution GNSS-RO water vapor density profiles and InSAR PWV maps. Those three datasets will allow an assessment of temporal, vertical and horizontal accuracy of the tomography against state-of-the-art NWP results.

Figure 2 compares the PWV values estimated by GNSS observations (taken as the reference) versus the PWV simulated by WRF, ERA5, and tomography (TOMO) at GNSS station locations. The WRF data show the worst RMSE value, reaching 5.18 kg/m², the ERA5 2.70 kg/m², and the TOMO 0.96 kg/m². The bias is positive for all three datasets, 2.19, 0.10, and 0.38 kg/m² for WRF, ERA5, and TOMO. The correlation coefficient is very high in all cases, in the same order: 0.947, 0.982, and 0.998. The PWV retrieved from tomography is very close to that observed by GNSS. These findings are consistent with those presented by [19]. The spread observed in WRF results (up to 40 kg/m² away from GNSS observations), which were driven by ERA5 boundary conditions, results from small-scale variability (more evident in later figures with InSAR data) that may be out of phase with the real world and is partially filtered out by GNSS processing.

A complementary analysis is presented in Figure 3a, which shows the average PWV over the tomographic domain. If we take TOMO results as the reference, as direct GNSS observations are only available at specific grid points, we find that ERA5 is much closer to tomography than WRF. Figure 3b displays a zoom over August when the maximum PWV value reached 76.4 kg/m² (observed by GNSS) during Hurricane Harvey’s passage (25 August to 1 September). The green stars represent the GNSS-RO PWV.

Figure 4 uses 30 profiles from GNSS-RO (stars in Figure 3) to assess the vertical distribution of water vapor density. Mean profiles for the all cases (Figure 4a) are close, with ERA5 being wetter than GNSS-RO, WRF, and tomography in the lower troposphere. Again, the profiles of tomography compare better against the reference profile (GNSS-RO) with an RMSE well below 1 g/m³ at all levels (Figure 4b). All retrieved profiles have a negative bias in the low troposphere (smaller in the case of tomography) and a less clear signal aloft. Tomography shows the best temporal correlation with GNSS-RO, exceeding 0.9 at all levels. Due to the sparse distribution of these observations, these results are not robust and show some noise. While GNSS-RO has been found to be a critical source of data for reanalysis [31], with unique high vertical resolution information, its reliability as an independent observation is not as strong as PWV from GNSS or InSAR. The good match with tomography is reassuring but does not substitute the need for in situ radiosonde observations, which are not available in Houston.

The final assessment uses 85 available InSAR observations. Figure 5 shows the InSAR PWV map obtained on 29 August at 00:26 UTC (at the end of Hurricane Harvey’s passage by Houston). With a temporal difference of approximately 5 min, we compared the corresponding PWV maps for WRF, ERA5 data, and the tomographic solution, the latter on a smaller domain. Notably, this specific InSAR map exhibits the highest PWV value among all InSAR cases, reaching 69.6 kg/m², and a pronounced southeast–northwest gradient with maximum values near the ocean, which was accurately simulated by ERA5 (Figure 5c). The WRF data also show similarities with InSAR but tend to depict a drier atmosphere in the northwest at higher altitudes. This difference may be attributed to WRF’s 4 km horizontal resolution, which explicitly represents deep convection processes and the heterogeneities of Earth’s topography in greater detail than traditional lower-resolution climate models like ERA5 [63]. Interestingly, the water concentration variability observed in the WRF map is less evident in the InSAR map, likely due to a certain degree of smoothing introduced during the computation of the absolute interferometric phase [51]. Taking InSAR as a reference, the ERA5 (Figure 5c) exhibits the best performance metrics, with an RMSE of 3.22 kg/m², followed closely by TOMO with 3.88 kg/m², and WRF with 5.83 kg/m². The BIAS shows similar values for all cases, in the same order: 3.08, 3.29, and 5.38 kg/m². The tomographic solution achieves a spatial correlation of 0.849, close to the value achieved by ERA5 of 0.924; the WRF only achieves a correlation of 0.596.

Figure 6 and Figure 7 compare the TOMO and model data with InSAR PWV maps. A spatial bilinear interpolation was performed to match the TOMO spatial resolution across all datasets before comparison. Figure 6 provides a temporal evaluation using InSAR maps as the “ground truth”. Averaging over all cases, the TOMO solution shows an RMSE of 1.41 kg/m² and a bias of 0.08 kg/m², representing a significant improvement compared to the WRF model, which has 3.61 and 1.35 kg/m², respectively, and ERA5, with 2.26 and −0.38 kg/m² (Figure 6a,b). The average TOMO correlation coefficient shows an improvement of about 0.4 over the WRF model. On average, the models’ PWV data do not exhibit a spatial correlation with the InSAR PWV data, with correlation coefficients of approximately 0.214 and −0.029 for WRF and ERA5, respectively. In only 10 cases, the TOMO solution reveals a lower correlation (<0.2), being greater than 0.6 for 50% of all cases. Upon closer inspection, it can be observed that the higher TOMO RMSE in the southern and northern sectors of the domain (Figure 6c) can be explained by the lower number of rays in the surface voxels, which significantly impacts the accuracy of the inversion. However, this is not easy to resolve, as the number of rays in the surface voxels, unlike the higher voxels, depends on the number and location of GNSS stations, as well as the region’s topography. In this study, we adhere to the initially proposed algorithm, incorporating the abovementioned changes without any additional modifications. However, we believe there is still room for improvement. Figure 7 presents a spatial evaluation, showing similar metric values as before (Figure 7a–c), except for the correlation values, which on average exceed 0.9 in all cases, with the highest values observed for TOMO. As previously stated, considering InSAR as the “ground truth”, we calculated the spatial skill score (SS) relative to WRF and ERA5 (Figure 7d,e). An average positive score of approximately 68%, ranging from 49% to 78%, was obtained for the WRF model, indicating a significant gain in accuracy for the TOMO solution across the entire domain. Relative to ERA5, an average positive SS of about 42%, ranging from −16% to 69%, was obtained. Notably, negative SS values in the southern part of the domain suggest ERA5’s superiority in that location.

The continuous evolution of the averaged tomographic vertical profile for WRF, ERA5, and TOMO for the month of August is shown in Figure 8. In the analysis above, both the WRF model and ERA5 exhibit a higher concentration of ${\rho }_v$ near the surface. Throughout this month, the integrated values range between 40 and 60 kg/m², except for a few days coinciding with the Hurricane Harvey episode (indicated by dashed vertical lines). The vertical distribution of ${\rho }_v$ is quite similar between the two models, but WRF tends to concentrate more moisture close to the surface, resulting in a more defined planetary boundary layer compared to ERA5 and TOMO. Remarkably, TOMO shows numerous similarities with the other two solutions. It captures significant oscillations, including short, dry, and wet episodes during extreme weather events. The vertical distribution of ${\rho }_v$ aligns with the two models, with approximately 50% of water concentrated below 2000 m, maintaining consistent top levels as simulated by the models. Notably, the tomographic solution exhibits a diurnal oscillation, also observed at higher altitudes between 6 August and 26 August—this behavior is not simulated by the models but is confirmed by three GNSS-RO vertical profiles. Analyzing the full dataset, we verified that TOMO could capture the annual cycle of water vapor density, showing a less-noisy behavior when the wet season is established.

An example of a vertical cross-section of the tomographic inversion for 22 August and 28 August at 00:00 UTC, corresponding to the southeast–northwest red dashed line in Figure 1, is shown in Figure 9. The tomography inversion shows some water vapor near the top of the domain during the Hurricane Harvey passage (see Figure 8 region between the vertical dashed lines and Figure 9c), only detected in this period, which could be an artifact in the tomographic inversion; however, validation remains challenging due to the absence of radiosondes or GNSS-RO observations during the hurricane event. Additionally, the TOMO results demonstrate striking similarities with WRF data, particularly on the 22nd. Disturbances of humid air at upper levels originating from lower levels are clearly captured by WRF and tomography (Figure 9a,c). The heterogeneous humid atmosphere observed on day 22 can be linked to turbulent wind field conditions, coinciding with the more vigorous effects of the hurricane. In contrast, on day 28, steadier wind fields prevail at high speeds. On day 28, the tomography is quite close to ERA5. Taking TOMO as reference, the RMSE for the two instants in Figure 9 is 0.8 kg/m², the BIAS is 0.3 kg/m², and the correlation is 0.995. Compared with the WRF data, the RMSE and BIAS increased only by 0.2 kg/m², and the correlation decreased by 0.008. Analyzing the entire 1-year dataset, ERA5 and WRF yielded mean RMSE values of 1.27 kg/m² and 1.34 kg/m², mean BIAS values of −0.03 kg/m² and 0.17 kg/m², and mean correlation coefficients of 0.958 and 0.953, respectively. These results show that the tomographic inversion is well within the observed range and captures most of the observed vertical water vapor density variability in different synoptic and seasonal time scales.

We employed the NNLS algorithm to derive the tomographic solution, avoiding the occurrence of negative values that can arise at higher altitudes when using the initially proposed pseudoinverse (pinv). Compared to the previous method, it takes on average 20–30% longer to find a positive solution. In our case study, inverting a tomographic domain with dimensions of 10 × 7 × 20 voxels requires approximately 20 s (CPU i7 13th with 32 GB RAM). As anticipated, significant differences between the two methods manifest at upper levels; however, they have a residual impact on the integrated water vapor values. Our hypothesis attributes the negative values to a poorly conditioned system at higher altitudes. Adding more satellites (from other constellations) may obviate the need to apply the NNLS algorithm, requiring confirmation through future studies.

4. Summary and Conclusions

The variability of water vapor in the troposphere presents a challenging task for weather models and observation techniques. It justifies investigating new methods to estimate local water vapor. In the present study, we tested the ability of the unconstrained GNSS tomography algorithm introduced by [17] with improvements added by [18,19] in direct estimates of the evolution of the water vapor density profile in Houston, U.S., based on GNSS observations from the Continuously Operating Reference Station (CORS) Network.

Results shown here indicate that GNSS tomography outperforms ERA5 reanalysis and convection permitting simulations with WRF in the description of the temporal variability of specific PWV observations (as given by direct 30 min GNSS retrievals), in the vertical distribution of water vapor density (as given by GNSS-RO), and in the spatial distribution of PWV (as given by InSAR). Since the two latter datasets are sparsely distributed, many studies have focused on the simple comparison with the always available ERA5 data, which is not an observation. Present results do indicate that ERA5 is very good and better than WRF in the representation of the temporal evolution of PWV and in the vertical distribution of water vapor density, although it compares less favorably with observations than tomographic results. In what concerns spatial variability, WRF shows small-scale features not unlike those observed in InSAR and that are absent from low-resolution ERA5 and appear mitigated in tomography. Due to inevitable phase errors (in space and time), those features do not lead to improved point statistics for WRF and are difficult to validate with current observations.

The tomographic solution was found to effectively maintain integrated water vapor values, not only at GNSS station locations but also in areas without stations, surpassing both weather models, as previously shown for the Amazon climate [19]. Remarkably, the vertical distribution of tomographic water vapor density closely aligns with GNSS-RO, exhibiting an RMSE below 1 g/m³ across all vertical levels and a correlation exceeding 0.9. The algorithm captures the annual water vapor density cycle, exhibiting reduced noise during the wet season and agreement with a vertical density distribution simulated by the two weather models. Results also show that the system responds sensibly to the passage of an extreme weather system, such as record-breaking Hurricane Harvey. However, validation through radiosonde profiles remains necessary. Since the tested tomographic algorithm does not require any external data to constrain the system inversion, it provides an independent observation of the atmospheric state, which can be implemented operationally at high temporal resolution.

The Houston GNSS network constitutes a rare example of a continuous high density GNSS network, freely available with several years of continuous data, and still operational. For the purposes of the present study, we selected the full year of 2017, because it has the advantage of containing a rare passage of a major Hurricane. The network was not designed for tomography and its horizontal distribution is very heterogeneous. That implied the choice of a relatively coarse horizontal resolution of 9 km, in a subdomain. The network is also limited to GPS in most stations, and there are no radiosonde observations in its vicinity. Present results, however, indicate the feasibility of tomographic inversions from such data and its added value against reanalysis and NWP data. Improved results are to be expected with better adapted networks operating with the full set of GNSS constellations (GPS, GLONASS, Galileo, and BeiDou).

Author Contributions

Conceptualization, P.M., J.C. and P.M.A.M.; methodology, P.M.; software, P.M. and R.F.; validation, P.M., J.C. and P.M.A.M.; formal analysis, P.M., J.C. and P.M.A.M.; investigation, P.M.; resources, P.M. and R.F.; data curation, P.M.; writing—original draft preparation, P.M.; writing—review and editing, J.C. and P.M.A.M.; visualization, J.C. and P.M.A.M.; supervision, P.M.; project administration, P.M.; funding acquisition, J.C. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors on request.

Acknowledgments

We thank the ECMWF for making the ERA5 hourly data available at Climate Data Store. We thank CDAAC for providing the GNSS-RO data level 2. The CORS/NOAA for providing the GNSS observations (rinex format). We also thank ESA for the Copernicus Open Access Hub, which provides the Sentinel-1 imagens acquisitions and makes the SNAP software (version 9) used to process the interferograms freely available.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1. GNSS tomographic 2D grid over Houston city (TX, USA). Blue triangles display the location of GPS stations used in the tomographic process, and red triangles are the stations used to obtain the absolute and calibrated InSAR PWV maps used to evaluate the tomographic solution spatial variability. The rectangles display the two footprints of Sentinel-1 used (ascending and descending orbits). The background colormap represents the elevation. The red dashed line indicates the vertical cross-section location used further on. The black line is the best track (determined by the National Hurricane Center) of Hurricane Harvey from 26 August to 29 August 2017 (hours in UTC).

Figure 2. TOMO-, WRF-, and ERA5-PWV derived values at the GNSS station’s location versus the PWV estimated from GNSS observations. The black line corresponds to the perfect fit. Legend indicates the correlation and the RMSE in kg/m2 (within []).

Figure 3. Temporal series of the average PWV over the topography domain. (a) Averaged WRF-, ERA5-, and TOMO-PWV; the green stars represent the PWV derived from RO. The vertical dashed lines delimit the month of August; (b) zoom over August. Date in month/day format.

Figure 4. Vertical profiles of TOMO, ERA5, and WRF assessed against 30 available GNSS-RO observations: (a) average water vapor density; (b) RMSE; (c) BIAS; and (d) temporal correlation coefficient.

Figure 5. PWV maps for 29 August around 00:30 UTC (at the end of Hurricane Harvey’s passage over the Houston area): (a) estimated using the interferometric phase (ascending orbit); (b) simulated by WRF model; (c) derived from ERA5 reanalysis; and (d) obtained via GNSS tomography (area corresponding to the TOMO footprint, black rectangle in (a–c)). The dotted rectangle corresponds to the InSAR footprint in ascending orbit.

Figure 6. TOMO and models temporal evaluation using InSAR PWV maps as the “ground truth”. (a) mean RMSE; (b) mean BIAS; and (c) mean correlation coefficient. Before applying the statistical metrics, a spatial bilinear resampling method (that uses the distance-weighted average of the four nearest pixel values to estimate a new pixel value) was performed to attain the TOMO spatial resolution (9 km).

Figure 7. Spatial evaluation using InSAR PWV maps. (a) Spatial RMSE for WRF; (b) for ERA5; and (c) for TOMO. (d) Skill score (SS) taking InSAR as “ground true”, WRF as first model, and TOMO as second; and (e) the same as before, but with ERA5 as the first model. The mean correlation and RMSE in kg/m2 are within [].

Figure 8. Hovmöller diagram of water vapor density profiles every 30 min. (a) WRF; (b) ERA5; and (c) TOMO solution. Averaged over the tomographic domain.

Figure 9. An example of a vertical cross-section of the tomographic inversion for 22 August at 00:00 UTC (a–c), corresponding to the beginning of the hurricane’s passage over Houston, and for 28th at 00:00 UTC (d–f) during Harvey near peak intensity over Houston. The cross-section corresponds to the southeast–northwest red dashed line in Figure 1. (a,d) WRF; (b,e) ERA5; and (c,f) TOMO solution.

Table 1

Specific processing strategies used in the GAMIT/GLOBK (GG) package.

Name of Strategy	Strategy Setting
Approach	DD (Double Difference)
Choice of observable	LC (Linear Combination of L1 and L2)
Cut-off elevation angle	5°
Sampling interval	30 s
Dry a priori model	SAAS (Saastamoinen model)
Dry and Wet mapping function	VMF3
Ocean tidal model	FES2004
Solid tide models	IERS10
Orbit and clock	IGS final products
Satellite and antenna phase center	IGS14 ANTEX file
Tropospheric Gradient Estimation	2 h
Tropospheric delay	$\mathrm{Random} \mathrm{walk} \mathrm{with} \mathrm{process} \mathrm{noise} \mathrm{of} 0.02 \mathrm{m}\sqrt{\mathrm{h}}$
Ionosphere effect	1-order (LC), 2- and 3-order correction

References

1. Georgiev, C.G.; Santurette, P.; Maynard, K. Diagnosis of Thermodynamic Environment of Deep Convection. Weather Analysis and Forecasting; Academic Press: Cambridge, MA, USA, 2016; pp. 157-222. [DOI: https://dx.doi.org/10.1016/B978-0-12-800194-3.00004-2]

2. Flores, A.; Ruffini, G.; Rius, A. 4D Tropospheric Tomography Using GPS Slant Wet Delays. Ann. Geophys.; 2000; 18, pp. 223-234. [DOI: https://dx.doi.org/10.1007/s00585-000-0223-7]

3. Rohm, W.; Bosy, J. The Verification of GNSS Tropospheric Tomography Model in a Mountainous Area. Adv. Space Res.; 2011; 47, pp. 1721-1730. [DOI: https://dx.doi.org/10.1016/j.asr.2010.04.017]

4. Zhao, Q.; Yao, Y.; Yao, W.; Xia, P. An Optimal Tropospheric Tomography Approach with the Support of an Auxiliary Area. Ann. Geophys.; 2018; 36, pp. 1037-1046. [DOI: https://dx.doi.org/10.5194/angeo-36-1037-2018]

5. Yang, F.; Guo, J.; Shi, J.; Meng, X.; Zhao, Y.; Zhou, L.; Zhang, D. A GPS Water Vapour Tomography Method Based on a Genetic Algorithm. Atmos. Meas. Tech.; 2020; 13, pp. 355-371. [DOI: https://dx.doi.org/10.5194/amt-13-355-2020]

6. Benevides, P.; Catalao, J.; Miranda, P.M.A. On the Inclusion of GPS Precipitable Water Vapour in the Nowcasting of Rainfall. Nat. Hazards Earth Syst. Sci.; 2015; 15, pp. 2605-2616. [DOI: https://dx.doi.org/10.5194/nhess-15-2605-2015]

7. Benevides, P.; Catalao, J.; Nico, G.; Miranda, P.M.A. 4D Wet Refractivity Estimation in the Atmosphere Using GNSS Tomography Initialized by Radiosonde and AIRS Measurements: Results from a 1-Week Intensive Campaign. GPS Solut.; 2018; 22, 91. [DOI: https://dx.doi.org/10.1007/s10291-018-0755-5]

8. Chen, B.; Liu, Z. Assessing the Performance of Troposphere Tomographic Modeling Using Multi-Source Water Vapor Data during Hong Kong’s Rainy Season from May to October 2013. Atmos. Meas. Tech.; 2016; 9, pp. 5249-5263. [DOI: https://dx.doi.org/10.5194/amt-9-5249-2016]

9. Xia, P.; Cai, C.; Liu, Z. GNSS Troposphere Tomography Based on Two-Step Reconstructions Using GPS Observations and COSMIC Profiles. Ann. Geophys.; 2013; 31, pp. 1805-1815. [DOI: https://dx.doi.org/10.5194/angeo-31-1805-2013]

10. Notarpietro, R.; Cucca, M.; Gabella, M.; Venuti, G.; Perona, G. Tomographic Reconstruction of Wet and Total Refractivity Fields from GNSS Receiver Networks. Adv. Space Res.; 2011; 47, pp. 898-912. [DOI: https://dx.doi.org/10.1016/j.asr.2010.12.025]

11. Trzcina, E.; Hanna, N.; Kryza, M.; Rohm, W. TOMOREF Operator for Assimilation of GNSS Tomography Wet Refractivity Fields in WRF DA System. J. Geophys. Res. Atmos.; 2020; 125, e2020JD032451. [DOI: https://dx.doi.org/10.1029/2020JD032451]

12. Haji-Aghajany, S.; Amerian, Y.; Verhagen, S. B-Spline Function-Based Approach for GPS Tropospheric Tomography. GPS Solut.; 2020; 24, 88. [DOI: https://dx.doi.org/10.1007/s10291-020-01005-x]

13. Zhao, Q.; Yao, Y. An Improved Troposphere Tomographic Approach Considering the Signals Coming from the Side Face of the Tomographic Area. Ann. Geophys.; 2017; 35, pp. 87-95. [DOI: https://dx.doi.org/10.5194/angeo-35-87-2017]

14. Yang, F.; Guo, J.; Shi, J.; Zhou, L.; Xu, Y.; Chen, M. A Method to Improve the Distribution of Observations in GNSS Water Vapor Tomography. Sensors; 2018; 18, 2526. [DOI: https://dx.doi.org/10.3390/s18082526] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30072630]

15. Adavi, Z.; Mashhadi-Hossainali, M. 4D Tomographic Reconstruction of the Tropospheric Wet Refractivity Using the Concept of Virtual Reference Station, Case Study: Northwest of Iran. Meteorol. Atmos. Phys.; 2014; 126, pp. 193-205. [DOI: https://dx.doi.org/10.1007/s00703-014-0342-4]

16. Yao, Y.; Liu, C.; Xu, C. A New GNSS-Derived Water Vapor Tomography Method Based on Optimized Voxel for Large GNSS Network. Remote Sens.; 2020; 12, 2306. [DOI: https://dx.doi.org/10.3390/rs12142306]

17. Miranda, P.M.A.; Mateus, P. A New Unconstrained Approach to GNSS Atmospheric Water Vapor Tomography. Geophys. Res. Lett.; 2021; 48, e2021GL094852. [DOI: https://dx.doi.org/10.1029/2021GL094852]

18. Miranda, P.M.A.; Mateus, P. Improved GNSS Water Vapor Tomography With Modified Mapping Functions. Geophys. Res. Lett.; 2022; 49, e2022GL100140. [DOI: https://dx.doi.org/10.1029/2022GL100140]

19. Miranda, P.M.A.; Adams, D.K.; Tomé, R.; Fernandes, R.; Mateus, P. Optimizing Boundary Conditions in GNSS Tomography: A Continuous 7-Month Case Study in the Amazon. Geophys. Res. Lett.; 2023; 50, e2023GL105030. [DOI: https://dx.doi.org/10.1029/2023GL105030]

20. Haji-Aghajany, S.; Amerian, Y. Hybrid Regularized GPS Tropospheric Sensing Using 3-D Ray Tracing Technique. IEEE Geosci. Remote Sens. Lett.; 2018; 15, pp. 1475-1479. [DOI: https://dx.doi.org/10.1109/LGRS.2018.2853183]

21. Flores, A.; De Arellano, J.V.G.; Gradinarsky, L.P.; Rius, A. Tomography of the Lower Troposphere Using a Small Dense Network of GPS Receivers. IEEE Trans. Geosci. Remote Sens.; 2001; 39, pp. 439-447. [DOI: https://dx.doi.org/10.1109/36.905252]

22. Champollion, C.; Masson, F.; Bouin, M.N.; Walpersdorf, A.; Doerflinger, E.; Bock, O.; Van Baelen, J. GPS Water Vapour Tomography: Preliminary Results from the ESCOMPTE Field Experiment. Atmos. Res.; 2005; 74, pp. 253-274. [DOI: https://dx.doi.org/10.1016/j.atmosres.2004.04.003]

23. Bender, M.; Dick, G.; Ge, M.; Deng, Z.; Wickert, J.; Kahle, H.G.; Raabe, A.; Tetzlaff, G. Development of a GNSS Water Vapour Tomography System Using Algebraic Reconstruction Techniques. Adv. Space Res.; 2011; 47, pp. 1704-1720. [DOI: https://dx.doi.org/10.1016/j.asr.2010.05.034]

24. Foken, T. Springer Handbook of Atmospheric Measurements; Springer Handbooks Springer International Publishing: Cham, Switzerland, 2021; ISBN 978-3-030-52170-7

25. Rohm, W.; Yuan, Y.; Biadeglgne, B.; Zhang, K.; Marshall, J. Le Ground-Based GNSS ZTD/IWV Estimation System for Numerical Weather Prediction in Challenging Weather Conditions. Atmos. Res.; 2014; 138, pp. 414-426. [DOI: https://dx.doi.org/10.1016/j.atmosres.2013.11.026]

26. Perler, D.; Geiger, A.; Hurter, F. 4D GPS Water Vapor Tomography: New Parameterized Approaches. J. Geod.; 2011; 85, pp. 539-550. [DOI: https://dx.doi.org/10.1007/s00190-011-0454-2]

27. Yang, F.; Gong, X.; Wang, Y.; Liu, M.; Li, J.; Xu, T.; Hao, R. GNSS Water Vapor Tomography Based on Kalman Filter with Optimized Noise Covariance. GPS Solut.; 2023; 27, 181. [DOI: https://dx.doi.org/10.1007/s10291-023-01517-2]

28. Bonafoni, S.; Biondi, R.; Brenot, H.; Anthes, R. Radio Occultation and Ground-Based GNSS Products for Observing, Understanding and Predicting Extreme Events: A Review. Atmos. Res.; 2019; 230, 104624. [DOI: https://dx.doi.org/10.1016/j.atmosres.2019.104624]

29. Saxena, S.; Dwivedi, R. GNSS Ground-Based Tomography: State-of-the-Art and Technological Challenges. Int. J. Remote Sens.; 2023; 44, pp. 5313-5343. [DOI: https://dx.doi.org/10.1080/01431161.2023.2247526]

30. Yang, F.; Sun, Y.; Meng, X.; Guo, J.; Gong, X. Assessment of Tomographic Window and Sampling Rate Effects on GNSS Water Vapor Tomography. Satell. Navig.; 2023; 4, 7. [DOI: https://dx.doi.org/10.1186/s43020-023-00096-4]

31. Penrose, R. A Generalized Inverse for Matrices. Math. Proc. Camb. Philos. Soc.; 1955; 51, pp. 406-413. [DOI: https://dx.doi.org/10.1017/S0305004100030401]

32. Lawson, C.L.; Hanson, R.J. Solving Least Squares Problems; Society for Industrial and Applied Mathematics Prentice-Hall, Inc.: Englewood Cliffs, NJ, USA, 1995; ISBN 978-0-89871-356-5 [DOI: https://dx.doi.org/10.1137/1.9781611971217.fm]

33. Hersbach, H.; Bell, B.; Berrisford, P.; Hirahara, S.; Horányi, A.; Muñoz-Sabater, J.; Nicolas, J.; Peubey, C.; Radu, R.; Schepers, D. et al. The ERA5 Global Reanalysis. Q. J. R. Meteorol. Soc.; 2020; 146, pp. 1999-2049. [DOI: https://dx.doi.org/10.1002/qj.3803]

34. Skamarock, W.C.; Klemp, J.B.; Dudhia, J.; Gill, D.O.; Liu, Z.; Berner, J.; Wang, W.; Powers, J.G.; Duda, M.G.; Barker, D.M. et al. A Description of the Advanced Research WRF Model Version 4.3 (No. NCAR/TN-556+STR). Available online: https://opensky.ucar.edu/islandora/object/opensky:2898 (accessed on 1 August 2024).

35. Anthes, R.A. Exploring Earth’s Atmosphere with Radio Occultation: Contributions to Weather, Climate and Space Weather. Atmos. Meas. Tech.; 2011; 4, pp. 1077-1103. [DOI: https://dx.doi.org/10.5194/amt-4-1077-2011]

36. Mateus, P.; Miranda, P.M.A. Using InSAR Data to Improve the Water Vapor Distribution Downstream of the Core of the South American Low-Level Jet. J. Geophys. Res. Atmos.; 2022; 127, e2021JD036111. [DOI: https://dx.doi.org/10.1029/2021JD036111]

37. Mateus, P.; Miranda, P.M.A.; Nico, G.; Catalao, J. Continuous Multitrack Assimilation of Sentinel-1 Precipitable Water Vapor Maps for Numerical Weather Prediction: How Far Can We Go With Current InSAR Data?. J. Geophys. Res. Atmos.; 2021; 126, e2020JD034171. [DOI: https://dx.doi.org/10.1029/2020JD034171]

38. Hoffmann, L.; Günther, G.; Li, D.; Stein, O.; Wu, X.; Griessbach, S.; Heng, Y.; Konopka, P.; Müller, R.; Vogel, B. et al. From ERA-Interim to ERA5: The Considerable Impact of ECMWF’s next-Generation Reanalysis on Lagrangian Transport Simulations. Atmos. Chem. Phys.; 2019; 19, pp. 3097-3124. [DOI: https://dx.doi.org/10.5194/acp-19-3097-2019]

39. Stull, R.B. Practical Meteorology: An Algebra-Based Survey of Atmospheric Science; 1.02b. Sundog Publishing, LLC: Vancouver, BC, Canada, 2017.

40. Neiman, P.J.; Ralph, F.M.; Wick, G.A.; Kuo, Y.-H.; Wee, T.-K.; Ma, Z.; Taylor, G.H.; Dettinger, M.D. Diagnosis of an Intense Atmospheric River Impacting the Pacific Northwest: Storm Summary and Offshore Vertical Structure Observed with COSMIC Satellite Retrievals. Mon. Weather Rev.; 2008; 136, pp. 4398-4420. [DOI: https://dx.doi.org/10.1175/2008MWR2550.1]

41. Palmer, P.I.; Barnett, J.J.; Eyre, J.R.; Healy, S.B. A Nonlinear Optimal Estimation Inverse Method for Radio Occultation Measurements of Temperature, Humidity, and Surface Pressure. J. Geophys. Res. Atmos.; 2000; 105, pp. 17513-17526. [DOI: https://dx.doi.org/10.1029/2000JD900151]

42. Hanssen, R.F. Radar Interferometry, Data Interpretation and Error Analysis; Remote Sensing and Digital Image Processing Springer: Dordrecht, The Netherlands, 2001; Volume 2, ISBN 978-0-7923-6945-5

43. Kinoshita, Y.; Shimada, M.; Furuya, M. InSAR Observation and Numerical Modeling of the Water Vapor Signal during a Heavy Rain: A Case Study of the 2008 Seino Event, Central Japan. Geophys. Res. Lett.; 2013; 40, pp. 4740-4744. [DOI: https://dx.doi.org/10.1002/grl.50891]

44. Brenot, H.; Ducrocq, V.; Walpersdorf, A.; Champollion, C.; Caumont, O. GPS Zenith Delay Sensitivity Evaluated from High-Resolution Numerical Weather Prediction Simulations of the 8–9 September 2002 Flash Flood over Southeastern France. J. Geophys. Res. Atmos.; 2006; 111, 15105. [DOI: https://dx.doi.org/10.1029/2004JD005726]

45. Goldstein, R.M.; Werner, C.L. Radar Interferogram Filtering for Geophysical Applications. Geophys. Res. Lett.; 1998; 25, pp. 4035-4038. [DOI: https://dx.doi.org/10.1029/1998GL900033]

46. Nico, G.; Fortuny, J. Using the Matrix Pencil Method to Solve Phase Unwrapping. IEEE Trans. Signal Process.; 2003; 51, pp. 886-888. [DOI: https://dx.doi.org/10.1109/TSP.2002.806990]

47. Bevis, M. GPS Meteorology: Mapping Zenith Wet Delays onto Precipitable Water. J. Appl. Meteorol.; 1994; 33, pp. 379-386. [DOI: https://dx.doi.org/10.1175/1520-0450(1994)033<0379:GMMZWD>2.0.CO;2]

48. Jean Rüeger, A.M. Refractive Index Formulae for Radio Waves. Proceedings of the FIG XXII International Congress; Washington, DC, USA, 19–26 April 2002; pp. 19-26.

49. Mateus, P.; Nico, G.; Catalão, J. Maps of PWV Temporal Changes by SAR Interferometry: A Study on the Properties of Atmosphere’s Temperature Profiles. IEEE Geosci. Remote Sens. Lett.; 2014; 11, pp. 2065-2069. [DOI: https://dx.doi.org/10.1109/LGRS.2014.2318993]

50. Mendes, V.B. Modeling the Neutral-Atmospheric Propagation Delay in Radiometric Space Techniques; University of New Brunswick: Fredericton, NB, Canada, 1999; 353.Available online: https://gge.ext.unb.ca/Pubs/TR199.pdf (accessed on 1 August 2024).

51. Mateus, P.; Catalão, J.; Nico, G.; Benevides, P. Mapping Precipitable Water Vapor Time Series from Sentinel-1 Interferometric SAR. IEEE Trans. Geosci. Remote Sens.; 2020; 58, pp. 1373-1379. [DOI: https://dx.doi.org/10.1109/TGRS.2019.2946077]

52. Mateus, P.; Catalão, J.; Nico, G. Improving the Accuracy and Spatial Resolution of ERA5 Precipitable Water Vapor Using InSAR Data. IEEE Geosci. Remote Sens. Lett.; 2024; 21, 1001905. [DOI: https://dx.doi.org/10.1109/LGRS.2024.3379249]

53. Meroni, A.N.; Montrasio, M.; Venuti, G.; Barindelli, S.; Mascitelli, A.; Manzoni, M.; Monti-Guarnieri, A.V.; Gatti, A.; Lagasio, M.; Parodi, A. et al. On the Definition of the Strategy to Obtain Absolute InSAR Zenith Total Delay Maps for Meteorological Applications. Front. Earth Sci.; 2020; 8, 359. [DOI: https://dx.doi.org/10.3389/feart.2020.00359]

54. Remy, D.; Chen, Y.; Froger, J.L.; Bonvalot, S.; Cordoba, L.; Fustos, J. Revised Interpretation of Recent InSAR Signals Observed at Llaima Volcano (Chile). Geophys. Res. Lett.; 2015; 42, pp. 3870-3879. [DOI: https://dx.doi.org/10.1002/2015GL063872]

55. Wadge, G.; Webley, P.W.; James, I.N.; Bingley, R.; Dodson, A.; Waugh, S.; Veneboer, T.; Puglisi, G.; Mattia, M.; Baker, D. et al. Atmospheric Models, GPS and InSAR Measurements of the Tropospheric Water Vapour Field over Mount Etna. Geophys. Res. Lett.; 2002; 29, 11. [DOI: https://dx.doi.org/10.1029/2002GL015159]

56. Zebker, H.A.; Rosen, P.A.; Hensley, S. Atmospheric Effects in Interferometric Synthetic Aperture Radar Surface Deformation and Topographic Maps. J. Geophys. Res. Solid Earth; 1997; 102, pp. 7547-7563. [DOI: https://dx.doi.org/10.1029/96JB03804]

57. Mateus, P.; Catalão, J.; Nico, G. Sentinel-1 Interferometric SAR Mapping of Precipitable Water Vapor over a Country-Spanning Area. IEEE Trans. Geosci. Remote Sens.; 2017; 55, pp. 2993-2999. [DOI: https://dx.doi.org/10.1109/TGRS.2017.2658342]

58. Mateus, P.; Miranda, P.M.A.; Nico, G.; Catalão, J.; Pinto, P.; Tomé, R. Assimilating InSAR Maps of Water Vapor to Improve Heavy Rainfall Forecasts: A Case Study with Two Successive Storms. J. Geophys. Res. Atmos.; 2018; 123, pp. 3341-3355. [DOI: https://dx.doi.org/10.1002/2017JD027472]

59. Herring, T.A.; King, R.W.; Floyd, M.A.; McClusky, S.C. GAMIT Reference Manual Reference Manual GPS Analysis at MIT; Department of Earth, Atmospheric, and Planetary Sciences Institute of Technology: Cambridge, MA, USA, 2018; Available online: http://geoweb.mit.edu/gg/docs.php (accessed on 1 August 2024).

60. Chen, G.; Herring, T.A. Effects of Atmospheric Azimuthal Asymmetry on the Analysis of Space Geodetic Data. J. Geophys. Res. Solid Earth; 1997; 102, pp. 20489-20502. [DOI: https://dx.doi.org/10.1029/97JB01739]

61. Gal-Chen, T.; Somerville, R.C.J. On the Use of a Coordinate Transformation for the Solution of the Navier-Stokes Equations. J. Comput. Phys.; 1975; 17, pp. 209-228. [DOI: https://dx.doi.org/10.1016/0021-9991(75)90037-6]

62. Tomasi, C. Vertical Distribution Features of Atmospheric Water Vapor in the Mediterranean, Red Sea, and Indian Ocean. J. Geophys. Res. Atmos.; 1984; 89, pp. 2563-2566. [DOI: https://dx.doi.org/10.1029/JD089iD02p02563]

63. Prein, A.F.; Rasmussen, R.; Castro, C.L.; Dai, A.; Minder, J. Special Issue: Advances in Convection-Permitting Climate Modeling. Clim. Dyn.; 2020; 55, pp. 1-2. [DOI: https://dx.doi.org/10.1007/s00382-020-05240-3]