2.1 The ZTD and tropospheric gradient operator

The ZTD is calculated through

The refractivity, Ψ, is a function of pressure, temperature, and humidity (Thayer, 1974), with z denoting the height above the station. The forward operator for the ZTD, along with the tangent linear and adjoint operators, is already integrated into the WRFDA system. For more details, please refer to the three respective routines in the WRFDA code.

1. In da_get_innov_vector_gpsztd.inc, the innovation (i.e., the difference between the observed and forward-modeled ZTD) is calculated.

2. The routine da_transform_xtoy_gpsztd.inc represents the corresponding tangent linear code.

3. On the other hand, the routine da_transform_xtoy_gpsztd_adj.inc represents the corresponding adjoint code.

Regarding tropospheric gradients, we have identified two possible approaches for the forward operator, which we refer to as the rigorous and fast methods. The rigorous approach involves the following steps: we compute numerous tropospheric delays, considering various elevation and azimuth angles for the given station location, and then perform a least-squares fit to obtain the north and east gradient components (as described in Zus et al., 2018). The rigorous approach aims at replicating how tropospheric gradients are estimated in GNSS analysis. However, it has a drawback – it can be challenging to implement it into DA systems due to the need for a ray-tracing algorithm (Zus et al., 2015).

Therefore, we prefer the fast approach, which works as follows: for the given station location, we utilize a closed-form expression that depends on the north–south and east–west horizontal gradients of refractivity (as outlined in Davis et al., 1993). This enables the calculation of the north and east gradient components through

Here, x, y, and z represent Cartesian coordinates and the subscripts denote partial derivatives. Therefore, similarly to ZTD, the tropospheric gradients are computed through numerical integration. In essence, we define a sequence of integration points, calculate the integrant at these points, and then sum them up along with their respective interpolation weights (as described in Zus et al., 2012). The critical aspect lies in the computation of the horizontal refractivity gradients, Ψ_y and Ψ_x, at the location of these integration points. To start, we express these gradients in terms of the station's longitude (λ) and latitude (ϕ) as follows:

where ρ represents the radial distance or the distance to the center of the osculating sphere. Our proposed method for calculating the horizontal refractivity gradients, Ψ_ϕ and Ψ_λ, involves a least-squares adjustment. Essentially, we establish a relationship between the refractivity, Ψ, and the horizontal refractivity gradients, Ψ_ϕ and Ψ_λ, at a specific (geometric) height, h, and the vertically adjusted refractivity at neighboring grid points of the weather model using Taylor’s series expansion.

Through vertical interpolation, the vertically adjusted refractivity, Ψ_i, can be determined by

Above the weather model top, the vertically adjusted refractivity, Ψ_i, is determined using the hydrostatic equation

where ${\mathrm{\Psi }}_i^n$, $T_i^n$, and $h_i^n$ represent the respective refractivity, temperature, and height at the top of the weather model, and G is the hydrostatic constant. In this equation, we approximate the geopotential height by the geometric height. To obtain the horizontal refractivity gradients, Ψ_ϕ and Ψ_λ, we invert the linear system of Eq. (8) using the least-squares method. Consequently, the horizontal refractivity gradient at a specific height can be expressed as a linear combination of the vertically adjusted refractivity, Ψ_i:

where w and v represent the interpolation weights. The least-squares adjustment incorporates grid points within a 35 km radius to compute the horizontal refractivity gradients. For example, assuming the weather model has a 10 km horizontal resolution, this computation involves 8×8 grid points. The selection of grid points included in the least-squares fit is not arbitrary; it is justified by comparing tropospheric gradients derived from the rigorous and fast approaches (as discussed in Zus et al., 2023).

The tangent linear operator is derived by applying the chain rule of differential calculus in forward mode. Similarly, the adjoint operator is derived by applying the same in the reverse mode (Giering and Kaminski, 1998). However, two simplifications are introduced. First, the geometric height, derived from the geopotential height (the natural height of the weather model), is not treated as a control variable. In other words, the partial derivatives in the construction of the tangent linear (adjoint) code are ignored. This simplification is also applied to other operators, such as the operator for radio occultation refractivity profiles. Second, the temperature we use to compute the refractivity above the weather model's top is not treated as a control variable. This is not problematic because tropospheric gradients are not very sensitive to refractivity at high altitudes (Zus et al., 2018).

The newly developed operator is not implemented as a stand-alone entity. Instead, we integrated the tropospheric gradient operator as an add-on to the existing ZTD operator. For each respective station, the code for tropospheric gradients is automatically executed alongside the ZTD code. To control the assimilation of tropospheric gradients, users can use a single variable in the namelist, called use_gpsgraobs, to switch it on or off. Therefore, for in-depth information about the forward operator, readers are directed to the routine da_get_innov_vector_gpsztd.inc. Details on the tangent linear and adjoint operators can be found in the routines da_transform_xtoy_gpsztd.inc and da_transform_xtoy_gpsztd_adj.inc, respectively.

We would like to note that the MPI implementation of the newly developed operator presents a unique challenge. While the existing ZTD operator is considered “local” because it utilizes variables from nearby (four surrounding) grid points, the tropospheric gradient operator is considered “non-local” due to the inclusion of numerous nearby grid points when calculating horizontal refractivity gradient components. This necessitates adjustments to the size of the “halo” region. Instead, we opted to implement an approach similar to the one used for the non-local excess-phase path operator (Chen et al., 2009). This approach ensures that the global refractivity (and temperature) field is available to individual processors. For more in-depth information, readers are encouraged to consult the WRFDA system documentation.

3.1 WRF model and configuration

We employed the non-hydrostatic WRF model version 4.4.1 for our impact study, implementing the gradient observation operator. WRF has a strong track record in both the research community and the operational forecasting agencies worldwide, making it an ideal platform to test our gradient operator. We specifically used the Advanced Research WRF core (Skamarock et al., 2008).

The model domain (Fig. 1) was configured to have a 0.1° (≈ 11 km) horizontal resolution, featuring a grid of 200×200 points. Vertically, the model includes 50 levels, extending to a model top of 50 hPa. We set the model time step for the WRF forecast simulation to 30 s. Our model forecast simulations were driven by the National Centers for Environmental Prediction (NCEP) Global Forecast System (GFS) operational analysis data, with a spatial resolution of 0.25° (approximately 27 km).

Figure 1The WRF model domain at a horizontal resolution of 0.1° (≈ 11 km) with the orography and the locations of GNSS stations.

For our model physics, we used the following parameters, as listed in Table 1. The radiation parameterization included the scheme of rapid radiative transfer model for general circulation models (RRTMG) (Iacono et al., 2008), which is known for accurately and efficiently calculating long-wave and short-wave fluxes and heating rates, particularly suited for general circulation model applications. The cloud microphysics scheme was the Thompson double-moment scheme (Thompson et al., 2008), capable of predicting mixing ratios for cloud water, rain, ice, snow, and graupel.

Table 1Model physics.

Download Print Version | Download XLSX

Our planetary boundary layer scheme of choice for this simulation was the Yonsei University (YSU) scheme (Hong et al., 2010; Hong and Lim, 2006). The YSU scheme, a non-local scheme with first-order closure, incorporates counter-gradient and explicit entrainment terms in the turbulence flux equation.

We utilized the unified Noah land surface model (Chen and Dudhia, 2001) for this study. This model comprises four layers and predicts soil temperature and moisture, canopy moisture, and snow cover. It considers various factors, including root zone dynamics, evapotranspiration, soil drainage, runoff variables, vegetation categories, and soil texture. This comprehensive approach provides valuable information on sensible and latent heat fluxes related to the boundary layer and incorporates an enhanced urban treatment.

To accurately simulate the model at a non-convective-scale resolution, it is essential to include convection parameterization. This helps represent the statistical impact of sub-grid-scale convective clouds. For this purpose, we employed the Grell–Freitas ensemble scheme (Grell and Freitas, 2014), which combines a probability density function with DA techniques.

3.2 Data assimilation system

The DA systems in WRF are typically classified into three categories: deterministic, probabilistic, and hybrid systems, which combine elements from both. In this study, we utilized the deterministic three-dimensional variational (3D-Var) DA system.

In the 3D-Var DA system, the objective is to iteratively minimize the cost function, J(x), where the independent or control variable is the analysis state vector, x. The cost function for the 3D-Var system is represented as follows:

The cost function, J(x), consists of two main terms: a background term and an observation term. The variables x, x_b, and y are column vectors representing the analysis state, the background (or first guess), and the observation state, respectively. The forward operator, denoted as H, is responsible for mapping the analysis state vector space to the observation vector space. In this study, we have implemented the forward operator, H, specifically for GNSS tropospheric gradients.

Alongside the column vectors, there are two square matrices that play a crucial role in minimizing the cost function: the background error covariance matrix, B, and the observation error covariance matrix, R. R is a diagonal matrix because we assume that observation errors from various sources are uncorrelated. On the other hand, B is a square, positive semi-definite, and symmetric matrix with positive eigenvalues. It includes variances in background forecast errors in the diagonal and covariances between them on the symmetric upper and lower triangular elements. After assimilating an observation, the variances and covariances in B significantly influence the analysis response. Therefore, accurately determining B is essential in a variational DA system.

In this research, we computed the B matrix using the National Meteorological Center (NMC) method (Parrish and Derber, 1992). The NMC method is widely used for generating B by estimating climatological background error covariances. We selected the NMC method due to its ability to yield physically sound results within regional model domains and its lower computational cost compared to ensemble methods. The NMC method involves calculating forecast difference statistics to obtain the forecast error covariance tailored to a specific domain. However, it does have limitations, such as overestimating covariances in large-scale simulations and regions with poor observation (Berre, 2000; Fischer, 2013; Berre et al., 2006).

For our regional simulations, forecast statistics were derived from analyzing forecast differences over a month, utilizing both 24 and 12 h predictions. These statistics were obtained from data in June 2021. We chose the CV5 option as it allows independent control of moisture levels without interference from other variables. CV5 is a version of the background error covariance matrix used in WRF, which used five control variables such as stream function (Ψ), unbalanced velocity potential (χ_u), unbalanced temperature (T_u), pseudo-relative humidity (RH_s), and unbalanced surface pressure (P_s,u). The pseudo-relative humidity is expressed as $Q/Q_{\mathrm{b},\mathrm{s}}$, where Q_b,s represents the saturated specific humidity of the background field.

3.3 Experimental setup

The assimilation experiment was conducted for 2 months, in June and July 2021, using a rapid update cycle (RUC) approach, as illustrated in Fig. 2. The RUC was configured for 6-hourly DA cycles. The datasets employed for assimilation included: (1) conventional data, comprising surface stations (SYNOP), radiosondes, and Tropospheric Airborne Meteorological Data Reporting (TAMDAR) observations; (2) ZTDs; and (3) tropospheric gradients.

Figure 2Schematic of the 3D-Var rapid update cycle initialized from the GFS analysis. A spin-up of 12 h was performed until 00:00 UTC on 1 June 2021. Four experiments with different setups are performed: a control run (black) assimilating conventional data, a ZTD run (purple) assimilating ZTDs on top of the control run, a ZTDGRA run (red) assimilating ZTD and gradients on top of the control run, and a GRA run (green) assimilating gradients on top of the control run. The WRFDA namelist switch use_gpsgraobs has a value of 0 for ZTD assimilation and 1 for ZTD and gradient assimilation.

Download

We performed four primary experiments for this study: (1) a control run, incorporating only conventional observations (SYNOP, radiosondes, and TAMDAR); (2) an impact run (ZTD run), assimilating ZTD observations in addition to the control run; (3) an impact gradient run (ZTDGRA run), where both ZTD and gradient observations were assimilated on top of the control run, using the newly developed gradient operator; and, finally, (4) a gradient run (GRA run), assimilating only gradients on top of the control run. The gradient-only assimilation was performed by down-weighting ZTDs by assigning a large observation error. The assimilation period spanned from 1 June 2021 at 00:00 UTC to 31 July 2021 at 18:00 UTC, with 6 h intervals, totaling 244 DA cycles. The initial DA cycle commenced after a 12 h spin-up run aiming to stabilize the model's initial and boundary conditions, ensuring reliable forecasts for subsequent DA.

3.4 Datasets for assimilation

3.4.1 Conventional datasets

To enhance the capabilities of the DA system, we established a comprehensive network of surface reports, SYNOP, across Europe. Radiosonde measurements, obtained through TEMP, provide a detailed view of the atmospheric thermodynamic structure at launch points. TEMP is a collection of alphanumerical codes established by the WMO. These codes represent upper-air soundings obtained through weather balloons launched from either land or sea level. These observations report on the weather conditions in the upper regions of the atmosphere. In WRFDA, FM-35 TEMP is the observation code used to identify radiosonde observations launched from land and FM-36 for those launched from ships. In order to compensate for the underrepresented radiosonde network during specific time periods, such as 06:00 and 18:00 UTC rather than 00:00 and 12:00 UTC, we utilized a series of TAMDAR observations. Table 2 summarizes the average number of observations assimilated per time step over 2 months. These valuable observations are conveniently accessible via the WMO's Global Telecommunication System data archive, which is housed at ECMWF.

Table 2Average number of observations assimilated per time step for June and July 2021.

Download Print Version | Download XLSX

To maintain simplicity within the DA system, we limited the use of conventional datasets to surface observations, radiosondes, and TAMDAR observations. It is worth noting that the primary objectives of this work are as follows: (1) to test the functionality of the newly developed operator (code) and (2) to analyze the difference between the experiment where we assimilate ZTDs only and the experiment where we assimilate both ZTDs and tropospheric gradients. In essence, our focus lies in assessing the relative impact rather than the absolute impact of GNSS data in variational DA.

4.1 Single observation tests

To assess the impact of assimilating a single observation from a station location, we conducted the SOT in WRF. Through the SOT, our aim was to gain an insight into the model's behavior in response to GNSS gradient observations. The DA system employed here is the 3D-Var, where the extent of the assimilation impact largely depends on the B matrix. We sought answers to the following questions:

1. How does the impact region differ when compared to assimilating ZTD observations?

2. What distinguishes the SOTs resulting from the assimilation of gradients from those involving ZTD observations?

3. How sensitive is the model to the assimilation of a single gradient observation?

To address these questions, we selected a point at the center of the model domain and conducted SOTs with varying pseudo-ZTD and pseudo-gradient observations and their associated errors.

For questions 1 and 2, we conducted separate SOTs for gradients and ZTDs. In the gradient SOT, we selected an increment of −1.0 mm for the north gradient observation value and 0 mm for the east gradient observation value, with an observation error of 1.0 mm. In the ZTD SOT, we used an increment of 1.0 cm for the ZTD value and an observation error of 1.0 cm. We opted for unit values as errors to facilitate the understanding of the increment's impact without a scaling factor.

Figure 3 provides a comparison of the spatial impact resulting from SOTs using gradient and ZTD observations. This comparison helps us visualize how the impact spreads. The gradient impact plot in Fig. 3a justifies the name gradient as it reveals increased moisture in the south (positive lobe) and decreased moisture in the north (negative lobe), indicating a redistribution of moisture in the analysis. Visualizing the gradient increment as a vector pointing south aligns with the input values of −1.0 mm for the north gradient and 0 mm for the east gradient. Tropospheric gradients can be considered moisture vectors, similarly to wind vectors, encompassing both north–south and east–west components. Typically, tropospheric gradients point from dry to moist areas.

Figure 3Single observation test spatial plot of gradient and ZTD. An N gradient increment of −1.0 mm and an error of 1 mm are applied for the gradient SOT on the left. A ZTD increment of 1 cm and an error of 1 cm are applied on the right.

Comparing the spatial impact of SOTs, we observe that the maximum impact response for a −1 mm north gradient increment with a 1 mm error is 0.062 g kg⁻¹. In contrast, the maximum response for the SOT due to a 1 cm magnitude ZTD increment with a 1 cm error is 0.2 g kg⁻¹. The impact radius reduces by 50 % from 0.062 to 0.032 g kg⁻¹ within a radius of 67 km for the gradient and around 80 km for the ZTD.

Figure 4 displays the water vapor mixing ratio (WVMR) and temperature profiles of the SOT analysis with respect to model levels. The profiles were drawn from the locations of the positive and negative lobes. In order to find out the locations of the lobes, the following steps were adopted:

1. We average the analysis over the vertical levels between 2 and 4 km, where the gradient observations have the maximum influence.

2. We determine the latitude and longitude of the point in the domain where the absolute value of the water vapor mixing ratio is maximum.

3. Through step 2, we get two maximum value points on both sides of the gradient observation location. The profile is derived from these two locations; hence, we get the lobes.

The gradient SOT profile exhibits both positive and negative lobes, explaining the positive and negative responses in the WVMR profile. The temperature response in both gradient and ZTD SOTs is negligible. It is evident that the impact, to a large extent, depends on the B matrix, particularly with the CV5 option, which assumes moisture independence from other variables.

Figure 4The vertical profile of the SOTs with the same settings as in Fig. 3. The gradient SOT on the left shows the positive and negative (dash) lobes of the impact.

Download

Figure 5The vertical cross section of the SOTS with the same settings as in Fig. 3.

Download

The profile plots highlight that gradient observations have a more significant impact on the lower troposphere than on the surface level, where the impact is minimal. In contrast, ZTD impacts are more pronounced in the lower troposphere but significant at the surface level as well. A cross section in Fig. 5 further emphasizes this difference, showing that ZTD impact extends from the lower troposphere to the surface with only a slight decrease in the WVMR value compared to the gradient impact. The gradient impact is most prominent between 1 and 5.5 km, with minimal influence at the surface level.

To gain a more realistic understanding of the impact of tropospheric gradients, we conducted an SOT using an actual gradient observation from a GNSS station. We selected a station near the center of the model domain with the coordinates $\mathrm{52}\mathit{°}{\mathrm{38}}^{\prime }{\mathrm{35}}^{\prime \prime }$ N and $\mathrm{9}\mathit{°}{\mathrm{12}}^{\prime }{\mathrm{22}}^{\prime \prime }$ E. The observation values were 0.497 mm for the east gradient, 0.099 mm for the north gradient, with an observation error of 0.65 mm assigned. Unlike previous SOTs that focused on one gradient direction, this scenario involved increments in both the east and the north gradient components, resembling real scenarios. In this case, the gradient dipole's direction is the vector sum of the two gradient components. Figure 6 displays the spatial and profile plots of this SOT.

Figure 6A real observation SOT with gradient observations: the east gradient component equals 0.497 mm, and the north gradient component equals 0.099 mm. The observation error value for the SOT is chosen to be 0.65 mm, which is the standard error used for this study.

We modified the WRFDA code to allow users to conduct SOTs by adding specific options to the WRF namelist.input file.

4.2 Qualitative analysis of the gradient assimilation impact

This section aims to assess the noticeable impact of gradient observations on assimilation, specifically analyzing the spatio-temporal features emerging in the model as a result of assimilating gradients. To illustrate the influence of gradients, we compared the analyses obtained from the ZTDGRA and ZTD experiments using difference plots. The analysis is averaged over vertical levels between 2 and 4 km to capture the comprehensive impact from the lower troposphere, where gradient observations have the most significant influence.

Assimilation was conducted every 6 h in a RUC DA environment. In a RUC DA environment, the initial assimilation impact is of utmost importance as it reflects the analysis increment resulting from the assimilation of fresh observations, localized to their respective locations. As the cycles progress, the features tend to expand on a larger scale, making visual analysis more challenging. Therefore, our focus is primarily on the first and second DA cycles.

The initial DA cycle occurred on 1 June 2021, at 00:00 UTC. Following this, a 6 h free forecast was initiated based on the analysis until the next assimilation at 06:00 UTC. Figure 7 illustrates the initial DA cycle, with Fig. 7a showing the analysis difference observed at 00:00 UTC, while Fig. 7b–f display the forecast differences from 01:00 to 05:00 UTC. By comparing Fig. 7a with the station network in Fig. 1, it is evident that data from GNSS receivers in the region have had a substantial impact. Approximately 100 stations located exclusively in Germany have contributed to this impact. The absolute magnitude of the impact ranges from −0.15 to 0.14 g kg⁻¹.

Figure 7Spatial analysis and forecast differences: ZTDGRA minus the ZTD experiment. The figure shows the exclusive impact of the gradient observation for the first DA cycle at 00:00 UTC (analysis) and the 5 h forecast from the analysis.

While the impact of assimilating gradients may seem relatively small in magnitude, it significantly affects the distribution of the moisture field in the initial model state (Fig. 7a). While these figures provide valuable insight, a quantitative analysis is necessary to determine if the assimilations have successfully corrected the model's moisture fields. Subsequent sections will delve into this topic, providing detailed comparisons for clarity.

Intercomparing Fig. 7b–f reveals that the impact of the gradients persists in the model for a significant period (6 h; until the next assimilation). This persistent impact serves as compelling evidence of the reliability of gradient data. Figure 8 reinforces this statement by displaying the analysis of the second DA cycle and the subsequent 6 h forecast difference. The structures developed during the forecast difference at 05:00 UTC are comparable to those observed in the analysis difference at 06:00 UTC during the second assimilation, as demonstrated by the comparison of Fig. 7f with Fig. 8a. The structures that emerged during the second assimilation were not significantly different from the forecast made before the assimilation, indicating that the gradient observations were accurate and that the assimilation did not degrade the model state nor did it tamper with the flow of the moisture fields.

Figure 8Same as Fig. 7 but for the second DA cycle.

4.3 Quantitative analysis of the gradient assimilation impact

To comprehensively assess the impact of assimilating GNSS observations, we conducted a quantitative analysis and compared the results with the original station data from the GNSS network. Our objective was to gauge the extent of improvement achieved through this assimilation process. Our observation network initially comprised just over 100 GNSS stations. However, for validation purposes, we deliberately excluded 18 stations, categorizing them as “excluded” stations.

These excluded stations were strategically selected to ensure a balanced spatial distribution, aligning them with the locations of the German Weather Service (DWD) radar stations. They were set aside solely for independent validation purposes. The data from the remaining GNSS stations, known as the allowed stations, were assimilated in the ZTD, GRA, and ZTDGRA experiments. Consequently, we could evaluate the progress by comparing the assimilated data with the independent observations from these allowed stations.

This approach allowed us to rigorously assess the impact of our assimilation efforts and validate the improvements achieved. We will go through the Root Mean Square Error (RMSE) impact with respect to the allowed and excluded stations and also look into the average impact (in terms of RMSE) as a function of the forecast length in between the assimilation intervals.

4.3.1 RMSE with respect to allowed stations

To assess the impact of assimilating ZTD and gradient observations on the analyses in observation space, we computed the RMSE for four experiments: (1) control, (2) ZTD, (3) ZTDGRA, and (4) GRA using data from the allowed GNSS stations. We compared the ZTD and gradient values derived from the analyses at the station coordinates with the observations collected by GNSS stations. This comparison was performed over the entire 2-month period, using hourly data and considering all allowed stations. Figure 9 presents the station-specific RMSEs for the control, ZTD, ZTDGRA, and GRA experiments. The mean RMSE is listed in Table 3.

Figure 9The station-specific RMSE of the ZTD, north and east components (allowed stations): control (black), GRA (green), ZTD (purple), and ZTDGRA (red).

Download

Table 3Comparison of the mean μ in mm of station-specific RMSE of ZTDs and gradients (W for allowed and B for excluded stations).

Download Print Version | Download XLSX

From the figure, it is evident that when assessing the improvement in the ZTD variable, the mean RMSE values for the ZTD parameter are the lowest for the ZTDGRA experiment compared to the other runs, providing clear evidence of the successful gradient assimilation impact. The mean RMSE of the ZTD variable for the control run was 14.4 mm, which was reduced to 9.7 mm in the ZTD run and further decreased to 9.3 mm in the ZTDGRA run. The GRA experiment also had an impact on the ZTD analysis. The GRA RMSE of 12.2 mm was in between that of the control and the ZTD experiments.

Regarding the impact on gradient components, both the north and east components exhibit similar enhancements. Notably, the ZTDGRA assimilation leads to the most significant improvement when compared to the other runs. The north (east) gradient RMSE decreased from 0.68 mm (0.69 mm) in the control run to 0.61 mm (0.62 mm) in the ZTD run and further decreased to 0.56 mm (0.56 mm) in the ZTDGRA run. From the GRA experiment, the impact of gradients was clearly evident. There was a significant reduction in RMSE from 0.68 mm (0.69 mm) in the control to 0.58 mm (0.57 mm) in the GRA for the north (east) gradients. The results are almost equal to those of the ZTDGRA run. The assimilation of gradients significantly reduced the RMSE values, indicating a substantial enhancement of the moisture field in the model state.

The results shown in Fig. 9 can be interpreted as follows: when assimilating only ZTDs, not only are the ZTDs adjusted, but the tropospheric gradients are also affected. Also, the assimilation of the gradients alone improves not only gradient analyses, but also the ZTDs. This suggests that tropospheric gradients contain valuable information. Such adjustment would not be possible if tropospheric gradients contained no useful data. However, as long as we do not assimilate tropospheric gradients, we do not make use of their information. This is addressed in the experiment, where both the ZTDs and the tropospheric gradients are assimilated. It can be observed that the tropospheric gradients are further adjusted, demonstrating the functionality of our implementation and that the ZTDs are also further adjusted, providing evidence that we extracted valuable information from the tropospheric gradients. Also, through the GRA run, it is clear that ZTDs do not provide weightage to improve the gradients. However, the ZTD plays a role in improving the gradient and vice versa.

4.3.2 RMSE with respect to excluded stations

To further evaluate the consistency of our findings from the allowed stations, we conducted a similar analysis using the 18 excluded stations. We measured the RMSE values as described previously but specifically for these 18 stations. This comparison was performed using hourly data over the entire 2-month period to provide a robust statistical assessment. Figure 10 presents the RMSE comparisons for the control, ZTD, ZTDGRA, and GRA experiments, using data from the excluded stations.

Figure 10The station specific RMSE of the ZTD, north and east components (excluded stations): control (black), GRA (green), ZTD (purple), and ZTDGRA (red).

Download

In the figure, it is evident that the RMSE values for the ZTD variable have decreased. The mean RMSE for the ZTD variable decreased from 14.2 mm in the control run to 10.2 mm in the ZTD run and was further reduced to 9.7 mm in the ZTDGRA run. The GRA run also contributed to a reduction in RMSE from 14.2 mm in the control to 12.2 mm. Table 3 lists the mean RMSE for all the runs.

Similar improvements in the north and east gradient components were observed, mirroring the results obtained with the allowed stations. The RMSE values for the north (east) gradient component dropped from 0.68 mm (0.68 mm) in the control run to 0.62 mm (0.61 mm) in the ZTD run and further decreased to 0.58 mm (0.57 mm) in the ZTDGRA run. In the GRA experiment, it was observed that gradients had a clear impact, similarly to the comparisons of allowed stations. There was a significant decrease in RMSE from 0.68 mm in the control to 0.58 mm in the north and 0.57 mm in the east gradients. These results were nearly equivalent to those of the ZTDGRA run.

These results clearly demonstrate that the assimilation of gradients has enhanced the model's analyses. The 2-month statistics using independent GNSS station data validate the improvements achieved through gradient assimilation.

4.3.3 Average impact as a function of forecast length

As our assimilation process occurs every 6 h, we inherently have a 6 h forecast lead time. To calculate the mean RMSE for both ZTD and gradients, while considering both allowed and excluded stations, we examined the mean RMSE values across forecast lead times ranging from 0 to 5 h, similarly to our previous analysis. Figure 11 illustrates the RMSE variation with respect to forecast length.

Figure 11Average impact with respect to forecast lead times from analyses for 2 months of simulation. The control run (black), GRA run (green), ZTD run (purple), and ZTDGRA run (red) are shown for the allowed stations (a, b, c) and for the excluded stations (d, e, f).

Download

In both allowed and excluded stations, the ZTDGRA experiment consistently outperforms the other runs, displaying the lowest RMSE values. As expected, the RMSE tends to increase as the lead time extends. However, what is noteworthy is the sustained improvement in RMSE across different forecast lead times in both allowed and excluded stations. This consistency of the improvement provides further evidence of the positive impact of gradients on the model state's moisture field. It is more evident in the GRA run, where the gradient forecast RMSE outperforms the ZTD. This is a clear sign that there is no weightage from the ZTD observations for improving the gradients. Additionally, the ZTDGRA experiment maintains a stable RMSE improvement throughout the forecast lead time, indicated by the small but consistent offset.

4.4 Comparison with ERA5

In addition to GNSS observation data, we conducted validation using the ERA5 dataset, which is considered one of the most comprehensive atmospheric reanalyses of the global climate to date, produced by the Copernicus Climate Change Service at ECMWF. ERA5 provides hourly estimates of climate variables globally on a 30 km grid, offering detailed atmospheric information with 137 levels up to 80 km in altitude. Given its high quality, ERA5 serves as a valuable reference dataset. To perform this validation, we selected five locations within the model domain, as shown in Fig. 12. With a total of 244 assimilation cycles, each generating five profiles, we had a substantial dataset of 1220 profiles for comparison with ERA5. This large number of profiles allowed us to obtain robust statistical results.

Figure 12ERA5 comparison domain. The station map and the selected points for the comparison of the ERA5 data with the analyses.

We compared the water vapor (WV) profiles for three different lead times after the assimilation cycles: 0, 3, and 5 h. Figure 13 presents the WV 1-sigma deviation between the assimilation experiments and ERA5 as a function of altitude for these three lead times. The findings from Fig. 13 indicate that assimilating ZTDs has a positive impact across the entire altitude range. This confirms that the assimilation of ZTDs results in an improvement in the WV field.

Figure 13ERA5 profile comparison. The statistics of 1220 profiles for the control run (black), GRA (green), ZTD run (purple), and ZTDGRA run (red) are shown for the analyses (00:00 UTC) in the leftmost column, 3 h time lead in the middle column, and 5 h time lead in the rightmost column. The bottom row shows the magnified lower-troposphere profiles.

Download

Furthermore, the assimilation of tropospheric gradients, in addition to ZTDs, leads to further improvements in the WV field. However, it is important to note that this improvement is primarily restricted to altitudes above 2.5 km. Below an altitude of 2.5 km, particularly near the Earth's surface, the assimilation of gradients, in addition to ZTDs, does not have a significant impact. This lack of impact close to the Earth's surface is not surprising, as gradients are not sensitive to surface variables (as indicated by Eqs. 4 and 5). Any impact near the surface can be attributed to the background error covariance matrix and its correlations between different altitudes, as observed in the SOT experiment. It is possible to gain more insight into the impact of observations by differentiating the ZTD and GRA profiles. GRA has more influence on the lower troposphere but not on the surface, whereas the ZTD has a greater impact on the surface level.

The positive impact of gradient assimilation remains consistent across different lead times. In essence, the reduction in the standard deviation above 2.5 km observed in the ZTDGRA experiment remains consistent even at the 5 h lead time.

4.5 Forecast comparison with radiosonde profiles

Forecast validation was carried out for the analysis forecasts generated in all four experiments. The analyses valid at 00:00 UTC from each of these experiments were subjected to a free forecast of 12 h in length. The 12 h forecasts were prepared starting from the DA cycle at 00:00 UTC on 1 June 2021, and continuing until 00:00 UTC on 30 July 2021. Figure 14 shows the average RMSE of the mixing ratio up through the radiosonde profile for 60 d. Since radiosonde profiles are usually available at irregular heights, we interpolated the model profiles to those of the radiosondes. Hence, the model profiles of the 12 h forecast from the analyses were compared with the radiosonde profiles at that time instant to calculate the RMSE. We compared the forecasts valid at 12:00 UTC (forecasted from analyses at 00:00 UTC) with the radiosondes at 12:00 UTC, since there were more radiosondes during that time, to get the statistics. There were a minimum of 38 radiosondes available for comparison.

Figure 14Water vapor mixing ratio comparison of the 12 h forecasts from the analyses valid at 00:00 UTC to the radiosonde profiles at 12:00 UTC (around 38 radiosondes per epoch). The mean RMSEs (g kg⁻¹) of the respective runs for the 2 months are specified in the legend.

Download

While the improvement observed in the ZTDGRA experiment is notable, it remains relatively modest. It is important to note that this improvement may not persist beyond the 12th hour of forecasting. Additionally, any non-uniformity among the radiosonde profiles could potentially impact the final assessment of these improvements. However, from the mean RMSE, we understand that incorporating the gradient observations improves the forecasts. The mean RMSEs in descending order are control, GRA, ZTD, and, finally, ZTDGRA.