1 Introduction

Urban heat island (UHI) effect arises as urban regions become warmer than their rural environments (Roth 2013). The UHI is mainly caused by the heat absorbed by built structures and anthropogenic heat sources in cities (Roth 2013). The UHI intensity (UHII) is related to many factors such as regional climate, urbanization, and topography.

The UHI effect in many metropolises are reported to be substantially worse associated with the rapid development of urbanization since the last a few decades (Zhai et al. 2016; Xu et al. 2018; Zou et al. 2019; Jiang et al. 2019). The UHI effect severely affects the human health and even the urban ecology system. For instance, it increases the heat stress of citizens and triggers cardiovascular, respiratory and mental diseases. Studies on the formation and spatiotemporal evolution of UHI have become the focus of attention of many scholars in recent years (Rizwan et al. 2018; Memon et al. 2009; Azevedo and Leal 2017; Lamarca et al. 2018).

Traditionally, the UHII can be estimated by using temperature measurements from various techniques, such as ground-based (Ramamurthy and Sangobanwo 2016), airborne-based (Peng et al. 2017), and satellited-based (Schwarz et al. 2011; Wu et al. 2014; Fang et al. 2016; Kayet et al. 2016) sensors. However, these approaches are limited by specific disadvantages, such as the low spatial resolution of the sparse-distributed ground-based sensors, the high cost of the airborne measurements, and the weather-dependent availability of the satellite measurements (Jorge et al. 2021). Global navigation satellites system (GNSS) is a novel atmospheric sounding technique characterized with advantages of high accuracy, (near) real time, all-weather availability, no need for human interference, and no need for instrument calibration (Kouba and Héroux 2001; Cai and Gao 2013; Cao et al. 2015; Yuan et al. 2023a). It has been widely used to retrieve atmospheric water vapor and related phenomena, such as extreme weather (Bonafoni et al. 2019), monsoon and drought (e.g., Jiang et al. 2017), and climate change (e.g., Yuan et al. 2021, 2023b).

Recently, the applications of GNSS in atmospheric environmental science have been extended to the UHII monitoring as presented by Jorge et al. (2021). This algorithm was based on the relationship between single-station GNSS-derived Zenith Tropospheric Delay (ZTD) and different meteorological variables, like pressure, water vapor partial pressure, and temperature. The UHII was calculated as the GNSS-inverted temperature differences between the urban and rural stations. However, due to limited number of GNSS stations in most cities, the site-specific GNSS atmospheric temperature inversion is difficult to provide a comprehensive monitoring of the heat island effect. Nevertheless, GNSS tomography has been an effective means to retrieve three-dimensional (3-D) distribution of wet refractivity (Troller et al. 2006; Bender et al. 2011; Rohm 2013; Chen and Liu 2014; Xia et al. 2018) and it is potential for the monitoring of UHII. A key advantage of GNSS tomography compared to the site-specific GNSS temperature inversion is that it is capable to measure the temperature changes in both horizontal and vertical directions in addition to the three-dimensional distribution of the temperature over the study region.

In order to monitor the UHII, the temperature can be estimated from GNSS-derived wet refractivity. It is easy to know that the quality of the temperature estimates is related to the GNSS-derived wet refractivity according to the error propagation equation. The wet refractivity obtained from GNSS tomography are influenced by various factors, such as the accuracy of Slant Wet Delay (SWD), the number of satellite rays, grid division, and tomography model. Previous studies have focused on how to improve the accuracy of GNSS tomography. Such as, combining multi-source water vapor data (Chen and Liu 2016) or reducing the error in signal propagation can also further improve the accuracy of wet delay information (Ye et al. 2016; Möller and Landskron 2019; Heublein et al. 2019). Besides, some studies tried to use more signals to establish the observation equation (Champollion et al. 2005; Yao and Zhao 2016; Xiong et al. 2021; Zhang et al. 2022a, b). In terms of the grid division, Chen and Liu (2014) proposed a novel method to establish the optimal horizontal distribution of voxels. Chen et al. (2020) developed an improved parameterized algorithm to refine the tropospheric tomographic model to enhance the performance of the wet refractivity reconstruction. Zhang et al. (2022a, b) developed a GNSS combining remote sensing (RS) tomography model to exploit the adding value of RS measurements to GNSS tomography. Compared to reference values from radiosonde, the root-mean-square (RMS) error of the water vapor profiles derived from the developed approach was reported to be reduced by 28% with respect to the GNSS-only results. Yang et al. (2023) investigated the tomographic window and sampling rate for the modeling of the water vapor tomography, and then, the authors recommended to set the tomographic window width and the sampling rate as 10 min and 300 s, respectively, for the GNSS water vapor tomography with dense stations.

In this paper, we proposed an approach for the evaluation of UHII by using GNSS 3-D troposphere tomography as well as GNSS radio occultation (RO) and external radiosonde data. We first developed an optimization of GNSS tomography technique with radiosonde and RO historical data for the retrieval of high-accuracy 3-D wet refractivity distributions. Then, we calculated the temperature from the tomography-derived wet refractivity. Feasibility of the approach was verified by using the GNSS observations in Hong Kong, China, in 2020. In addition, temperature measurements at five synoptic stations in Hong Kong were taken as references for the validation of the GNSS-inverted temperatures. The paper is organized as follows: Section 2 presents the method and Sect. 3 describes the processing of the data. The discussions are given in Sect. 4, followed by the conclusions in Sect. 5.

2 Methodology

This section first introduces the principles of GNSS tomography model, and then describes the algorithms for tomography grid division and temperature calculation from wet refractivity, and at last presents the calculations of UHII from temperature.

2.1 Tomography model

The SWD along the ray paths of dual-frequency GNSS data traversing the imaging region should first be derived to reconstruct the 3-D images of the atmospheric wet refractivity distributions. The SWD is defined by the line integral of atmospheric wet refractivity (Nw) along the ray path from the satellite to the ground-based GNSS receiver (Perler et al. 2011; Jiang et al. 2014) as follows:

$$ {{\rm{SWD}}}=10^{-6} \cdot \int\limits^{\infty}_{{{h_{0}}}} N_{w} \cdot dh = 10^{-6} \cdot \sum \limits_{i=1}^{m} S_{i} \cdot N_{w} ({S_{i} }), $$
(1)

in which h0 is the height of the station above Mean Sea Level (MSL), m is the number of the grids which the signals have passed, and Si is the length that GNSS rays span ith voxel; \({N}_{w}\left({S}_{i}\right)\) is the wet refractivity corresponding to the ith grid. Alternatively, Eq. (1) can be expressed in a concise matrix form as follows:

$$ {\mathbf{SWD}} = {\varvec{S}} \cdot {\varvec{N}}_{{\varvec{w}}} + \Delta_{{{\mathbf{SWD}}}} , $$
(2)

in which S represents the distance of the GNSS signals spanning the voxel, \({{\varvec{N}}}_{{\varvec{w}}}\) represents the wet refractivity of voxel nodes, and ΔSWD is the noise. In the parameterized tomographic model, the wet refractivity of each voxel is no more regarded as invariable but varies with the position. The wet refractivity of a generic point is expressed by a weighted mean of the wet refractivity values at the eight voxel nodes based on Newton–Cotes quadrature, where the point is located (Perler et al. 2011; Chen et al. 2020). In this study, we used the fourth-order Newton–Cotes quadrature to obtain the wet refractivity of each voxel from the wet refractivity values at the eight voxel nodes.

GNSS signals cannot pass through all voxels over the tomography area because of the nearly cone geometry of GNSS signals received at a specific ground-based station. Accordingly, there are many null values in the design matrix S (Bender and Raabe 2007; Benevides et al. 2016). Therefore, it is necessary to apply appropriate constraints so that the tomographic system can be inverted. The Gauss weighting method is commonly used for horizontal direction constraint (Song 2004; Jiang et al. 2014; Chen et al. 2014; Xia et al. 2018), whereas the vertical direction constraint is modeled by an exponential equation, taking account of water vapor in the vertical direction (which usually decreases with increasing height) as follows:

$$ N_{w} \left( h \right) = N_{C} \cdot e^{{\left( { - \frac{{h - h_{0} }}{{H_{z} }}} \right)}} , $$
(3)

in which \({N}_{w}\left(h\right)\) denotes the atmospheric wet refractivity at the height of h, \({H}_{Z}\) represents the height index of \({N}_{w}\), and \({N}_{c}\) is the constant value of wet refractivity. As can be known from Eq. (3), the relationship between atmospheric wet refractivity in adjacent vertical layers is as follows:

$$ \frac{{N_{w}^{i,j,k + 1} }}{{N_{w}^{i,j,k} }} = e^{{\left( {\frac{{h_{k} - h_{k + 1} }}{{H_{z} }}} \right)}} , $$
(4)

in which the superscripts I, j, and k are the indexes of the voxels in the east–west, north–south, and vertical directions, respectively. The \({h}_{k}\) is the height of the \(k\) th voxel. The tomography equation (Eq. 2) can be solved by adopting Kalman filtering based on horizontal and vertical constraints.

2.2 Tomography grid division

Generally, both the lower and upper limits of the tomographic grid refer to the height from the ground to tropopause (Flores et al. 2000; Chen et al. 2014). However, the wet refractivity is mostly clustered at a height that significantly below the tropopause. If tropopause is set as the top of the grid for the tomography inversion, the results could be unrealistically negative because the wet refractivity is very sparse near the tropopause (Flores et al. 2000). In order to define the upper limit appropriately, we attempted two schemes according to the zenith wet delay (ZWD) variations obtained from radiosonde and GNSS RO data. The first grid top is the upper limit of the tomography grid, but when the Nw value between the first and the second grid tops cannot be calculated based on the tomography equation, the radiosonde- and RO-derived Nw are used to take its place. As a result, when the height of the grid top decreases, the effective number of satellite rays increases. In the tomography, only the rays that penetrate into the grid from the top boundary are used for the tomography processing.

The GNSS tomography aims to reconstruct the vertical distribution of Nw. The division of the vertical grid severely affects the tomography solutions. Conventionally, two approaches have been used for dividing the grid: the uniform division (Flores et al. 2000; Xia et al. 2013, 2018) and non-uniform division (Perler et al. 2011; Rohm 2012; Jiang et al. 2014). Considering the practical distribution characteristics of Nw are sparse in high layers and dense in low layers, the non-uniform division is used here.

2.3 Calculation of temperature from Nw

The wet refractivity Nw of the troposphere is defined as:

$$ N_{w} = \left( {k_{2} - k_{1} \frac{{R_{d} }}{{R_{w} }}} \right) \cdot \frac{{{\text{P}}_{{\text{w}}} }}{T} + k_{3} \cdot \frac{{{\text{P}}_{{\text{w}}} }}{{T^{2} }}, $$
(5)

in which the constants are \({\text{k}}_{1}\text{=}\text{77.674 K}\bullet {\mathrm{hPa}}^{-1}\), \({\text{k}}_{2}= \text{7} \text{1.97 K}\bullet {\mathrm{hPa}}^{-1}\), and \(k_{3} = 3.75406 \times 10^{5} {\text{K}}^{2} {\text{hPa}}^{ - 1}\) (Rüeger 2002), \({\text{P}}_{\text{w}}\) is water vapor partial pressure in hPa, T is temperature in Kelvin, and \({\text{R}}_{\text{d}}\) and \({\text{R}}_{\text{w}}\) are the mean specific gas constant for dry air and water vapor, respectively.

Equation  (5) can be expressed as follows:

$$ T^{2} \cdot N_{w} - T \cdot \left( {k_{2} - k_{1} \frac{{R_{d} }}{{R_{w} }}} \right) \cdot P_{w} - k_{3} \cdot P_{w} = 0, $$
(6)

The \({\text{P}}_{\text{w}}\) and T are the variables to be estimated in Eq. (6), whereas the Nw has already been estimated by using Eq. (1). It is easy to know that Eq. (6) is underdetermined, and therefore, additional conditions are needed to obtain a proper solution of \({\text{P}}_{\text{w}}\) and T.

To this end, the vertical changes of temperature and water vapor partial pressure are modeled by linear and exponential (Callahan 1973) functions, respectively

$$ T_{i + 1} = T_{i} - \beta \cdot \left( {h_{i + 1} - h_{i} } \right) $$
(7)
$$ P_{w}^{i + 1} = P_{w}^{i} \cdot \exp \left( {a \cdot \left( {h_{i + 1} - h_{i} } \right) - b \cdot \left( {h_{i + 1} - h_{i} } \right)^{2} } \right) $$
(8)

in which \(\beta \) denotes the lapse rate of temperature, \({T}_{i}\), \({P}_{w}^{i}\), and \({h}_{i}\) denote the temperature, water vapor partial pressure, and height at the grid nodes at the ith levels, respectively. The \({T}_{i+1}\), \({P}_{w}^{i+1}\), and \({h}_{i+1}\) denote the corresponding terms at the (i + 1)th levels, respectively. The a and b are constants which can be obtained from numerical weather prediction products. We used Eq. (8) to fit a and b with the ‘lsqcurvefit’ function provided by MATLAB. The initial values of a and b are 0.248 km−1 and 0.048 km−2, respectively (Callahan 1973).

We first calculated the β, a and b using the fifth-generation reanalysis model (ERA5) products. We then took the temperature and water vapor partial pressure from ERA5 as the initial values. At last, we calculated the search ranges for water vapor partial pressure and temperature based on the bias between radiosonde data and ERA5. When the water vapor partial pressure and temperature satisfy Eqs. (6), (7), and (8) and reach the optimal set of values, they will be used as the best values in our search. Details on the calculation is provided in Supplement 3.

2.4 Calculation of the UHII

The UHII can be calculated as the temperature difference between urban and rural grid points inverted by the GNSS tomography:

$$ {\text{UHII}}_{{{\text{GNSS}}}} = T_{{{\text{GNSS}}}} \left( {{\text{urban}}} \right) - T_{{{\text{GNSS}}}} \left( {{\text{rural}}} \right), $$
(9)

in which \({T}_{\mathrm{GNSS}}\) are the temperature in Kelvins obtained from GNSS tomography.

The UHII obtained from the synoptic (SYN) temperature measurements were taken as references:

$$ {\text{UHII}}_{{{\text{SYN}}}} = T_{{{\text{SYN}}}} \left( {{\text{urban}}} \right) - T_{{{\text{SYN}}}} \left( {{\text{rural}}} \right), $$
(10)

3 Results and analyses

3.1 Data description

Since 2001, satellite-based GNSS RO has being providing substantial continuous temperature and pressure measurements with high accuracy, high vertical resolution, and global coverage. In this paper, we used the latest Wegener Center (WEGC) multi-satellite GNSS RO data, OPSv5.6, from May 2001 to December 2020. The WEGC OPSv5.6 dataset has been widely used for weather, climate, space weather, and geodetic studies as it provides quality-controlled global upper-air satellite measurements from multiple RO satellite missions, including CHAMP, GRACE, SAC-C, Formosat-3/COSMIC, and Metop (Schreiner et al. 2007; Anthes et al. 2008).

In addition, the radiosonde technique is a traditional tool for the meteorological measurement from the ground to the lower stratosphere. The ground-based GNSS observations were obtained from the Hong Kong Satellite Navigation System (HKSN) network, which is composed of 18 continuously operating stations with inter-station distances of 10–15 km as shown in Fig. 1. All the 18 stations are provided by “LEICA GRX1200 + GNS” receivers with a sampling rate of 5 s.

Fig. 1
figure 1

Geographical distribution of the 18 GNSS stations in Hong Kong and the horizontal grids (Xia et al. 2013, 2018; Jiang et al. 2014; Chen and Liu 2014) for tomography. The horizontal and temporal resolutions are set as 6 km and 30 min, respectively

Full size image

In this work, we collected the GNSS observations of all the 18 stations of the HKSN network in the year of 2020. Moreover, the GNSS RO profiles over Hong Kong and the radiosonde profiles observed at the WMO45004 station from 2010 to 2019 are used as historical data for the optimization of tomography. Supplement 1 provides the details of the RO events.

The GAMIT software (version 10.71, Robert 2023) was used in this study for the calculation of GNSS ZTD. For the simulation of near real-time GNSS tomography, the GNSS data were processed with a sliding session window with a width of 6 h and a step of 1 h.

A sliding time window strategy is a commonly used approach for the simulation of near real-time GNSS tomographic experiment (Foster et al. 2005; Herring et al. 2010). Moreover, providing a logical time interval is meaningful in the framework of the rain now-casting. Furthermore, a 6-h time interval for the minimum broken line length is recommended empirically to allow the linear fitting algorithm can conduct a better discretization of the PWV signal features without being affected by the noisy features (Benevides et al. 2015). Therefore, a 6-h interval’s time window is used, moving forward by an hour each time.

3.2 Definition of the grid top and optimization of vertical resolutions

In the zenith direction, the wet tropospheric delay can be expressed as:

$$ {\text{ZWD}} = 10^{ - 6} \mathop \smallint \limits_{{h_{s} }}^{\infty } N_{w} \cdot dh, $$
(11)

in which \(\mathrm{ZWD}\) is the wet tropospheric delay in meter, \({h}_{s}\) is the height of the observation station above mean sea level in meter, and \({N}_{w}\) is the atmospheric wet refractivity (unitless) which can be calculated from Eq. (5).

In this study, the ZWD was obtained from radiosonde data and RO profiles were used to define the grid top. The radiosonde sensors can measure several meteorological variables such as pressure, temperature, and relative humidity. Similarly, RO can also provide the profiles of temperature, water vapor pressure, etc. Taking the characteristics of exponential decreasing of the atmospheric refractivity into account, the formulas for wet delays from the radiosonde measurements and RO profiles can be derived:

$$\begin{aligned} ZWD &= 10^{ - 6} \sum\limits_{i} {\left[ {(h_{i} - h_{(i + 1)} )(N_{w}^{(i + 1)} - N_{w}^{i} )} \right]}\\ & \quad \mathord{\left/ {\vphantom {{\left[ {(h_{i} - h_{(i + 1)} )(N_{w}^{(i + 1)} - N_{w}^{i} )} \right]} {({\text{ln}}N_{w}^{i} - lnN_{w}^{(i + 1)} )}}} \right. \kern-0pt} {({\text{ln}}N_{w}^{i} - lnN_{w}^{(i + 1)} )}\end{aligned} $$
(12)

Afterward, SWD can be obtained from ZWD based on the wet Niell mapping function (Niell 1996; Chen and Herring 1997).

$$\begin{aligned} {\text{SWD}}\left( {\varepsilon ,{ }\alpha } \right) &= {\text{ZWD}} \cdot M_{w}^{{{\text{Niell}}}} \left( \varepsilon \right) + \frac{{1}}{{{\text{sin}}\left( \varepsilon \right) \cdot {\text{tan}}\left( \varepsilon \right){ + 0}{\text{.0007}}}}\\ &\quad \left[ {{\text{G}}_{{{\text{NS}}}} \cdot {\text{cos}}\left( \alpha \right){\text{ + G}}_{{{\text{EW}}}} \cdot {\text{sin}}\left( \alpha \right)} \right],\end{aligned} $$
(13)

in which \({M}_{w}^{\mathrm{Niell}}\) is the wet Niell mapping function; \(\varepsilon \) means the satellite elevation cut-off angle, and the satellite cut-off angle is set to 10° in this study, GNS is the north–south wet gradient component, and GEW the wet gradient component relative to the east–west direction, GNS and GEW were determined based on the numerical weather models (NWMs) using the ray-tracing method (Daniel and Johannes 2017), and α is the satellite azimuth. When the \({\mathrm{SWD}}_{h}\) is less than or equal to 1 mm, the corresponding height is defined as the first grid top (FGT). The difference between ZWD and \({N}_{w}\) between two adjacent time periods can be calculated:

$$ \Delta {\text{ZWD}}_{h}^{t} = \left| {{\text{ZWD}}_{h}^{t} - {\text{ZWD}}_{h}^{t + 1} } \right|, $$
(14)
$$ {\text{RMSE}}\left( {N_{w} } \right) = \sqrt {\frac{{\mathop \sum \nolimits_{i = 1}^{n} \left( {N_{w}^{{h_{i} ,t}} - N_{w}^{{h_{i} ,t + 1}} } \right)^{2} }}{n},} $$
(14)

in which \({\mathrm{ZWD}}_{h}^{t}\) and \({\mathrm{ZWD}}_{h}^{t+1}\) are the ZWD at height h at epoch t and t + 1, respectively. \({N}_{w}^{{h}_{i},t}\) and \({N}_{w}^{{h}_{i},t+1}\) are the height h at time t and t + 1; n is the number of the layers of radiosonde and RO data. When the \(\Delta {\mathrm{ZWD}}_{h}^{t}\) is less than or equal to 0.5 mm or \(\mathrm{RMSE}({N}_{w})\) is less than or equal to 0.5 N, the corresponding height is defined as the second grid top (SGT).

RO wet profile and radiosonde products at the WMO45004 station from 2010 to 2019 are used to determine the grid top based on Eqs. (13), (14) and (15). Individual radiosonde data, which have different vertical resolutions, are linearly interpolated to a 100-m vertical grid before the grid top identification. Then, we compute the daily mean for the first grid top and the second grid top as shown in Fig. 2.

Fig. 2
figure 2

The first grid top and the second grid top obtained by radiosonde and RO products

Full size image

It can be seen from Fig. 2 that FGT and SGT are about 1–2 km higher in summer than in winter. In addition, the fluctuations of the SGT are more significant than the FGT. With the considering of the uneven vertical distribution of water vapor in the addition tree, we developed an elaborate division scheme for the vertical layers as shown in Fig. 3. For the height range from surface to 1000 m, there are three layers with height of 300, 300, and 400 m. The height range from 1000 m to \({H}_{\mathrm{SGT}}\) is divided into layers with an identical height \(\Delta {H}_{\mathrm{S}}\) which is ranging between 400 and 600 m. In addition, the height range from \({H}_{\mathrm{SGT}}\) to \({H}_{\mathrm{FGT}}\) is also divided into layers with an identical height \(\Delta {H}_{\mathrm{F}}\) which is ranging between 600 and 1000 m. Usually, there are 11 layers between 1000 m to \({H}_{\mathrm{SGT}}\), whereas 6 layers between \({H}_{\mathrm{SGT}}\) and \({H}_{\mathrm{FGT}}\).

Fig. 3
figure 3

The elaborate division scheme for the vertical layers

Full size image

3.3 Obtaining the \({N}_{w}\) between FGT and SGT

The RO wet profile and Radiosonde product have been quality controlled, and Nw can be obtained from the water vapor pressure and temperature provided by the two meteorological data in Eq. (5). We add daily and semidiurnal terms to the annual and semiannual cycle variation characteristics of Nw, and the Nw time series obtained by the following equation which is layered for periodic fitting (from 5 to 11 km, it is divided into 12 layers on average, that is, a layer of 500 m).

$$\begin{aligned} N_{w}^{j} &= a_{0}^{j} + \mathop \sum \limits_{n = 1}^{2} a_{n}^{j} \cos \left( {2n\pi \frac{{{\text{doy}} - b_{n}^{j} }}{365.25}} \right)\\ &\quad + \mathop \sum \limits_{n = 1}^{2} d_{n}^{j} \cos \left( {2n\pi \frac{{{\text{hod}} - c_{n}^{j} }}{24}} \right),\end{aligned} $$
(16)

in which j is the number of layers;\({N}_{w}^{j}\) is the wet refractivity of the jth layer; doy is the annual cumulative day; hod is the UTC time; a0, a1 and a2 are the annual mean; annual cycle variation amplitude and semiannual cycle variation amplitude of Nw, respectively; b1 and b2 are the annual cycle variation initial phase and semiannual cycle variation initial phase, respectively; d1 and d2 are the daily cycle variation amplitude and semiannual cycle variation amplitude, respectively; c1 and c2 are the daily cycle variation initial phase and semiannual cycle variation initial phase, respectively. The values of an and bn in Eq. (16) at different altitude levels can be fitted by selecting the RO and radiosonde products in Hong Kong from 2010 to 2019.

SWD is the input value of the GNSS tomography technique, and its accuracy directly affects the accuracy of tomography-derived Nw. To evaluate the fitting accuracy of Nw between SGT and FGT, we selected the 2020 Hong Kong area radiosonde and RO products as the benchmark values. The difference of SWD can be obtained between benchmark value and Nw for different seasons derived from Eq. (16), as shown in Eq. (17). The statistical results are displayed in Table 1.

$$ \Delta {\text{SWD}}_{{{\text{SGP}}}}^{{{\text{FGP}}}} = \left( {{\text{ZWD}}_{{{\text{model}}}} - {\text{ZWD}}_{T} } \right) \cdot M_{w}^{{{\text{Niell}}}} \left( \varepsilon \right), $$
(17)

in which \({\mathrm{ZWD}}_{\mathrm{model}}\) denotes the ZWD obtained from Eq. (5) and Eq. (16); \({\mathrm{ZWD}}_{T}\) denotes the ZWD estimated from RO and radiosonde products using Eq. (5); ε is the elevation angle of satellite; \({M}_{w}^{\mathrm{Niell}}\) is the Niell wet projection function.

Table 1 Results of the average differences of SWD between model-derived and benchmark values at different cut-off angles (Unit mm)
Full size table

Table 1 shows that the value of \(\Delta {\mathrm{SWD}}_{\mathrm{SGP}}^{\mathrm{FGP}}\) is significantly larger at lower satellite cut-off angles. In addition, \(\Delta {\mathrm{SWD}}_{\mathrm{SGP}}^{\mathrm{FGP}}\) is significantly larger in summer than in other seasons, with a deviation of more than 5 mm. At satellite cut-off angles above 45°, the effect of \(\Delta {\mathrm{SWD}}_{\mathrm{SGP}}^{\mathrm{FGP}}\) is less than 1 mm, while at satellite cut-off angles below 30°, the effect of \(\Delta {\mathrm{SWD}}_{\mathrm{SGP}}^{\mathrm{FGP}}\) is greater than 1 mm.

4 Discussions

The radiosonde station WMO45004 was carried aloft once every 12 h in Hong Kong, and it was equipped with a configured sensor that collects information about temperature, pressure, relative humidity and so on. Here, the WMO45004 radiosonde products in 2020 were used to evaluate the accuracy of the tomography results. These products were further used to validate the accuracy of the temperature obtained from Nw.

4.1 GNSS tomographic results

We utilized GAMIT software to obtain the ZTD based on GNSS data from the Hong Kong SatRef network with IGS (International GNSS Service) ultra-rapid products orbit file. The Saastamoinen model and synoptic observation data were used to obtain the ZHD and then the ZWD was obtained by deducting the ZHD from the ZTD. The SWD was computed from ZWD using Niell wet mapping function. Finally, we estimated the 3-D distribution of atmospheric wet refractivity using parameterized approaches, in which Eq. (16) was used for deriving the Nw between SGT and FGT. Besides, Kalman filtering algorithm was used for tomography solutions.

To evaluate the new method, the tomography results were compared with those derived from the traditional tomography method, in which the atmospheric wet refractivity from SGT to FGT was estimated as unknown. The differences between the traditional and optimized tomography models are shown in Table 2. In the tomography, only the rays that penetrate into the grid from the top boundary were used for the tomography processing. The height of the grid top was decreased, conversely increasing the effective number of the satellite rays. First, the numbers of signals passing through the voxel (NSV) were compared when the optimized method and traditional method were used to invert the Nw given in Fig. 4. Then, tomography solutions were compared with external results derived from the radiosonde. The results are shown in Fig. 5.

Table 2 The processing strategy for traditional and optimized tomography models
Full size table
Fig. 4
figure 4

Comparison results of the number of signals passing through voxel between the traditional method and the optimized method. Opti refers to the optimized method; Trad refers to the traditional method

Full size image
Fig. 5
figure 5

Wet refractivity obtained from tomography-derived and radiosonde-derived data. RS is the wet refractivity derived using radiosonde products, Trad is the wet refractivity derived using the traditional tomography method, and Opti is the wet refractivity derived using the optimized method

Full size image

In Fig. 5, Trad is the wet refractivity derived using the traditional tomography method, and Opti is the wet refractivity derived using the optimized method. As shown in Fig. 4, the average number of NSVs per month of the optimized method is higher than that of the traditional method. From August to December, the average monthly NSVs of the optimized tomography are more than a thousand signals than the traditional tomography. Statistics show that the NSVs in the optimized technique are 5.8% better than that of the traditional technique. As can be seen in Fig. 5, there is a good agreement between the changing trends of wet refractivity with height across the tomography-obtained and data from radiosonde. However, in the case when the “inversion layer” is present, GNSS tomography fails to accurately represent in this situation. The wet refractivity derived from our optimized method is better than that from the traditional method since the blue curve is closer to the red curve. In Table 3, we present the deviation statistics for GNSS tomography-obtained and radiosonde-obtained wet refractivity over the whole year 2020.

Table 3 The statistical results between tomography-derived and radiosonde-derived wet refractivity (Unit: N)
Full size table

Table 3 provides statistical values of the differences between GNSS tomography-obtained and radiosonde-obtained results. As seen from the statistical results, the RMS and mean values of troposphere tomography using the optimized technique is less than that of the traditional method. Especially in summer, the optimized method is slightly better than other seasons. In addition, compared with the radiosonde data, the test results show that the wet refractivity quality obtained by the optimized technique is 16.5% better than that of the traditional technique.

4.2 Validation of temperature results

After obtaining the wet refractivity profile based on the GNSS tomography method, the temperature was estimated by the optimal search method using Eqs. (6), (7) and (8). The fifth-generation reanalysis model (ERA5) could provide temperature and water vapor partial pressure, which were selected as the initial values in this study. Since the temperature and water vapor pressure provided by ERA5 are inconsistent with the spatial and temporal resolution of the tomographic results, the Gaussian distance weighting function in the horizontal direction and the exponential function in the vertical direction are used to interpolate ERA5 to be consistent with them. In terms of time, the temperature and water vapor partial pressure of ERA5 can be interpolated by the Chebyshev function of order 9 (Press et al. 1992), which can achieve a time resolution consistent with the tomography results. Supplement 2 give the flowchart of data processing. Since our research area is Hong Kong, China, and the tallest building in this area is not more than 600 m, we only calculated the temperature at the vertices of each grid layer below 600 m. If determining the appropriate search range, it is crucial to find the range of percentage deviation between benchmark value and ERA5 product. Then, using the radiosonde product as the benchmark value, calculate the difference between the temperature and water vapor partial pressure provided by ERA5 and the radiosonde product below 600 m. This deviation can be formulated as follows:

$$ {\text{DT}} = \frac{{T_{{{\text{RS}}}} - T_{{{\text{ERA}}}} }}{{T_{{{\text{RS}}}} }} \cdot 100\% , $$
(18)
$$ {\text{DWP}} = \frac{{{\text{WP}}_{{{\text{RS}}}} - {\text{WP}}_{{{\text{ERA}}}} }}{{{\text{WP}}_{{{\text{RS}}}} }} \cdot 100\% , $$
(19)

in which the \({T}_{\mathrm{RS}}\) and \({\mathrm{WP}}_{\mathrm{RS}}\) are the temperature and pressure provided by radiosonde, respectively, and the \({\mathrm{T}}_{\mathrm{ERA}}\) and \({\mathrm{WP}}_{\mathrm{ERA}}\) are the temperature and water vapor partial pressure provided by ERA5, respectively. To study the range of percentage deviation of DT and DWP, we computed the situation in Hong Kong from 2010 to 2019 based on Equations (18) and (19).

Table 4 provides the statistics on the scope of DT and DWP in Hong Kong. If the ranges of DT and DWP are too large, some of the temperature and water vapor partial pressure may be over-corrected, but if the range of DN is too small, the temperature and water vapor partial pressure may be under-corrected. In this study, [− 0.75%,0.75%] in temperature and [− 7.5%,7.5%] in water vapor partial pressure are selected as the range of theoretical retrieval. Then, the theoretically retrieved value of atmospheric temperature is obtained at each layer as CT + CT·DT, where the search step size of DN is 0.25%. The theoretically retrieved value of atmospheric water vapor partial pressure is obtained at each layer as CWP + CWP·DWP, where the search step size of DN is 2.5%. Finally, the optical CT + CT·DT values are derived based on Eqs. (6), (7) and (8). Figure 6 gives the 3-D temperature distribution on Hong Kong below 600 m on April 2 and 3, 2020.

Table 4 Statistics of the change intervals of temperature and water vapor pressure between ERA5 and Radiosonde from 2010 and 2019. DT denotes the difference in temperature; DWP denotes the difference in water vapor partial pressure
Full size table
Fig. 6
figure 6

3-D distribution of atmospheric temperature below 600 m on April 2 and 3, 2020

Full size image

Figure 6 describes the atmospheric temperature changes at different heights. It shows that the atmospheric temperature tends to decrease significantly with elevation. In the horizontal direction, the temperature of the first layer does not change significantly over time, while the temperature of the second and third layers changes more obviously over time. In order to verify the accuracy of the inversion results of the temperature and water vapor pressure, we selected the radiosonde products in 2020 as the true value, and compared them with the inversion results corresponding to time and space. The statistical results are shown in Table 5

Table 5 Statistical results between GNSS-inverted and radiosonde-derived temperature and water vapor partial pressure below 600 m
Full size table

Table 5 provides the different maxima, minima, means and RMSs of GNSS-inverted and radiosonde-derived temperatures and water vapor partial pressures. In terms of the statistical results, the accuracy of GNSS-inverted temperature and water vapor partial pressure in autumn is better than other seasons. In addition, the best statistical accuracy of GNSS-inverted water vapor partial pressure is in winter while the worst is in summer. This can be attributed to summer and winter usually being the most and least humid seasons of the year, respectively. Supplement 4 provides a comprehensive evaluation of GNSS-inverted temperature (Fig. 7).

Fig. 7
figure 7

Distribution of selected GNSS stations and weather stations in Hong Kong. The blue circle indicates the GNSS station, and red circle indicates the weather station

Full size image

4.3 The urban heat island

The UHII is defined by the difference between the temperature in urban areas and surrounding rural areas. In urban areas, anthropogenic sources of heat are present, such as transportation and air conditioning equipment. In contract, the quantity and variety of anthropogenic heat sources are less in rural areas because there are few existing buildings and most of them are occupied by nature. It is common for rural and urban areas to be interdependent, with rural areas located outside of urban or city areas (Memon et al. 2009). In order to monitor the intensity of the UHI in Hong Kong, we selected several GNSS stations in the urban area (equipped with meteorological observation) as urban stations, and a weather station as a rural station which is located on a surrounding independent island. The distribution of the stations is shown as follows:

The daily maximum, minimum and average values have been obtained with meteorological data. We fitted these values into a second-order polynomial function separately. Thus, the maximum, minimum and average values of the UHII in meteorological data were calculated using Eq. (10). In addition, in order to validate the UHII in GNSS data, the temperature obtained by GNSS was interpolated to a same spatial and temporal resolution as the meteorological data using a second-order polynomial function. Similarly, the maximum, minimum and average of the UHII in GNSS data were calculated using Eq. (9). The results of one of the meteorological stations and the GNSS results that matched with meteorological stations are shown in Figs. 8 and 9, respectively.

Fig. 8
figure 8

UHII estimated with meteorological data between HKKT station and CC station. CC refers to Cheung Chau station

Full size image
Fig. 9
figure 9

UHII estimated with ground-based on GNSS observation data. GHKKT refers to the GNSS-inverted temperature matched with HKKT meteorological stations. GCC refers to the GNSS-inverted temperature matched with Cheung Chau meteorological stations

Full size image

In Fig. 9, GCC refers to the GNSS-inverted temperature matched with Cheung Chau meteorological stations. The range of the UHII obtained with meteorological data and GNSS data between HKKT and CC is shown in Figs. 8 and 9. The shape of the graphs obtained using both data is very similar. In summer, the UHII increases compared to winter. In addition, compared with meteorological data-derived UHI, the UHI obtained from GNSS data is smaller. Beyond that, the 5 pairs of meteorological and GNSS data were used for validation purposes, and the root mean square of the differences between rural and urban areas from meteorological data and GNSS data in different seasons are shown in Tables 6 and 7. Finally, the validation of the algorithm had been carried out by comparing the UHII which determined from GNSS data (UHIIGNSS) with the UHII which calculated from temperature sensors at weather stations (UHIISYN). The difference in intensity on a given day of the year (Diff_UHII(DOY)) had been compared using the following simple calculation:

$$ Diff\_{\text{UHII}}\left( {DOY} \right) = {\text{UHII}}_{{{\text{GNSS}}}} \left( {DOY} \right) - {\text{UHII}}_{{{\text{SYN}}}} \left( {DOY} \right) $$
(20)
Table 6 Average UHII from each season obtained using meteorological data in 2020 (Unit: K)
Full size table
Table 7 Average UHII from each season obtained using GNSS data in 2020 (Unit: K)
Full size table

The RMS values of the differences of the results obtained from GNSS and from meteorological which using both all data and seasonal data are shown in Table 8. The RMS values of the differences were used to validate the algorithm. The 5 pairs of meteorological and GNSS data used for validation purposes are clearly related to location, as described in Table 9.

Table 8 Relation of meteorological and GNSS pairs
Full size table
Table 9 RMS of the differences between UHII obtained with meteorological data and GNSS data in Hong Kong in 2020 (Unit: K)
Full size table

Tables 6 and 7 show the mean UHII of meteorological data and GNSS data in 2020 using 1-year data and data for each season. In all cases, the UHI intensity is the highest during spring, and the lowest during autumn. The mean UHII in different seasons is less than 0.6 K at the same station while the mean UHII of one-year data is less than 0.4 K. As shown in Fig. 8, all RMS of the differences between the UHII obtained with GNSS data and meteorological data are below 1.5 K. In addition, compared with meteorological data, the accuracy of the UHII is 1.20 K at a 95% confidence level using a full year of GNSS data.

5 Conclusion

Global navigation satellite system radio occultation provides high-precision middle and upper atmospheric parameter profiles (pressure, water vapor partial pressure and temperature). In this paper, historical radiosonde data and radio occultation data were used to optimize the ground-based GNSS tomography model to improve the accuracy of tomography-derived wet refractivity. After obtaining the wet refractivity, the ERA5 product was used as the initial value, and the search method was used to obtain the best temperature for the wet refractivity. The developed algorithm demonstrated the possibility of using GNSS data to monitor the UHII. Ground-based GNSS data can be used for micro- and meso-scale urban climate studies and has the following advantages: (1). the ground-based GNSS tomography technique works in all weather conditions, and its data are widely available as GNSS constellations are designed to cover the earth at all times; (2). GNSS data have a very high temporal resolution and can be processed in real time or near real time.

This study overcomes two major challenges in the algorithm development. The first challenge is the determination of the GNSS tomographic top grid height. Here, we obtained the SGT and FGT based on the RO data and radiosonde products in Hong Kong, and fitted the wet refractivity between FGT and SGT to a multi-order spherical harmonic function based on historical radiosonde and RO products. The height between the earth’s surface and SGT was divided into several layers, and the wet refractivity at the vertex of the voxels was used as an unknown parameter for GNSS tomography. While several voxels are also divided between SGT and FGT, and the wet refractivity at the vertex of voxels was directly obtained based on Eq. (16). Thus, the height of the grid top is decreased, conversely increasing the effective number of the GNSS satellite rays. Moreover, the number of unknowns in GNSS tomography can be reduced, and the accuracy of the tomography results can be improved.

The second challenge is the estimation of temperature from wet refractivity. Based on the relationship between wet refractivity and temperature and water vapor partial pressure, as well as the linear variation of temperature with elevation and the approximate exponential change of water vapor partial pressure with elevation, the optimal search method was used to obtain water vapor partial pressure and temperature from wet refractivity. After selecting five meteorological observing stations inside the city of Hong Kong as urban stations, and a station on an island in Hong Kong as a rural station, we used Eq. (10) to estimate the UHII as the benchmark value of the UHII obtained from GNSS data.

By using the data of 18 stations in Hong Kong in 2020 for a case study, the following conclusions are obtained:

  1. (1)

    Compared with the radiosonde data, the test results show that the wet refractivity quality obtained by the optimized technique is 16.5% better than that obtained by the traditional technique.

  2. (2)

    Using the radiosonde product as the benchmark value, the accuracy of the temperature obtained by GNSS data below 600 m is better than 1.35 K.

  3. (3)

    By solving the RMS of the differences between UHII obtained from GNSS data and meteorological data on the 5 selected locations, it has been shown that the difference of the UHII obtained from GNSS data and the measured UHII using temperature data in spring and summer is higher than other seasons, because the water vapor content is more abundant in these two seasons. Therefore, the water vapor partial pressure is not accurately calculated in spring and summer. The discrepancies between the HUI estimated by the algorithm and the UHII obtained from meteorological stations can be attributed to the lack of water vapor partial pressure data and GNSS processing. Another major reason is that we fitted the daily maximum, minimum and average temperature values obtained with meteorological data and GNSS-inverted temperature into a second-order polynomial function, respectively. Due to the limited fitting ability of the second-order polynomial function, this can result in the loss of some of the accuracy of UHII between GNSS-inverted and meteorological data-derived during the fitting process. The new algorithm can be used to monitor the diurnal cycle of the UHI.