A new global grid-based weighted mean temperature model considering vertical nonlinear variation

Global Navigation Satellite Systems (GNSS) have been proved to be an excellent technology for retrieving precipitable water vapor (PWV). In GNSS meteorology, PWV at a station is obtained from a conversion of the zenith wet delay (ZWD) of GNSS signals received at the station using a conversion factor which is a function of weighted mean temperature (Tm) along the vertical direction in the atmosphere over the site. Thus, the accuracy of Tm directly affects the 10 quality of the GNSS-derived PWV. Currently, the Tm value at a target height level is commonly modelled using the Tm value at a specific height and a simple linear decay function, whilst the vertical nonlinear variation in Tm is neglected. This may result in large errors in the Tm result for the target height level, as the variation trend in the vertical direction of Tm may not be linear. In this research, a new global grid-based Tm empirical model with a horizontal resolution of 1°× 1°, named GGNTm, was constructed using ECMWF ERA5 monthly mean reanalysis data over the 10-year period from 2008 to 2017. A 15 three-order polynomial function was utilized to fit the vertical nonlinear variation in Tm at the grid points, and the temporal variation in each of the four coefficients in the Tm fitting function was also modelled with the variables of the mean, annual and semi-annual amplitudes of the 10-year time series coefficients. The performance of the new model was evaluated using its predicted Tm values in 2018 to compare with the following two references in the same year 1) Tm from ERA5 hourly reanalysis with the horizontal resolution of 5°× 5°; 2) Tm from atmospheric profiles from 428 globally distributed radiosonde 20 stations. Compared to the first reference, the mean RMSEs of the model predicted Tm values over all global grid points at the 950hPa and 500hPa pressure levels were 3.35K and 3.94K respectively. Compared to the second reference, the mean bias and mean RMSE of the model predicted Tm values over the 428 radiosonde stations at the surface level were 0.34K and 3.89K respectively; the mean bias and mean RMSE of the model’s Tm values over all pressure levels in the height range from the surface to 10 km altitude were −0.16K and 4.20K respectively. The new model results were also compared with that of the 25 GPT3, GTrop and GWMT_D models in which different height correction methods were also applied. Results indicated that significant improvements made by the new model were at high-altitude pressure levels; in all five height ranges, GGNTm results were generally unbiased, and their accuracy varied little with height. The impact of Tm on GNSS-PWV was evaluated in terms of relative error, and significant improvement was found compared to the widely used GPT3 model. These results suggest that considering the vertical nonlinear variation in Tm and the temporal variation in the coefficients of the Tm model 30 can significantly improve the accuracy of model-predicted Tm for a GNSS receiver that is located in anywhere below the https://doi.org/10.5194/amt-2020-274 Preprint. Discussion started: 7 October 2020 c © Author(s) 2020. CC BY 4.0 License.

tropopause (assumed to be 10 km), which has significance for applications needing real-time or near real-time PWV converted from GNSS signals.

Introduction
Water vapor, as an important greenhouse gas, is tightly related to weather variations, hence, it is crucial to monitor the 35 water vapor content in the atmosphere for reliable weather forecast. The meteorological parameter that is closely related to water vapor is precipitable water vapor (PWV) and it can be measured by various technologies such as radiosondes, remote sensing satellites and water vapor radiometers. Global Navigation Satellite Systems (GNSS), which were initially designed for positioning, navigation, and timing, can be used to retrieve the zenith tropospheric delay (ZTD) of the GNSS signal over an observation station. The ZTD can be divided into zenith hydrostatic delay (ZHD) and zenith wet delay (ZWD). The ZHD can 40 usually be obtained at a high accuracy from a standard empirical model together with some surface meteorological data at the station. The atmospheric water vapor information is contained in the GNSS-ZTD, more precisely, in the GNSS-ZWD, which can be converted into PWV. The GNSS were first applied to meteorological research in the 1990s (Bevis et al., 1992). Near real-time GNSS-ZTD products estimated from GNSS data processing have been routinely assimilated into numerical weather models (NWM) for improving the performance of weather forecast (Bennitt and Jupp, 2012;Dousa and Vaclavovic, 2014; 45 Guerova et al., 2016;Le Marshall et al., 2012. Research has demonstrated that the accuracy of GNSS-PWV can meet the accuracy requirements for most meteorological applications, and the applicability of GNSS-PWV for studying extreme weather conditions and climatic events has been investigated by many researchers (Bonafoni et al., 2013;Calori et al., 2016;Chen et al., 2018;Choy et al., 2013;Eugenia Bianchi et al., 2016;Shi et al., 2015;Wang et al., 2016Wang et al., , 2018Zhang et al., 2015).
To obtain GNSS-PWV over a station, the first step is to estimate the ZTD of the station from GNSS data processing, and the 50 two most common data processing strategies are the network approach and precise point positioning (PPP) approach (Ding et al., 2017;Douša et al., 2016;Guerova et al., 2016;Li et al., 2015;Lu et al., 2015;Rohm et al., 2014;Yuan et al., 2014;Zhou et al., 2020). The former uses double-differenced observations, while the latter uses un-differenced observations in the observation equation system. The ZWD can be obtained from subtracting the ZHD from the GNSS-ZTD, or directly estimated if the ZHD has been corrected in the GNSS observation equation system, depending on the processing strategies adopted. Then 55 the GNSS-PWV can be converted by: where is the conversion factor (Askne and Nordius, 1987;Bevis et al., 1992), which is given by: where is the density of liquid water; is the specific gas constant for water vapor; 2 ′ and 3 are atmospheric refractivity constants; is the weighted mean temperature over the GNSS site, which is defined and approximated through where and are the water vapor pressure (hPa) and absolute temperature (K) respectively; is the number of the layers; , and ∆ℎ are the mean water vapor pressure, mean temperature, and thickness of the th layer respectively.
From Eq.
(2), one can see that is a crucial variable for the determination of the conversion factor , which in turn affects the determination of PWV expressed by Eq.
(1). The significance of obtaining accurate values has been demonstrated by previous researches (Bevis, 1994;Jiang et al., 2019;Ning et al., 2016;Wang et al., 2005Wang et al., , 2016. can be calculated from an observed atmospheric profile or NWM data (Wang et al., 2005. This observed atmospheric profile can be acquired from a radiosonde station, which is valid only for the sounding site. In fact, for GNSS stations, they are usually not co-located with any regional radiosonde stations, i.e. observed atmospheric profiles are unavailable, as a result, equation (3) is not applicable for GNSS stations. Moreover, even a GNSS station is co-located with a radiosonde station, due to the low temporal and spatial resolution of radiosonde data, the temporal resolution of its resultant is also low, which cannot meet 70 the requirements of GNSS near real-time/real-time (NRT/RT) applications such as the conversion of GNSS-ZWD time series into PWV time series. Although atmospheric profiles can be obtained from NWM data, it is still difficult for users to obtain predicted results from the NWM data for NRT/RT GNSS-PWV sensing. Thus, it is of great importance to develop empirical models for time-critical applications. Some models have been developed with a focus of improving the accuracy of the , and these empirical models can be classified into two categories. One category is such a model that depends on in situ 75 surface temperature observation , like the Bevis model, which is a simple linear function expressed as: = + (Bevis et al., 1992). The two coefficients of such a linear function can be determined from the linear regression method based on long-term regional radiosonde data. However, the deployment of radiosonde stations is geographically sparse due to their high cost, and even worse is that there are no radiosonde stations at all in some areas. Yao et al. (2014a) developed a global latitude-dependent − linear model using data from the global geodetic observing system (GGOS) and data 80 from the European center for medium-range weather forecasts (ECMWF). Jiang developed a time-varying global gridded − model using both and derived from ERA-Interim (Jiang et al., 2019). Ding (2018Ding ( , 2020 developed two generations of global models using the neural network algorithm, in which temperature observations were required for the input and the models performed well. The models mentioned above need in situ meteorological observations as the model's input. However, for GNSS stations, again, not all stations are equipped with meteorological sensors. In this case, the type of empirical models that are independent of meteorological observations had to be constructed. Yao et al. (2014bYao et al. ( , 2012Yao et al. ( , 2013 used spherical harmonics to develop the GWMT, GTM-Ⅱ and GTM-Ⅲ models, in which both the height and the periodicity of were taken into account. and altitude. The widely used GPT2w model (Böhm et al., 2015) and its successor, GPT3 (Landskron and Böhm, 2018),

90
provided gridded results with both 1°×1°and 5°×5° horizontal resolutions and the models also contain a few terms related to temporal variations in including the mean, annual and semiannual amplitudes. However, the height differences between the user site, e.g. a GNSS station, and its nearest four surrounding grid points were not considered.  (Huang et al., 2019b), GTm_R  and GPT2wh (Yang et al., 2020). The GTrop model (Sun et al., 2019), developed for predicting both ZTD and , also took into account the lapse-rate, and it 100 outperforms GPT2w obviously at altitudes under 10 km.
As previously discussed, considering the lapse-rate in a model can improve the model's accuracy. However, the assumption that linearly varies with height, which many recently developed models were based on, may not agree well with the truth. In this research, a new global grid-based empirical model, named GGNTm, in which the vertical nonlinear variation of was taken into account, was developed using a three-order polynomial function and ERA5 105 monthly mean reanalysis data over the 10-year period from 2008 to 2017, and the temporal variation in each of the four coefficients in the fitting function was also modelled with the variables of the mean, annual and semi-annual amplitudes of the 10-year time series coefficient.
The outline of the paper is as follows. The features of the vertical nonlinear variation in were investigated in Sect.
2.2, then a three-order polynomial function fitting the 10-year profiles obtained from ERA-5 monthly mean reanalysis 110 data was developed for the GGNTm model. In Sect.3, the performance of GGNTm was validated using the values from ERA5 hourly reanalysis and globally distributed radiosonde profiles in 2018 as the references. Conclusions are summarized in the final section.
2 Methodology for new model construction

115
ERA5 reanalysis data were the latest reanalysis data developed by the ECMWF. In this research, ERA5 monthly mean reanalysis data in the 10-year period from 2008 to 2017containing geopotential heights, temperatures, and specific humidity at 37 pressure levels with a horizontal resolution of 1°× 1° were downloaded from the web server of the Copernicus Climate Change Service (C3S, https://climate.copernicus.eu/climate-reanalysis). The geopotential heights, which are often used in meteorology, were then converted to WGS-84 ellipsoidal heights. Water vapor pressure was calculated by (Nafisi et al., 2012): where is the specific humidity, which can be obtained from NWM data; is the atmospheric pressure.

Vertical variation of
The ERA5 monthly mean products were used to analyze the vertical variation of . As defined in Eq. (3), is a function of water vaper pressure and temperature. The variation in water vapor pressure in the vertical direction has been 125 known nonlinear, while the vertical variation in temperature is often assumed to be a linear decay function (Dousa and Elias, 2014). In fact, there is such a phenomenon that temperature increases with the increase in height, the so-called temperature inversion, which occurs in both the upper atmosphere and near ground surface, meaning that the vertical variation in temperature is complex. As a result, in the vertical direction varies nonlinearly due to the irregular variations in both water vapor pressure and temperature in the vertical direction. Fig. 1 shows four vertical profiles of water vapor pressure, 130 temperature, and at the pressure levels that were under a 10 km ellipsoidal height at four grid points obtained from ERA5 monthly mean reanalysis in December 2017. It should be noted that the surface heights of the four grid points were different, and they were 0 m, 301 m, 13 m, and 180 m respectively. Sub-figures (a) and (b) show that, in the height range near the surface, temperature increases with the increase in height. In addition, all the four profiles (the black curves with dots) in these sub-figures show a nonlinear variation trend. This implies that using a constant lapse-rate to model the vertical variation trend will result in large errors, i.e. the profiles cannot be accurately modelled through a constant lapse-rate. This finding aligns well with other researchers (e.g. (Yao et al., 2018)).

Three-order polynomial function for vertical fitting
A linear decay function with a constant lapse-rate can be expressed as: where is the value at the reference height ℎ 0 ; is the lapse-rate and is the ellipsoidal height (km) of the user 140 site. An equivalent expression of Eq. (5) is: where , , , are the four unknown coefficient parameters of the fitting function.
named scheme-1 and scheme-2, were for the two functions fitting the sample data of profiles of the 120 monthly mean reanalysis data over the 10-year period from 2008 to 2017 at each grid point. It should be noted that, only those values from the heights under 10 km were selected for the sample data. For measuring how well the fitting function fits the sample 150 data, the root mean square (RMS) of the differences between the values resulting from the fitting function and the sample data was calculated by: where is the residual of at the th pressure level over the grid point. Fig. 2 shows the map for the mean of the RMSs of the fitting residuals of the from the aforementioned 120 monthly mean profiles (the samples) at each of the grid points. The mean of the mean RMSs at all global grid points for scheme-1 and scheme-2 were 1.26 K and 0.30 K respectively.

155
In addition, the RMS results in the left sub-figure (for linear function) were latitude-dependent, and small RMSs (blue) were in mid-latitude regions; large RMS values in both sub-figures were in Antarctica. Comparing the two subfigures, we could find that the RMS values shown in the right sub-figure were all very small and significantly smaller than that of the left sub-figure, meaning that the three-order polynomial fitting function superior to the linear fitting function.

160
In the previous section, the 10-year time series of coefficients in the three-order polynomial function expressed in Eq. (7)  165 are shown in Fig. 3, which presented noticeable annual and semi-annual amplitudes. Similar periodicities were also found at other grid points. According to these characteristics, the fitting model for GGNTm containing three terms including mean, annual and semi-annual amplitudes for each coefficient time series at each grid point expressed by the following was adopted in this study: where 0 , 1 and 2 are the mean, annual and semi-annual amplitudes respectively; denotes "day of year"; 1 and 170 2 are the initial phases of the annual and semi-annual periodicities, which are estimated together with the mean and amplitudes.
Then, the mean, annual and semi-annual amplitudes, and initial phases for each coefficient at each of the grid points over the globe (with the resolution of 1°× 1°) were determined using the least squares estimation method and the 10-year time series of the coefficient. To calculate for a specific site and time, e.g. for a GNSS station at an observing time, the following 3-step procedure needs to be carried out: 1) using Eq. (9) to calculate each of the four coefficients at each of the four grid points surrounding the user site; 2) using Eq. (7) to calculate the values at the height of the user site at each of the above four grid points (which is for the height dimension); 3) using an interpolation method, such as the inverse distance weighting or bilinear interpolation, on the four values 180 from step 2) to obtain the value for the user site (which is for the horizontal dimension, as is shown in Fig. 4).
Till now the new model has been developed based on the 10-year sample data from 2008 to 2017. This model will be validated using the model predicted results in 2018 compared against the same year's (i.e. out-of-sample) reference data.
Results will be discussed in the next section.

185
Our new model was developed for obtaining predicted Tm values over a site that is located in anywhere below the tropopause (assumed to be 10 km), which has significance for applications needing real-time or near real-time PWV converted from GNSS signals received by a GNSS receiver located in a flat area, an ocean area, a high mountainous area, or even a flight vehicle. For the performance assessment of the newly developed model, values over different pressure levels obtained from both ERA-5 reanalysis and radiosonde profiles in 2018 were selected as the references due to the fact that GNSS 190 receivers may be located in any reasonable altitudes. The two statistics, bias and RMSE, were utilized to measure the systematic discrepancy and the accuracy of the model results. Their formulas are: where i is the index of the data element; denotes the model resultant value; denotes the reference value; is the number of the values in the statistics.

195
As the first set of the reference selected for the evaluation of the new model, ERA-5 hourly data (with the resolution of 5°× 5°) at 12:00 UTC on each day in 2018, which were out-of-sample data, were downloaded from the C3S. Then they were converted to profiles, and values at each of five pressure levels: 950hPa, 800hPa, 650hPa, 500hPa and 350hPa were used to calculate the bias and RMSE of the new model's results at the pressure level. In addition to the GGNTm model, other three empirical models developed in recent years including GPT3, GTrop and GWMT_D, in which different vertical correction methods were also applied, were also evaluated for performance comparisons of GGTNm and these three models. Table 1 shows the mean bias and mean RMSE of the values over all global grid points resulting from each of the above four models. As we can see, on a global scale, GGNTm outperformed all the other three models, especially at high pressure levels. The mean bias and RMSE of GPT3 varied significantly due to its lack of height refinement. GTrop was considerably better than GPT3, owing to its use of the lapse-rate, although its results were still had large errors at 205 high pressure levels, which is most likely to be resulted from the neglecting of the nonlinear vertical variation in . The large bias and RMSE of the GWMT_D results were possibly because its modelling was based on NCEP reanalysis data, and there may exist systematic differences between the reanalysis data from ECMWF and NCEP. However, GWMT_D still significantly outperformed GPT3, due to its voxel-based structure. Compared to GTrop and GWMT_D, GGNTm performed very well at all pressure levels. This is because the model accounted for the vertical nonlinear variation in .

210
The results shown in Table 1 were the statistics of all global grid points at each of the five pressure levels selected. For more refined results, Fig. 5 shows the map for the RMSE of at each grid point at either the 950hPa or 500hPa pressure levels resulting from three models, and the reason for not selecting GPT3 here was due to its much poor performance. GGNTm, over the other two models were at high-altitude pressure levels.

Comparison with radiosonde data
In this section, from radiosonde profiles were used as the reference for the performance assessment of the models selected. The original radiosonde data at all globally distributed stations in 2018 were downloaded from the website of the University of Wyoming (http://weather.uwyo.edu/upperair/). Different from the use of reanalysis data as the reference, water vapor pressure at each pressure level from a radiosonde profile was calculated through a mixing ratio: where denotes the mixing ratio (g/kg).
An additional data pre-processing procedure needs to be conducted for data quality control. Those poor radiosonde profiles needed to be identified and excluded from their use for the reference. The first check was for a valid mixing ratio value: if a pressure level lacks a valid mixing ratio value, then it is regarded invalid and thus to be excluded. After this initial checking was performed, further identifications were also carried out. A profile would be excluded if it met any one of the following 230 four conditions: (1) the profile lacks surface meteorological observations; https://doi.org/10.5194/amt-2020-274 Preprint. Discussion started: 7 October 2020 c Author(s) 2020. CC BY 4.0 License.
(2) the pressure value of the top pressure level is greater than 100 hPa; (3) the difference in the pressure values at two successive levels is under 200hPa; (4) the profile consists of a few pressure levels, e.g. if / ≤ 30 ℎ (where is the difference of the pressure 235 values at the surface and the 100 hPa pressure levels, and n is the number of all pressure levels from the surface to the 100hPa pressure levels), then the profile was regarded to have sufficient number of pressure levels, otherwise it would be excluded from the use in the testing.
Sounding balloons are commonly launched twice a day (at 00:00 and 12:00 UTC). In this research, only those stations that had at least 300 profiles in 2018 were selected in the model performance assessment. After the above 5-step quality control 240 procedure was performed, a total of 260140 profiles from 428 global radiosonde stations were finally used in the performance evaluation of four selected models. close; GWMT_D was the worst, with the largest bias and RMSE values, which may be due to its low horizontal resolution (5°×5°). The other set of results, the RMSE of under 10 km, was calculated using the differences between model-predicted values and the reference Tm values over all pressure levels that with a height less than 10 km. A small RMSE of Tm under 10 km indicates that the model performs well at any altitudes below the tropopause. As we can see, GPT3 performed the worst and significantly worse than the other three models, mainly due to lack of an appropriate modelling for the vertical 250 variation; GWMT_D was slightly better than GTrop, possibly because the value from the former was interpolated from the values at four height levels; the mean bias of from the new model, GGNTm, was the least, with the value of −0.16, which was close to 0, meaning nearly unbiased; the RMSE of the new model was also the least, among the four models, which suggest that the vertical nonlinear variation of was modelled in the new model more accurately than the other existing models.
255 Similar to Fig. 5, Fig.6 shows the map for the RMSE of values at each of the 428 radiosonde stations in 2018 at the surface pressure level (the left column) and all pressure levels that with a height less than 10 km (the right column) resulting from GTrop, GWMT_D, and GGNTm. It can be found that the RMSEs of all models were latitude-dependent, and those stations that had a large RMSE value were most located in north Africa and north-east America. At the four stations located in Antarctic, their surface values were accurately modelled by these models. However, in terms of the RMSE of all pressure 260 levels under 10 km, the GTrop results were relatively large at the four stations, whilst both GWMT_D and GGNTm performed well at three of the stations.
To further evaluate the performance of the three models at different height ranges under 10 km, the models' values from the aforementioned radiosonde profiles at the 428 global stations were divided into five height ranges, and Fig. 7 shows each height range's bias and RMSE. We can see the following results 1) in the height ranges above 4 km, the GTrop results had the largest bias and largest RMSE, and GWMT_D was considerably better than GTrop; 2) in low height ranges the GWMT_D results were the worst; 3) in all height ranges the GGNTm results were nearly unbiased and their accuracy varied little with height. The GGNTm model's consistent high accuracy in all height ranges suggest that the characteristics of the vertical nonlinear variation in is modelled by the proposed model more accurately than the other models.

270
The accuracy of GNSS-PWV over a GNSS site at an observing time is dependent upon the accuracies of the ZWD and the conversion factor. Uncertainty analysis has been conducted by some researches to study the uncertainty of the GNSS-derived PWV, including the uncertainty in the conversion factor (Bevis, 1994;Jiang et al., 2019;Ning et al., 2016). This research mainly focuses on the impact of the newly developed model on PWV, however, it is difficult to evaluate the impact of on the GNSS-PWV directly, as GNSS stations and radiosonde stations are usually not collocated.

275
Thus, a theoretical model was adopted by several studies to study the error in PWV resulting from (He et al., 2017;Huang et al., 2019a;Li et al., 2020;Wang et al., 2005Wang et al., , 2016. The relationship of the relative error between and PWV can be expressed as: where ∆ is the absolute error in PWV resulting from ; ∆ / is the relative error of PWV. Due to the fact that 2 ′ / 3 is quite small (≈ 6 −5 −1 ), the relative error of PWV resulting from is approximately equal to that of .

280
In this research, surface and PWV were obtained from the abovementioned radiosonde profiles downloaded from the University of Wyoming. The mean absolute error and mean relative error of PWV resulting from GGNTm model at the abovementioned radiosonde stations are shown in Fig. 8. As is shown in the figure, the distribution of both mean absolute error and mean relative error of PWV resulting from GGNTm are latitude-dependent, and stations at high latitudes tended to have smaller absolute errors but larger relative errors compared to stations located in low-latitude regions. The mean of mean absolute error 285 and mean relative error of PWV resulting from derived from GGNTm on surface level were 0.19 mm and 1.13%, respectively. values over different pressure levels of radiosonde profiles at the 428 radiosonde stations in 2018 were then divided into five height ranges to calculate the relative error of PWV resulting from selected empirical models. As is shown in

4 Conclusions
In GNSS meteorology, is an essential parameter for converting GNSS-ZWD to PWV over the GNSS observing station. In practice, the value over a GNSS station at an observing time is commonly obtained from an empirical respectively; and the mean bias and mean RMSE of resulting from GGNTm at all pressure levels from surface to 10 km 305 height were -0.16K and 4.20K respectively, which were significantly smaller than that of all the other three models. In all five height ranges from surface to 10 km in altitude, the GGNTm results were nearly unbiased, and their accuracy varied little with height. This result suggests that the characteristics of the vertical nonlinear variation in is modelled by the approach proposed in this study more accurately than the existing models. In addition, the impact of GGNTm on GNSS-PWV was analyzed using a theoretical function. The results showed that the relative error of PWV resulting from GGNTm outperformed 310 GPT3, GTrop and GWMT model.
The improvement in the accuracy of the new model has significance for GNSS meteorology. Our future work will be focussing on using high temporal resolution atmospheric data such as ERA5 hourly reanalysis data, instead of monthly mean data used in this study, to model the temporal variation of the coefficents in the fitting function for further improving the accuracy of the GGNTm model.