Tropospheric parameters from multi-GNSS and numerical weather models: case study of severe precipitation and flooding in Germany in July 2021

Storms and floods are the most frequent natural disaster in western and central Europe. Due to climate change, intensive storms with prolonged rain episodes will continue to cause even more destructive flooding. The good understanding and forecasting of such events are thus of uttermost importance. One of the ways to improve weather forecasts is the assimilation of external data, such as the global navigation satellite systems (GNSS). In this study, a preparation of the multi-GNSS tropospheric products—zenith total delays, tropospheric gradients and slant total delays—for future operational assimilation is shown. For a severe precipitation event in July 2021 in Germany, the GNSS parameters from three systems—GPS, GLONASS and Galileo—are compared to three Numerical Weather Models (NWMs)—ERA5 reanalysis of ECMWF, ICON run by the German Weather Service and GFS run by the US Weather Service. The flood that followed the rainfall was the deadliest natural disaster in Germany since 1962. The results show that all considered GNSS solutions have a similar level of agreement with the NWMs. However, for the flood region in western Germany, the biases from the multi-GNSS solutions are smaller compared to the GPS-only solution. From the models, ICON has the highest agreement with the GNSS data for all considered tropospheric parameters. The best agreement with the ICON is probably due to its high horizontal resolution and, thus, low representative errors and the fine tuning of DWD’s regional model for the specific region (Germany).


Introduction
Global navigation satellite systems (GNSS) meteorology is a way of observing the meteorological parameters (mostly water vapor) by exploiting the GNSS. For GNSS, the tropospheric delay used to be just a source of errors, but currently, it is also considered an important meteorological signal. It originated as GPS meteorology in a paper by Bevis et al. (1992) and has been gaining growing interest over the past 30 years.
The interest in GNSS meteorology is demonstrated in two major European Cooperation in Science and Technology (COST) actions and many other scientific projects. The COST action 716 Exploitation of ground-based GPS for operational weather prediction and climate applications (period: 1999-2004) was the first large, coordinated activity in the field of GNSS meteorology. It dealt with the assessment of the operational potential on an international scale of the exploitation of a ground-based network of GPS receivers to provide near real-time observations for Numerical Weather Models (NWMs) and climate applications (Elgered et al. 2005). This action resulted also in the establishment of several large projects in GNSS meteorology, such as the project Towards Operational Use of GPS Humidity Measurements in Meteorology (TOUGH, http:// tough. dmi. dk/).
The success of the COST action 716 was extended by another COST action ES1206 Advanced GNSS Tropospheric Products for Monitoring Severe Weather Events and Climate (Jones et al. 2020). The action duration was 2013-2017 and it focused on the usage of GNSS meteorology for severe 49 Page 2 of 17 weather events. One of the main outcomes was establishing a benchmark campaign using different techniques (i.e., GNSS, NWM, radiosondes, water vapor radiometer and radar images) of monitoring heavy precipitation events that occurred in May-June 2013 causing severe flooding of the Danube, Moldau and Elba rivers (Dousa et al. 2016). However, there were also several other studies that involved the GNSS monitoring of severe weather events, especially precipitation. Most of the studies focused on the relations between the Integrated Water Vapor (IWV) calculated from GPS observations and the precipitation events (Seco et al. 2012, Choy et al. 2013, Zhang et al. 2015, Bonafoni and Biondi 2016, Priego et al. 2017. Brenot et al. (2013) additionally presented the usefulness of the GPS tropospheric gradients as preliminary signs of deep convection. There are also some more recent studies showing that GNSS time series, not only IWV but also zenith total delays (ZTDs), can be used as indicators of precipitation (Zhao et al. 2018;Li et al. 2020Li et al. , 2021. This study is conducted within another GNSS meteorology project Advanced MUlti-GNSS Array for Monitoring Severe Weather Events (AMUSE, https:// www. gfz-potsd am. de/ en/ secti on/ space-geode tic-techn iques/ proje cts/ amuse/) of the German Research Foundation (DFG). The main objectives of this project are the developments to provide high-quality, multi-GNSS advanced tropospheric products that could be used in the future for the assimilation into NWMs. In this paper, we compare the tropospheric parameters, ZTDs, tropospheric gradients and slant total delays (STDs), from three GNSS solutions with three NWMs for a case of severe rainfall that lead to deadly flooding in northwestern Germany in July 2021. This paper is a continuation of a study by Wilgan et al. (2022), where the tropospheric parameters from three GNSS systems were compared with ERA5 reanalyzes for the year 2020. Wilgan et al. (2022) found (1) ZTDs average bias for Germany of 3 mm and a standard deviation (SD) of 7 mm; (2) GE bias of 0.01 mm with 0.4 mm SD; (3) GN bias of 0.03 mm with 0.4 mm SD; (4) STD bias of 4 mm and 26 mm SDs. All three GNSS solutions showed very similar agreement with ERA5. In this study, we put the focus on a specific extreme weather event as well as on the data from three different NWMs. The next section presents the meteorological conditions during the event, then, the data are described and the results of the comparisons are shown. The last section summarizes the results.

Meteorological conditions
The July 2021 flood affected Germany, Belgium and France causing more than 180 casualties and large material damage. The flood mechanism was as follows: from July 6 to July 12, there was a series of precipitation events that along with a slow moving pressure system led to an accumulation of rain on July 13-15 and following heavy and long lasting rainfall events (Puca et al. 2021). In Germany, the largest flooding was in the state of North Rhine-Westphalia and Rhineland-Palatinate, especially in the district of Ahrweiler. The precipitation and flooding were actually well predicted. The European Centre for Medium Weather Forecasts (ECMWF) indicated the presence of a very unusual precipitation event over Western Germany. On July 13, the German Weather Service (DWD) issued heavy precipitation warnings. The indications of the flood were also found in the soil wetness. The anomalies of the H SAF Soil Wetness Index showed very wet conditions in the flooded regions (10-30% above average) for July 14. A few days before the flooding and on the day after the heaviest rainfall, the soil wetness index increased towards saturation (between 10 and 15 July). The very wet antecedent soil moisture conditions and the significant precipitation amounts were key drivers of the catastrophic flood occurred on July 14-15 (Puca et al. 2021). In this study, we focus on one week, July 12-18, which includes the severe precipitation and flooding events. Figure 1 shows the accumulated weekly rainfall, and Fig. 2 shows the accumulated daily rainfall measured by a network of rain gauges in Germany provided by the DWD.
As shown in Fig. 1, the majority of rainfall was observed in western Germany, in the two federal states mostly affected by the flood. Figure 2 shows that the highest precipitation was observed on July 14, with a record of 154 mm at station Köln-Stammheim (Cologne). There was also severe precipitation in the south of Germany (Bavaria), especially for July 17; however, in this study, we focus on the western region (longitude between 6° and 7.5° E and latitude between 49.5° and 52° N). The region was chosen based on the precipitation and the availability of the GNSS station data.

Data
We compare the tropospheric parameters: ZTDs, tropospheric gradients and STDs derived from multi-GNSS observations and three NWMs: global reanalysis ERA5, global forecast model GFS and local analysis model ICON. Two areas are considered: entire Germany and the region with the highest rainfall that includes the GRE-capable stations.

Global navigation satellite systems
The GNSS data are processed using the Precise Point Positioning (PPP) technique (Zumberge et al. 1997) with the GFZ-developed software EPOS.P8 (Dick et al. 2001;Gendt et al. 2004). In the preprocessing step, the GFZ high-quality orbits and clocks are estimated using a base of approx. 100 stations located uniformly around the world (Ramatschi et al. 2019). Then, the tropospheric parameters are adjusted using the 24 h data intervals with a sampling rate of 15 min for ZTD and tropospheric gradients and 2.5 min for STD. The parameters are calculated for 199 stations in Germany. Figure 3 shows the stations for Germany capable of receiving GPS only (G), GPS/GLONASS (GR) and GPS/GLO-NASS/Galileo (GRE) solutions. Not all stations are GRE capable yet, but the number of such stations is increasing compared to Wilgan et al. (2022) for the year 2020. In this study, we use only the GRE-capable stations from the German national SAPOS network. The red rectangle denotes the region of severe rainfall and flooding. Figure 4 shows the zoom at the stations for the flood region (13 stations).
In the GNSS analysis, the tropospheric delay is approximated according to: where ZHD and ZWD are the hydrostatic and wet parts of the ZTD, respectively, G N and G E denote north-south and east-west gradient components, MF h , MF w and MF g are the mapping functions for the hydrostatic, wet part (here, the Global Mapping Function (GMF) from Böhm et al. (2006)) and gradients (here, Bar-Sever et al. (1998)); el is the elevation angle, A is the azimuth angle and res are the postfit phase residuals. The a priori ZHDs are taken from the Global Pressure Temperature 2 (GPT2) model (Böhm et al. 2007;Lagler et al. 2013). The cut-off elevation angle is 7°. More details about the processing can be found in Wilgan et al. (2022).

Numerical weather models
Three NWMs are used for the comparisons in this study: ERA5, ICON and GFS. ERA5 is the 5th generation reanalysis from ECMWF (https:// www. ecmwf. int/ en/ forec asts/ datas ets/ reana lysis-datas ets/ era5). The horizontal resolution of the model is 0.25° × 0.25° with 31 pressure levels up to 80 km. The temporal resolution is 1 h. There is no groundbased GNSS data assimilation in the model, but GNSS Radio Occultations (RO) are assimilated (Healy et al. 2005).
The Icosahedral Nonhydrostatic (ICON) model is run and provided by DWD (https:// www. dwd. de/ EN/ resea rch/ weath erfor ecast ing/ num_ model ling/ 01_ num_ weath er_ predi ction_ model ls/ icon_ descr iption. html). Here, we use the regional ICON-D2 model, which is nested from the global ICON model. The horizontal resolution used in this study is 0.02° × 0.02° with 65 vertical model layers up to 20 km. The temporal resolution is 1 h. The GNSS ZTDs and ROs are assimilated into the ICON global model, but not (yet) into the nested, regional model. We use the ICON analysis, not the forecast data to obtain the best possible ICON parameters.
The Global Forecast System (GFS) is run by the US National Centers for Environmental Prediction (NCEP) (https:// www. ncei. noaa. gov/ produ cts/ weath er-clima te-models/ global-forec ast). The horizontal resolution is 0.25° × 0.25° with 34 vertical pressure levels up to 0.01 hPa which corresponds to about 40 km. The temporal resolution is 1 h. The GNSS data are not assimilated into GFS. We use short-term forecasts (between 6 and 11 h lead time).
The NWMs provide gridded pressure, temperature and humidity fields. To calculate the tropospheric parameters, at first, the gridded refractivity field N is computed using the following equation (Thayer 1974): where p is the atmospheric air pressure [hPa], e is the water vapor partial pressure [hPa], T is the temperature [K], k 1 = 77.60 [K · hPa −1 ], k 2 = 70.4 [K · hPa −1 ] and k 3 = 373,900 [K 2 · hPa −1 ] are the refractivity coefficients from Bevis et al. (1994); Z −1 d and Z −1 v are the inverse compressibility factors for dry air and water vapor, respectively, usually approximated as 1.
To interpolate the refractivity at any point, four surrounding refractivity profiles are identified; for each profile, a logarithmic interpolation adjusts the refractivity vertically, and then, a bilinear interpolation including the vertically adjusted refractivity values is performed (Zus et al. 2012). Global NWMs, such as ERA5 and GFS, typically provide refractivity values up to high altitudes (e.g., 80 km); thus, no special routine is required to calculate refractivity above the model top. In essence, the interpolation routine turns into an extrapolation routine for points above the model top. The situation is different for limited area NWMs, such as ICON-D2, where the model top is at an altitude of 20 km. Therefore, to mitigate errors caused by the extrapolation routine, a different approach is chosen. If we assume that the humidity is negligible and the temperature is constant above the model top, the pressure can be derived from the hydrostatic equation. Subsequently, the refractivity N at the geopotential height H is given by where N 0 , T 0 and H 0 denote the refractivity, temperature and geopotential height at the model top, respectively, and GMR denotes the hydrostatic constant. This approach above the model top is justified by the following observation: when we compare the numerically integrated ZHDs with the ZHDs derived from the semi-analytical formula proposed by Davis et al. (1985), we obtain a sub-mm agreement. On the other hand, when we simply extrapolate the refractivity above the model top, the numerically integrated ZHDs deviate from the ZHDs derived from the semi-analytical formula by about 3 mm.

Fig. 4 GRE-capable stations in the flood region
Page 5 of 17 49

Results
We compare the tropospheric parameters from the GNSS and NWM data for entire Germany and the flood region for July 12-18, 2021. The following subsections contain the comparisons for ZTDs, tropospheric gradients and STDs.

Comparisons of zenith total delay
At first, we compare all three GNSS solutions with the three NWMs. Table 1 shows the biases and SDs between the estimates from NWMs and the GNSS solutions (G, GR, GRE). Figure 5 shows a map of biases and SDs between the models and the GRE solution. For Germany, and especially for the flood region, there is a small reduction of bias for the residuals between each NWM and GNSS when we take GRE solutions versus GR or G; however, for ERA5 and GFS, the SDs are slightly increased. Only for ICON, the SDs remain at the same level, while biases decrease. Visually, the differences between the particular GNSS solutions and NWMs are very small and thus, are not shown.
As shown in Fig. 5, the residuals between the ERA5 and GFS models are more similar to each other than to the ICON model. For ICON, the biases are smaller and below 5 mm for most of the country, except for eastern Germany. The biases for ERA5 are not uniformly distributed-they are negative for central Germany and positive for the east, west and south of the country. For GFS, the biases in the south are the highest and reach even 10 mm. The SDs are generally the lowest for ICON, which we also can see in Table 1; the SDs for ICON are around 2 mm smaller than for the other models. However, for all NWMs, the SDs follow the same pattern: they are the highest in the east of Germany. There are two reasons why the SDs are the lowest for ICON. First, we take the analysis of the model and not the forecasts to have the best possible representation of the model. Second, the resolution of the ICON model is one order of magnitude higher than of the global models, and thus, ICON can resolve small-scale processes better than the other models. Figure 6 and Table 2 show the statistics between the NWMs, calculated only from the interpolated points of the GNSS stations.  Table 2 shows that the biases between the models are smaller than between NWMs and GNSS. Figure 6 shows that the biases between ERA5 and GFS are uniform and negative for most of the country, except for the east. The largest differences between ICON and GFS are in the south, where the biases reach even 10 mm for some stations. The differences between ERA5 and ICON are positive for the south and the west, while negative for central Germany; thus, the average bias equals only 0.5 mm. The biases are not evenly distributed, but the SDs are as follows: the patterns also look very similar to the patterns in Fig. 5. The SDs between all models are of around 13 mm and they are in general the highest for eastern Germany. For the flood region, the SDs are in general smaller than for Germany. Moreover, the bias between ICON and GFS is close to zero, meaning that for this region these two models are more similar than ERA5. To take a closer look at the flood region, Fig. 7 shows the statistics between the NWMs and GRE solution. Figure 7 shows that the residuals between the models and GRE are different for each NWM. For ERA5, we have a large, positive bias, of around 6 mm. For ICON, the differences vary mostly between 2 and 4 mm and for GFS they are positive for the southern stations of the region, but close to zero or even negative for the northern. The SDs are the smallest and very uniform for ICON (around 8 mm), around 9-10 mm for GFS and positive and even reaching 12 mm for ERA5. Thus, for this region, similarly to the entire Germany, ICON has the best agreement with GNSS; however,  here GFS performs much better than ERA5. The majority of precipitation in the flood area happened on July 13-15, as shown in Fig. 2. Table 1 includes also statistics only for the period of precipitation. For this period, we can observe an increase in ERA5 bias and SDs compared to the whole week. For ICON, the biases are smaller and SDs are at a similar level, while for the GFS model, the biases are smaller, but the SDs are larger. However, the statistics for different GNSS solutions follow the same pattern as for the entire week: the biases for GR and GRE solutions are smaller than for G, while the SDs slightly increase. It shows that we can draw similar conclusions from the statistics for the whole week as from the days with increased precipitation. Next, we plot the ZTD time series for sample stations. Figure 8 shows the values for station 2615, located in the middle of the flood region and Fig. 9 for station 2623, located in the north. Figure 8 shows the evolution of the ZTD values for the week that includes the heavy rainfall for station 2615. The ZTD values for the days of rainfall are steadily rising, while on the night between July 14 and 15 we see a sudden drop in the values, due to the drop of humidity values after rainfall. Then, the values rise again on the 15th and are very stable until July 17, when they drop by around 10 cm. For station 2623 (Fig. 9), we can see similar behavior in the GNSS data: the values are increasing until the night between July 14-15. However, the highest rainfall for this station is on July 15, when the ZTD values start to drop. The large drop in values is also visible for July 17. Figure 10 shows the differences between the NWMs and the GRE solution for station 2615 and Fig. 11 for station 2623.
Even though stations 2615 and 2623 are located relatively close, the statistics between the particular NWMs and GRE differ significantly. For station 2615, as shown in Fig. 10, there is a big overestimation of ERA5, especially on the night of July 14-15, right after the large precipitation on the 14th. For station 2623, the large precipitation only occurs on July 15 and for this day we have a large misestimation of ERA5 values (however, there are both positive and negative values; thus, the average bias is not so large). It is also obvious that the SDs are the largest from the ERA5 model; however, the values also differ for these two stations. For this region, the smallest biases are from both ICON and GFS models; however, GFS has here smaller biases, which are again different for the two stations. Figures 10 and 11 show the differences between the NWMs and the GRE solution; however, it is worth mentioning that there are small differences between the particular GNSS solutions. For example, for station 2615, the bias between the ERA5 and G solution is 7.9 mm, while between GRE it is 6.4 mm. Such a 1.5 mm difference between the G and GRE solutions is visible for all three NWMs. The SDs remain at the same level for all the solutions. For station 2623, the biases between ERA5 and G, GR and GRE are 2.9 mm, 2.7 mm and 2.2 mm, respectively, so there is a small improvement for the GRE solution. The SDs are around 0.1 mm larger for GRE than G. Such improvement is visible for around half of the stations in the region, and for the other half, there is no difference between the GNSS solutions. However, it shows that using GRE instead of GPS only can improve slightly the agreement with NWMs for severe weather events.

Comparisons of tropospheric gradients
In the next step, we compare the tropospheric gradients. The gradients are the measure of the anisotropy in the troposphere in the east-west (GE) and the north-south (GN) directions. Their magnitude is small, usually close to 1 mm and can reach up to a few millimeters in extreme cases.   Figure 12 shows the biases and SDs for all stations, and Table 3 shows the averaged statistics for Germany and the flood region. Figure 12 shows the differences between the three NWMs and the GRE solution for both gradient components. The biases are rather small and uniform across the models and the country, contrary to the ZTD case. For the SDs, we can see a geographical pattern. For GE, the SDs are the highest for eastern Germany and in the flood region for ERA5 and GFS. In the middle of the country, the values are the lowest and are also lower for ICON than for the global models. For GN, the values are the highest in the east and the north, and the smallest for the south of the country, like in the GE case. The values in the east are, however, larger than in the GE case.
We do not show the analogical figure for the flood region, but we mention that for GE, the biases are very small and similar for all models. A similar situation is observed for GN except for a few stations in the middle of the region (it is also visible as red/yellow dots in Fig. 12). The SDs for both GE and GN are larger for ERA5 and GFS than for ICON and the values are larger in the west of the region (also visible in Fig. 12). Table 3 shows the statistics between the NWMs and the three GNSS solutions. We can see smaller differences between the models and GNSS than in the case of ZTDs. The GE biases for most of the cases are smaller than 0.01 mm and the differences between the GNSS solutions are really small. For GN, there is a small reduction of bias for GR and GRE compared to the G solution for the ICON model, but for the GFS model, we observe the opposite. For ERA5, the biases are at a similar level. For the SDs, the GRE values are the largest for both gradients and all models. For the flood region and GE, the biases are slightly reduced for GR and more for GRE compared to G for all models. For GN, we can see the reduction for ERA5 and ICON (but GRE is here slightly larger than GR), but for GFS, the biases are increased for GRE. Here, the SDs have similar behavior as for entire Germany. The smaller difference between the GNSS solutions compared to the ZTD values is because the gradients are more tightly constrained between epochs. Similarly to the comparisons for ZTD, if we consider only the period of precipitation, i.e., July 13-15, we obtain worse statistics than for the entire week. The SDs from all models are higher than for the entire week, which is understandable as gradients are sensitive to extreme weather events. Moreover, the biases are larger than for the entire week, which could also be explained by the fact that the values of gradients, especially the GE, are very large for the rain period (cf. Figure 13). The general patterns in the GNSS solutions are also similar for the entire week and the chosen period: all biases are smaller for GR and GRE compared to G, while the SDs are slightly larger. Next, we compare the gradient values between the NWMs. Table 4 shows the results. Table 4 shows the differences between the particular NWM solutions. The biases are rather small; for the entire Germany, they do not exceed 0.02 mm for GE and 0.04 mm for GN, while for the flood region the largest differences exist between ERA5 and GFS of − 0.05 mm for GE and between GFS and ICON of − 0.07 mm for GN. For GN, also the biases between ERA5 and GFS are high and equal to 0.06 mm. For SDs, the values are very similar for all the models. We can conclude that for the gradients, the differences between the particular NWMs are smaller than in the case of ZTDs, which means that the NWMs are estimating gradients at the same level.
Next, we take a look into the time series of gradients for particular stations. Figure 13 shows the values of the gradients for station 2615, and Fig. 14 shows the differences between the three NWMs and the GRE solution for this station. The same station was previously used for the ZTD comparisons.
As seen in Fig. 13, the GE values dropped very significantly on the night of July 14-15, which coincides with the    Figure 14 shows that the differences between NWMs and GNSS for this night are higher than for the other days. Only for ERA5, there are several points where the model overestimates the values for GE. For GN, we observe such a peak in the gradient values, but on July 17 (when also a drop in the ZTDs is observed). Moreover, for both gradient components, we can see some point-wise overestimation of ERA5 values (e.g. for July 14). The histograms in Fig. 14 show the values of average biases and SDs from the three NWMs, which, unlike in the ZTD case, are very close to each other. Other stations in the region (figures not shown) exhibit a similar gradient behavior, especially a drop in the GE values on the night of July 14-15 and a drop in GN values for July 17. For the gradients, all GNSS solutions behave similarly; however, we can see some improvements of using GRE solution versus G or GR for station 2615. For GN, we observe a bias reduction for all three NWMs of about 0.1 mm, which is a few percent of the total gradient value. The SDs remain at a similar level, as in the case of ZTD. For GE, the reduction in bias is much smaller; however, for this station, the GRE solution seems to perform the best for the chosen period. To look closer at the situation on July 14-15, we plot the two components of the tropospheric gradients as vectors. Figure 15 shows the values of the gradients from GNSS and NWM for selected dates: July 14, 7:00 UTC; 12:00 UTC; 18:00 UTC; and July 15, 0:00 UTC covering the main precipitation event. Along with the gradients, we also show the hourly rainfall information for the selected dates.
As shown in Fig. 15, for July 14, 7:00 UTC, the gradients in the south and east of Germany are pointing to the same north direction. The values from GNSS are higher than from NWMs but very similar to each other. Also, all NWMs are quite similar to each other, just the values from ICON are slightly larger. We can also see that the precipitation in the flood region has started. For 12:00 UTC, all gradients except for ERA5 are very similar. The ones in the middle of Germany are pointing to the flood region. For ERA5, the gradients have a similar direction, but the values are much higher, which can also be seen in Fig. 14. At 18:00 UTC, all gradients are pointing in the same direction (from the middle of Germany outwards), but the values from NWMs are larger. There is less precipitation in the flood region. Finally, for July 15, 0:00 UTC, after the main precipitation event, the gradients are generally smaller and noisier. The patterns are similar for the three GNSS solutions and the NWMs for most of the country, except for the east of Germany, where even the NWMs are not consistent. For this date, the precipitation  in the flood region has mostly stopped, but there is some in Baden-Württemberg and the local gradients are also pointing toward it. Such behavior of gradients pointing towards higher water vapor concentration can also be found in Kacmarik et al. (2019).

Comparisons of slant total delay
We compare the STDs from the GNSS solutions and two NWMs, ERA5 and ICON. ERA5 and GFS have similar statistics for the previous comparisons and thus, GFS is not shown here. Table 5 shows the statistics between the two NWMs and three GNSS solutions for Germany and the flood region, while Fig. 16 shows the differences between the NWMs and GRE solutions for Germany. Table 5 shows the statistics between the NWMs and GNSS. It consists of absolute values, relative values and mapped ZTD values. For Germany, the biases from ERA5 are of about 7 mm and from ICON around 8 mm. The SDs are, however, much larger from ERA5 (around 40 mm) than ICON (33 mm). For ERA5, there is a very small reduction in biases for GR and GRE compared to the G solution; however, the SDs are slightly higher. For ICON, the GR solution shows better statistics than G, and GRE is further better than GR. For the flood region, the biases for ICON are much smaller (6-7 mm) than for ERA5 (around 12 mm). The biases for ICON are larger for GRE than for G, but the SDs are still slightly smaller, while for ERA5 it is vice versa. Table 5 also shows the statistics for the mapped ZTDs (using the 1/sin(el) mapping function). For both Germany and the flood region, the biases correspond well with the values in Table 1, while the SDs are around 1 mm larger, probably due to using the simple mapping function. Figure 16 shows the differences between the NWMs and GRE solution for Germany, and as already shown in Table 5, it demonstrates that the biases from both models are similar, slightly smaller for ERA5, while the SDs are much smaller for ICON. Here, also all the differences between NWMs and GRE are shown. The shape of the figures is similar for both models; however, ICON has fewer extreme values and hence the smaller SD. The differences are very dependent on the elevation angle. Thus, to remove this dependence, we plot the relative differences (divided by the absolute STD GNSS value) w.r.t. the elevation (Fig. 17) and azimuth angle (Fig. 18). The relative statistics are also included in Table 5.
As shown in Figs. 17 and 18, also the relative values are smaller from the ICON model than from ERA5. Figure 17 shows that for most of the angles, the STDs do not depend on the elevation angle except for close to the zenith. For such angles, for ERA5, the relative biases are becoming closer to zero, while for ICON, the biases close to the zenith are the largest, probably due to the limited height of the model. Figure 18 shows that there is no dependence on the azimuth angle. However, only the GRE solution is shown here. For the G solution (not shown), the dependence on the angle is still visible. This also corresponds with the findings in Wilgan et al. (2022).

Discussion and conclusions
We presented the comparisons of the three tropospheric parameters (ZTDs, tropospheric gradients and STDs) obtained from three NWMs (ERA5, ICON, GFS) and three GNSS solutions (G, GR, GRE) for a week in July 2021 that contained strong precipitation leading to deadly flooding. We considered two regions: entire Germany and the 'flood region', which was the most affected by the flooding. The comparisons of ZTDs showed that all GNSS solutions were similar when compared to the NWMs. However, the biases from GR and GRE were slightly smaller than for the G solution and the SDs were slightly larger. For the flood region, the reduction of the biases was more pronounced. Of the three NWMs and Germany, ERA5 had the smallest biases, of around 3 mm, then ICON with 3.5 mm and GFS with around 4.5 mm. The SDs were the smallest from the ICON model, of less than 10 mm and the highest from ERA5 and GFS-around 11-12 mm. For the flood region, both the biases and SDs were the smallest from the ICON model. From the comparisons for particular stations, we could see that the rise and drop in the ZTDs values were connected to high precipitation. For tropospheric gradients, for Germany, the biases were very small, below 0.01 mm for most of the cases, and very similar for all GNSS solutions. The SDs were of around 0.5 mm and were usually slightly higher for the GR/GRE solution compared to G. The SDs from ICON were smaller than from the other models. For the flood region, where high gradient values are correlated with high precipitation, the biases were an order of magnitude larger. For the STDs, we compared only the ICON and ERA5 models. For Germany, the biases from the ERA5 model were the smallest (around 8 mm), while from ICON were slightly larger (around 9 mm). However, the SDs were much smaller from ICON than from ERA5. For ICON, there was even an improvement for GRE SDs. For the flood region, the biases from ICON were much smaller (around 6-7 mm) than from ERA5 (around 13 mm). Here, there was an increase in biases and SDs for the GRE models, but a decrease in ICON SDs.
Compared to Wilgan et al. (2022), we obtained a worse agreement with ERA5; however, our comparisons were only for one week in the summer, when the water vapor variability is the highest and the values of tropospheric parameters are the largest. For this special weather event, we can also see that using the GRE solution is more advantageous for more cases than it was shown in Wilgan et al (2022). The biases from ZTDs and STDs were slightly smaller for GRE compared to the other solutions. The best agreement of GNSS data was achieved with the ICON model. We believe that this best agreement depends on few factors: (1) the horizontal resolution of ICON is much higher than of the other models, (2) we use the analysis data and not the forecast, (3) the tuning of the DWD for a specific region (here, Germany). On the other side, the GFS and ERA5 have low horizontal resolution and additionally, GFS is a free forecast, unlike ERA5, which is a reanalysis.