A Novel Approach for the Determination of the Height of the Tropopause from Ground-Based GNSS Observations

: In this paper, we present a new method to calculate the height of the second lapse-rate tropopause (LRT2) using GNSS high-precision data. The use of GNSS data for monitoring the atmosphere is possible because as the radio signals propagate through the troposphere, they are delayed according to the refractive index of the path of the signal. We show that by integrating the vertical proﬁle of the refractive index in the troposphere, we are able to determine the altitude of LTR2. Furthermore, as GNSS data is available from many stations around all latitudes of the globe and make up a network with high spatial and temporal resolution, we can monitor the diurnal cycle of the variables related to the refractive index of the path of the signal. A comparison between the heights of the LRT2 obtained with radiosonde data and with this novel method is presented in the paper, and it shows good agreement. The average di ﬀ erence found is ≤ 1 km for stations between the latitudes of 30 ◦ S and 30 ◦ N.


Introduction
Accurate weather data is required to correctly assess weather issues and to produce accurate models that help mitigate climate issues such as global warming. One atmospheric parameter that is important to monitor is the height of the tropopause. The tropopause is defined by the second and first lapse-rate tropopauses (LRT2 and LRT1, respectively). Changes in the height of the global tropopause can be an indicator of anthropogenic climate change [1][2][3]. Also, it is relevant to study the tropopause layer because it separates the troposphere and stratosphere, which have different characteristics with respect to chemical composition and static stability, which affect the tropopause layer [4]. Furthermore, the tropopause region plays an important role in the stratosphere-troposphere exchange [5].
The troposphere is the layer of the atmosphere closest to the earth's surface which extends to about 18 km above sea level. Between the troposphere and the stratosphere, there is a middle layer called the tropopause [6]. This layer starts with the LRT1 and ends at the LRT2 [7]. The next layer is called the stratosphere. Figure 1 illustrates the approximate height of each layer, as well as the temperature profile through the layers. The LRT1 can be defined in terms of the gradient of temperature against altitude. For example, Seidel and Randel [8], based on the World Meteorological Organization (WMO), defined the LRT1 as being at the altitude where the lapse rate of the temperature decreases to −2 °C/km, provided that within 2 km above this point, the average lapse rate does not exceed 2 °C/km. The LRT2 is defined as the point above the first tropopause, 1 km beyond which the average lapse exceeds 3 °C/km; any point higher than this 1-km layer, the lapse exceeds 3 °C/km.
In the profile of temperature showed in Rieckh et al. [9], as the altitude increases to 10 km, the temperature decreases up to 200 K (−73 °C). Then, there is an inflection point where the temperature remains almost constant as the altitude further increases; this point is defined as LRT1. It has been found that in the northern hemisphere between latitudes 20° and 40°, the height of the LRT1 is approximately 15 km above sea level [9]. According to the results obtained by Feng et al. [6], the LRT2 can be found at altitudes of around 14 and 25 km, depending on the latitude where the measurement is done. The tropopause layer is, therefore, defined as the distance between LRT2 and LRT1. The tropical tropopause is defined as being between the latitudes of 30°N and 30°S.
There are three main methods to obtain the height of the tropopause: the use of radiosonde data [4,6,8,10], the use of GPS Radio occultation (GPS-RO) [9], and the use of Very High Frequency (VHF) radars [11]. Radiosonde is the data source most commonly used. Other sources that have been used to study the tropopause are Atmospheric Infrared Sounder (AIRS) [12], Infrared Atmospheric Sounding Interferometer (IASI) [13], and the High-Resolution Infrared Radiation Sounder (HIRS).
Radiosonde data is readily available from stations around the globe; they provide environmental measurements such as temperature, humidity, pressure, etc. at different heights in one sounding. Radiosondes take measurements for the whole troposphere up to the geopotential height where the pressure is 100 hPa; however, radiosondes have a poor temporal resolution, because the measurements are only taken every 12 h. Its vertical resolution of 2 km requires the use of complex mathematical models to interpolate the temperature transition to determine the height of the tropopauses. In addition, the unique temperature profile along the altitude also increases the complexity of the mathematical models required to process the radiosonde data [8]. Because of a few potential inflection points in the temperature changes along the altitude, the mathematical model needs to describe a curve, rather than a linear relationship. The spatial resolution of the stations remains poor, even though developed algorithms are used to increase it.
The GPS-RO method provides near-vertical profiles of atmospheric thermodynamic variables with high vertical resolution (better than 1 km) and global coverage [9]. Using AIRS data, an error of approximately 4 K is found for temperatures around the tropopause [12]. In a study by Lerner et al. [13], the authors used IASI data and found that in the troposphere (i.e., below 200 hPa), the retrieved profiles exhibited temperature errors of less than 1 K and specific humidity errors of minus 10% at most heights, associated with a vertical resolution of 1.5-2 km. HIRS data provides a vertical resolution of 3-5 km in the troposphere. The LRT1 can be defined in terms of the gradient of temperature against altitude. For example, Seidel and Randel [8], based on the World Meteorological Organization (WMO), defined the LRT1 as being at the altitude where the lapse rate of the temperature decreases to −2 • C/km, provided that within 2 km above this point, the average lapse rate does not exceed 2 • C/km. The LRT2 is defined as the point above the first tropopause, 1 km beyond which the average lapse exceeds 3 • C/km; any point higher than this 1-km layer, the lapse exceeds 3 • C/km.
In the profile of temperature showed in Rieckh et al. [9], as the altitude increases to 10 km, the temperature decreases up to 200 K (−73 • C). Then, there is an inflection point where the temperature remains almost constant as the altitude further increases; this point is defined as LRT1. It has been found that in the northern hemisphere between latitudes 20 • and 40 • , the height of the LRT1 is approximately 15 km above sea level [9]. According to the results obtained by Feng et al. [6], the LRT2 can be found at altitudes of around 14 and 25 km, depending on the latitude where the measurement is done. The tropopause layer is, therefore, defined as the distance between LRT2 and LRT1. The tropical tropopause is defined as being between the latitudes of 30 • N and 30 • S.
There are three main methods to obtain the height of the tropopause: the use of radiosonde data [4,6,8,10], the use of GPS Radio occultation (GPS-RO) [9], and the use of Very High Frequency (VHF) radars [11]. Radiosonde is the data source most commonly used. Other sources that have been used to study the tropopause are Atmospheric Infrared Sounder (AIRS) [12], Infrared Atmospheric Sounding Interferometer (IASI) [13], and the High-Resolution Infrared Radiation Sounder (HIRS).
Radiosonde data is readily available from stations around the globe; they provide environmental measurements such as temperature, humidity, pressure, etc. at different heights in one sounding. Radiosondes take measurements for the whole troposphere up to the geopotential height where the pressure is 100 hPa; however, radiosondes have a poor temporal resolution, because the measurements are only taken every 12 h. Its vertical resolution of 2 km requires the use of complex mathematical models to interpolate the temperature transition to determine the height of the tropopauses. In addition, the unique temperature profile along the altitude also increases the complexity of the mathematical models required to process the radiosonde data [8]. Because of a few potential inflection points in the temperature changes along the altitude, the mathematical model needs to describe a curve, rather than a linear relationship. The spatial resolution of the stations remains poor, even though developed algorithms are used to increase it.
The GPS-RO method provides near-vertical profiles of atmospheric thermodynamic variables with high vertical resolution (better than 1 km) and global coverage [9]. Using AIRS data, an error of approximately 4 K is found for temperatures around the tropopause [12]. In a study by Lerner et al. [13], the authors used IASI data and found that in the troposphere (i.e., below 200 hPa), the retrieved profiles exhibited temperature errors of less than 1 K and specific humidity errors of minus 10% at most Remote Sens. 2020, 12, 293 3 of 18 heights, associated with a vertical resolution of 1.5-2 km. HIRS data provides a vertical resolution of 3-5 km in the troposphere.
The numerical weather prediction (NWP) and climate monitoring communities require a retrieval accuracy of less than 1 K in 1-km layers for temperature and better than 15% in 2-km layers for water vapor in the troposphere [13]. Thus, in order to improve the spatial resolution of the estimated thickness of the troposphere without increasing the complexity of mathematical models, an alternative method of measuring the height of the first and second lapse rate tropopauses is proposed based in Global Navigation Satellite System data (GNSS).In this paper, a new method is proposed to obtain a fine spatial resolution of both LRT1 and LRT2 at different latitudes from GNSS data.
Global Navigation Satellite System (GNSS) has been widely used in positioning [14,15] and attitude determination [16,17]. Also, GNSS data has been used to monitor the atmosphere [18]. The usage of GNSS signals for climate studies is gaining more interest as global warming becomes a serious challenge.
The method is novel because it uses of the Zenith Tropospheric Delay (ZTD), a by-product of the Precise Point Positioning (PPP) technique using Global Navigation Satellite System (GNSS) data. The ZTD can be used to determine the height of the troposphere because radio signals are delayed depending on the radio refractive index of the travelling path. Thus, by mathematically modeling the refractivity profile of the troposphere, it is possible to solve for the path of the signal which corresponds to the troposphere thickness. The main advantage of using the ZTD is that the temporal resolution is increased, as there is GNSS data every 30 s. Also, any GNSS receiver can be turned into a measuring station which makes it possible to increase the spatial resolution. There are GNSS receivers embedded in devices such as smartphones, car-navigation systems, etc. If the device supports carrier-phase tracking, the ZTD can be estimated from the PPP technique. The PPP positioning technique is widely implemented in many software and online services. A previous study [19] found that most online software available for PPP obtained very close values for the estimation of the ZTD. Thus, we can expect consistent results using any software.
The method presented in this paper makes it possible to estimate the height of the LRT1 with a precision of 500 m from GNSS data once an adjustment parameter (derived with radiosonde data) is known. In this study, only one year of data is used to test the method because in previous studies, such as the one by Mathew and Kumar [20], it was argued that once the height of the LRT had been calculated, it could be applied globally to data from different years.
The paper is structured as follows. Section 2 describes the data used to measure the height of the LRT1 and LRT2, i.e., radiosonde data and GNSS data, and the new method proposed to measure the height of the LRT1 from GNSS data is presented. Section 3 shows the results of the presented technique. Section 4 comprises discussion of the results. Section 5 presents the conclusion drawn from this analysis.

Description of GNSS Data Sources
The International GNSS Service (IGS) operates 505 GNSS stations across the world. Different products and observation data are provided by each station. Collectively, they provide data every day in a year (some stations do not provide data daily). They are made publicly available through their ftp site (ftp://cddis.gsfc.nasa.gov/gnss). Observation data with a frequency of 30 s is available for all stations managed by the IGS. The IGS also provides navigation data and ephemeris used in the PPP technique. The latitudes, longitudes, and heights of 15 selected GNSS stations covering the latitudes of 35 • N to 35 • S are given in Table 1 (values provided by the IGS). This data is used to determine LRT1 and LRT2.

Description of Radiosonde Data
The Integrated Global Radiosonde Archive (IGRA) run by the National Ocean and Atmosphere Administration (NOAA) (https://www.ncdc.noaa.gov/data-access/weather-balloon/integrated-globalradiosonde-archive) contains radiosonde data from more than 2500 stations around the world. Most stations contain data from two soundings every day: one at 00:00 Coordinated Universal Time (UTC) and one at 12:00 UTC. For the validation of the values of LRT1 and LRT2 obtained with GNSS data, 15 stations were chosen with the same latitude (because the height of the tropopause is highly associated with latitude) as the GNSS stations previously described. The latitude, longitude, and height of the radiosonde stations used in this study as reported by the IGRA are given in Table 2.  Figure 2 shows in a block diagram of the methodology used to calculate the LRT1 from GNSS data. Radiosonde data collected from stations described in Table 2 were used to construct the vertical profile of the refractivity of the troposphere and to calculate the height of LRT1 and LRT2 applying the Word Meteorological Organization (WMO)'s definition. The Zenith Tropospheric Delay ZTD is obtained by processing GNSS data in the Precise Point Positioning Technique. The ZTD is the input to the algorithm here presented together with the altitude of the receiver found using the PPP technique. Validation of the technique here presented was done by comparing the height of LRT2 obtained with GNSS data and with radiosonde data. Each of the blocks shown in Figure 2 are described in detail in the following subsection.

LRT1 and LRT2 from Radiosonde Data
According to the WMO, the tropopause layer can be defined by two points, the LRT1 and the LRT2, which can be obtained from the vertical profile of temperature and height obtained from radiosonde data [21]. The first tropopause is defined as the lowest level at which the lapse rate decreases to 2 °C/km [7]. If above the first tropopause, the average lapse rate between any level of all higher levels within 1 km exceeds 3 °C/km, then a second tropopause is defined by the same criterion as before [7]. This troposphere may be either within or above the 10 km layer. The height of the LRT1 is relevant for climate studies and studies where the depth of the troposphere closest to the earth ground is needed.

Profile of the Refractivity of the Troposphere
According to the literature [22][23][24][25], the refractivity of the troposphere (N) at height z depends on environmental variables at the point of measurement. The empirical relation used to describe the refractive index in terms of environmental variables is shown in Equation (1).
The constants k1, k2, and k3 have been empirically determined: k1 = 7.76 × 10 −1 K/Pa, k2 = 7.04 × 10 −1 K/Pa, and k3 = 3.739 × 10 3 K 2 /Pa [24]. p(z) is the air pressure in hPa at height, z, T(z) is the absolute temperature in Kelvin at height z, and e(z) is the water vapor partial pressure in hPa at height z [24]. Radiosonde data collected from stations described in Table 2 were used to construct the vertical profile of the refractivity of the troposphere and to calculate the height of LRT1 and LRT2 applying the Word Meteorological Organization (WMO)'s definition. The Zenith Tropospheric Delay ZTD is obtained by processing GNSS data in the Precise Point Positioning Technique. The ZTD is the input to the algorithm here presented together with the altitude of the receiver found using the PPP technique. Validation of the technique here presented was done by comparing the height of LRT2 obtained with GNSS data and with radiosonde data. Each of the blocks shown in Figure 2 are described in detail in the following subsection.

LRT1 and LRT2 from Radiosonde Data
According to the WMO, the tropopause layer can be defined by two points, the LRT1 and the LRT2, which can be obtained from the vertical profile of temperature and height obtained from radiosonde data [21]. The first tropopause is defined as the lowest level at which the lapse rate decreases to 2 • C/km [7]. If above the first tropopause, the average lapse rate between any level of all higher levels within 1 km exceeds 3 • C/km, then a second tropopause is defined by the same criterion as before [7]. This troposphere may be either within or above the 10 km layer. The height of the LRT1 is relevant for climate studies and studies where the depth of the troposphere closest to the earth ground is needed.

Profile of the Refractivity of the Troposphere
According to the literature [22][23][24][25], the refractivity of the troposphere (N) at height z depends on environmental variables at the point of measurement. The empirical relation used to describe the refractive index in terms of environmental variables is shown in Equation (1).
The constants k 1 , k 2 , and k 3 have been empirically determined: k 1 = 7.76 × 10 −1 K/Pa, k 2 = 7.04 × 10 −1 K/Pa, and k 3 = 3.739 × 10 3 K 2 /Pa [24]. p(z) is the air pressure in hPa at height, z, T(z) is the absolute temperature in Kelvin at height z, and e(z) is the water vapor partial pressure in hPa at height z [24]. e(z) is calculated using the following model where A, B, and C are Antoine's constants (A = 8.071 B = 1730.63 C = 233.43), and T(z) is the temperature in Kelvin [26].
Equation (2) is valid for temperatures greater than 0 • C. Lower temperatures produce ice fog instead of water vapor. In this research, two values of e are found, i.e., one at 00:00 UTC and another at 12:00 UTC, because radiosonde data is available at those times. Only temperatures above 0 • C are processed. Therefore, if a temperature is lower than 0 • C, the measurement is discarded, and the value of e is not updated. The profile of the water vapor partial pressure obtained using Antoine's model for temperatures between 1 • C and 50 • C (typical surface temperatures on Earth) is shown in Figure 3. e(z) is calculated using the following model where A, B, and C are Antoine's constants (A = 8.071 B = 1730.63 C = 233.43), and T(z) is the temperature in Kelvin [26].
Equation (2) is valid for temperatures greater than 0 °C. Lower temperatures produce ice fog instead of water vapor. In this research, two values of e are found, i.e., one at 00:00 UTC and another at 12:00 UTC, because radiosonde data is available at those times. Only temperatures above 0 °C are processed. Therefore, if a temperature is lower than 0 °C, the measurement is discarded, and the value of e is not updated. The profile of the water vapor partial pressure obtained using Antoine's model for temperatures between 1 °C and 50 °C (typical surface temperatures on Earth) is shown in Figure 3.  Figure 3 shows that at 1 °C, the e = 5 hPa and at 50 °C e = 120 hPa. The variation of the air temperature at locations near the equator is very small. Therefore, the change of e is very small. In this research, it is assumed that the water vapor partial pressure can be fitted into the curve shown in Figure 3.

Profile of the Refractivity
An example of the profile of the refractivity vs height estimated based on radiosonde data and Equation (1) (radiosonde USM00072376 located in the coordinates lat = 35.05°, long = −111°, height = 2179 m) is shown in Figure 4. The curve is fitted to an exponential function in order to find a mathematical relationship between the two variables.
The decrease of the refractivity as the geodetic height increases can be fitted to an exponential function. Thus, Equation (1) can be written as the empirically determined equation: where N0 is the refractivity at height = 0 and Nh is a factor that indicates the exponential decay of the refractivity with respect to the height.  Figure 3 shows that at 1 • C, the e = 5 hPa and at 50 • C e = 120 hPa. The variation of the air temperature at locations near the equator is very small. Therefore, the change of e is very small. In this research, it is assumed that the water vapor partial pressure can be fitted into the curve shown in  Figure 4. The curve is fitted to an exponential function in order to find a mathematical relationship between the two variables.

Estimation of ZTD
The troposphere causes a delay to the signal travelling in the zenith direction ZTD (in meters) which can be expressed as an integral of the total refractivity N along a propagation path from the satellite to the receiver. Therefore, the ZTD can be defined as [26]: The decrease of the refractivity as the geodetic height increases can be fitted to an exponential function. Thus, Equation (1) can be written as the empirically determined equation: where N 0 is the refractivity at height = 0 and N h is a factor that indicates the exponential decay of the refractivity with respect to the height.

Estimation of ZTD
The troposphere causes a delay to the signal travelling in the zenith direction ZTD (in meters) which can be expressed as an integral of the total refractivity N along a propagation path from the satellite to the receiver. Therefore, the ZTD can be defined as [26]: where z site is the height of the receiver with respect to sea level and z trop is the height of the tropopause. The integral is calculated from the geopotential height where the receiver is located (z site ) to the top of the troposphere (z trop ). z trop defines the layer of the atmosphere that causes the tropospheric delay to the travelling GNSS signal. It corresponds to the height of the LRT2. In terms of Equation (3): where N 0 and N h are the parameters of the refractivity of the tropopause defined in Equation (3). z is the height of the LRT2. The path of the signal causing the tropospheric delay is defined between z site and the height of LRT2, as shown in Figure 1. The Precise Point Positioning (PPP) technique is a method that performs positioning determination using a single GNSS receiver [27]. Double frequency GNSS pseudorange and carrier-phase data is used as an input to mostly eliminate the effect of the ionosphere on the GNSS signal [28]. The technique requires precise orbit and clock products which are available from the International GNSS Service (IGS) or other sources at different resolutions [28]. The PPP algorithm requires an estimator which is typically the Least Squares [29] or the Extended Kalman Filter [31]. In both cases, raw GNSS data is the input, and the output is the estimated location coordinates and other parameters such as tropospheric delay and ambiguities.
In this research, the ZTD has been obtained as a by-product of the Precise Point Positioning (PPP) technique [27,29] implemented with the open-source software RTKLIB [31]. The PPP implementation in RTKLIB uses the Extended Kalman Filter EKF [30] to estimate the position of the receiver (x, y, and z in the Cartesian coordinates system, which can be converted to latitude, longitude, and height in the Geodetic coordinates system) and other parameters such as ZTD and ambiguities [30]. The inputs are GNSS observation data and precise ephemeris obtained from thIGS. RTKLIB developed by T. Takatsu [31], has been chosen as the processing software because in a previous study [14], it was found that the estimations from it were very near those obtained using other software. Furthermore, open-source RTKLIB allowed us modify the code for our development, and it has the flexibility of using different orbit and clock products. In this research, GNSS observation data from the whole year 2018 from the stations described in Table 1 were used to test the proposed algorithm.

Computing LRT2 from GNSS Data
The height of LRT2 is the whole path through which the GNSS signals propagate the troposphere; it is equivalent to the height of the tropopause. Therefore, from Equation (5), the height of tropopause z trop can be calculated from the ZTD by using the derivative of Equation (5). Reorganizing Equation (5), z trop is obtained in terms of ZTD and the profile of the refractivity of the troposphere: N h and N 0 are determined from the profile of refractivity at the site of measurement from radiosonde data. The ZTD and z site are the Zenith Tropospheric Delay estimated from GNSS data and the height of the station, respectively.

Results
In this section, the outcomes of each of the steps of the algorithm described in Figure 2 are shown. The first outcome shown here is the annual cycle of the refractivity profile, in order to show the diurnal cycle of that profile (N against height). Also, the diurnal cycle of the N0 and Nh is shown. Then, in Section 3.2, the definition of a parameter called D1 is defined. D1 is used to calculate LRT1 from LRT2. In Section 3.3, the LRT2 calculated from GNSS data is shown in two different time frames: one week and the whole year, 2018. In Section 3.4, the results are validated by comparing them to the LRT2 values obtained with radiosonde data and the WMO definition of LRT1 and LRT2. A statistical analysis is presented in order to validate the algorithm.

Annual Cycle of the Refractivity Profile
The annual cycle of the refractivity profile was calculated using 10 years of radiosonde data, as described in Table 2. In Figure 5, the values of N 0 and N h for station USM00072376 are shown. A clear periodicity is seen. We compared all years and found that the variations between years is small (a variation of the height of the tropical tropopause of 0.5-1.2 Km has been found by Seidel and Randel [8]). This confirms that we can use radiosonde data for the slow-varying cycle, and GNSS data for the fast-varying cycle, (diurnal cycle). The numerical values are unique for each station, but an annual cycle is found for all stations. So, in order to compute Equation (6), we use the closest N 0 and N h to the DOY when the ZTD was determined. Diurnal Cycle of N 0 and N h Figure 5 indicates a yearly periodicity of the profile of N 0 and N h with data from radiosonde station USM0006727376 (Northern Hemisphere). During the summer months, the values of N 0 are the highest, and those of N h are the lowest. In order to analyze the diurnal cycle of N 0 and N h , the same radiosonde data was used. The data was fitted to an exponential in order to visually display the behavior of N 0 and N h during a single day. Figure 6 shows the profile of the fitted N 0 and N h from radiosonde data over four days in the summer of 2018.  The diurnal cycle of the profile of the refractivity shows two inflection points: the maximum is found every day at 00:00 Coordinated Universal Time (UTC) and the minimum at 12:00 UTC. From 00:00 UTC to 12:00 UTC, the values of N0 decrease an average of 40 units for N0. The behavior of Nh is the inverse to the behavior of N0, i.e., the minimum value is found at 00:00 UTC and the maximum value is found at 12:00 UTC. Most radiosonde stations described in Table 2 provide data at both times.

D1 Parameter
The use of Equation (5) yields the height of LRT2 directly because the whole troposphere causes a tropospheric delay. Therefore, in order to calculate the LRT1 with GNSS data, an adjustment parameter D1 is needed which corresponds to the thickness of the tropopause obtained with radiosonde data [6]  The diurnal cycle of the profile of the refractivity shows two inflection points: the maximum is found every day at 00:00 Coordinated Universal Time (UTC) and the minimum at 12:00 UTC. From 00:00 UTC to 12:00 UTC, the values of N 0 decrease an average of 40 units for N 0 . The behavior of N h is the inverse to the behavior of N 0 , i.e., the minimum value is found at 00:00 UTC and the maximum value is found at 12:00 UTC. Most radiosonde stations described in Table 2 provide data at both times.

D1 Parameter
The use of Equation (5) yields the height of LRT2 directly because the whole troposphere causes a tropospheric delay. Therefore, in order to calculate the LRT1 with GNSS data, an adjustment parameter D1 is needed which corresponds to the thickness of the tropopause obtained with radiosonde data [6] (i.e., D1 = LRT2 RS − LRT1 RS ). Therefore, LRT1 GNSS = LRT2 GNSS − D1.

Results of Proposed Algorithm
GNSS data is available every 30 seconds every day of the year for most of the stations presented in Table 1. However, in some cases, some data is missing. All available GNSS data of the selected stations and periods was processed to test the algorithm. Daily, there are 2880 epochs of data, which were processed with RTKLIB. Each processed epoch produced a value for ZTD and for z site , which were used to test the proposed algorithm. The results are presented at three different time frames; first, one week in the month of July for two stations are presented as an example of how the results can be used to study the diurnal cycle of the height of the tropopause. Second, the yearly profile of the height of the tropopause obtained with the algorithm is shown for two stations as an example. Finally, the height of LRT1 and LRT2 obtained from radiosonde data and from the algorithm are shown for all latitudes.
As an example, the weekly results of stations JPLM, GMSD, DJIG, NAUR, ZAMB, and NNOR during the week of the 9 to 15 September (autumn), 2018 are presented in Figure 7. In the same figure, the LRT2 obtained from radiosonde data for the same time period is presented. Regarding the GPS data, the week is GPS week 2018, days of the year (DOY): 252-258. The diurnal cycle can only be obtained with GNSS data because radiosonde data is available only twice a day.
were processed with RTKLIB. Each processed epoch produced a value for ZTD and for zsite, which were used to test the proposed algorithm. The results are presented at three different time frames; first, one week in the month of July for two stations are presented as an example of how the results can be used to study the diurnal cycle of the height of the tropopause. Second, the yearly profile of the height of the tropopause obtained with the algorithm is shown for two stations as an example. Finally, the height of LRT1 and LRT2 obtained from radiosonde data and from the algorithm are shown for all latitudes.
As an example, the weekly results of stations JPLM, GMSD, DJIG, NAUR, ZAMB, and NNOR during the week of the 9 to 15 September (autumn), 2018 are presented in Figure 7. In the same figure, the LRT2 obtained from radiosonde data for the same time period is presented. Regarding the GPS data, the week is GPS week 2018, days of the year (DOY): 252-258. The diurnal cycle can only be obtained with GNSS data because radiosonde data is available only twice a day.  Figure 7 shows the diurnal cycle of the height of the tropical tropopause obtained with GNSS data and the algorithm presented in the previous section. In all stations, during the week, rapid variations of the height of the tropopause can be seen. For some days, there is a clear diurnal cycle: at the beginning of the day, the height is higher than at the end of the day. However, some days exhibit more chaotic variation, and the diurnal cycle is less obvious.  Figure 7 shows the diurnal cycle of the height of the tropical tropopause obtained with GNSS data and the algorithm presented in the previous section. In all stations, during the week, rapid variations of the height of the tropopause can be seen. For some days, there is a clear diurnal cycle: at the beginning of the day, the height is higher than at the end of the day. However, some days exhibit more chaotic variation, and the diurnal cycle is less obvious. Figure 8 shows the annual profile of the height of LRT2 for stations JPLM, GMSD, DJIG, NAUR, ZAMB, and NNOR obtained with GNSS data and the presented algorithm.
(e) (f)  Figure 7 shows the diurnal cycle of the height of the tropical tropopause obtained with GNSS data and the algorithm presented in the previous section. In all stations, during the week, rapid variations of the height of the tropopause can be seen. For some days, there is a clear diurnal cycle: at the beginning of the day, the height is higher than at the end of the day. However, some days exhibit more chaotic variation, and the diurnal cycle is less obvious. Figure 8 shows the annual profile of the height of LRT2 for stations JPLM, GMSD, DJIG, NAUR, ZAMB, and NNOR obtained with GNSS data and the presented algorithm. According to Figure 8a-f, Stations JPLM and GMSD follow the same pattern because they are further away from the equator. During the summer months, the height of LRT2 is minimum. Station NAUR shows a very small variation, i.e., the height of LRT2 seems to be constant except for three peaks. Station NAUR is located in Nauru, which is on the equator; therefore, the height of LRT2 is expected to be higher than other locations and the weather to be almost constant. Station NNOR is south of the equator; therefore, the maximum values for the height of LRT2 are found during the According to Figure 8a-f, Stations JPLM and GMSD follow the same pattern because they are further away from the equator. During the summer months, the height of LRT2 is minimum. Station NAUR shows a very small variation, i.e., the height of LRT2 seems to be constant except for three peaks. Station NAUR is located in Nauru, which is on the equator; therefore, the height of LRT2 is expected to be higher than other locations and the weather to be almost constant. Station NNOR is south of the equator; therefore, the maximum values for the height of LRT2 are found during the summer. The behavior is the opposite to that of stations north of the equator.
Equation (6) is applied to estimate the height of the tropical tropopause using GNSS data; however, radiosonde data is also used for validation purposes. The height of the tropical tropopause measured using annual GNSS and radiosonde data is shown in Figure 9.
(e) (f) According to Figure 8a-f, Stations JPLM and GMSD follow the same pattern because they are further away from the equator. During the summer months, the height of LRT2 is minimum. Station NAUR shows a very small variation, i.e., the height of LRT2 seems to be constant except for three peaks. Station NAUR is located in Nauru, which is on the equator; therefore, the height of LRT2 is expected to be higher than other locations and the weather to be almost constant. Station NNOR is south of the equator; therefore, the maximum values for the height of LRT2 are found during the summer. The behavior is the opposite to that of stations north of the equator.
Equation (6) is applied to estimate the height of the tropical tropopause using GNSS data; however, radiosonde data is also used for validation purposes. The height of the tropical tropopause measured using annual GNSS and radiosonde data is shown in Figure 9.  Figure 9 shows the height of LRT2 obtained with radiosonde and with GNSS data. The LRT2 was calculated using the algorithm previously described with GNSS data and applying the WMO's definition of the LRT2 to the radiosonde soundings in the stations described in Table 2. Every day, there are 2880 epochs of GNSS data. Therefore, a daily average is calculated; these daily averages are themselves averaged to find a single value at a given station in a year, which is plotted in Figure 9.  Figure 9 shows the height of LRT2 obtained with radiosonde and with GNSS data. The LRT2 was calculated using the algorithm previously described with GNSS data and applying the WMO's definition of the LRT2 to the radiosonde soundings in the stations described in Table 2. Every day, there are 2880 epochs of GNSS data. Therefore, a daily average is calculated; these daily averages are themselves averaged to find a single value at a given station in a year, which is plotted in Figure 9. There are two radiosonde soundings every day. Therefore, a daily average of the obtained height of LRT2 is calculated and plotted per station in Figure 9.

Validation of Results
Validation of the algorithm was done by comparing the results of LRT2 obtained with GNSS and radiosonde data. The results are shown in Figure 9. The differences of 1.4 km between 30 • S and 30 • N with a 95% confidence interval were obtained by doing a P-value test. Data from radiosondes were averaged and a single value per day was used. Figure 10 shows the residuals of the comparison of LRT2 obtained with GNSS and radiosonde data.
Validation of the algorithm was done by comparing the results of LRT2 obtained with GNSS and radiosonde data. The results are shown in Figure 9. The differences of 1.4 km between 30°S and 30°N with a 95% confidence interval were obtained by doing a P-value test. Data from radiosondes were averaged and a single value per day was used. Figure 10 shows the residuals of the comparison of LRT2 obtained with GNSS and radiosonde data.  Table 3 shows the height of LRT2 measured with radiosonde and GNSS data from the year 2018 (LRT2_RS and LRT2_GNSS). Also, the mean of the differences, the standard deviation, the maximum and minimum values, the median, and the RMSE of the differences are shown in Table 3 for validation purposes.  Table 3 shows the height of LRT2 measured with radiosonde and GNSS data from the year 2018 (LRT2_RS and LRT2_GNSS). Also, the mean of the differences, the standard deviation, the maximum and minimum values, the median, and the RMSE of the differences are shown in Table 3 for validation purposes.

Discussion of Results
GNSS data can be used to calculate the height of LRT2 and LRT1. Furthermore, it is possible to detect the diurnal cycle of the height of LRT2. The results in Figure 7 show the diurnal cycle of the height of the tropopause found with GNSS data and the algorithm presented in this paper. Stations JPLM and GSMD are located at latitudes 35 • and 30 • north of the equator. For both stations, during the morning, the values were higher than during the afternoon. That's the reason why a peak was found at the beginning of the day. The average height found for JPLM was 13.4 km, while for GSMD, it was 16.5 km, as expected, because latitudes closer to the equator have higher values for the height of LRT2.
Stations DJIG and NAUR are closer to the equator: DJIG is located at latitude 10 • and NAUR is on the equator. Figure 7c,d show that the values for the height of LRT2 are bigger for NAUR than for DJIG because, according to previous research, the highest values of height of LRT2 are found at the equator. The profile of diurnal cycle shows a little change between heights in one day. Stations ZAMB and NNOR are represented in Figure 7e,f respectively. Both stations are south of the equator, at latitudes of −15 • and −30 • respectively. The maximum values were found during the afternoon rather than in the morning. As expected, the heights were lower than those for the stations near the equator.
The heights of the LRT2 obtained with radiosonde and GNSS data for six stations are shown in Figure 8. Figure 8a−f show the annual profile of the height of the tropopause found with GNSS data. Figure 8a,b show the profile for stations JPLM and GSMD at latitude 35 • and 30 • north of the equator. In both cases, the minimum values were found during the summer months. The maximum values were found during the spring. GSMD lacks data from July and August 2018, but the tendency is clear for lower values. The values for the height were lower at station JPLM than GSMD. Figure 8c,d show the annual profile of the height of LRT2 at stations DJIG and NAUR, located at latitudes 10 • and on the equator. As expected, the values for NAUR were the greater than those for all other stations because it is on the equator. Furthermore, the profile is almost flat. DJIG has bigger changes than NAUR; the maximum values for this latitude are found during the months of August and September, which correspond to the summertime in the northern hemisphere. Figure 8e,f show the annual profile of the height of LRT2 obtained with GNSS data at stations in the southern hemisphere ZAMB (lat: −15 • ) and NNOR (lat: −30 • ). The profile at ZAMB has the opposite behavior to the profile of stations north of the equator. The maximum values were found during the months of February and March, which are the summer months in that hemisphere, and the minimum values were found in the months of August and September, which are the winter months in the southern hemisphere. The values of the height of LRT2 at NNOR were smaller than for ZAMB, and the profile showed small changes between the maximum and minimum values. The maximum values were found during June and July. In contrast, the minimum values were found during the months of April and May. The averaged values of LRT2 at different latitudes yielded the results presented in Figure 9, which show that the maximum values were found for stations closer to the equator with both GNSS and radiosonde data.
Validation of the results was done by comparing the results obtained with radiosonde data and with GNSS data. Both outcomes are shown in Figure 9. The quality of the results was evaluated with statistical indicators. In Table 3, the following statistical indicators are shown: the average of the differences, the maximum and minimum values, the standard deviation of the differences, and the median and the Root Mean Square Error of the differences. Furthermore, a P-value test found that the difference in the 95% confidence interval for stations between −35 • and 35 • latitude is 2.7 km, which is close to the resolution obtained using other techniques and data sources. Seven stations had an average difference of less than 1 km, five stations had an average difference between 1 and 2 km, and the remaining three stations had an average difference between 2 km and 3 km. However, none of the stations had an average greater than 3 km. The discrepancies between the results are due to the model of the e, and to errors in the estimation of ZTD, which are also linked to environmental variables not being accurately measured or estimated. Furthermore, missing radiosonde or GNSS data also contributed to the discrepancies.
In order to reduce discrepancies, one year of data was processed so that different weather conditions would be present, especially in the stations further away from the equator which have different weather conditions through the year. Ten years of radiosonde data were processed to obtain the equation here presented. The profile of the refractivity through the years changes very little; therefore, we used only one year of data to test the algorithm and assumed that in the other years, the results would be similar. A limitation of the algorithm is that special weather conditions can increase the discrepancies between the obtained height of the LRT2 with GNSS data and the measured LRT2 from radiosonde data.
The advantage of the proposed technique is that the temporal resolution is increased. The ZTD can be estimated with the PPP technique as frequently as the ephemeris and the GNSS data are available. In this case, there is a solution available every 30 seconds. Final precise ephemeris was used in our test; therefore, the solutions of the proposed algorithm is postprocessing with the latency of the availability of the final precise ephemeris product from IGS. Rapid precise ephemeris and improvements of the PPP technique towards real time would make the computation of ZTD possible in real time. Therefore, the height of LRT2 could be measured in real time. In contrast, the use of radiosonde data is limited by the number of soundings, i.e., two per day in most radiosonde stations. Another advantage of the proposed algorithm is that it is possible to observe the diurnal cycle of the height of LRT2 at a higher resolution, because there are GNSS measurements throughout whole the day.

Conclusions
The Zenith Tropospheric Delay (ZTD) can be used to monitor the environment because of two relationships, i.e., the relationship between the ZTD and the refractivity index of the troposphere (N), as defined using Equation (4), as well as the relationship between N and environmental variables such as temperature, water vapor partial pressure, and atmospheric pressure. Previously, algorithms have been developed to monitor the water vapor and the precipitable water using ZTD. However, this research showed that monitoring the height of the tropopause layer using the ZTD as input is also possible because of the relationship of the ZTD and the path of the signal which passes the tropopause.
The profile of refractivity was obtained with the measurements of pressure (p) and temperature (T) at different levels. It was assumed that the water vapor is distributed in the troposphere, so it was calculated with the first temperature reported by the radiosonde; the same value was used at different heights. The profile was fitted into an exponential function, and two values were found: N 0 and N h . Ten years of data were processed, and it was found that the profiles of N 0 and N h were periodical at each latitude. Therefore, a single average value of N 0 and N h per latitude was obtained.
In this paper, we have presented an algorithm to calculate the height of the tropical tropopause using high-resolution GNSS data. The algorithm was tested using high-resolution GNSS data from 15 stations between latitudes of 30 • S and 30 • N. The algorithm was validated using radiosonde data from stations in the same latitudes. It was found that the profile of refractivity against height is periodical through the years. Also, it was found that the annual change of the values N 0 and N h are negligible; therefore, those parameters can be calculated with data from a single year and applied to different years.
The profile of the height of the tropopause (LRT1 and LRT2) can be obtained with GNSS data, as shown in Figure 9. Validation was done by comparing the results obtained with GNSS and radiosonde data, and the definition of the WMO. It was found that the average difference between both results is less than 3 km for all stations; five stations yielded averaged differences of less than 1 km. Furthermore, it is shown in Figure 8 that it is possible to monitor the diurnal cycle of the height of the tropopause using GNSS data because it has a higher temporal resolution. The ZTD can be estimated from PPP at a frequency of 30 s. Therefore, the height of LRT2 can be estimated at the same frequency as the ZTD.