Assimilation of GNSS tomography products into WRF using radio occultation data assimilation operator

From Global Navigation Satellite Systems (GNSS) signals, accurate and high-frequency atmospheric parameters can be determined in all-weather conditions. GNSS tomography is a novel technique that takes advantage of these parameters, especially of slant troposphere observations between GNSS receivers and satellites, traces these signals through 10 a 3D grid of voxels and estimates by an inversion process the refractivity of the water vapour content within each voxel. In the last years, the GNSS tomography development focused on numerical methods to stabilize the solution, which has been achieved to a great extent. Currently, we are facing new challenges and possibilities in the application of GNSS tomography in numerical weather forecasting the main research objective of this paper. In the first instance, refractivity fields were estimated using two different GNSS tomography models (TUW, WUELS), which cover the area of Central Europe during 15 the period of 29 May 14 June 2013, when heavy precipitation events were observed. For both models, Slant Wet Delays (SWD) were calculated based on estimates of Zenith Total Delay (ZTD) and horizontal gradients, provided for 72 GNSS sites by Geodetic Observatory Pecny (GOP). In total, three sets of SWD observations were tested (set0 without compensation for hydrostatic anisotropic effects, set1 with compensation of this effect, set2 cleaned by wet delays outside the inner voxel model). The GNSS tomography outputs have been assimilated into the nested (12and 36-km horizontal 20 resolution) Weather Research and Forecasting (WRF) model, using its three-dimensional variational data assimilation (WRFDA 3DVar) system, in particular its radio occultation observations operator (GPSREF). As only total refractivity is assimilated in GPSREF, it was calculated as the sum of the hydrostatic part derived from the ALADIN-CZ model and the wet part from the GNSS tomography. We compared the results of the GNSS tomography data assimilation to the radiosonde (RS) observations. The validation shows the improvement in the weather forecasting of relative humidity (bias, standard 25 deviation) and temperature (standard deviation) during heavy precipitation events. Future improvements to the assimilation method are also discussed.

from GNSS tropospheric parametersalso known as GNSS tropospheric tomography.This technique is one of the most promising since it uses Slant Wet Delay (SWD) observations together with the principles of tomography to obtain not only the amount of water vapour, but also its distribution in space and time.For the inversion of the equation system Flores et al. (2000) were using singular value decomposition (SVD) methods together with constraints imposed on the system of equations (smoothing, boundary and vertical constraints).In the following years, a number of refinements in the technique have been implemented.Inclusion of supplementary wet refractivity data from external data sources into a functional model (Bender et al., 2011;Rohm et al., 2014;Benevides et al., 2015;Möller, 2017) resulted in stabilization of the solution without need of smoothing constraints, which is especially important during severe weather phenomena with high variability of water vapour in the troposphere.A new approach of parametrization of the domain, using trilinear and spline functions instead of constant values inside each volume element (voxel) showed the substantially smaller maximum error of the solution (Perler et al., 2011;Ding et al., 2018).Another approach, using the Kalman Filter instead of Least Squares (Rohm et al., 2014), led to better responsiveness to the severe weather.A number of experiments have been carried out showing positive results for detecting heavy precipitation events (Troller et al., 2006 Chen et al., 2017).Besides the development of technical aspects, which was widely performed in the last years, GNSS tomography should be also tested on its ability to fulfil the requirements of up-to-date weather prediction methods, i.e. on its suitability for assimilation into the NWP models (Innes and Dorling, 2012).The first research on this topic was carried out by Möller et al. (2015).Its aim was to transform wet refractivities derived by GNSS tomography into humidity and temperature profiles, and then assimilate them into the AROME NWP model for a selected test period.Verification was made using surface station measurements, radio sounding data, and precipitation analysis.Tomography data assimilation results compared with the results of ZTD data assimilation, indicated a significantly larger impact of the new technique, especially within the first 6 to 12 hours of the forecast (a drying effect in the AROME forecast field).As the results clearly showed the potential of 3D refractivity observations, another case study has been performed by Trzcina and Rohm (2018), concerning assimilation of a Near-Real-Time (NRT) tomographic solution into WRF model using an operator dedicated to radio occultation (RO) observations of total refractivity.
Comparison with several external observations has shown the positive impact on 6 to 12 hours forecast of humidity in autumn, but not in summer.Such assimilation-oriented tomography data analyses are essential for further research on the utilization of GNSS troposphere tomography output as a valuable data source for NWP forecasts.
In this work, we present the assimilation of a 3D wet refractivity field derived by GNSS tomography technique, into the WRF model using the GPSREF observation operator.The experiment was performed during a heavy precipitation event in Central Europe in order to investigate the impact on the weather forecast.Since two tomographic models (TU Wien, WUELS) and three different approaches of SWD data calculation were used, the results of several tomographic approaches were discussed in terms of assimilation impact.The structure of the paper is as follows: Section 2 describes the derivation of slant wet delays from GNSS measurements.Section 3 is dedicated to the methodology of GNSS troposphere tomography and provides a detailed description of the models used for tomography processing.Section 4 introduces the methodology of

GNSS slant wet delays
GNSS tropospheric estimates (1h ZTDs and 6h horizontal gradients), as provided by Geodetic Observatory Pecny (GOP) for 72 GNSS sites in Central Europe (see Fig. 3), were utilised for the computation of Slant Wet Delays (SWDs).For more details about the GNSS processing strategy at GOP, the reader is referred to (Dousa et al., 2016, Section 5.1).
In a first step, the tropospheric parameters were linearly interpolated over time to 0, 6, 12, 18 UTC.Afterwards, for each epoch the Zenith Hydrostatic Delay (ZHD) was computed by means of Saastamoinen model (Saastamoinen, 1972) and removed from the ZTD estimates.The required air pressure values have been derived from meso-scale ALADIN-CZ 6 hour forecast data, as provided by the Czech Hydrometeorological Institute (CHMI), and spatially interpolated to each GNSS reference site.

𝑍𝑊𝐷 = 𝑍𝑇𝐷 − 𝑍𝐻𝐷.
( The Zenith Wet Delays (ZWDs) and the horizontal gradients (GN, GE) were mapped into direction of the GPS and GLONASS satellites in view above 3° elevation angles.Therefore and for highest consistency with GNSS data processing, the VMF1 mapping function (Böhm et al., 2006) was used for mapping (  ) of the ZWDs and the (Chen and Herring, 1997) gradient mapping function (  ) was applied for mapping of the horizontal gradient parameters as follows: The elevation () and azimuth angles () of each satellite in view have been computed from broadcast ephemerides.The resulting dataset of SWDs is hereafter referred to as 'set0'.
Based on set0 two additional, more refined datasets were derived, whereby 'set1' compensates for hydrostatic anisotropic effects and 'set2' was cleaned by wet delays 'outside' the inner voxel model (see Table 3 for further details about the voxel model).In both cases, 2D ray-tracing through ALADIN-CZ model level data (6 hour forecast data) was carried for the determination of the delay corrections.For more details about the applied ray-tracer, the reader is referred to (Möller, 2017).
For set1, i.e. for compensation of hydrostatic asymmetric effects, ray-traced delays were determined for each GNSS satellite in view, and for equidistant azimuth angles (separated by 30°).From the obtained set of ray-traced hydrostatic delays, the mean hydrostatic delay was computed and removed from the ray-traced hydrostatic delay in direction of the satellite in view.
Figure 1 shows the resulting hydrostatic asymmetric delay, obtained for all observations within the study period (29 May - 14 June, 2013).In a final step, the hydrostatic asymmetric delays were removed from set0.Assuming that the ALADIN-CZ model pressure error is small, the obtained SWDs (set1) are widely free from hydrostatic effects.Furthermore, for set2 the ray-tracer was applied for determining intersection points between the GNSS signal paths and the inner voxel model boundaries.In case, the signal enters the inner voxel model not through the top layer but through any lateral surface, the outer tropospheric wet delay (between intersection point and the upper rim of the troposphere at about 13.5 km height above mean sea level) is determined by ray-tracing and removed from set0.

GNSS tomography
For conversion of precise integral measurements (like SWDs) into three-dimensional structures, a technique called tomography has been invented.According to Iyer & Hirahara (1993) the general principle of tomography is described as follows: where   is the projection function, () is the object property function and  is a small element of the ray path  along which the integration takes place.In order to adapt the tomography principle for GNSS tropospheric delay parameters and to find a solution for Eq.(3), first the troposphere is discretised in volume elements (short: voxels) in which the index of refraction is assumed as constant and the ray path is assumed as straight line.Further, by replacing   with  and () by wet refractivity (  ), the basic function of GNSS tomography reads: whereby  , is the constant wet refractivity and   is the ray length in volume element k.In matrix notation Eq. ( 4) reads: where  is the observation vector of size (l, 1),  is the vector of unknowns of size (m, 1) and  is a matrix of size (l,m) which contains the partial derivatives of the slant wet delays with respect to the unknowns.Since Eq. ( 5) is linear, the partial derivatives of SWD are the ray lengths (  ) in each voxel k.A solution for  is obtained by inversion of matrix .
Unfortunately in GNSS tomography, matrix  is usually not of full rank (i.e. has zero singular values).In consequence Eq. ( 6) is ill-posed, i.e. not uniquely solvable.In order to find a solution for Eq. ( 6), truncated singular value decomposition (TSVD) methods are applied (Xu, 1998).Therefore, matrix  is split into three components: where  and  are orthogonal matrices and matrix  (l, m) is a diagonal matrix.The first m diagonal elements of  are the singular values ( , ), all other elements are zero.The ratio between largest and smallest singular value defines the condition number () of matrix .
A condition number near '1' indicates a well-conditioned matrix, a larger condition number indicates an ill-conditioned problem.As consequence the tomography solution is sensitive to observation errors, changes in the observation geometry but also to the solving strategy and its parameters defined within the analysis.In the following, the two processing strategies as used for parameter estimation are described more in detail.The general, underlying tomography settings and input data are summarised in Table 3.

Least-squares solution (TUW)
The TU Wien tomography solution was calculated using the ATom software package (Möller, 2017).It is based on a least-squares adjustment, whereby each observation type is processed individually as partial solution (see Möller, 2017 Eq. 7.46 ff.).This approach allows for proper Truncated Singular Value Decomposition (TSVD), whereby an eigenvalue threshold of 0.006 km², found by l-curve technique, was set for singular value selection.
Weighting of the GNSS slant observations was based on the elevation angle (sin 2 ), with a standard deviation of 2.8 mm for the zenithal observations (ZTD).A height-dependent weighting model was applied to the a priori wet refractivity information, whereby the height-dependent standard deviations were computed by comparison of the a priori data with radiosonde data and interpolated to the 15 height levels of the voxel model (see Table 1).
The first two TU Wien solutions are based on SWD set0 and set1.Therefore, the voxel model boundaries were defined according to Table 3 and wet refractivity was estimated in both, the inner and outer domain.In consequence, all available SWDs could be used for estimation of the wet refractivity fields (sol0 and sol1).The third TU Wien solution is based on the SWDs of set2.For processing, the outer voxel model was removed since the outer tropospheric delay was already considered in preparation of set2.In consequence, all available SWDs could be used and wet refractivity was estimated for the inner voxel model (sol2).Quality control was based on analysis of the SWD residuals ().
=   −   •   (9) Those residuals larger than 120 times the RMS of the residuals were discarded.The processing stopped after 10 iterations or before, if the change in RMS was lower than 3 percent or if the RMS was lower than 0.5 mm.

Kalman filter solution (WUELS)
Estimation of the wet refractivity distribution in the WUELS solution was performed using TOMO2 software where Φ denotes a state transition matrix, indicating predicted changes of the wet refractivity between consecutive epochs;  is a dynamic disturbance noise matrix.The state transition matrix is a diagonal matrix, set based on a priori information for both epochs where: The dynamic disturbance noise matrix was set based on mean wet refractivity changes between epochs, calculated from ALADIN-CZ model data for the whole period.The noise was calculated separately for each height level of the tomographic domain (see Table 3).After the prediction step, corrections for the state vector and its covariance matrix were computed from the observations.Hereby, the Kalman gain matrix K was defined as follows: where  stands for a vector of observations and  for its covariance matrix.The vector of observations consists of two parts, the SWD observations and the a priori information derived from ALADIN-CZ model data.Errors of the slant delays   are dependent on the elevation angle  and were calculated using the following mapping function: whereby   = 10 is the assumed ZWD error (Kacmarík et al., 2017).Errors of the a priori information were computed as for the TUW solution by comparison with radiosonde observations (stations 10548 and 10771), separately for each height of tomographic domain (Table 1).
For tomography processing, all available SWD observations above 3° elevation angle were taken into account and the signal paths were assumed as straight-lines.The WUELS tomography solutions (sol0 and sol1) are based on the SWD 5 observations from set0 and set1, respectively.The voxel model boundaries were defined according to Table 3 (both inner and outer voxel model taken into account).A main quality indicator in the WUELS model was based on a comparison of the SWD residuals (see Eq. 9).Observations with residual values   larger than the standard deviation multiplied by 3 were removed from the solution.The filter process was stopped after five iterations.Thereby, stations with a number of removed observations higher than 150 (for the whole period) were rejected.part of the model.Nesting enables also to run the model in a high-resolution very small domain, avoiding a mismatch time and spatial lateral boundary conditions (Skamarock et al., 2008).A coarse grid and the fine grid can interact, depending on the chosen feedback option.Nested grid simulations can be produced using either 1-way nesting or 2-way nesting.In both, the 1-way and 2-way simulation modes, the fine grid boundary conditions (i.e., the lateral boundaries) are interpolated from the coarse grid forecast.In a 1-way nest, the only information exchange between the grids comes from the coarse grid to the fine grid.In the 2-way nest integration, the fine grid solution replaces the coarse grid solution for coarse grid points that lie inside the fine grid.This information exchange between the grids in both directions (Skamarock et al., 2008).
In this research, the WRF model includes two domains, i.e. a parent domain, which spans the area of almost the whole Europe, with 36 km horizontal spacing, and a nested (local) domain within the area of Central Europe, with 12 km horizontal spacing (Fig. 3).The information between the domains flow in both directions   (Dudhia, 1989), and Kain-Fritch scheme for the cumulus parametrization (Kain, 2004).

Assimilation operator
Data assimilation is the technique by which observations are combined with an NWP product (the first guess or background forecast) and their respective error statistics to provide an improved estimate (the analysis) of the atmospheric (or oceanic, Jovian, etc.) state.Variational (Var) data assimilation achieves this through the iterative minimization of a prescribed cost (or penalty) function.The aim of a variational data assimilation scheme is to find the best least-square fit between an analysis field  and observations () with an iterative minimization of a cost function (): In Eq. 17,  is a vector of analysis,   is the forecast or background vector (first guess),  is an observation vector,  is the background error covariance matrix,  is an observation operator, and  is the observation covariance matrix.The observation operator , which can be non-linear, converts model variables to observation space.Differences between the analysis and observations/first guess are penalized (damped) according to their perceived error.
In order to assimilate observations into the WRF data assimilation system, a special function, i.e. data assimilation operator (), should be used.Although the WRF DA system does not provide a direct assimilation operator for the GNSS tomography wet refractivity fields, it is possible to apply a radio occultation operator (GPSREF) for the GNSS tomography output data.The GPSREF assimilation operator is used to map model variables to refractivity space (Cucurull et al., 2007): where  is the total atmospheric pressure (hPa),  is the atmospheric temperature (K), and   is the partial pressure of water vapor (hPa).The GPSREF operator enables assimilation of total refractivity ()at each observation height.Since GNSS tomography provides   only, the hydrostatic (dry) component of refractivity ( ℎ ) has to be modeled.Therefore we use meteorological parameters (total atmospheric pressure and atmospheric temperature) provided by ALADIN-CZ model.
The ALADIN-CZ model is a local area, hydrostatic model provided by Czech Hydrometeorological Institute (CHMI).
Model simulations run at 4.

Assimilation settings
The 3DVar assimilation of five sets of tomography data (TUW set0, TUW set 1, TUW set2, WUELS set0, WUELS, set1) and the control run of the base weather forecast (without assimilation) has been performed.In order to allow the NWP model to approach its own climatology after being started from the initial conditions (NCEP final data), the assimilation was performed six hours after model run (spin-up time).The analysis has been performed at 00, 06, 12 and 18 UTC.Each analysis includes 18 hours forecast with one hour resolution.Table 2 shows the data assimilation settings.In the process of assimilation, only the data from the voxels crossed by the GNSS signal were used.

Case study
Our studies are based on the complex dataset collected within the European COST Action ES1206, as described in detail by Dousa et al. (2016).The study area comprises central parts of Europe (Fig. 2).The study period (29 May to 14 June, 2013) covers the events of heavy precipitation, which finally led to severe river flooding in all catchments on the Alpine north-side (Switzerland, Austria, and southern Germany) and of the mountain ranges in Southern and Eastern Germany as well as in the Czech Republic.
As reported by Grams et al. (2014), the core period of the heavy precipitation occurred from 31 May 2013, 00:00 UTC to 3 June 2013, 00:00 UTC.Within these 72 hours, accumulated precipitation of 50 mm was reported in regions of Eastern Switzerland, the Austrian Alps and Czech Republic, and exceeded even 100 mm in several regions of the Northern Alps.The events of heavy precipitation are connected to the presence of three cyclones: "Dominik", "Frederik", and "Günther".These cyclones formed over (South-) Eastern Europe as a "cluster" with very similar tracks, and were rather shallow systems with relatively high values of minimum sea level pressure.These systems unusually track westward, maintaining a northerly flow against the west-east oriented mountain ranges in Central Europe.However, within this period their movement was blocked by an anticyclone located over the North Atlantic, which forced the cyclones to form an equatorward conveyor belt.
The GNSS tomography solutions covers the region of East Germany and western parts of Czech Republic, including the Erzgebirge (see Fig. 2).Within this area, the equatorward flowing warm air masses started to lift up, which goes along with a local maximum of precipitation of more than 75 mm within the core period of heavy precipitation.For tomography processing, the study area was divided into an inner and an outer voxel model.The parameters of the tomography model and model settings are summarized in Table 3. From the benchmark dataset, 72 GNSS sites could be identified within lay within the inner model, with an average distance of 48 km and altitude ranging from 70 m to 885 m.

Tomography results
In the following, the TUW and WUELS tomography solutions are compared with each other and validated against GNSS ZWDs and radiosonde measurements.

Intra-technique comparison
Over the study period of 17 days, in total 68 tomography epochs (four solutions per day) have been processed.
Based on the wet refractivity differences (TUW minus WUELS) for each voxel, bias and standard deviation were computed.A similar analysis has been carried out for TUW solution 1 and 2 (Fig. 5).The major differences between both solutions are caused by the different approaches for compensation of the outer voxel delay.While in solution 1 the tropospheric delay in the outer voxel model was estimated, in solution 2 the outer delay was removed beforehand by raytracing through ALADIN-CZ six hour forecast data.Largest offsets between both TUW solutions appear in Northern parts of the study area (voxels columns #55-63).In this part, the tropospheric delay is systematically overestimated (compared to ray-traced delays) in the outer voxel model and in consequence, underestimated in the inner voxel model.This leads to the positive bias as visible in Fig. 5.In all other parts of the voxel model, the differences are widely averaged out over the study period of 17 days.However, the standard deviation shows that variations over time cannot be avoided.Especially between the 31 st of May 18 UTC and the 4 th of June 18 UTC, larger standard deviations were detected (not shown).These variations are caused by changes in the observation geometry but also by changes in atmospheric conditions not described by the forecast data (error source in set2).

Time series of integrated zenith wet delays
From the obtained wet refractivity fields, time series of ZWDs were computed for the 72 GNSS sites within the 5 study area by vertical integration using Eq. 19,  = 10 −6 ∑   • ℎ    0 (19) where  0 is the height of the GNSS site and   is the height of the voxel top (at about 13.5 km height above mean sea level).
Beforehand,   was horizontally interpolated from adjacent voxel centre points to GNSS site and the vertical resolution was further increased to 20 m. Figure 6 shows the derived ZWD time series for the entire study period of 17 days with 6 h 10 temporal resolution, exemplary for GNSS site WTZR.Both tomography solutions are more consistent with the GNSS derived ZWDs than the ALADIN-CZ data.The ZWDs derived from ALADIN-CZ 6 hour forecasts are occasionally biased, a few hours 'ahead' or 'behind' the GNSS derived 5 ZWDs.The WUELS solution is negatively biased, while discrepancies between GNSS and TUW are smaller and equally distributed in positive and negative values.Table 4 shows the statistic of the ZWD differences for GNSS site WTZR but also the mean values over all 72 GNSS sites within the study domain.forecast.Especially the standard deviation of the ZWD differences can be reduced if GNSS tomography is applied.
However, significant differences are visible between both tomography solutions.The WUELS solution tends to overdetermine the water vapour content in the atmosphere slightly, this goes along with a higher variability of the ZWD differences (about 2 times larger standard deviation than for the TUW solution).This shows that weighting of a priori data and applied quality control methods have a great impact on the results of tomography solutions.

Validation of tomography results with radiosonde data
Radiosonde measurements of temperature (  ) and dew point temperature (   ) were downloaded from  (20) and wet refractivity with  2 ′ = 22.9744/ℎ and  3 = 375463 2 /ℎ , see Rüeger, 2002, best average.For comparison, the obtained profiles of wet refractivity (  ) were linearly interpolated to the voxel centre heights (see Table 1).
Figure 7 shows wet refractivity profiles as obtained at radiosonde site RS10548 on the 6 th of June 2013, 00 UTC and radiosonde site RS10771 on the 13 th of June 2013, 00 UTC.Both plots were selected to highlight typical characteristics of the tomography solution.Due to proper weighting, the tomography solution can correct deficits in the a priori model (see the lower 2 km of TUW solution 1 in Fig. 7 left and the WUELS solution between 2 and 5 km height in Fig. 7   Fig. 8 shows the corresponding bias and standard deviation over all analysed epochs (34 for RS10548 and 68 for RS10771), separately for each radiosonde site.While RS10548 is located nearby a GNSS site, RS10771 lies within a voxel 5 in which no GNSS site is located.Nonetheless, the quality of the tomography solution does not vary significantly according to the locations of RS stations.Up to 2 km height, the statistics for both tomographic models are on the same level.In the upper parts of the troposphere, the RS -TUW comparison shows similar accuracy as the RS -ALADIN-CZ data, whereas the standard deviation of the WUELS model is noticeably higher.Figure 9 shows the differences between solution 1 and 2 of the TUW model, separately for each radiosonde site.In fact, the strategy for compensation of the outer delay has only a small impact on the tomography solution.Thereby, the impact is independent from the location of the radiosonde site within the inner voxel model but mostly related to the tomography settings.The largest differences are visible in the lower 2 km of the atmosphere.However, the overall impression is that solution 1 (estimation of the tropospheric delay in the outer voxel model) provides slightly better results than solution 2 (compensation of the outer delay by ray-tracing through ALADIN-CZ 6 hour forecast data).This is most likely related to the quality of the weather model forecast data during the period of extreme precipitation.Nevertheless, since both solutions are rather close to each other, especially with respect to standard deviation, only small improvement is expected by ray-tracing through more reliable weather forecast data.7 Assimilation results

Diagnosis output
The WRFDA system and the GPSREF operator are equipped with a quality control diagnostic tool, which allow for verification of all input data before assimilation.As a result, not all of the refractivity observations are actually assimilated into the model.Figure 10 presents a percentage of successfully assimilated observations for each tomography set, as a function of height.The height range with the largest number of assimilated observations is between 4 and 10 km.For these heights, more than 90% observations of the TUW solutions and about 50-80% of those form the WUELS solution were assimilated.Below 4 km, the percentage of assimilated observations grows systematically with height, from 0% at surface to 70% for TUW and 40% for WUELS.Above 10 km, no observations were assimilated.There is a noticeable difference between the numbers of assimilated observations for the two models.At all heights, about 5%-40% more observations were assimilated by the TUW than the WUELS model.Since the comparison of tomographic observations with radiosonde data showed that in general the TUW solution has smaller errors than the WUELS, the number of observations that passed the quality control is connected to the quality of the tomographic data.Apart from the number of successfully assimilated observations, the reason of the rejection was also studied.In the quality control diagnostics of the GPSREF operator, each observation is assigned to one group, based on the information of acceptance or rejection (and its cause).Figure 11 presents the results of the quality control check, separately for each group and for each tomography solution as function of height in colour lines.For both tomography solutions, a large number of observations is assigned to a group 0 (orange line) that denotes successfully assimilated data.The number of observations from this group corresponds well with the one already presented in Fig. 10.The quality flags with numbers -88, -77, and -3 are the general flags of the WRFDA system, whereas the numbers from -31 to -36 are the flags assigned by the GPSREF operator.
The general flags with numbers -88, -77, and -3 denote respectively: missing data, data laying outside of horizontal domain, data failing maximum error check.Observations assigned to the first group (blue line) occur in the surface layer only, in number of about 2000 observations for the TUW model and 1200 for the WUELS.The second group (yellow line) includes observations from the two highest layers (11.260 km and 12.814 km).The number of observations assigned to the third group (green line) is about 0 at all heights, only for the WUELS model at heights 4 -8 km the number slightly grows up to about 50 observations.

Assimilation results at analysis time
The WRF model outputs with assimilated tomographic data were compared to the radiosonde measurements of relative humidity, temperature, and wind speed (Fig. 12).Systematic errors and standard deviations are presented as a function of pressure.Relative humidity slightly improved in terms of systematic error in the boundary layer (900 -700 hPa) by 1% -2%.In higher layers, bias of relative humidity increased (by 1% -5%) but standard deviation decreased (by 1% -2%) in the case of assimilation of the TUW data (especially within the pressure range of 700 to 400 hPa).Bias of temperature is growing after assimilation in the pressure range of 900 -400 hPa, especially for the WUELS data assimilation (by 0.2 -0.5°C) but also for the TUW model (by 0.1°C).In the higher troposphere, the bias is decreased when compared with the base run, by 0.5°C for WUELS and 0.2°C for TUW.Impact of assimilation on the standard deviation of temperature is neutral in the lower troposphere and negative for the higher parts (600 -250 hPa for WUELS and 400 -250 hPa for TUW).For the wind speed, there is a small impact of data assimilation on the values of bias, but in terms of standard deviation there is a noticeable improvement for pressure range of 600 -400 hPa (0.2 m s -1 ).In the lower troposphere, standard deviation of wind speed decreases after assimilation of the TUW data (by 0.1 m s -1 ) but increases after assimilation of the WUELS output (by 0.2 m s -1 ).

Assimilation output results at simulation time
Based on the radiosonde observations (10548, 10771), we calculated the bias and standard deviation between the radiosonde data and the model run.The results of the statistics for the weather forecast with assimilated tomographic output data have been compared to the statistics for the base run (raw run) of the model.This comparison has been conducted 6 hours after assimilation (i.e., 12 hours since model start) for the whole tested period (29 May 2013 00:00 UTC -13 June 2013 18:00 UTC), in the location of radiosonde stations (10548, 10771).We analysed three meteorological parameters: relative humidity, temperature, and wind speed (Fig. 13).
In order to assess the impact of the assimilation of the GNSS tomography outputs on the weather forecast, we subtracted the absolute values of statistics for the case of assimilation and the base model run.If the difference between the statistics is negative, it means that the assimilation of the GNSS tomography output decreases the value of bias and/or standard deviations between the model and radiosonde observation.It implies that the model after assimilation becomes closer to the real observations.In the beginning of the tested period, i.e., during the first two days, the weather in the Central Europe area was rather calm, thus for all assimilation cases the differences between the statistics for the weather forecast after tomography data assimilation and the base run are on similar level.In case of a relative humidity assimilation, in time of the heavy precipitation events (01-03 June 2013), the absolute differences of bias between the models after assimilation of the tomography data and the base run have negative values.The highest positive impact (i.e., decrease of the value of bias) occurs in case of assimilation of the WUELS tomography outputs and exceeds -4.5 % for both bias and the standard deviation differences.In case of the TUW GNSS tomography data assimilation, the differences in the statistical numbers reach around -1 % for set2 and around -2.5 % for set0 and set1.5 We observe that the assimilation of the tomography outputs does not have any significant influence on the other examined meteorological parameters.In case of the temperature, the differences hold the level up to 0.1 ℃ for both bias and standard deviations.Since 04 June 2013, the tomography wet refractivity assimilation gives mainly a positive effect for the differences in case of WUELS data (set0, set1), whereas the differences for the TUW data (set0, set1, set2) are close to zero.
The statistics for the wind speed look very similar; the assimilation of the tomography data modifies barely the difference 10 between results of statistics, up to 0.5 m/s for both bias and standard deviation differences.
Atmos.Meas.Tech.Discuss., https://doi.org/10.5194/amt-2018-419Manuscript under review for journal Atmos.Meas.Tech.Discussion started: 12 February 2019 c Author(s) 2019.CC BY 4.0 License.tomographic data assimilation into the WRF model; followed by a description of the meteorological situation under study in Section 5.In section 6 and 7 the results of tomography and of the assimilation experiments are analysed.The main conclusions are presented in section 8.

Figure 1 :
Figure 1: Hydrostatic component of azimuthal asymmetry, derived from the ALADIN-CZ model level data by ray-tracing.
(2-way nested run).The vertical resolution of model includes 35 pressure levels, whereby the model top is defined at 50 hPa air pressure.The model background (initial and boundary conditions) is provided by the NCEP FNL (Final) operational global analysis data with 1° x 1° horizontal resolution and 26 vertical layers from the National Centers for Environmental Prediction (NCEP).These data are available on http://rda.ucar.edu/datasets/ds083.2/.

Figure 2 :
Figure 2: The location of the WRF domains (d01 = coarse grid with 36 km horizontal spacing, d02 = fine grid with 12 km horizontal spacing).
7 km x 4.7 km horizontal resolution and 87 hybrid-model levels.NWM 3D analysis fields are provided in GRIB format for 00, 06, 12 and 18 UTC together with 6 hours forecast fields with one hour resolution Atmos.Meas.Tech.Discuss., https://doi.org/10.5194/amt-2018-419Manuscript under review for journal Atmos.Meas.Tech.Discussion started: 12 February 2019 c Author(s) 2019.CC BY 4.0 License.(Kacmarik et al, 2017).In this study, the analysis data of the ALADIN-CZ model is used to: 1) Compute pressure values at the height of GNSS station to reduce the hydrostatic part from the Slant Total Delays (STDs); 2) Compute temperature and pressure at the height of each   observation in order to model the hydrostatic part of refractivity.In order to calculate the hydrostatic component of refractivity, the meteorological parameters provided by ALADIN-CZ are interpolated to each observation height.Finally, total refractivity field, i.e. the sum of  ℎ and   , is assimilated.

Figure 4
Figure 4 shows the results for the 63x15 voxels of the inner voxel model, whereby voxel number 1 indicates the lower South-West corner and the voxel number 63 the North-East corner of the domain.Since both tomography solutions (TUW and WUELS) are based on same input data (set 1), similar wet refractivity fields are obtained.The largest discrepancies are visible within the boundary layer, i.e. the lowest 1-2 km of the atmosphere, and for the layer of 3-6 km height.The TUW solution tends to produce larger wet refractivity fields in specific voxels near surface (#19, 20, 24, 25, 35, 38, 41, 48, 50,…),in return slightly lower values in the 2-3 layers above, and then slightly higher values for the heights of 3-6 km (up to 5 ppm).Above the 6 km altitude, the absolute bias and the standard deviation between TUW and WUELS solutions decrease significantly (below 3 ppm).The discrepancies between both models are caused by different approaches of the tomography settings.The tomography models differ in application of quality indicators.In the WUELS approach, the GNSS stations that

Figure 4 :
Figure 4: Bias (left) and standard deviation (right) of the differences in wet refractivity [ppm] between TUW and WUELS solution 1. Analysed period: 29 May -14 June 2013 (68 epochs).Grey triangles denote GNSS stations removed from the specific voxels of the WUELS solution (one triangle stands for one rejected GNSS station).

Figure 6 :
Figure 6: GNSS derived ZWD and integrated ZWD time series at GNSS site WTZR.The top plot shows the absolute values and the bottom plot highlights the ZWD differences with respect to the GNSS derived ZWDs.

Table 4 :
Statistics of the differences in ZWD for GNSS site WTZR and all 72 GNSS sites within the study domain.The reference solution (GNSS) was obtained by parameter estimation using GNSS phase measurements.The ALADIN-CZ, the TU Wien and 10 WUELS solutions were computed by vertical integration through the wet refractivity fields: Analysed period: 29 th of May until 14 th of June 2013.GNSS minus ALADIN-CZ GNSS minus TU Wien GNSS minus WUELS bias(dZWD) stddev(dZWD) bias(dZWD) stddev(dZWD) bias(dZWD) stddev(dZWD) WTZR -0.5 mm 7.2 mm -0.9 mm 2.6 mm -2.5 mm 3.9 mm Atmos.Meas.Tech.Discuss., https://doi.org/10.5194/amt-2018-419Manuscript under review for journal Atmos.Meas.Tech.Discussion started: 12 February 2019 c Author(s) 2019.CC BY 4.0 License.With respect to ZWD, both tomography solutions are closer to the GNSS solution than the a priori ALADIN-CZ 6 hour right).Sharp changes in wet refractivity can be recovered from TUW tomography solution in the boundary layers and from the WUELS solution in the mid-troposphere (related to weighting model).The tomography solutions are still affected by reconstruction errors (artefacts in the Nw profiles as visible in Fig. 7 left between 2 and 5 km height).Atmos.Meas.Tech.Discuss., https://doi.org/10.5194/amt-2018-419Manuscript under review for journal Atmos.Meas.Tech.Discussion started: 12 February 2019 c Author(s) 2019.CC BY 4.0 License.

Figure 8 .
Figure 8. Statistics of the differences in wet refractivity at radiosonde site RS10548 (left) and RS10771 (right).

Figure 9 :
Figure 9: Statistics of the differences in wet refractivity between radiosonde data and various tomography solutions.Differences between solution 1 and 2 are related to the strategy for compensation of the outer tropospheric delay.

Figure 10 :
Figure 10: Number of successfully assimilated observations for TUW (green) and WUELS (blue) models in a function of height.
Atmos.Meas.Tech.Discuss., https://doi.org/10.5194/amt-2018-419Manuscript under review for journal Atmos.Meas.Tech.Discussion started: 12 February 2019 c Author(s) 2019.CC BY 4.0 License.The quality flags assigned by the GPSREF operator (from -31 to -36) are connected with the values of assimilated refractivity data.The first type of diagnostics is based on the comparison between the assimilated observations and the background values of refractivity (Cucurull et al. 2007).The discrepancies between the two refractivity values should not be larger than 5% (below 7 km) or 4% (7 -25 km) of the mean value.Observations that do not meet this requirement are assigned to the group -31 (below 7 km, red line) or -32 (above 7 km, purple line).Additionally, if an observation gets the flag -31, all observations in the same vertical profile (same latitude and longitude) below that observation, are also assigned to the flag -31 and they are not assimilated into the model.For the TUW model, the largest number of observations with flag -31 is at height of 0.675 km (about 3000 observations); this number systematically decreases with height, to about 0 observations at 4 km.For the WUELS model, the largest discrepancies between the observations and the background occur between the heights of 0.675 km and 2.972 km (2000 -3000 observations); no significant discrepancies are noticed above 6 km.The second type of quality check inside the GPSREF operator is based on a refractivity lapse rate (Poli et al., 2009).The -34 flag (brown line) is assigned to the observations where   is smaller than -50 km -1 , whereas the -35 flag (pink line) indicates observations where an absolute value of  2   2 is larger than 100 km -2 .For the TUW model there are no observations rejected from the assimilation based on the lapse rate of refractivity.It indicates internal coherency of the model's output.In the WUELS model, for each layer above 6 km, more than 1000 observations were rejected from the assimilation process based on the refractivity lapse rate.The last type of quality check (flag -36, grey line) is based on the discrepancies between observations and background, as proposed by Cucurull (2010).In the TUW model, these flags occur mainly in the bottom parts of the troposphere (2 km), whereas in the WUELS model inconsistency between model and background is noticed also for the higher parts (5 -8 km).Atmos.Meas.Tech.Discuss., https://doi.org/10.5194/amt-2018-419Manuscript under review for journal Atmos.Meas.Tech.Discussion started: 12 February 2019 c Author(s) 2019.CC BY 4.0 License.

Figure 11 :
Figure 11: Number of flagged observations as function of height (flags explanation in the text) -separately for TUW sol1 (left) and WUELS sol1 (right).

WRF Assimilation operator and settings 4.1 WRF model
The Weather Research and Forecasting (WRF) model is a mesoscale numerical weather prediction system designed for both, atmospheric research and operational forecasting needs.It provides 1) two different numerical coresthe Non-In the WRF model, a nest simulation can be performed.The nest is a finer-resolution model run, which can be embedded simultaneously within a coarser-resolution (parent) model run, or run independently as a separate model forecast.The nested solution allows for the low-resolution model run within the parent domain and the high-resolution model run in the nested domain, what helps enhance the efficiency of the model and reduces the time of calculation of the high-resolution Atmos.Meas.Tech.Discuss., https://doi.org/10.5194/amt-2018-419Manuscript under review for journal Atmos.Meas.Tech.Discussion started: 12 February 2019 c Author(s) 2019.CC BY 4.0 License.

Table 3 : Applied tomography model settings.
TUW: outlier tests based on the post-fit residual RMS WUELS: filtration of dSWD based on synthetic data for detecting outlier observations