Urban Heat Island monitoring with Global Navigation Satellite System (GNSS) data

The Urban Heat Island (UHI) effect occurs when the temperature in an urban area is higher than the temperature at a rural area. UHIs are monitored using remote sensing techniques such as satellite imagery or using temperature sensors deployed in a metropolitan area. In this chapter we propose a methodology to monitor the UHI intensity from Global Navigation Satellite Systems (GNSS) data. As the GNSS signal travels from the satellite to the receiver it propagates through the troposphere. A delay (Tropospheric delay) affects the signal. The delay is proportional to environmental variables. Also, the tropospheric delay in zenith direction (ZTD) is estimated as part of the Precise Point Positioning (PPP) technique. Therefore, in this chapter it is shown how to use process GNSS data to obtain ZTD and obtain temperature at an urban and a rural site simultaneously from the ZTD. The advantages of using GNSS data is its availability and many GNSS networks have been deployed in different cities so no need to deploy sensor networks. Furthermore, GNSS signal is less affected by bad weather conditions than satellite imagery.


The troposphere
The atmosphere is divided into several layers, including the ionosphere and the troposphere. The ionosphere is the ionized part of earth's upper atmosphere and is found from 60-1000 km above earth's surface. The troposphere is the layer closest to the earth's surface and is where all weather phenomena are located. with altitudes reaching up to more than 20 km at the equator, and 7 km or more at the poles. The troposphere contains approximately 75% of Earth's total atmospheric mass and 99% of its total water vapour [1]. Apart from water vapour, the troposphere comprises nitrogen (78%), and oxygen (20%) with the remaining 2% made up of other gases [1,2]. The composition of gases within the troposphere is essentially homogenous, except for water vapour which shows high temporal and spatial variability. There are no ionized gases in the troposphere; therefore, for electromagnetic waves in the radio-frequency spectrum (up to 15 GHz of fundamental frequency), the troposphere is a non-dispersive medium [3]. The environmental variables present in the troposphere are: pressure, water vapor and temperature.

Pressure in the troposphere
Pressure (P) is defined as force per unit of area [1]. According to the International Standard Atmosphere [4] at Mean Sea Level, the pressure is 1013.25 hPa. By definition 1013.25 hPa, equals a pressure of 1 kg/cm 2 of surface area. Due to gravity and decreased air density, air pressure decreases exponentially with increased height above the surface. The pressure at the surface of the earth can be measured using a barometer, and at different altitudes, using radiosondes.

Relative humidity (RH) in the troposphere:
Humidity is defined as the amount of water vapour in the air. It's temporal and spatial distribution is highly heterogeneous (highly variable) horizontally and vertically. For these reasons, it is hard to model or simulate with simple mathematical models. Furthermore, it has been discovered that humidity is found almost entirely within the first 10 km above the surface of the earth. Several measures are used to characterize water vapour: • Water vapour pressure-expressed in hPa or mbar.
• Absolute humidity-the amount of water vapour in the air, expressed in g/m 3 .
• Specific humidity-the ratio between the density of water vapour and the density of wet air. 3 • Relative humidity-the ratio of water vapour pressure to saturation vapour pressure, expressed in %. • Mixing ratio-the ratio of the density of water vapour to dry air.
Humidity is mostly measured as relative humidity, by sensors at surface level, such as weather stations. Another measurement typically available with humidity sensors is the Dew Point, which is the temperature at which enough water vapour is in the air for saturation to occur.

Temperature in the troposphere
The temperature (T) varies depending on several factors, including altitude, latitude, and time of the day. In the tropospheric region, temperature decreases with increasing altitude, at a rate of about -5 to -7 Kelvin / km. The decrease of temperature is due mainly to greater heat absorption by the sun-heated earth's surface, which via conduction then heats up the air closer to the ground. At a certain altitude-the first boundary of the tropopause-the temperature increases at a different rate to its decrease. After that, when the stratosphere layer is reached, the temperature remains constant. Temperature varies depending on the latitude and the day of the year. The yearly profile of temperatures at different latitudes in the northern hemisphere are shown in figure 1 and for latitudes in the southern hemisphere are shown in figure 2. As shown in Figure 1, the temperature has its maximum in the summertime, between the sixth and eighth months (between June and August) and the minimum is found between the second and third months (February and March). The stations in the arctic region (80 ºN and 70 ºN) exhibit extreme temperatures; during winter, the average temperature is around -30 ºC, while in the summer the temperatures reach an average of 0 ºC. On the other hand, stations near the equator show smaller temperature fluctuation between summer and winter. It can be seen in figure 1 that their temperatures hardly fluctuate through the year. Figure 2 shows the behavior of the temperatures throughout the year for stations in the southern hemisphere.  Figure 2 shows the yearly temperature profiles for southern hemisphere stations, which are the opposite of those for the northern hemisphere. The coldest temperatures are found between the 6 th and 8 th month. While the hottest temperature is found during the 12 th and first month. The behavior is the opposite in the southern hemisphere to the northern hemisphere.

Profile of P, T and RH vs Height
The profile of the variables in the troposphere, P, T and RH vs Height is shown in Figure 3.  Figure 3 shows that temperature decreases as altitude increases up to the tropopause. Then, the temperature increases at a different rate, until the second tropopause. After which the temperature remains constant. According to the International Standard Atmosphere [4] the average sea level temperature is 27 ºC, while the minimum at the tropopause is -64 ºC. Pressure decreases in an exponential fashion-as an inverse to altitude-and reaches almost 0 at a point in the tropopause layer. Humidity also behaves exponentially, decreasing up to 10 km altitude (this 10 km layer is where most of the moisture of air is concentrated) and then behaving in an irregular way through higher altitudes.

Effect of the Troposphere to GNSS signals
GNSS signals are susceptible to experiencing delays during their transmission from the satellite to the receiver. The troposphere induces a delay to the GNSS traveling signal called the tropospheric delay. In positioning applications, the tropospheric delay is translated into positioning errors. Therefore, it is important to understand this delay.
Mathematically, the propagation delay, (in meters), is defined as an integral of the refractive index, (n), of the media along the ray path, s, between the satellite W and the receiver R [3] [1], as shown in the next equation.
The refractive index (n) is defined as the ratio of the propagation velocity of the signal in a respective medium (v) and the propagation velocity of the signal in the vacuum (c). The tropospheric delay can be written in terms of the refractivity of the troposphere N as [5]: The refractivity of the troposphere (N) depends on environmental variables at the point of measurement. In the next equation. , and are empirically determined coefficients, is the air pressure in hPa, is the absolute temperature in Kelvin and is the water vapour partial pressure in hPa. Zd and Zw are unit-less compressibility factors for dry air and water vapour, respectively [6].
The compressibility factors are corrections to account for the deviation of atmospheric constituents from an ideal gas. For an ideal gas, the compressibility factors equal 1. For simplicity, it is assumed that the troposphere constitutes an ideal gas. Values for the empirically determined refractivity constants k1, k2 and k3 have been investigated by several authors, including Thayer [5], Smith and Weintraub [6] and Bevis et, al [7] -and their published values are listed in Table 1. Since each of the environmental parameters that are needed to calculate the refractivity of the troposphere using the previous equation depend on altitude, the refractivity profile will also depend on altitude.

Profile of the Refractivity vs Altitude.
In order to compute the profile of N vs altitude, 10 years of radiosonde data have been collected and processed for 12 stations at different latitudes [10]. Values at the same height have been averaged, and a single value per height has been used.
The values of p, T and e, at different altitudes were input into the definition of N to compute the refractivity at different heights. It has been found that all the profiles of the refractivity of the troposphere have a similar shape and values. The mathematical expression of the relationship between N and altitude is described in the following equation: a ) b ) 9 Where is the refractivity at height 0 m and is the ratio at which the refractivity decreases with altitude. A value for N0 and Nh can be obtained for different latitudes.

Tropospheric delay
The first term in the equation defining N is the refractivity caused by the induced dipole moment of the dry constituents of the atmosphere. The second term is the induced dipole moment of water vapour, and the third term shows the effect of the permanent dipole of the water vapour molecules [11]. Therefore, the tropospheric delay can also be separated into two components, the dry component, related to temperature and air pressure, and the wet component related to the amount of water vapour available in the troposphere. The dry component is relatively stable, while the wet component fluctuates and varies a lot [12]. Therefore, the tropospheric delay can be written as:

Zenith tropospheric delay
According to McCarthy and Pétit [13] the total delay in the line of sight between a GNSS receiver and a GNSS satellite is derived as the sum of the hydrostatic (or dry) and wet delays, in the zenith direction, multiplied by respective mapping functions and a gradient correction. Next equation represents the tropospheric delay, where and are the hydrostatic and wet zenith delays, respectively, with associated hydrostatic and wet mapping functions and . Symbol represents elevation angle of the satellite. 10 The term is called the tropospheric gradient correction and accounts for the azimuthal dependence of the tropospheric path delay, stands for the gradient mapping function, with respect to elevation angle . and denote the horizontal delay gradients in the north and east directions, respectively. is the azimuth angle of the received signal, measured from east to north. The tropospheric delay in the zenith direction is called the Zenith Tropospheric Delay (ZTD). It can be determined as an integral of refractivity N, in the zenith direction [14].
The ZTD is defined as the sum of the Zenith Hydrostatic Delay (ZHD) and the Zenith Wet Delay (ZWD) [5].
The tropospheric delay and the ZTD are related to each other using the following relation [5]: Where and are the hydrostatic and wet mapping functions, depending on the elevation angle. According to Vedel et al. [15], the ZTD can also be defined as the integral of the refractivity over a vertical column of the neutral atmosphere, as represented in following eqution. is the density of air, z is the geometric height, R is the gas constant, T is the temperature, zsite is the height of the receiver with respect to the ground and ztop is the height of the troposphere.

Estimation of ZTD
ZTD is a by-product of the Precise Point Positioning (PPP) technique [16] which estimates the positioning coordinates and other parameters. There are two main estimators used, the Extended Kalman Filter (EKF) [17] and the least squares adjustment [18]. A study of different software and online-services implementing the PPP technique is presented in Mendez Astudillo, et al. [19]. The ZTD is required as input for the algorithm, therefore, it can be estimated with a chosen software or online-service.

Urban heat island monitoring using GNSS data
The Urban Heat Island effect happens when an urban area is warmer than a rural area [20] because of the built structures releasing heat to the atmosphere and because of anthropogenic heat sources such as air conditioning units or transportation networks. The metric of the UHI effect is the Intensity of the UHI (UHII). The UHII can be monitored using GNSS data doing the following 4 steps: • Data collection simultaneously in an urban and a rural area. Therefore, it is necessary to classify the stations as being in an urban or a rural area. • Processing of GNSS data to obtain the Zenith Tropospheric Delay (ZTD) and location coordinates of the station. The Precise Point Positioning technique is used (PPP). • Calculation of temperature using ZTD in the urban and the rural station • Computation of the UHII.

Classification of the type of environment around the station
There are many definitions for urban and rural areas. Some definitions are in terms of economical activity, land cover or geographic location. Typically, urban areas are densely populated and densely constructed [21]. Also, in urban areas, there are anthropogenic sources of heat, such as the means of transportation and air conditioning units, among others. In contrast, a rural area has few built structures, is mostly covered by nature and the numbers and types of anthropogenic heat sources are reduced. Usually rural and urban areas are attached to each other, with rural areas found outside the city or urban area [21].
Stewart and Oke [22,23] have developed a classification to be used in urban climatology studies depending on land cover surrounding a station. Their technique is called Local Climate Zones, it defines 17 zones summarized in Table 2

Algorithm to calculate temperature from GNSS data
The algorithm developed to calculate temperature from GNSS data is shown in the block diagram in Figure 5. The algorithm requires values for the height of the troposphere (S), air pressure at the place of measurement (P) and the water vapor partial pressure (e) at the site of measurement. These values are obtained from radiosonde data and are called the universal parameters.

Calculation of Universal parameters:
Height of the troposphere (S) The height of the troposphere S can be calculated in two ways, using radiosonde data only or using GNSS and radiosonde data. As shown in Figure 1, the height of the troposphere can be defined in terms of the First and Second Lapse Rate Tropopause (LRT1 and LRT2). According to the World Meteorological Organization (WMO), the LRT1 is defined as the lowest level at which the lapse rate (change of temperature with height) decreases to 2 ºC / km or less, provided also that the averaged lapse rate between this level and levels within the next 2 km vertically does not exceed 2 ºC / km. If above the LRT1 the average lapse rate between any level and all higher levels with 1 km exceeds 3 ºC / km, then a second tropopause (LRT2) is defined. LRT2 may be either within or above the 1 km layer [24]. Radiosonde data is used to compute the vertical profile of the temperature and with the WMO's definition, the height of LRT2 is found.
A second technique used to calculate the height of the troposphere S is described in Mendez Astudillo et al. (2020). It requires the profile of the refractivity previously obtained with radiosonde data, N0 and Nh are needed as inputs (defined in Section 1). The ZTD is also an input into the following equation [25].

Pressure (P)
The pressure is obtained either from radiosonde sounding at its lowest level or from weather stations. In this algorithm, the pressure is considered a constant at a given latitude because its annual fluctuation is negligible.

Water vapor partial pressure at different heights e(z)
Water vapor partial pressure is the partial pressure of water vapor in any gas mixture in equilibrium with solid or liquid water. This variable cannot be obtained with direct measurements such as weather sensors or radiosondes. Therefore, many models have been developed to estimate this parameter. Models in term of ZWD can be found in Younes [26]. Another model based on the Clausius-Clapeyron relation for gases is Antoine's model. According to Tolman [17] Antoine's model is defined as: Where A, B and C are Antoine's constants (A=8.071, B=1730.63 and C=233.43). T(z) is the temperature in Kelvin [17] at heigth Z. Equation 11 is valid for temperatures greater than 0 ºC and lower than 100 ºC. It is assumed that urban temperatures in the metropolitan areas under study are between that range most of the time. Negative temperatures are assumed to yield 0 hPa of water vapor partial pressure.

Calculation of local parameters
The previous equations defined the refractivity of the troposphere in terms of the environmental variables surrounding the GNSS receiver. Next equation defines the profile of the refractivity in the troposphere. The troposphere causes a delay to the signal in the zenith direction (in meters) which can be expressed as an integral of the total refractivity N along propagation path s from receiver r to satellite w [5]: A local N0 (N0_local) is found using Nh, Ztrop (height of the tropopause), the altitude of the receiver (Zrec) and the ZTD in the following equation: N0_local is the value of the refractive index near the surface of the place of measurement.

Calculation of temperature
The temperature of the site is calculated equating to zero the equation of refractivity using local parameters and solving the quadratic equation:

Calculation of the UHII intensity (UHII) from GNSS data
The UHI intensity is measured by subtracting the temperature at an urban station with the temperature at a rural station. The UHII can only be calculated by comparing temperatures obtained at the same epoch. In some applications, hourly, daily or monthly averages are needed to be compared. Mathematically, the UHII obtained from GNSS data is defined as: where Tgnss (Location) is the temperature in ºC obtained with the algorithm at the chosen location. The validation of the algorithm is done by comparing the UHII obtained with GNSS data with the UHII obtained from meteorological data. A good indicator of the quality of the calculation of UHII is the Root Mean Square of the difference UHIIGNSS and UHIIMET. The algorithm here presented has been validated using Radiosonde, GNSS and meteorological data from the metropolitan area of Los Angeles, California. The Special Administrative Region of Hong Kong, China and Ningbo, China. It has been concluded that it is possible to estimate the UHII with GNSS data with an accuracy of 3ºC.