A comparison between the GNSS tomography technique and the 1 WRF model in retrieving 3 D wet refractivity field in Hong Kong

Water vapor plays an important role in various scales of weather processes. However, there are 7 limited means to monitor its 3-dimensional (3D) dynamical changes. The Numerical Weather Prediction 8 (NWP) model and the Global Navigation Satellite System (GNSS) tomography technique are two of the 9 limited means. Here, we conduct an interesting comparison between the GNSS tomography technique and 10 the Weather Research and Forecasting (WRF) model (a representative of the NWP models) in retrieving Wet 11 Refractivity (WR) in Hong Kong area during a rainy period and a rainless period. The GNSS tomography 12 technique is used to retrieve WR from the GNSS slant wet delay. The WRF Data Assimilation (WRFDA) 13 model is used to assimilate GNSS Zenith Tropospheric Delay (ZTD) to improve the background data. The 14 WRF model is used to generate reanalysis data using the WRFDA output as the initial values. The radiosonde 15 data are used to validate the WR derived from the GNSS tomography and the reanalysis data. The Root Mean 16 Square (RMS) of the tomographic WR, the reanalysis WR that assimilate GNSS ZTD, and the reanalysis 17 WR that without assimilating GNSS ZTD are 6.50 mm/km, 4.31 mm/km and 4.15 mm/km in the rainy period. 18 The RMS becomes 7.02 mm/km, 7.26 mm/km and 6.35 mm/km in the rainless period. The lower accuracy 19 in the rainless period is mainy due to the sharp variation of WR in the vertical direction. The results also 20 show that assimilating GNSS ZTD into the WRFDA model only slightly improves the accuracy of the 21 reanalysis WR and that the reanalysis WR is better than the tomographic WR in most cases. However, in a 22 special experimental period when the water vapor is highly concentrated in the lower troposphere, the 23 tomographic WR outperforms the reanalysis WR in the lower troposphere. When we assimilate the 24 tomographic WR in the lower troposphere into the WRFDA model, the reanalysis WR is improved. 25


Introduction
Water vapor (WV), mostly contained in the troposphere, plays an important role in various scales of atmospheric processes.But due to its active nature, there are limited models and techniques that can accurately describe or monitor its 3-dimensional (3D) dynamical changes (Rocken et al., 1993).
The development of Global Navigation Satellite System (GNSS) technique and the densely deployed GNSS receivers provide us the opportunity to monitor the WV field in near real time.When GNSS signal travels through the neutral atmosphere, it undergoes time delay and bending due to the atmospheric refractivity.This effect is usually called the tropospheric delay in the GNSS field (Altshuler, 2002).The tropospheric delay is usually considered as the product of the zenith delay and the mapping function (Lanyi, 1984;Niell, 1996).
The Zenith Tropospheric Delay (ZTD) consists of two parts: the hydrostatic part and the wet part.The wet delay is mainly associated with the WV and reflects WV content in the troposphere.Bevis et al. (Bevis et al., 1992) introduced the principle of using GNSS Zenith Wet Delay (ZWD) to retrieve the Precipitable Water Vapor (PWV).Since then, many scientists carried out the GNSS PWV experiments (Askne and Nordius,1987;Bokoye et al., 2003;Yao et al., 2014;Lu et al., 2015;Shoji and Sato, 2016).Now, the GNSS PWV can be retrieved with an uncertainty of 1-2 mm in post-processing (Tregoning et al., 1998;Adams et al., 2011;Grejner-Brzezinska, 2013) or real-time modes (Yuan et al., 2014;Li et al., 2014;Li et al., 2015).
The GNSS WV tomography technique was first proposed to monitor the 3D or 4D WV in 2000 (Flores et al., 2000;Seko et al., 2000;Hirahara et al., 2000).Since then, many scientists have proposed refined methods to improve the GNSS WV tomography (Flores et al., 2001;Nilsson and Gradinarsky, 2006;Rohm and Bosy, 2011;Wang et al., 2014;Wang et al., 2014;Zhao and Yao, 2017).The GNSS WV tomography methods can be roughly categorized into two groups.One group solves the tomography equation in the least squares scheme or in the Kalman filter scheme with additional constraints or using a priori information (Flores et al., 2000;Rohm and Bosy, 2011;Cao et al., 2006;Zhang et al., 2017).The other group uses the algebraic reconstruction algorithm or similar methods (Bender et al., 2011;Wang et al., 2014;Zhao and Yao, 2017).Some scientists also use different methods from the above to solve the GNSS WR tomography (Nilsson and Gradinarsky, 2006;Perler et al., 2011;Altuntac, 2015).Though the tomography technique has the advantages of (1) free of weather conditions and (2) retrieve 3D WR field in near real time, it still suffers some problems.The sparse distribution of the GNSS receivers and the bad satellite-receiver geometry lead to serious ill-posed and illconditioned problems, and also limit the WR retrieve resolution in both vertical and horizontal domains.
Besides the GNSS tomography technique, the WR can also be retrieved by numerical weather prediction models (Gutman and Bock, 2007;Perler et al., 2011).The Weather Research and Forecasting (WRF) model is a state-of-the-art atmospheric modeling system that is used to simulate the dynamic processes of the atmosphere (Jankov et al., 2005;Carvalhoaabc et al., 2012).It is mainly developed and supported by Mesoscale and Microscale Meteorology (MMM) Laboratory of the National Center for Atmospheric Research.
Both the GNSS tomography technique and the WRF model could retrieve 3D WR field.It will be interesting to compare their capabilities in retrieving WR field under different weather conditions and to explore the feasibility to combine them.For this purpose, we conduct GNSS tomography and data assimilation experiments in Hong Kong area using the Hong Kong SatRef Network in a rainy period and a rainless period.
WR fields retrieved from GNSS tomography and WRF reanalysis are validated by the radiosonde data.We also explore the feasibility of assimilating the GNSS tomographic WR into the WRF model to further improve the WR field.

Study Area
The study area is between 113.75°E-114.5°E and 22°N-22.6°N

Data Analysis
Two periods of GNSS observation data are processed to generate ZTD and Slant Wet Delay (SWD).One is a rainy period from July 20 to 26, 2015 when Hong Kong suffered the heaviest daily rainfall of 2015 (191.3 mm rainfall on July 22).The other is a rainless period from August 1 to 7, 2015.
We adopt the same settings as detailed in Zhang et al. (2017) to process the GNSS data.The precise point positioning module in Bernese 5.0 software is used to process the GNSS observations and generate GNSS ZTD data.The International GNSS Service final orbit and clock products are used.The differential code Biases (DCB) is corrected by products from the Center for Orbit Determination in Europe.Antenna phase center offsets and variations, phase wind-up, Earth tides, Earth rotation, ocean tides and relativistic effects are corrected by conventional methods detailed in (Kouba and Hé roux, 2001).We use the ionosphere-free combination of double frequencies to eliminate the first order ionospheric delay and the higher-order terms are ignored.The tropospheric delay models are Saastamoinen model (Saastamoinen, 1972) and Neill mapping functions (Niell, 1996).The zenith hydrostatic delay is estimated and removed from the ZTD by using the insitu pressure observations and Saastamoinen model.The cut-off elevation angle is 10°.The SWD is reconstructed by mapping the ZWD and horizontal gradients onto the ray direction.The phase residuals are added to SWD to consider the inhomogeneity of the troposphere.

Method
The Rapid Radiative Transfer Model (Mlawer et al., 1997) and Dudhia's scheme (Dudhia, 1989) were used for longwave radiation and shortwave radiation, respectively.The nested mode is not used.
We use the reanalysis data from European Center for Medium-Range Weather Forecasts (ECMWF) ERA-Interim as the background data, whose nominal spatial resolution is 0.125°× 0.125°.The procedures to do the assimilation experiments are shown in Figure 2.  The background data are processed by WRF preprocessing system (WPS).The WRFDA is run with the generic CV3 option, and the default background error is adopted in this study.The GNSS ZTDs are the input observations for WRFDA.We run WRFDA and obtain the output which is then used to initialize the WRF model.For comparison's sake, we also run the WRF model using the WPS output as the initial conditions.
When we obtain the reanalysis from the WRF model, we use Equation (1) to calculate WR (Vedel and Huang, 2004).where P is the pressure in Pascal, q is the specific humidity in kg/kg.

GNSS tomography
The limited number of stations, flat vertical distribution of stations, and bad satellite-station geometry impose serious ill-posed problem in the GNSS WR tomography.To well handle this problem, we use the tomography method proposed by Zhang et al. (2017).This method is based on the adaptive Laplacian smoothing and Helmert Variance Component Estimation.It also uses the meteorological data from each GNSS station to constrain the WR near the ground.This tomography strategy is free of a priori information, which makes it an independent technique and thus ensures the fairness when the tomography technique is compared with the WRF model.The WR can be retrieved directly by this tomography strategy when the SWDs are used as observations.The troposphere is vertically divided into 13 layers with a constant thickness of 800 meters, and horizontally divided into grids whose resolution is ~10 km in longitudinal direction and ~8 km in latitudinal direction.

Results
The radiosonde data are used to validate the WR derived from GNSS tomography and reanalysis.Since the radiosonde launches at 0:00 and 12:00 UTC daily, the WR at these time points are validated.Equation ( 1) is also used to calculate WR from radiosonde data.In the horizontal direction, we use WR from reanalysis at the nearest four grids to interpolate the WR at the radiosonde location by the bilinear interpolation method.In the vertical direction, we interpolate the above results and the radiosonde observation to the centers of the associated tomography voxels by the linear interpolation method.Finally, the WR from reanalysis and GNSS tomography are validated by the radiosonde data.For simplicity, WR from radiosonde data, reanalysis that assimilate ZTD, reanalysis that without assimilating GNSS ZTD, and GNSS tomography are denoted as "Radiosonde", "Reanalysis1", "Reanalysis2", and "Tomography" hereinafter.
Figures 3 and 4 show the vertical profiles of the Reanalysis1, the Reanalysis2, the Tomography, and the Radiosonde in the July and August periods, respectively.The Reanalysis1, Reanalysis2, and Tomography agree well with the Radiosonde, which indicates that these three methods successfully retrieved the vertical profile of the WR.It is also observed that the Reanalysis1, the Reanalysis2, and the Tomography agree better with the Radiosonde in the July period than in the August period.This difference should be due to the vertical distribution of WR.Though Hong Kong suffered heavy rain in the July period, the WR was more evenly distributed from 0 to 10 km height than that in the August period.In the rainless August period, the WR was highly concentrated in the lower troposphere (< 6 km) and its vertical changes were very sharp.This situation decreased the performance of the tomography technique and the WRF model.This may also indicate that both methods have decreased capabilities in retrieving WR in highly changing troposphere.Compared with Reanalysis2, the Reanalysis1 is slightly improved, but the improvement is not significant.The difference   Figure 5 shows the statistics of the bias, standard deviation (STD), and Root Mean Square (RMS) of the Tomography, the Reanalysis1, and the Reanalysis2 validated by the Radiosonde at different heights.In the rainy period, bias of Reanalysis1 is smaller than that of Reanalysis2, but the differences are not obvious in terms of STD and RMS.In the rainless period, the bias of Reanalysis1 in the lower troposphere is slightly greater than that of the Reanalysis2.Overall, the differences between Reanalysis1 and Reanalysis2 are not significant.
In the rainy period, the bias, STD, and RMS of the Tomography are greater than that of the Reanalysis1 in most of the time.But in the rainless period, the STD and RMS of the Tomography tend to be smaller than that of the Reanalysis1 in the lower troposphere, but its bias is still greater.In general, the WRF model performs better than the tomography technique in most of the cases, but in the lower troposphere in the rainless period the tomography may perform better than the WRF model.Table 1 shows the bias, STD, and RMS of Tomography, Reanalysis1, and Reanalysis2 validated by the Radiosonde.In the whole troposphere in the rainy period, Tomography has the smallest bias but the largest STD and RMS.The Reanalysis1 and Reanalysis2 have the similar STD and RMS that are much smaller than that of the Tomography.But the Reanalysis2 has the largest bias than Reanalysis1 and the Tomography.In the lower troposphere in the rainy period, Reanalysis1 has the smallest STD and RMS while the Tomography has the largest ones.The bias of Tomography is positive in the low troposphere but negative in the upper troposphere, this should be due to the vertical smoothing constraints imposed on the WR.In the upper troposphere in the rainy period, Tomography has the largest bias, STD, and RMS while Reanlysis1 has the smallest ones.Overall, both the tomography and the reanalysis results have larger bias, STD, and RMS in the lower troposphere than in the upper troposphere, indicating both the tomography technique and the WRF model has deceased capabilities in the lower troposphere.
In the whole troposphere in the rainless period, Reanalysis2 has the smallest bias but the largest STD and RMS.The STD and RMS of Tomography are larger than Reanalysis1.In the lower troposphere in the rainless period, Reanalysis2 has the largest RMS and STD while Reanalysis1 has the smallest ones.In the low troposphere in the rainless period, the performance of Tomography is not as good as Reanalysis1 in terms of RMS.However, in the upper troposphere in the rainless period, the Tomography has relatively larger bias, STD and RMS than the reanalysis results.In general, assimilating GNSS ZTD into the WRF model could slightly improve the WR but the improvement is not significant.The reanalysis WR overall outperforms the tomography WR but the tomography WR may be better in the lower troposphere in the rainless period.The results are better in the rainy period than in the rainless period, which is mainly due to the sharp vertical variation of WR in the rainless period.

Discussion
In the rainless period, due to the sharp vertical variations of WR, the Tomography, the Reanalysis have decreased performance in retrieving the WR, especially in the lower troposphere.Compared with the results in the rainy period, the RMS of the Tomography and the Reanalysis1 increases by 0.94 mm/km, 3.24 mm/km in the rainless period, respectively.The accuracy decrease is more significant in the Reanalysis1 than in the Tomography, resulting in that the tomographic WR becomes better than the reanalysis WR (Figures 5d and     5f) in the low troposphere.
When assimilating ZTD into the WRFDA, we only use the total column water vapor and do not know its vertical structure.This leads to that the assimilation of ZTD has limited improvement in retrieving the vertical structure of the WR.Therefore, it is natural to consider assimilating good tomographic WR into the WRFDA model to improve the retrieval of the vertical structure of WR.At present, WRFDA could not assimilate WR directly, but can assimilate meteorological parameters such as relative humidity, temperature and pressure.To assimilate the tomographic WR, we convert WR to relative humidity.
The relationship between relative humidity (RH) and w P is shown as Equation (3).(3) where s P is the saturated water vapor pressure which is related to temperature and can be calculated by Wexler formula (Wexler, 1976(Wexler, , 1977)).The w P is calculated by Equation ( 1).The needed temperature and pressure data are from the reanalysis data.After converting the tomographic WR to RH, we assimilate the RH together with the corresponding temperature and pressure to the WRFDA.Then, the similar procedures as described in Section 3.1 are performed to generate reanalysis.
The Tomography agrees better with the Radiosonde than the Reanalysis1 and Reanalysis2 in the lower troposphere below 3 km at 12:00 on August 6 (Figure 4l) and at 12:00 on August 7 (Figure 4n).So, we assimilate the tomographic WR below 3 km into the WRFDA model at these two time points.The generated reanalysis data are denoted as "Reanalysis3".The difference between Reanalysis3 and Radiosonde is denoted as "DA-Tomo".The difference between Reanalysis1 and Radiosonde is denoted as "DA-ZTD".The difference between Tomography and Radiosonde is denoted as "Tomo".The differences at different heights at 12:00 on August 6 and 7 are shown in Figure 6.In the rainy period in the whole troposphere, the RMS of Tomography, Reanalysis1 and Reanalysis2 are 6.50 mm/km, 4.31 mm/km, and 4.15 mm/km.The RMS becomes 7.02 mm/km, 6.35 mm/km, and 7.26 mm/km in the rainless period.Both methods obtained better WR in the rainy period than in the rainless period.We infer that the sharp vertical variations of WR reduced the WR retrieving accuracy in the rainless period.In most of the cases, the reanalysis WR outperforms the tomographic WR but the tomographic WR may be better than the reanalysis WR in the lower troposphere in the rainless period.By assimilating better tomographic WR in the lower troposphere into the WRFDA model, we slightly improve the reanalysis WR.
The above results suggest that both the WRF model and the tomography technique can retrieve good WR but also have drawbacks.If we combine the two by assimilating good tomographic WR into the WRFDA, we may further improve the performance of the WRF model in retrieving the water vapor field.
as shown in Figure 1.There are 15 continues GNSS stations belonging to the Hong Kong SatRef Network deployed in the study area.They are all equipped with Leica GNSS receivers and antennas to receive the GNSS signals and automatic meteorological devices Ann.Geophys.Discuss., https://doi.org/10.5194/angeo-2018-84Manuscript under review for journal Ann.Geophys.Discussion started: 16 July 2018 c Author(s) 2018.CC BY 4.0 License. to record the temperature, pressure, and relative humidity.The average inter-distance between stations is about 10 km.The altitudes of the highest station (HKNP) and the lowest station (HKLM) are 354 m and 10 m.Therefore, this network is vertically flat.

Figure 1 .
Figure 1.Research area of the experiment.The red triangles indicate the GNSS stations and the blue star indicates the radiosonde station in Hong Kong.

Figure 2 .
Figure 2. Flowchart of data assimilation using the WRF model.
. Discuss., https://doi.org/10.5194/angeo-2018-84Manuscript under review for journal Ann.Geophys.Discussion started: 16 July 2018 c Author(s) 2018.CC BY 4.0 License.where w P is the water vapor pressure in Pascal, T is the temperature in Kelvin.

Figure 3 .
Figure 3. Comparisons among WR derived from reanalysis, tomography, and radiosonde in the rainy158

Figure 5 .
Figure 5. Statistics of bias, STD, and RMS of Tomography, Reanalysis1, and Reanalysis2 validated by the Radiosonde.

Figure 6 .
Figure 6.Differences between WR obtained by various methods and radiosonde WR.

Figure 6
Figure6shows that the DA-ZTD is very close to the DA-Tomo.The average DA-ZTD is 6.04 mm/km and the average DA-Tomo is 5.92 mm/km.This indicates that assimilating tomographic WR into the WRF model can slightly improve the WR retrieve.The large uncertainty (8.35 mm/km) of the tomography WR in the lower troposphere may limit the improvement.6.ConclusionsGNSS WR tomography and data assimilation experiments are conducted in Hong Kong during a rainy and a rainless period to test the capabilities of the tomography technique and the WRF model in retrieving WR.The results show that both the tomography technique and the WRF model can retrieve WR that agrees well with the radiosonde data.

Table 1 .
Statistics of bias, RMS and STD of Tomography, Reanalysis validated by the radiosonde WR.