4DVAR assimilation of GNSS zenith path delays and precipitable water into a numerical weather prediction model WRF

The GNSS data assimilation is currently widely discussed in the literature with respect to the various applications for meteorology and numerical weather models. Data assimilation combines atmospheric measurements with knowledge of atmospheric behavior as codified in computer models. With this approach, the “best” estimate of current conditions consistent with both information sources is produced. Some approaches also allow assimilating the non-prognostic variables, including remote sensing data from radar or GNSS (global navigation satellite system). These techniques are named variational data assimilation schemes and are based on a minimization of the cost function, which contains the differences between the model state (background) and the observations. The variational assimilation is the first choice for data assimilation in the weather forecast centers, however, current research is consequently looking into use of an iterative, filtering approach such as an extended Kalman filter (EKF).

This paper shows the results of assimilation of the GNSS data into numerical weather prediction (NWP) model WRF (Weather Research and Forecasting). The WRF model offers two different variational approaches: 3DVAR and 4DVAR, both available through the WRF data assimilation (WRFDA) package. The WRFDA assimilation procedure was modified to correct for bias and observation errors. We assimilated the zenith total delay (ZTD), precipitable water (PW), radiosonde (RS) and surface synoptic observations (SYNOP) using a 4DVAR assimilation scheme. Three experiments have been performed: (1) assimilation of PW and ZTD for May and June 2013, (2) assimilation of PW alone; PW, with RS and SYNOP; ZTD alone; and finally ZTD, with RS and SYNOP for 5–23 May 2013, and (3) assimilation of PW or ZTD during severe weather events in June 2013. Once the initial conditions were established, the forecast was run for 24 h.

The major conclusion of this study is that for all analyzed cases, there are two parameters significantly changed once GNSS data are assimilated in the WRF model using GPSPW operator and these are moisture fields and rain. The GNSS observations improves forecast in the first 24 h, with the strongest impact starting from a 9 h lead time. The relative humidity forecast in a vertical profile after assimilation of ZTD shows an over 20 % decrease of mean error starting from 2.5 km upward. Assimilation of PW alone does not bring such a spectacular improvement. However, combination of PW, SYNOP and radiosonde improves distribution of humidity in the vertical profile by maximum of 12 %. In the three analyzed severe weather cases PW always improved the rain forecast and ZTD always reduced the humidity field bias. Binary rain analysis shows that GNSS parameters have significant impact on the rain forecast in the class above 1 mm h⁻¹.

1 Introduction

The data assimilation in weather forecasts is one of the key components in all prediction systems as it is an initial value problem and the quality of the initial field has a large impact on the forecasts. Currently, the leading weather agencies assimilate operationally dozens of observation data types such as radiosonde (RS) profiles, radiances from satellite observations, SYNOPs, refractivities from radio occultations, pilot reports and many others (Barker et al., 2004). With the advent of European Cooperation in Science & Technology (COST) actions 716 (1999–2004), 1206 (2013–2017), as well as the project funded in the fifth framework program “Targeting Optimal Use of GPS Humidity Measurements in Meteorology” (TOUGH), the adoption of the ground based global navigation satellite systems (GNSS) observations to the operational forecasts by most of the weather services in Europe become a fact. In this study the term GNSS covers all navigation systems used world-wide, whereas the term Global Positioning System (GPS) is related to only one source of observations – the US based GPS. There are many publications related to either (1) performance of large-scale weather forecast systems augmented with many observations including GNSS, (2) added value of GNSS observations in nowcasting services, or (3) case-based studies showing the impact of GNSS data in particular cases. The following three approaches are discussed below.

A very comprehensive study done by Poli et al. (2007) on the global forecast model Arpage using four dimensional variational assimilation (4DVAR) shows that the impact of GPS zenith total delay (ZTD) on forecasts is different in winter (improving pressure), spring (reducing surface humidity root mean square error) and summer (positive impact on wind, geopotential and precipitation, negative on humidity). A similar, very detailed study was done by Bennitt and Jupp (2012), where the authors discussed the operational assimilation of GPS ZTD in MetOffice into the North Atlantic and European 12 and 24 km model in spring, summer and autumn. The results were mixed: for all cases the introduction of GPS ZTD increased the humidity bias, however the improvements of clouds forecasts were observed. Bennitt and Jupp (2012) also identify no clear benefit of 4DVAR against 3DVAR. Lindskog et al. (2017) in their study in the Nordic countries of GNSS ZTD impact on forecasts, confirmed that the forecasts are sensitive to thinning distance. Shorter distance between stations (below 100 km) leads to a larger humidity bias in the lower troposphere, which may explain the humidity bias in the Bennitt and Jupp (2012) solution. Lindskog et al. (2017) showed that the humidity forecast is better when the GNSS ZTD is assimilated with other meteorological observations such as the Advanced Microwave Sounding Unit (AMSU) or Infrared Atmospheric Sounding Interferometer (IASI) radiances. The authors also showed that the adopted bias correction strategy and GNSS ZTD estimation procedure have marginal impact on the forecasts. All studies run in large weather forecasting systems suggests that the assimilation of GNSS ZTD, either 3D or 4DVAR, on average has a mostly neutral impact on the forecast if the system is already saturated with meteorological observations.

Another branch of weather models are those used in nowcasting, such as the (legacy) rapid update cycle RUC 20 km and RUC 40 km (Benjamin et al., 2002) or currently operational rapid refresh (RAP) (Benjamin et al., 2016) that are targeting short 12 and 24 h predictions for decision making and safety operations with a large number of observations assimilated every hour. One of the first experiments using GPS precipitable water (PW) in nowcasting service in the USA (RUC model 60 km) (Smith et al., 2000), showed a 1 % improvement of the relative humidity forecast in the bottom part of the atmosphere. However, in specific cases related to active frontal weather, the improvement was much larger: 14 % in the moistening and 24 % in the drying stage of the advection. The increased spatial resolution to 20 km of RUC20 (Smith et al., 2007) shows stronger improvements in humidity field and convective available potential energy, CAPE, than with RUC40. The 850 hPa relative humidity (RH) forecasts improve more in the nighttime, and in the colder season than that in the warmer season. The current RAP model, running on a 13km grid, continuously assimilates GPS PW every hour from 300 stations across the US (Benjamin et al., 2016). It shows that there is clear benefit in using GPS observations, especially for short term (nowcasting) predictions.

The third type of studies that are appearing in the literature are case based, showing the impact of GNSS on particular weather events. One of the first to test the impact of GPS based ZTD observations in Europe was Cucurull et al. (2004). The authors used the NCAR/Penn State Mesoscale Model 5 (MM5) ZTD 3DVAR assimilation for a case of a snow storm on 14–15 December 2001 over the western Mediterranean Sea. They found that there are reductions of root mean square error (RMSE) of wind by 1.7 %, temperature by 4.1 % and surface humidity by 17.8 %. The authors also noted that the forecasts work better if the ground based automatic weather stations were used in the same assimilation run. Another example of an early stage case-based research is the assimilation of GPS PW by Nakamura et al. (2004) with a 4DVAR scheme into a mesoscale Japanese Meterological Agency (JMA) model for summer intensive rain cases. The assimilation of GPS data improved the precipitation location, but the statistics did not show large improvement. One of the first GPS 4DVAR ZTD studies in the US was by De Pondeca and Zou (2001), who ran assimilation of GPS observations in MM5 together with the wind profiler data and radio acoustic sounding system (RASS) virtual temperature. Five 12 h experiments for California's December frontal system passage were performed. It was found that the ZTD assimilation corrects the underestimation of accumulated rain by 33.15 % and 25.08 % for 6 and 12 h respectively. Adding the wind profiler improves the forecast by 88.26 % and 32.53 % and adding further RASS observations increases the performance to 93.21 % and 50.58 %, respectively. In a more recent study by Boniface et al. (2009), the GPS ZTD was assimilated (3DVAR) for 280 stations over 15 days into the high-resolution (2.5 km) AROME model. The results were positive for poorly predicted precipitation and neutral for well predicted ones. More recently, Tilev-Tanriover and Kahraman (2014) studied the impact of the GPS PW assimilation in the Weather Research and Forecasting (WRF) model in a 2 day case of intense snowfall forecasts in the central Anatolia. The authors performed three experiments: base run, cold start and cycling all with PW 3DVAR operator. Results show that the cycling assimilation mode decreases the temperature and humidity biases, whereas the cold start performs worse than the control run. Saito et al. (2016) studied the impact of ensemble prediction that did not produce enough precipitation. They found that even downscaling from 10 to 2 km still does not improve locations of precipitation's cores. Finally, the 4DVAR PW assimilation into a non-hydrostatic model improved the location of scattered intensive rain. In summary, most of the literature reported a substantial increase in the quality of the forecast of humidity, rain location and sometimes also the rain accumulated total amount. Less significant improvement was achieved for wind speed and temperature. All studies that used additional observations, especially these resolving vertical structure of the water vapor and temperature, complemented GNSS observations and improved the forecast even more.

The literature review shows that the impact of the GNSS ZTD/PW assimilation depends on the number of already assimilated observations and applied preprocessing (Bennitt and Jupp, 2012; Lindskog et al., 2017; Poli et al., 2007) as well as on the type of weather conditions. The main aim of this paper is to quantify the impact of the GNSS data, both ZTD and PW, gathered operationally in Poland, in weather forecasting. The study is based on the WRF model with a high spatial resolution of 4 km × 4 km supported with the WRF data assimilation (WRFDA) package. We show the importance of the GNSS data assimilation for cases of various meteorological conditions observed in May and June 2013, which is a benchmark period for COST Action ES1206. To our best knowledge, no GNSS ZTD/PW assimilation experiment was carried out in Poland yet. Moreover, we found only one publication (Tilev-Tanriover and Kahraman, 2014) dealing with the assimilation of GPS ZTD and PW into widely adopted WRF model using the WRFDA package. We present a study showing the impact of GNSS ZTD and PW observations on the forecasts for a longer time period – 2 months (May 2013 – calm weather conditions and June 2013 – active stormy weather), followed by quantifying improvements of adding RS and SYNOP data into the assimilation system already run with GNSS observations, and finally we verified impact of GNSS observations on prediction for specific cases.

The paper has following structure: after the introductory section a short overview of used data and methodology is presented in Sects. 2 and 3, respectively. These sections are followed by the experiment setup description and results (Sect. 4). The paper is closed with a conclusion section.

2 Data

The GNSS PW/ZTD data are assimilated into the NWP WRF model. The chosen period covers May and June 2013, with a special focus on 5–23 May 2013 and three shorter cases: (a) 29–31 May, (b) 17–19 June and (c) 24–26 June 2013. The period is chosen in accordance with the COST Action ES1206 GNSS meteorology benchmark (Douša et al., 2016).

2.1 WRF model

In this study the WRF model is used, it is numerical weather prediction system designated for simulation of multiscale, spatial and temporal atmosphere flows. The WRF configuration (Kryza et al., 2013) is based on two nested model domains. The first domain covers the European area with a 12 km × 12 km grid spacing. The second, nested domain covers Poland and Central Europe with a 4 km × 4 km grid spacing (Fig. 1).

Figure 1 The WRF model configuration (the inner square represents the nested domain with 4 km × 4 km resolution). Colors denote the orography of the terrain (m a.s.l.).

Table 1 The WRF configuration used in the experiment.

Initial and boundary conditions are taken from the National Center for Environmental Prediction Final Analysis, Operational Model Global Tropospheric Analyses (NCEP FNL) database (National Centers for Environmental Prediction, 2000). The data are available with 1^∘ × 1^∘ horizontal and 6 h temporal resolution and with 26 vertical levels from 1000 to 10 hPa. The WRF model for Poland is calculated and provided by the Department of Climatology and Atmosphere Protection of the University of Wrocław. The details of the WRF configuration are presented in Table 1. Data assimilation was run using a 4DVAR WRF DA system, only for inner domain (d02). Prediction model was started once a day, at 00:00 UTC. Assimilation window was centered around 00:00 UTC. Background error covariance (BE) was selected for regional application (cv_options =5) (BE depends on the WRF domain). BE was constructed based on a forecast for convection events in the first week of May 2013.

Quality control was selected for SYNOP and RS data in observation processor (obsproc) and in WRFDA. For ZTD and PW data, quality control was conducted before processing, followed by obsproc verification and the last step in 4DVAR assimilation. The WRF configuration was based on the best ensemble members (ens1) with a small modification from ensemble system dedicated to Poland during summertime (Guzikowski et al., 2015).

2.2 GNSS data

The GNSS data are calculated by the GNSS and Meteo working group from the Institute of Geodesy and Geoinformatics, Wrocław University of Environmental and Life Sciences (http://www.igig.up.wroc.pl/igg, last access: 17 August 2018). The PW and ZTD values are calculated at 106 stations of the European Position Determination System Active Geodetic Network (ASG-EUPOS, http://www.asgeupos.pl, last access: 17 August 2018) in Poland and adjacent areas (Fig. 2).

Figure 2 Location of the GNSS, SYNOP and radiosonde stations in Poland.

The GNSS parameters are calculated from GPS data only, using the Bernese GNSS Software version 5.0 (Dach et al., 2007). The parameters (coordinates and troposphere) are estimated in a near-real time (NRT) regime, 30 min after each full hour, without the gradient estimation. The dry troposphere a priori model is taken from Saastamoinen (1972), mapped with dry Niell MF (Niell, 2000) and the ZTD relative constraining of 3 mm is applied. International GNSS Service (IGS) ultra-rapid orbits, clocks and Earth rotation parameters are used. These parameters are now altered to fit the more recent version of Bernese (5.2) (Dach et al., 2015; Dymarska et al., 2017), but this study uses NRT data, originally processed in 2013. In this way, our impact study will show the minimum potential of GNSS data assimilation in weather model exactly. More details on the GNSS data processing and quality monitoring of the data can be found in Bosy et al. (2012). Fifteen of the stations are a part of the EUREF Permanent Network (EPN) and provide the tropospheric parameters with the accuracy required by NWP data assimilation (Dymarska et al., 2017), i.e., the standard deviation between GNSS ZTD and WRF ZTD is 10 mm and the standard deviation between radiosonde ZTD and WRF ZTD is 14 mm. In the inter-comparison study using multiple techniques (Wilgan et al., 2015), the discrepancy between GNSS observations and radiosonde was found to be 10 mm. According to the EGVAP requirements (Offiler, 2012), this accuracy of the GNSS data is sufficient for the assimilation in NWP models.

2.3 Model evaluation

The WRF model runs are compared with surface meteorological measurements of air temperature, relative humidity, wind speed and precipitation, radiosonde temperature and relative humidity profiles in three locations across Poland, GNSS observations in 106 locations.

For the SYNOP, measurements are available every hour at 66 SYNOP stations, evenly distributed over the area of Poland, operated by the Institute of Meteorology and Water Management – National Research Institute. Model evaluation is performed only for the nested domain. Four error metrics are calculated to assess the forecast performance.

• Mean error (ME), which describes the model tendency of overestimation (ME > 0) or underestimation (ME < 0) of the given meteorological parameter. The ME (bias) is calculated as a mean difference between the modeled and observed values for all stations (domain wide). The units are the same as for the analyzed meteorological parameters.

• Root mean square error (RMSE), which takes only non-negative values. The RMSE (scatter) is calculated as a root of the squared differences between the modeled and observed values for all stations. The units are the same as for the analyzed meteorological parameters.

• Pearson correlation coefficient (corr), which takes values from −1 to +1, and the expected value is 1. Corr is unitless.

• Index of agreement (IOA), developed by Willmott (1981) as a standardized measure of the degree of model prediction error. IOA varies between 0 and 1, and 1 indicates a perfect match.

Model performance verification is done using observations for rainfall, wind speed, relative humidity and air temperature at 2 m. The model evaluation is done for each simulation considering the entire period and for different lead times and for selected days during which severe weather was observed (case studies). For rainfall forecasts, binary evaluation is also presented using performance diagrams (Roebber, 2009), separately for five different precipitation intensity thresholds. The closer data set is towards the upper-right corner of the plot the better performance of the forecasts.

Additionally, the model performance for air temperature and relative humidity was compared with radiosonde data, available with a high vertical resolution for three stations located in Poland: Łeba, Warsaw and Wrocław, For the May and June case (large data set) we provide (Fig. 4) profile mean error with standard deviation multiplied by 1.96 (p=0.05), for other cases (much fewer observations) only mean errors are provided (Figs. 6, 10, 12). A similar comparison was done by Guerova et al. (2005). The model based PW and ZTD are also compared to the GNSS based retrievals. Bias and standard deviation of the residual WRF based PW and ZTD minus GNSS observed PW and ZTD are calculated for 106 stations.

3 Methodology

The variational assimilation is based on the Bayesian probability theory and it states that the model analysis is inferred from two probabilities: background and observations. These can also be expressed as a minimization of a cost function, with two major components: background B and observations R error covariance (Ide et al., 1997; Lorenc, 1986; De Pondeca and Zou, 2001) in the 4DVAR implementation:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }J\left[x\left(t_o\right)\right]={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left[x\left(t_o\right)-x^{\mathrm{b}}\left(t_o\right)\right]}^T{B}_{\mathrm{0}}^{-\mathrm{1}}\left[x\left(t_o\right)-x^{\mathrm{b}}\left(t_o\right)\right]\\ \text{(1)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\sum }_{i=\mathrm{0}}^n\left(H_i\left[x\left(t_i\right)\right]-y_i^o\right)R_i^{-\mathrm{1}}\left(H_i\left[x\left(t_i\right)\right]-y_i^o\right),\end{array}$$

where x(t_i), x^b, x(t_o), are model state vectors at the time t_i, and background vector and model initial conditions t₀, respectively. In general cases, there are N kinds of observations y defined at discrete times t_i from t_o to t_n, where the assimilation window spans from the lowest to the highest t_i. The H_i[x(t_i)] is a forward operator that transforms parameters from the model space to the observation space. The 3DVAR differs to 4DVAR by taking t_i equal to observation time and analysis time. Minimization of the Eq. (1) requires also finding adjoint (ADJ) and tangent linear (TLM) operators, each related to the observation type and forward operator H_i[x(t_i)]. For more details, the readers are referred to, e.g., Barker et al. (2004) or Huang et al. (2009).

3.1 GPSPW operator

The WRFDA system employed in this study hosts variational (VAR) 3D/4D as well as hybrid variational ensemble algorithms (Barker et al., 2012). Currently, the system supports the assimilation of surface, radiosonde, aircraft, wind profile observations, as well as atmospheric motion vectors, radar reflectivities, spectrometric, GPS radio occultation and GPS ground-based data. The latter is linked directly to the GPSPW operator (the National Center for Atmospheric Research and WRF Model Users' Page, 2017). The operator defines the forward, tangent linear and adjoint of H for the 4DVAR and 3DVAR case for both ZTD and PW. The operator also defines the observation covariance R; in here diagonal matrix is assumed, with no correlation between observations, which requires spatial and temporal thinning (Bennitt et al., 2017; Bennitt and Jupp, 2012). The ZTD forward operator H reads as follows (Vedel and Huang, 2004) with further corrections made by Yong-Run Guo (from da_transform_xtoztd module of GPSPW):

$$\begin{array}{ll}{\displaystyle }\mathrm{ZTD}\left(i,j\right)& {\displaystyle }=\mathrm{ZHD}\left(i,j\right)+{\sum }_{k=kts}^{k=kte}\left({\displaystyle \frac{wd{k}_{\mathrm{1}}p(i,j,k)q(i,j,k)}{t(i,j,k)}}\right.\\ \text{(2)}& {\displaystyle }& {\displaystyle }\left.+{\displaystyle \frac{wd{k}_{\mathrm{3}}p(i,j,k)q(i,j,k)}{t^{\mathrm{2}}(i,j,k)}}\right){\displaystyle \frac{\mathrm{\Delta }h}{a_{ew}}},\end{array}$$

where ijk are indices of model nodes, p is a pressure, q is specific humidity, t is temperature, Δh is a height difference between two consecutive model layers, a_ew=0.622 is a constant, $wd{k}_{\mathrm{1}}={\mathrm{2.2110}}^{-\mathrm{7}}$, $wd{k}_{\mathrm{3}}=\mathrm{3.73}{\mathrm{10}}^{-\mathrm{3}}$ are compressibility constants, ZHD is a zenith hydrostatic delay computed according to the (Saastamoinen, 1972) explicitly given in Eq. (6).

The PW forward operator is formed similarly to the ZTD operator (following da_integrat_dz module of GPSPW operator):

$$\begin{array}{}\text{(3)}& \mathrm{PW}\left(i,j\right)={\sum }_{k=kts}^{k=kte}\left(\mathit{\rho }(i,j,k)q(i,j,k)\mathrm{\Delta }h\right),\end{array}$$

where ρ is air density.

3.2 GNSS data preprocessing

Two kinds of GNSS data are accepted by WRFDA package: ZTD and PW. In order to prepare the GNSS estimates for GPSPW, a preprocessing is required. The ZTD data are processed according to the following steps.

1. Calculation of the GNSS ZTD using Bernese software for all the stations.

2. Assimilation of the GNSS ZTD obtained in step (1) using the 3DVAR scheme.

3. Calculation of average `background' corrections from the WRF model for each station to reduce the systematic error between WRF and GNSS data and subtracting the corrections from the GNSS ZTD obtained in the step (1).

In the first step, we adjust the formal errors of GNSS ZTD by multiplying them by a factor of 10.5 mm, in which is the standard deviation of differences between WRF ZTD and GNSS ZTD according to Dymarska et al. (2017). Next, we remove the observations, which errors exceed 20 mm, which is the standard procedure in GNSS data assimilation (Bennitt and Jupp, 2012).

In the second step, the GNSS data are assimilated in the 3DVAR procedure in order to calculate the corrections that come from the “background”, which is the WRF model. The corrections for each day are expressed as O−B, where O is the “observation” ZTD (in this case same as ZTD_GNSS from the first step), and B is “background” ZTD, i.e., the WRF ZTD. The corrections are then averaged over the entire considered period to obtain one value (O−B)_av for each station.

In the third step, the corrected ZTDs are calculated as follows:

$$\begin{array}{}\text{(4)}& {\mathrm{ZTD}}_{\mathrm{corr}}={\mathrm{ZTD}}_{\mathrm{GNSS}}-{\left(O-B\right)}_{\mathrm{av}}.\end{array}$$

The PW data are processed in a similar way.

1. Calculation of the GNSS PW from GNSS data.

2. Assimilation of the GNSS PW obtained in the step (1) using the 3DVAR scheme.

3. Calculation of “background'` corrections and subtracting them from the GNSS PW obtained in the step  (1).

From GNSS processing, we can only estimate ZTDs. The PWs in step (1) are calculated in a standard way from GNSS and WRF data as follows:

$$\begin{array}{}\text{(5)}& \mathrm{PW}=Q\cdot (ZT{D}_{\mathrm{GNSS}}-{\mathrm{ZHD}}_{\mathrm{WRF}}),\end{array}$$

where ZHD_WRF is the hydrostatic delay calculated using Saastamoinen (1972) formula from pressure from WRF model p_WRF, height h and latitude φ of a GNSS station:

$$\begin{array}{}\text{(6)}& {\mathrm{ZHD}}_{\mathrm{WRF}}={\displaystyle \frac{\mathrm{0.0022767}{p}_{\mathrm{WRF}}}{\mathrm{1}-\mathrm{0.00266}\cos \left(\mathrm{2}\mathit{\phi }\right)-\mathrm{0.00000029}h}}.\end{array}$$

The proportionality factor Q is calculated as follows:

$$\begin{array}{}\text{(7)}& Q={\displaystyle \frac{{\mathrm{10}}^{\mathrm{6}}}{R_{\mathrm{w}}\left(k_{\mathrm{3}}/T_{\mathrm{m}}+k_{\mathrm{2}}{}^{\prime }\right)}}\end{array}$$

where R_w=461.525 [J (K kg)⁻¹] is the gas constant of a wet air, $k_{\mathrm{2}}{}^{\prime }=\mathrm{22.9726}$ [K hPa⁻¹] and k₃=375463 [K² hPa⁻¹] are the “best average” refractivity constants from Rueger (2002) and T_m is the mean temperature calculated from T_WRF as follows:

$$\begin{array}{}\text{(8)}& T_{\mathrm{m}}=\mathrm{70.2}+\mathrm{0.72}\cdot T_{\mathrm{WRF}}.\end{array}$$

After calculation of GNSS PW, the processing in steps (2) end (3) is analogical to GNSS ZTD.

4 Case studies

All cases presented in this study are selected from the period of May–June 2013 in Central Europe covering the benchmark campaign of COST Action ES1206 (Douša et al., 2016). The following experiments are considered (Fig. 3): (1) assimilation of ZTD or PW for all of May and June 2013, (2) assimilation of ZTD or PW and ZTD or PW with support of RS and SYNOP for 5–23 May 2013, (3) case studies a, b and c, showing impact of assimilation of ZTD or PW in severe weather cases which took place during May and June 2013. The mean and standard deviation of ZTD from GNSS and from WRF are presented in the Fig. 3. The overall agreement between GNSS and WRF traces are high; however, the WRF model is negatively biased with respect to the observations and shows fewer variations. Moreover, few cases of significant departure of WRF ZTD from GNSS ZTD are visible in June, two are collocated with case b and case c investigated in this study.

Figure 3 Time evolution of mean ZTD in the study domain for WRF (red) and GNSS (blue) data (solid lines). Filled area marks standard deviation spread around mean values. Arrows represent the time and duration of analyzed cases.

According to synoptic analysis presented in Douša et al. (2016), the beginning of May 2013 was characterized by a cyclonic field over 500 hPa, which in turn resulted in precipitation and convection development, moving from west to east. The mid-May weather in the region was developing under the influence of an upper-level cyclone (500 hPa) that brought the cold advection from west. Towards the end of May, a series of Atlantic cyclones approached Europe. The end of the month brought a stop to the advection of cold air by the upper-east ridge, which pushed the cyclones more to the south and brought humid and warm air to Central Europe. In June, three flooding events were recorded in the Czech Republic, which were associated with baroclinic instability developing over the area of interest, with a first one (1–3 June) event unexpected and of disastrous nature while two latter (9–11; 23–26 June) less severe and better predicted. As in this work, we use Poland as a study region, there is a time shift between the events recorded in Czech Republic (described by Douša et al., 2016) and Poland and also the precipitation effects were not as disastrous.

The first severe weather case study (case a) was observed in 29–31 May 2013. The weather event is related to an unusual, low-pressure regions: (1) that developed over Hungary and moving towards Czechia, and (2) that developed over Moldova and moving towards the east of Poland. In these two lows, in the presence of stratified clouds, the cumulonimbus clouds develop and form a supercell. It brought intensive rain and hail, however the precise location of such supercells is not easy to predict (http://www.meteo.pl, last access: 17 August 2018).

The second analyzed case (case b) occurred in 17–18 June 2013 and is related to two weather systems: (1) high pressure system with a center in Belarus affecting northern part of Poland, (2) low pressure system over the Bay of Biscay. Cold weather is observed in the north (20 ^∘C) and hot and humid in the south (above 30 ^∘C). The thermic contrasts and warm unstable air result in the occurrence of convective cells located southeast of the region. These cells merged in the late afternoon and formed a supercell storm that moved southward to the Moravia region (http://www.meteo.pl).

The third case analyzed was 24–25 June 2013. The weather in Europe was driven by high pressure system located over the Atlantic Ocean, as well as large and shallow trough extending north to Norway from a weak low centered over the northern part of the Adriatic Sea (with the atmospheric pressure of around 1010 hPa). Secondary cyclogenesis is organized over central Poland in the form of a thermal asymmetric low pressure system. A quasi-stationary anabatic cold front spread along this trough changing its position very slowly and bringing cold air from the north (in the western part of Poland), and warm, humid and unstable air masses from the south (in the eastern part of Poland). These conditions are prone to developing strong precipitation, thunderstorms and hail in central Poland (http://www.meteo.pl).

4.1 Assimilation of GNSS observations

The full period of 2 months is used as a first approach to validate the impact of PW or ZTD data on weather forecasts using radiosonde, GNSS and SYNOP observations.

4.1.1 Comparison to radiosonde profiles

It is expected that the ZTD and PW assimilation first will affect the 3D distribution of humidity and temperature. This is summarized in Fig. 4. and Table 2. In the case of relative humidity, assimilation of ZTD significantly improves the model performance in the layer from 2.5 to 10 km, with a 22 % improvement over the base simulation. At the same time, this assimilation increases the model error for air temperature, but only in the range 2.5–5.0 km and above 10 km. Assimilation of PW has a small impact on both relative humidity and air temperature, and the errors, both in terms of the value and vertical distribution is very similar if PW and base runs are compared (Fig. 4 and Table 2). The assimilation of PW result in significant increase of errors (by 20 %), only in the 2.5–5.0 km section of troposphere.

Figure 4 Vertical distribution of mean error for base, PW and ZTD for May and June 2013 for relative humidity (a) and air temperature (b). The errors are calculated for the lead time 12 h.

4.1.2 Comparison to GNSS observations

The assimilation of the GNSS products should also have an impact on the PW and ZTD estimates from the WRF model. We compare the PW and ZTDs calculated from WRF using formulas (2) and (3) respectively using NCAR command language (NCL) (NCAR, 2018), before the assimilation (“base”) and after the assimilation of both PW and ZTDs. The comparisons are performed for the GNSS data before the processing described in Sect. 3.1. Thus, the presented biases for the base run are removed from the GNSS to better fit the observations to the model. Station by station comparison (not shown) produces similar results across all stations, hence mean statistics is calculated (Table 3).

Table 2 Mean error for RH and T calculated for selected height classes for May and June 2013. Lead time is 12 h. Font denotes improvement (bold), deterioration (bold/italic) or no impact (italic) on the forecast with the use of GNSS data.

Table 3 Bias and standard deviation for PW and ZTD calculated for all GNSS stations in the experiment for May and June 2013. Lead time < 24 h. Font denotes improvement (bold), deterioration (bold/italic) or no impact (italic) on the forecast with the use of GNSS data.

In general, the assimilation of the GNSS does not bring a huge improvement in the WRF estimates. For PW, the bias, averaged from all stations equals to 2.6 mm for the base run, and is slightly improved by 0.1 mm the PW assimilation and remains same for ZTD assimilation. The average standard deviations for the estimates after assimilation improves by 0.2 mm for both PW and ZTDs. For the ZTD assimilation, there is degradation of the ZTD biases after the assimilation of PW (by 0.5 mm) but also improvement in the case of ZTD assimilation (by 0.2 mm). The assimilation of both PW and ZTD brings an improvement of the ZTD standard deviations for almost all of the stations, therefore the average standard deviations decrease for both PW and ZTD assimilation.

4.1.3 Comparison to SYNOP

Further forecast verification is based on 66 SYNOP stations distributed evenly across Poland (Fig. 2). Table 4 summarizes model forecast performance using ME, RMSE, corr and IOA. The statistics are calculated for lead times from 1 to 24 h for the following parameters: rain intensity (rain), wind speed (wspd), relative humidity (rh2) and temperature (T2).

Table 4 Impact of assimilation of PW and ZTD using 4DVAR operators, validated against SYNOP observations for June and May (lead time < 24 h). Font denotes improvement (bold), deterioration (bold/italic) or no impact (italic) on the forecast with the use of GNSS data.

The overall accuracy of the rain forecast is low, i.e., the base run prediction correlates with the observations in less than 15 %, while the same parameter for wind speed is close to 60 %, whereas corr for relative humidity is 82 % and for temperature 95 %. As the assimilation changes the initial conditions of parameters directly linked with the adjoint operator (a transpose of forward operator), the impact while using ZTD should be visible in pressure p (Eq. 2) (and thus also wind speed), specific humidity q (and thus relative humidity) and temperature t. Whereas PW should have impact mostly on specific humidity q (Eq. 3) and thus on the rh2 parameter. Rain as a parameter linked with physical parameterization and many other variables such as humidity, vertical and horizontal motion, or temperature profile are also sensitive to the GNSS data assimilation.

The results (Table 4) confirm that the assimilation of PW over the whole period of 2 months affects the forecasts only slightly, the assimilation increases the relative humidity scatter and has negative or neutral impact on the rain ME and RMSE, neutral or positive impact on wind speed and negative or neutral impact on temperature. Similarly, there is no gain for rainfall forecasts if ZTD is assimilated for the entire period. It has a negative impact on wind speed, but there are considerable improvements for the relative humidity forecast (15 % reduction of ME).

The negative impact of ZTD on the wind speed forecast could be linked to the representation of ZHD as a parameter related only to ground-based observations of temperature and pressure (Eq. 2), whereas in reality the ZHD is an integral of pressure and temperature across the whole troposphere (Vedel and Huang, 2004).

Figure 5 Performance diagram for assimilation of ZTD and PW for May and June 2013 (lead time < 24 h). The various colors represent rain intensity classes, whereas the shapes represent data sets.

If the rainfall forecasts are analyzed more closely using the binary verification with data stratification according to rainfall intensity (Fig. 5), it is clear that the PW run is very similar to the base run, regardless of the rainfall intensity. The ZTD assimilation leads to an overall decrease of the probability of detection.

4.2 Assimilation of GNSS, RS and SYNOP observations

In the second experiment, we focus on a short time span covering May with moderate precipitation and standard cyclonic weather, as opposed to June, with the occurrence of major severe weather events (analyzed in Sect. 4.3). This experiment is prepared to assess the impact of using RS and SYNOP together with GNSS data in 4DVAR assimilation.

Table 5 Mean error for RH and T calculated for selected height classes for 5–23 May 2013. Lead time is 12 h. Font denotes improvement (bold), deterioration (bold/italic) or no impact (italic) on the forecast with the use of GNSS data.

4.2.1 Comparison to radiosonde profiles

Comparison to radiosonde data (Fig. 6, Table 5) shows that largest impact on RH is visible in relatively high levels of troposphere 7–10 km, however improvements to this parameter are also present for 2.5–5.0 km range. Assimilation of PW+SYNOP+RS result in highest gain of quality (in the section 5–10 km by 10 %). In the lower part this impact is less visible (for any observation type), however below 2.5 km assimilation of PW and ZTD alone increase ME (8 % for PW and 2 % for ZTD). The impact on temperature is noted only in the 12.5–15 km sector. Use of SYNOP and RS data improves forecast in all cases but largest increase is visible in the PW observations and in the 2.5–10 km section.

Figure 6 Vertical distribution of mean error for the base, PW, PW+SYNOP+RS, ZTD and ZTD+SYNOP+RS for 5–23 May for relative humidity (a) and air temperature (b). The errors are calculated for the lead time 12 h.

4.2.2 Comparison to GNSS observations

Comparing all active GNSS stations ZTDs and PWs to the WRF based ZTDs and PWs (Table 6), one notice no impact on the PW for ZTD and PW assimilation (both bias and std), whereas for the same experiments ZTD bias increased but scatter decreased considerably (reduction of std by 1.2 mm). The combination of PW+SYNOP+RS and ZTD+SYNOP+RS if PW field is considered no different to the base run, whereas in the case of ZTD being used as a diagnostic variable, the PW+SYNOP+RS provides best solution from all assimilation cases, in terms of bias (only by 0.1 mm). Also mean station deviation is lower for assimilated cases than for base simulation.

Table 6 Bias and standard deviation for PW and ZTD calculated for all GNSS stations in the experiment for 5–23 May 2013. Lead time < 24 h. Font denotes improvement (bold), deterioration (bold/italic) or no impact (italic) on the forecast with the use of GNSS data.

4.2.3 Comparison to SYNOP

The MEs for base run forecasts in May (Table 7) are lower than for May and June, e.g., rain ME is approx. −0.5 whereas May and June is approx. −0.7, May relative humidity ME is approx. 0.4 % and May and June is approx. −2.1 %. Wind speed errors are similar or slightly higher in May than in May and June. A similar statement is correct for temperature errors. The overall correlation between observations and forecasts is in the range from 13 % to over 17 % for rain, 56 % for wspd, 82 % for rh2 and 91 % for T2, which is a few percent lower than in the May and June run (except for rain).

Table 7 Impact of assimilation of PW, ZTD, RS and SYNOP using 4DVAR operators, validated against SYNOP observations, for 5–23 May 2013 (lead time < 24 h). Font denotes improvement (bold), deterioration (bold/italic) or no impact (italic) on the forecast with the use of GNSS data.

Table 7 shows comparison to SYNOP stations and it confirms that the assimilation of either PW or ZTD has a negative impact on rain forecast in terms of mean error, however positive on all other statistics. The largest improvements are for relative humidity, where all the statistics, except the ME, are improved if compared to base. Surprisingly, adding more observations, i.e., SYNOP and RS data, does not improve rain or relative humidity forecast in the case of ZTD assimilation, but rather decreases the forecast's quality. Assimilation of both PW and ZTD deteriorates the model performance for wind speed. There is small but positive impact on the T2 forecast in terms of correlation coefficient, but, similarly to RH, mean errors are increased if compared to base run.

As two forecasted parameters are improved, that is relative humidity and rain (see Table 7), we investigate the lead time differences between the base run and four assimilation setups: namely PW, PW+SYNOP+RS, ZTD, ZTD+SYNOP+RS (Fig. 9).

Figure 7 Performance of rain intensity forecast with respect to lead time (positive – improvement with respect to the base run, negative – deterioration with respect to the base run): (a) ME, (b) RMSE.

The Fig. 7 ME (left panel) of rain forecast varies significantly over 24 h, especially in the lead time 10 to 24 h and is relatively stable between 1 to 9 h (nighttime). In the scattered section of Fig. 9, the ZTD+SYNOP+RS solution seems positive most of the time, while it is negative in the first 9 h of forecast. In the first 9 h of the forecast PW+SYNOP+RS reduces the forecast bias. PW and ZTD alone are rarely observed to improve ME of rain forecast. The RMSE pictured on the right panel of Fig. 7 shows, similar to ME, scattered and compacted sections, however there is clear positive impact of assimilating GNSS observations, especially ZTD+SYNOP+RS in the short run (until 24 h). Overall, the RS and SYNOP data helps to improve RMSE of rain forecast.

Figure 8 Performance of rh2 forecast with respect to lead time: (a) ME, (b) RMSE.

Less variation between the four scenarios is observed for relative humidity errors (Fig. 8). Both ME and RMSE are reduced while assimilating each data type, with an exception of lead time 15 h, when RMSE increases (more when ZTD is used, less when PW is used). Moreover, the highest reduction of ME is noticed between 12 and 18 h of lead time for PW+SYNOP+RS scenario, but other scenarios also shows a positive impact. It is also worth mentioning that the rh2 forecast RMSE is reduced after 12 h lead time whereas bias is constantly reduced starting from the first hour of forecast.

Figure 9 Performance diagram for assimilation of PW, PW+SYNOP+RS, ZTD+SYNOP+RS and ZTD for May 2013 (lead time < 24 h). The various colors represent rain intensity classes, whereas the shapes represent data sets.

Binary rain analysis shows that the impact of data assimilation on rainfall forecasts changes with rainfall intensity (Fig. 9). For the rainfall intensity above 0.1 mm h⁻¹, there are small improvements for all the model runs if compared to the base run in terms of success ratio, but the probability of detection is smaller. The positive impact of data assimilation is much stronger for higher rainfall intensities. For the thresholds exceeding 1.0 mm h⁻¹, both probability of detection and success ratio are improved if compared to the base run. The improvement is especially large for PW data assimilation and threshold > 5.0 mm h⁻¹.

4.3 Severe weather cases

The final test is performed using selected three cases with strong instabilities and supercell storms. The overall impact of GNSS data in all cases is similar: if there is any reduction in uncertainty it is visible mostly in the rain and relative humidity forecast, with a small negative or neutral impact on the wind speed and temperature forecasts.

4.3.1 Case (a) 29–31 May 2013

The comparison of WRF-based RH and T with radiosonde shows that for case (a) (Fig. 10) the impact of ZTD assimilation introduces large change into the initial conditions and forecast. This impact is RH positive in the first 2 km and between 6 and 10 km and might be negative close to 2.5 km. Temperature show mixed results, with stronger influence of ZTD and weaker for PW, but across the whole profile this impact is small but positive. Moreover ZTD increase agreement between model and RS in the bottom 2 km of temperature profile.

Figure 10 Vertical distribution of mean error for the base, PW, and ZTD for case (a) in 29–31 May 2013 for relative humidity (a) and air temperature (b). The errors are calculated for the lead time 12 h.

Table 8 Bias and standard deviation for PW and ZTD calculated for all GNSS stations in the experiment for selected cases (case a) 29–31 May 2013; (case b) 17–19 June 2013; (case c) 24–26 June 2013. Lead time < 24 h. Font denotes improvement (bold), deterioration (bold/italic) or no impact (italic) on the forecast with the use of GNSS data.

Comparing the forecast with assimilation of ZTD and PW with GNSS observations (Table 8), shows that both data reduce scatter of the observations but ZTD is also reducing systematic effects, while PW is increasing bias.

According to Table 9, which presents a comparison to SYNOP data, mixed results are observed for case (a). The rain forecast shows better performance if PW is compared to base run. Assimilation of ZTD for this case deteriorates model performance except for mean error. Humidity forecast (rh2) is improved, in terms of RMSE and corr, when ZTD is assimilated. For PW assimilation, model performance is worse if compared to base.

Figure 11 Performance diagram for assimilation of PW and ZTD for case (a) (lead time < 24 h). The various colors represent rain intensity classes, whereas the shapes represent data sets.

Binary analysis depicted in Fig. 11 shows the positive impact of assimilating of PW in a rainfall rate above 0.5 mm h⁻¹, and ZTD above 2.5 mm h⁻¹.

4.3.2 Case (b) 17–18 June 2013

The comparison of WRF-based RH and T with radiosonde shows that for case (b) (Fig. 12) the impact of ZTD assimilation introduces large change to the initial conditions and forecast, while PW has almost no impact. The ZTD impact is negative in RH in the first 5 km and positive from 7 to 10 km. The temperature profile shows positive results, with stronger influence of ZTD and almost none for PW. It is visibly positive in the bottom part of the troposphere (first 1 km).

Figure 12 Vertical distribution of mean error for the base, PW and ZTD for case (b) 17–18 June, for relative humidity (a) and air temperature (b). The errors are calculated for the lead time 12 h.

Comparing forecast with assimilation of ZTD and PW with GNSS observations (Table 8), shows that both data reduce the scatter of the observations. The bias for PW parameter slightly increases for both types of assimilated observations, however ZTD bias is reduced.

Table 9 Impact of assimilation of PW, ZTD, RS and SYNOP using 4DVAR operators, validated against SYNOP observations, for selected cases (case a) 29–31 May 2013; (case b) 17–19 June 2013; (case c) 24–26 June 2013 (lead time < 24h). Font denotes improvement (bold), deterioration (bold/italic) or no impact (italic) on the forecast with the use of GNSS data.

Comparison to SYNOP data (Table 9) shows that, the assimilation of PW does not change model performance for rainfall, except for ME reduction. The model performance for other parameters is similar or slightly worse than base. ZTD assimilation has small positive impact on RMSE, corr and IOA of rain intensity forecast, but negative on ME. Relative humidity ME is reduced by assimilation of ZTD by 43 %, all other measures are better for base run. In the local type of rain in southeastern Poland, as in this case, it is impossible to present statistically sound results for 5 rainfall classes, hence we did not provide binary rain analysis for this case.

4.3.3 Case (c) 24–25 June 2013

The comparison of WRF-based RH and T with radiosonde shows that for case (c) (Fig. 13) the impact of ZTD assimilation introduces large change to the initial conditions and consequently to the forecast. This impact is RH positive in the first 1 km, negative close to 2.5 km and clearly positive at 5 km. Temperature shows mixed results, however it is worth pointing out that from 2.5 to 5 km there is positive bias introduced to the forecast by both PW and ZTD. Temperature is improved once ZTD is assimilated in the bottom part of troposphere (first 1 km).

Figure 13 Vertical distribution of mean error for the base and PW for case (c) 24–26 June, for relative humidity (a) and air temperature (b). The errors are calculated for the lead time 12 h.

Comparing forecast with assimilation of ZTD and PW with GNSS observations (Table 8), shows that both data reduce scatter of the observations. The bias for PW parameters slightly increases (ZTD assimilation) or is neutral (PW assimilation). The bias for ZTD decreases for PW assimilation and slightly increases for ZTD assimilation.

The third case, once compared to the SYNOP, also shows a small but positive impact of ZTD, especially in PW data assimilation on rainfall forecasts, except for mean error. Errors statistics are also improved for wind speed. In the case of the PW run, there is also a gain for relative humidity, while for ZTD error statistics are worse if compared to the base run.

Figure 14 Performance diagram for assimilation of PW and ZTD for case (c). The various colors represent rain intensity classes, whereas the shapes represent data sets.

In the performance diagram in Fig. 14, the rain rate forecasts are improved with PW with respect to the base forecast, but worse when ZTD is assimilated. This effect is visible for all rain rates lower than 1 mm h⁻¹ and this discrepancy disappears for rain rates in the 2.5 mm h⁻¹ class, where both ZTD and PW have positive impact, whereas no impact is noticed for rainfall rates above 5 mm h⁻¹.

5 Summary and conclusions

In this study, we have analyzed 2 months (May and June 2013) of 4DVAR assimilation of GNSS ground-based observations in WRF model from over 100 stations in Poland. Two major approaches were investigated using GPSPW operator: assimilation of PW and ZTD. For a shorter time period of 18 days in May additional data were assimilated: namely RS and SYNOP observations across Poland. Moreover, three different case studies related to severe weather occurrence were investigated. All were linked to a supercell development and intense rain.

The major conclusion of this study is that for the analyzed time period, with more than 100 stations involved in the experiment, there are two parameters significantly changed once GNSS data are assimilated in the WRF model using GPSPW operator and these are the moisture field and rain. Other parameters such as pressure or temperature field are not changed from initial conditions significantly. The GNSS observations improves forecast in the first 24 h but with strongest impact starting from 9 h lead time. It is worth noticing that even moderate quality NRT estimates used in this study (ZTD discrepancy ∼10–15 mm) are improving relative humidity forecast, moreover the impact of ZTD is positive in the vertical profile (over 20 % decrease of mean error) starting from 2.5 km upward. Even if the humidity forecast in lower part of troposphere (below 2.5 km) after GNSS data assimilation deteriorates, the SYNOP observations confirms that ZTD has positive impact on the rh2 parameter. Assimilation of PW has less significant impact on both humidity and rain forecasts in the vertical profile and on the ground, it improves slightly rain forecast and deteriorates the humidity prediction by a few percent. The best results for use of integrated water vapor information are produced by adding radiosonde and SYNOP data to the assimilation system as it improves distribution of humidity in the vertical profile. In the three analyzed severe weather cases, PW always improved the rain forecast and ZTD always reduced the humidity field bias. Binary rain analysis shows that GNSS parameters have a significant impact on the rain forecast in the class above 1 mm h⁻¹, while adding more observations such as radiosonde and SYNOP result in a larger increase of quality. In each case however the impact of GNSS data was different, therefore below we provide a brief summary for each case.

The May and June case results showed that in the lead time of 24 h, the assimilation of GNSS data (both ZTD and PW) had a varying impact on a number of parameters. The assimilation of PW over the whole period of 2 months marginally diverted the humidity and temperature 3D performance, which is seen on the radiosonde profiles. It also improves the agreement between WRF model based and GNSS based PW observations by 0.1 mm. Assimilating PW slightly increased correlation of daily rainfall rates, increased relative humidity scatter and had positive impact on the rain RMSE, neutral impact on wind speed and negative or neutral on temperature. However, the binary analysis of rain rate in five intensity classes revealed that the forecasts with assimilation of PW improves forecast scores in high intensity rain above 2.5 mm h⁻¹.

Assimilation of ZTD had a large impact on the vertical profile of both temperature and humidity, as retrieved from comparison to radiosonde data especially in the range 5 to 10 km. It reduced bias by 0.1 mm and std by 0.5 mm between the base run and assimilation run for WRF derived ZTD and GNSS derived ZTD. Assimilation of ZTD reduced the ME for humidity by 16 %, while it had a slight negative impact on rain, temperature and wind speed bias. The binary analysis of rain rate in five intensity classes revealed that the forecasts with assimilation of ZTD improve forecast scores only in the highest rain rate class.

The more detailed study focused on 5–23 May (with non-severe weather events) and showed that the assimilation of either PW or ZTD reduced mean error in the humidity forecast in the vertical direction, most successfully for PW+SYNOP+RS data. It had a positive impact on the rain forecast and relative humidity forecast. The assimilation of PW improved rain rate RMSE by 5 % and had negative impact on bias (7 % increase). The relative humidity forecast bias was doubled with assimilation, however RMSE was reduced by minimum of 1 %. The assimilation of ZTD improved rain rate RMSE by 6 % and had a negative impact on bias (14 % increase). The relative humidity forecast bias was doubled with assimilation and RMSE was reduced by less than 1 %. Adding SYNOP stations and radiosonde did not bring any further improvements in forecasting humidity or rain but reduced the errors in wind speed and temperature data. Furthermore, the analysis of lead time with respect to the errors revealed that for rain rate ME error varies in time (both negative and positive impacts are present), whereas the RMSE for data with assimilation is considerably improving with time. In the case of relative humidity, both ME and RMSE are reduced when GNSS data are assimilated, the largest gain in quality is observed for PW+SYNOP+RS data set. The binary analysis show positive impact of GNSS data assimilation especially for rain rates above 1 mm h⁻¹

In the analyzed severe rain cases, the assimilation of GNSS in case (a), (b) and (c) brings reduction of ME and RMSE or at least RMSE for key sensitive parameters such as rain rate, relative humidity. Binary rain rate forecast performance analysis shows that the intensive rain is better predicted once GNSS data are assimilated. Further research, based on a larger number of cases, is required to investigate what the reasons of different impacts of GNSS data are on model forecasts.