Comparison of Radiosonde Measurements of Meteorological Variables with Drone, Satellite Products, and WRF Simulations in the Tropical Andes: The Case of Quito, Ecuador

Abstract

Radiosondes are the most widely used method for studies of vertical atmospheric behavior, but the high costs associated, and the logistic limitations have forced researchers to look for alternative methods for atmospheric profiling, such as lidar and satellite measurements, or modeling. However, the assessment of the accuracy of alternative methods is recommended, especially in complex terrain, such as the tropical Andes. In this research, the atmospheric profiling of satellite data from AIRS and MODIS products, simulations of the Weather Research and Forecasting model, WRF, and drone measurements are evaluated for a campaign of 10 radio soundings, between August 2021 and January 2022. Additionally, the capability to capture the planetary boundary layer height, hPBL, is studied. The measurements were conducted at Izobamba station near Quito, Ecuador. Temperature, T, Dew Point Temperature, TD, Mixing Ratio, Q, and Potential Temperature, PT, were evaluated from 0 to 300 m above ground level (magl.) for satellite, WRF, and drone data, and from 0 km to 15 km for satellite and WRF data. Additionally, the capability to capture the planetary boundary layer height, HPBL, was assessed. The results show that drone profiles best represented the magnitude of the analyzed variables showing mean RMSE of 0.79 for T, but the noise of the measurements caused a low correlation with radio sounding profiles, which was partially corrected with a quadratic fit on the profile. The WRF results achieved a positive representation in terms of correlation, but error metrics show that there are remarkable differences in magnitude in the first 300 magl., up to the tropopause height, which surpasses satellite representations for all variables. The MODIS profiles do not generally perform well due to their low vertical resolution and limitations with cloud coverage. However, AIRS data, despite its low resolution, show a better representation of vertical profiles than MODIS, for T and TD, surpassing WRF simulations in some dates. For the HPBL, the WRF results show that physical and atmospheric conditions limit its determination, and the methods and conditioning factors should be further analyzed.

1. Introduction

Meteorological data are the foundation of our knowledge of the atmospheric system. There cannot be a feasible pathway to understanding weather, climate processes, variability, extremes, and climate change, without long records of observations [1]. Increasingly, observations and measurements of atmospheric, geophysical, and environmental variables are necessary for atmospheric analysis at various scales, forecasting, air quality models, and many other applications in the study of the atmosphere at various scales. Meteorological data are also useful for hydrological and agricultural studies, and research in meteorological and climatological processes [2,3]. The main source of this data comes from surface meteorological stations, which provide information from several calibrated instruments, so it is reliable and serves as a reference to the development of new methodologies [4]. Worldwide, the number of meteorological stations has increased because of the need to provide constant information in real-time and the development of new data measurement systems [3]. Contrarily, in Ecuador since the 1990s, the number of ground stations has decreased mainly due to economic limitations and budget cuts to public institutions, despite this, the threats of climate extremes under a changing climate.

There is a global trend in the implementation of automated meteorological stations, and the data obtained is widely used. However, there are still several limitations in the use of automatic stations especially in developing countries, for instance, power supply to the equipment, communication with the data collection center that relies on network facilities, the high cost of operation, maintenance, and data quality assurance through regular calibration, among others [5]. Additionally, the difficulty of data collection in inhospitable locations, limits operational capabilities, and the need for extensive data coverage have accelerated the search for alternative measurement methods [6]. Radiosondes data are another valuable source of information as much as ground weather station data. This technique allows obtaining data of variables in the vertical profile up to 30–35 km, with good vertical resolution and high accuracy, but its high costs have limited its use in countries with low economic resources [7,8]. The scarcity of vertical information, particularly in places with complex topography, is a limiting factor for the development of meteorological research [9]. The existing limitations using radiosondes and meteorological stations have led to the pursuit of new methodologies for the determination of vertical profiles of atmospheric variables, including satellite data, modeling, ground-based remote sensing methods, or the use of drones for low vertical measurements [10,11,12,13,14].

Numerous studies have attempted to test alternative methodologies to the use of radiosondes and weather stations for vertical profiling of the atmosphere. For instance, Adamo et al., (2007) tested the effectiveness of MODIS vertical satellite data for two stations in Europe by comparing it with radiosonde data. In that study, it was determined that on average MODIS data have acceptable representation in relation to radiosondes, but it fails to represent the sharp variations in the profile [15]. Pérez-Planells et al., (2015) also demonstrated that the vertical resolution of MODIS data is not enough and there are data sources, such as the National Centers for Environmental Prediction (NCEP), reanalysis data that present better results for atmospheric variables. MODIS was also used for the calculation of the height of the Planetary Boundary Layer (hPBL) [16]. Feng et al., (2015) developed a method to derive the hPBL by the minimum mixing ratio gradient method, with a data filling technique to avoid gaps due to cloud cover, having an average RMSE of 370 m [17]. Onyango et al., (2020) also tested MODIS data to determine the hPBL, but it was concluded that it is not feasible because of the large differences in the values obtained due to the low vertical resolution of the profiles [18].

The Atmospheric Infrared Sounder (AIRS) is another satellite product that has been tested in the study of the atmospheric vertical profiling. In fact, in the study conducted by Feng et al., (2021) in northwest China, the temperature gradient method was used with AIRS data to derive the hPBL [19]. The results showed that for the structure of the Potential Temperature (PT) AIRS profile can successfully reflect the stability of the PBL, thus improving its detection. AIRS is also applicable during cloudy days for both convective and stable air conditions, demonstrating that the PT method is more reliable than the Mixing Ratio method with MODIS data, used by Feng et al., (2015), although the error caused by the vertical resolution of AIRS data is unavoidable [19]. Ding et al., (2021) compared the HPBL of the AIRS products, the global positioning system radio occultation (GPS RO) experiment, and the NASA Modern-Era Retrospective analysis for Research and Applications-2 (MERRA-2), from June 2006 to December 2015 across the globe. The improved water vapor retrieval in AIRS version 7 allowed the PBLH to better match that of GPS RO and MERRA-2, especially near the equator and at low latitudes [20]. Finally, Martins et al., (2010) used data from AIRS and the Rain in Cumulus over the Ocean (RICO) campaign to verify the accuracy and precision of the AIRS product. There was good agreement between AIRS and RICO data in a shallow oceanic cumulus regime that is known to be difficult to analyze with other remotely sensed data [21].

Vertical profiling has also been tested through the Weather Research and Forecasting (WRF) model. Caneo et al., (2011) tested the effectiveness of five configurations of the WRF model for vertical profiling with positive results in Mexico [22]. The WRF model has also been tested for simulating the vertical structure and diurnal cycle of the PBL with various configurations in Algeria [23], and in air quality studies in the Ecuadorian cities of Quito [24], and Cuenca [25]. The study conducted by Parra, (2017) tested five options of hPBL parametrizations in the WRF model in Quito, demonstrating that the parametrizations that best represent the conditions in the air quality models are the Quasi-Normal Scale Elimination (QNSE) and Mellor-Yamada-Janjic (MYJ) parametrizations [24]. Finally, technological advances have allowed different methods for atmospheric profiling to be tested such as lidar combined with artificial intelligence methods [26,27] or validation with drones [11]. The use of drones has become a research focus because they could replace radiosondes in studies of the lower troposphere and have the advantage that the sonde can be recovered after launch [10]. There are several studies worldwide that validate the use of drones for the determination of atmospheric variables at different altitudes [28,29,30], these studies explore the use of drones for vertical profiling, and they show that the results have low biases, but there are several conditions to be taken into account, such as the ascending rate, the influence of the drone components in the measures, and the flight time.

The search for alternative methods for vertical profiling of the atmosphere has received great attention for several years [30,31,32,33,34]. However, to the authors’ knowledge, no study has been conducted to test various methods of vertical profiling in Ecuador. Therefore, the aim of this study is to evaluate, for a campaign of 10 radio soundings, the representation of T, TD, Q, and PT with drone, MODIS, and AIRS satellite data, and WRF simulations from 0 magl. to 300 magl. to test the effectiveness of the methods in the first few meters from the surface, and to validate the accuracy of the drone data up to its maximum height (300 magl.). The vertical analysis was also developed from 0 magl. To tropopause height (15,000 magl.) for WRF and satellite data. Moreover, the capability to represent the HPBL is assessed. Thus, the article presents the following sections: Section 2 shows the study area describing the climatic characteristics of the zone. Section 3 details the data obtained for the profile comparison and summarizes the data preprocessing for the study and the selected profile comparison methods and metrics, in this case, vertical profiles, Taylor diagrams, and three evaluation metrics were selected. This section also describes the different methods used to determine the hPBL with the radiosonde data. Section 4 shows the results obtained in the study and their discussion, and Section 5 shows the main conclusions of the paper and future work in this area.

2. Study Area

Ecuador is located on the equatorial line, in the northwest of South America. The country has a strong geographic and topographic contrasts between the highlands, which correspond to the Andean Mountain range in the north–south direction, the coastal plains in the west, and the Amazon jungle in the east. The topographic gradients are sharp, making it possible to ascend from sea level to peaks above 5000 m above sea level (masl) in less than 300 km [35]. The climate in Ecuador shows high spatiotemporal variability. For instance, the Andes Mountain range acts as a weather divide, which also regulates precipitation regimes and seasonality [36]. In addition to the geographic influences, the climate variability is high due to the effects of the Intertropical Convergence Zone (ITCZ), El Niño Southern Oscillation (ENSO), and other low frequency phenomena, such as the Pacific Decadal Oscillation (PDO) [37]. Over the inter-Andean valleys, there is a bimodal regime related to the ITCZ, causing rainy periods from March–April and October–November [36,38].

The Metropolitan District of Quito (DMQ) is located in the Ecuadorian highlands between 1900 and 4800 masl., with valleys, elevations, and plains that affect wind circulation, generating topo-climates, e.g., accentuations or attenuations of the regional climate due to the relief [39]. The urban area of the DMQ has a mean altitude of 2850 masl. Currently, the DMQ is the most populated region of Ecuador with 3 million inhabitants approximately [40]. The observations of weather parameters in Quito are determined through governmental and municipal meteorological and air quality stations to obtain a high coverage of data due to the high irregularity of the area. The National Institute of Meteorology and Hydrology of Ecuador (INAMHI) operates three stations in the city, across a 50 km distance in the north–south direction [40]. Figure 1 shows the urban and the complete area of the DMQ and the location of the Izobamba (IZO) weather station, where the study is conducted.

IZ station is located in a district adjacent to Quito. However, due to the data availability and its proximity to the southern part of the capital of Ecuador, this station was selected for the present study. IZ has an altitude of 3048 m above sea level (masl.), and its minimum and maximum annual temperatures are 10 °C and 16 °C, respectively. Its location is classified as an urban area ecosystem and its annual precipitation is 1500 mm approximately [41]. IZ station was selected for this study because INAMHI has in place the necessary equipment to conduct radiosonde measurements.

3. Data and Methods

3.1. Radiosonde and Drone Data

Data from different sources were used to determine the effectiveness of alternative methods used to measure the vertical profile of different atmospheric variables in IZO station. The baseline information in this study was obtained from a sounding campaign of 10 radiosondes conducted between September 2021 and January 2022, which launching dates and time are shown in

Table 1. Weather Research and Forecasting (WRF) model configuration.

No	Radiosonde Launching Date and Local Time	Drone Launching Date and Local Time
1	6 August 2021—11:50	6 August 2021—10:56
2	17 September 2021—12:48	17 September 2021—11:46
3	22 October 2021—10:53	22 October 2021—10:10
4	29 October 2021—10:58	29 October 2021—10:16
5	12 November 2021—10:48	12 November 2021—09:50
6	19 November 2021—10:54	19 November 2021—10:18
7	26 November 2021—10:20	26 November 2021—09:36
8	10 December 2021—10:45	10 December 2021—09:38
9	17 December 2021—11:11	17 December 2021—10:09
10	21 January 2022—11:14	21 January 2022—10:06

. The radiosondes were conducted with Vaisala RS92-SGP sensors. The ground check of the radiosonde was performed with the Vaisala GC25 device and the sounding system, which automatically reads the calibration coefficients via a telemetry connection and removes any chemical contaminants to ensure the accuracy of humidity measurements [42]. The sounding system was set up with the surface variables data of the IZO station for all radiosondes when launched.

The drone campaign was conducted along with measurements captured by the same sensors used in the radiosonde but tethered to an ATyges FV8 drone instead of the weather balloon, as shown in Figure 2. The drone was launched before the radiosondes to use the same sensor and ground station; the launch time is shown in Table 1. The drone measurements were taken up to 300 magl., due to flight restrictions from the Ecuadorian Civil Aviation Authority, DAC. The advantage of using the drone is that it can carry the same sonde, so the format of the data obtained is the same as that of radiosondes; it is also an economical technique due to the potential for measurement in remote areas and it is easily conducted. In addition, its use allows the retrieval and reuse of sondes, which are generally lost with the balloon in radiosondes, thus reducing monitoring costs. Despite this, its use is limited mainly because of the flight altitude it can reach, and the interference that could exist in the measurements due to the operation of the drone, the effect of the propellers, or the temperature of the batteries [10].

The Vaisala sensor used for both weather balloon and drone measurements, captures information vertically at a preset time interval. For instance, the radiosonde system used in the project is configured to return data for all target variables every 20 s during flight. Drone and radiosondes were anchored to the same sensor, so the measures did not occur at the same time, the difference between the launches varied between 36 min and 68 min, so it is assumed that the weather conditions during this time do not vary too much. Another limitation of this type of data collection is that the ascending speed of the weather balloon and the drone is variable due to the technical conditions of the sonde and the atmospheric conditions in the place of the measurement. It is important to point out that the drone measurements needed to be taken with a special flight configuration to capture data correctly, the drone had an ascending rate of 3 m/s and stopped for two minutes every 100 m, so the ground station recognized that the sonde was launched. Due to the sonde configuration, the data was only taken while the drone was ascending, and not in the descending process.

The connection between the sensor and the surface station may present intermittencies, which could affect the data capture in the vertical profile. As the objective of this study is to make a comparison of vertical profiles, it is preferable that the information of the profiles have the same vertical resolution to be able to make the comparison. For this reason, the programming language R, version 4.1.2. was used for the interpolation of the measurements at uniform heights for all the soundings in all variables. As the information captured by the sonde is of high resolution, it was decided to linearly interpolate the profile values with a vertical resolution of 15 m due to the maximum height that the drone reaches, in this case, up to 300 magl. To perform the interpolation, the R based function called “approxExtrap” was used. Additionally, to determine the hPBL through three different gradient methods, which will be described in Section 3.5. below, a 40 m vertical resolution was applied to radiosonde data to compare it to the satellite and WRF data up to 15,000 magl. Drone data were not considered for the hPBL determination due to height limitation. For the drone data, a quadratic fit was tested to remove the noise of the data, the fit was performed in R and it was compared to the raw data to see the level of improvement in its representativeness up to 300 magl.

3.2. Satellite Data

Teledetection has become an essential tool for the analysis and evaluation of natural hazards and the effect of anthropogenic activities on the environment. The availability of global images and data at near real-time helps to warn of possible consequences of natural hazards [12]. The use of remote sensing not only reduces the time invested in research work but also offers other advantages, such as total coverage of the earth’s surface, panoramic vision, and homogeneity in data acquisition, among others [43]. The National Aeronautics and Space Administration (NASA) has been a pioneer in space research. One of NASA’s most important Earth Observation systems (EOS) includes the Moderate Resolution Imaging Spectroradiometer (MODIS) sensors aboard the TERRA and AQUA satellites. In addition to multispectral images, MODIS has 44 standard data products that have improved the understanding of the dynamics and processes taking place on the surface of the earth, the oceans, and the atmosphere [44]. In this study, the MOD07_L2 product of the TERRA satellite was evaluated, this product provides atmospheric measurements in the vertical profile. MODIS products have a grid resolution of 1 km, and a vertical resolution of 20 vertical pressure levels [17,44].

The satellite data contains gaps of information. While some gaps are internal to the satellite, mainly related to its orbits or to instrument failures, others are external, such as cloud coverage [17]. In the DMQ, cloud coverage is the main factor that leads to spatial gaps; in fact, the study area has approximately 8 months with high cloud coverage from October to May, being March the month with more cloud coverage of ca. 90% [45]. It is important to highlight that the satellite TERRA provides one measurement per day, and in Ecuador the data is taken at approximately 15:40 UTC, 10:40 local time. To provide a complete spatial coverage, it was decided to apply a gap-filling methodology adapted from the one proposed by Feng, (2015), which in this case, is only performed spatially, not temporally. The spatial component uses a 5 × 5 pixels area, which is a “quality control” parameter included in MODIS MOD07_L2 product data. If at least three data pixels are available in the area surrounding the point of interest, the inverse distance-weighted mean method (IDW) is used to fill the empty pixel. In this study, the pixel available data was set to three because of the limitations related to cloud coverage during the months of the study. A bigger area and the temporally gap-filling techniques were not considered due to the high spatiotemporal variability of the study area, which would have highly increased the uncertainty of this study [17,41].

The Atmospheric Infrared Sounder (AIRS) Version 7, Level 2 product was also used in this study. AIRS, as MODIS, is one of several satellites onboard the EOS Aqua spacecraft, being the most accurate and stable set of hyperspectral infrared radiance spectra measurements in space to date [19]. AIRS was designed to measure the Earth’s atmospheric water vapor and temperature profiles on a global scale, with a grid resolution of 50 km, and a vertical resolution of 28 pressure levels between 1100 and 0.1 hPa. The download format is HDF-EOS [19,46] and it was converted to GeoTIFF using the HEG conversion tool developed by NASA [47]. The download was performed from NASA’s “Goddard Earth Sciences Data and Information Services Center” (GES DISC) page [48]). The downloaded variables were geopotential height, atmospheric temperature, relative humidity, and water vapor mixing ratio. For quality control (QC) each atmospheric variable and pressure level has a field with three options: 0 indicates best quality, 1 indicates good quality, and 2 indicates no quality assurance [46].

Both MODIS and AIRS have information gaps that were not filled due to high cloud coverage, internal satellite failures, or low grid resolution. For this study, MODIS vertical profiles were available on 7 of the 10 dates of the study period, and for AIRS vertical profiles 6 of 10 were available and usable.

3.3. WRF Data, Parametrizations, and Processing

The Weather Research and Forecasting model, WRF, is a meteorological model widely used in both research and numerical weather forecasting. The WRF model solves the primitive equations, is non-hydrostatic, and is fully compressible. This model has been extensively used in research at universities and government laboratories, for operational forecasting by governments and private entities, and for commercial applications by industry [49,50]. Its broad use is due to its free access, and the permanent maintenance and updating to which it is subjected [51,52]. Studies that have been performed using the WRF model in Ecuador are still scarce but varied with respect to the use of different applications [49,53,54,55]. One of the advantages of the WRF model in mountainous areas such as the Andes is the use of sigma coordinates that allows increasing the vertical resolution of the simulation [49,54] and that was used in this study for the determination of vertical profiles.

A real-time WRF (Version 4.2) modeling was developed by performing an operational dynamical downscaling using three inner domains of four. Thus, the mesoscale synoptic characteristics are gradually reduced to local resolution, which is necessary for accurate forecasting required for local and low tropospheric studies [56].

Table 2. Weather Research and Forecasting (WRF) model configuration.

Configuration/Domain	27 km	9 km	3 km	1 km
Time Interval (min)	180	60	60	30
Model Data	Type: GRIB2 data Resolution: 1deg global data Output frequency: 6, hourly 27 pressure levels (1000–10 hPa)
Grid points	80 × 80	82 × 82
Vertical levels	60
Nesting	No	Yes
Microphysics	2—Lin et al. scheme
Radiation (longwave)	1—RRTM scheme
Radiation (shortwave)	2—Goddard Shortwave scheme
Surface layer	1—Monin-Obukhov Similarity scheme
Land surface	1—Thermal Diffusion scheme
Planetary Boundary Layer (PBL)	1—YSU scheme
Cumulus	10—KF-CuP scheme

summarizes the model configuration and parametrizations used in this study, which come from the published and unpublished work of the authors of this study [24]. The model was run using a two-way nesting with two input files configuration [57,58]. Three nested domains were set up with grid resolutions of 9 km, 3 km, and 1 km, while the outer 27 km domain was initialized with data from the NCEP Final (FNL) Operational Model Global Tropospheric Analyses (NCEP-FNL) with a resolution of 1° × 1° [59]. The boundary conditions of the outer domain are updated every 6 h from NCEP-FNL data; the outer domain provides the boundary conditions to the inner domains. In this study, a spin-up time of two days was evaluated for each of the modeled dates. The results of the higher resolution inner domain were generated every half hour to select the data closest to the time of the radiosonde launch to reduce possible variations in the data. Finally, the resulting data (wrfout) was processed with the library WRF-Python of the programming language Python version 2.7.5.

3.4. Results Evaluation

All new measurement methods must undergo a validation process to ensure that the information has been generated appropriately, and to allow acting in case of erroneous detection [3]. Validation is performed by comparing the measured data with the best quality available data from meteorological stations and radiosondes [2]. Validation and evaluation of results is an essential component of any scientific research, this is mainly because it allows showing, in a clear and simplified way, the performance of the techniques evaluated in scientific research and development studies [60,61]. There are several evaluation statistical metrics, such as the Root Mean Square Error (RMSE), and some graphical methods, such as Taylor Diagrams, which allows to compare the performance of various methods in a single plot [6,62,63,64,65]. The major challenge in selecting metrics to measure the performance of the methods studied for the profiling of climate variables is to determine which variables or phenomena are important to accurately assess and, therefore, which metrics to measure [63].

The performance of the variable measurement methods analyzed in this study was evaluated using three statistical metrics, two of error, and one of correlation, and a graphical method that relates different statistical indices. As error metrics the Root Mean Square Error (RMSE) (Equation (1)) and the Mean Absolute Error (MAE) (Equation (2)), and for correlation, Kendall’s coefficient was selected. The RMSE is one of the best-known and widely used metrics to determine the difference that exists between real and simulated data [66]. On the other hand, the MAE measures the average magnitude of errors such that all individual differences have equal weight [67]. There has been suggested that RMSE might be a misleading indicator of average error [68], and therefore the MAE would be a better error metric. However, Chai and Draxler, (2014) demonstrated that a combination of RMSE and MAE, among other metrics, are often required to successfully evaluate a model performance. For the two metrics, the closer the calculated value is to zero, the better the model represents the real values.

$$RMSE=\sqrt{\frac{1}{n}\sum _{i=1}^n{\left(\widehat{y_i}-y_i\right)}^2}$$

(1)

$$MAE=\frac{1}{n}\sum _{i=1}^n\left|\widehat{y_i}-y_i\right|$$

(2)

where: $\widehat{y_1}$, $\widehat{y_2}$, …, $\widehat{y_n}$ are the forecasted values.

$y_1$, $y_2$, …, $y_n$ are the observed values.

$n$ is the number of observations.

Kendall’s rank correlation coefficient or Kendall’s Tau is a non-parametric test used primarily in the analysis of climate time series data [69,70]. Kendall’s coefficient was used because the conditions to use Pearson, such as both variables being normally distributed, or linear relationship between both variables may not be meet by our variables. This coefficient returns values between −1 and +1, a positive correlation indicates a similar trend in the ranges of the analyzed data, while a negative correlation denotes that just as one data series increase, the other decreases [71].

Taylor diagrams provide a visual statistical representation that allows comparing a set of variables from one or more estimated datasets with benchmark observations of the same variables [64,72,73]. The relationship between the compared data is quantified in terms of the three most representative statistical values of a model/method: (1) Root Mean Square Error (RMSE), which represents the accuracy of the model; (2) Standard Deviation, which shows the variability between the compared data; and (3) Pearson Correlation Coefficient (R), indicates the linear relationship existing between the compared data [74,75]. These diagrams are commonly used with values obtained from atmospheric models, this is due to the large number of variables derived from these models, such as precipitation, temperature, and humidity, among others [72,73]. The potential of this graphical evaluation method motivated its use in this work for the comparison of the profiles, taking radiosondes as a reference, and determining the effectiveness of MODIS and WRF measurements up to 15 km above surface level and up to 300 m in drone measurements.

3.5. hPBL Calculation Methods

In addition to the profiles of the meteorological variables detailed above, we also considered determining the height of the Planetary Boundary Layer, because of its importance in the lower troposphere studies. The hPBL is a turbulent layer that extends from the surface to a variable height between 200 m and 4 km. The formation of the PBL is due to the daily cycle of radiative heating and cooling of the Earth’s surface, mainly under favorable meteorological conditions (clear skies), reaching minimum values at sunrise and maximum values around midday. This layer is characterized by presenting a cycle of variations in temperature, humidity, wind, and pollutants concentration in response to the variation of the surface cycles that produce it [76].

The study of the PBL is fundamental for atmospheric sciences due to the great influence it has on the physics of the atmosphere and the climate system in general [77]. There are two basic possibilities for the practical determination of the hPBL, its derivation from profile data (measurements, or simple calculations) and its parametrization by equations or models [7]. In this study, both approaches were compared. There are several methods to determine the hPBL from profile data, and every method has its own potential and some considerations to apply it, the selection of the used methods was made based on the available variables data and the characteristics of the study area. For the radiosonde and satellite data, the vertical gradients of potential temperature, relative humidity, and mixing ratio were used [7,31,78,79]. The results of these hPBL methods were also compared with the hPBL values extracted directly from the WRF model output files [56].

4. Results and Discussion

4.1. Lower Troposphere Analysis

The evaluation of the different sources of vertical measurements was conducted through various methods as was settled in the methodology section. Vertical profiling was the first method used to evaluate whether the measured atmospheric variables had a good graphical representation of the trend and variation when compared to radiosonde data. Figure 2 shows the temperature profiles, for 10 measurement events, up to 300 magl. The temperature profiles (Figure 3) show that the measurements of temperature have variations of no more than 4 °C, being the drone measurements the closer ones (less than ~2 °C). It can be seen recurrently that the MODIS and AIRS measurements (sky-blue and yellow), both for their vertical resolution and their surface value, do not match well with the reference values. However, AIRS shows better accuracy than MODIS. It can be observed that for the profiles of some dates there is no satellite data due to high cloudiness in MODIS measurements, which in the DMQ tends to have high percentages, especially in the study period of this research, and for AIRS data there are gaps of information that also limited its analysis. The profiles extracted from WRF model simulation and captured by the drone show better results, particularly those from the drone, although they present a noisy signal about 1 °C for three of the measurements. The noise of drone measurements was corrected through a quadratic adjustment (shown in purple) which helped the data to be the closest to the reference (black line).

The temperature (T) profiles measured by the drone and modeled with the WRF, in general, provided a good representation when compared to the reference values. However, when evaluating Figure 4, which contains the dewpoint temperature (TD) profiles, larger errors are observed. This may be related to the TD dependence on humidity, which is more difficult to estimate. Drone measurements are still the closest to the reference in terms of magnitude, but a strong noisy signal is shown, especially in the first few meters. Interpolation was not applied because the quadratic fit did not show good results for the humidity dependent data. This causes lower values of correlation than in the temperature profiles. In this case, MODIS presents much larger differences, this is mainly attributed to the vertical resolution that limits the analysis on the 300 m scale. AIRS, despite its low vertical resolution, has a better performance than MODIS in most of the profiles, with errors relatively low, and a better performance than WRF in some profiles. WRF has a moderately good representation, but it does not capture well the variations that occur mainly in the first few meters. For the rest of the variables: mixing ratio (Q) and potential temperature (PT), the profiles can be found in Appendix A.1 and Appendix A.2, respectively. Q and PT profiles show a similar result than TD with WRF and drone data closer to the reference.

Despite vertical profiles give a graphical overview of the performance of the different methods, they do not offer a quantitative measure of accuracy. Thus, Taylor plots were developed to evaluate the performance of each method, which can be observed based on the distance between the reference point and each plotted point. Figure 5 shows the normalized Taylor Plots for T, TD, Q, and PT up to 300 magl. for 10 measurement events, analyzing four different profiling methods (Drone, WRF, MODIS, and AIRS) with radiosonde measurements as a reference. In the Taylor Plots shown in Figure 5, it is notorious that the variable with the best representation in the different methods is the temperature (a), with correlations close to or higher than 0.9 and RMSE values lower than 1. In terms of standard deviation, it can be observed that MODIS (squares) and AIRS (stars) measurements have the highest standard deviation so it can be determined that up to 300 magl., they are the least accurate methods for measuring temperature in the vertical profile. The drone quadratic adjustment does improve some profiles, but the differences are hardly noticeable.

Regarding the other three variables analyzed in Figure 5, it can be observed that the representation is weak up to 300 magl. In general, correlation values are lower than 0.9, and even in some cases they have negative correlations, as in the MODIS (squares) and AIRS (stars) measurements of PT, and drone measurements (triangles) in the TD. In the TD Taylor plot (b), the four profile measurement methods have different results, but in terms of correlation, the drone profiles have results farthest from the reference due to the interferences in the flight. In the Q plot (c), several drone measurements have low correlation, while the WRF presents a normalized standard deviation close to zero in most of the events analyzed. Finally, regarding the PT, the MODIS and AIRS data have large distances from the reference values, and the drone maintains a low average correlation. The low correlation of drone data is related to the noisy signal (see Figure 3 and Figure 4) although the error is lower than WRF and satellite data. WRF has an overall good representation of the PT.

Although Taylor plots show a clearer idea of the measurements than the profiles presented before, it is difficult to determine which method is more effective in the graph if there are several points or if they are close to each other. For this reason, the analysis of the Taylor plots was complemented with three statistical indices shown in

Table 3. Average values of the statistical metrics of drone, MODIS, and WRF profiles for the four variables analyzed up to 300 magl.

Drone			MODIS			WRF			AIRS
RMSE	MAE	KENDALL	RMSE	MAE	KENDALL	RMSE	MAE	KENDALL	RMSE	MAE	KENDALL
Temperature (°C)
0.79 0.75	0.70 0.67	0.92 0.98	2.50	2.38	0.98	1.29	1.24	0.98	3.58	3.50	0.97
Dewpoint Temperature (°C)
0.83	0.65	0.48	6.06	5.99	0.88	1.96	1.85	0.88	2.58	2.51	0.89
Mixing Ratio (g/kg)
0.53	0.41	0.33	3.72	3.67	0.80	1.21	1.13	0.79	1.75	1.71	0.83
Potential Temperature (°C)
0.92	0.84	0.37	3.44	3.29	−0.22	1.47	1.40	0.30	4.66	4.54	−0.20

. The results show that in temperature measurements, even though the correlation values are higher in the MODIS, AIRS, and WRF profiles, the error metrics are lower in the drone measurements, so the magnitude of the variables is closer to the reference. Below the values of the metrics for the raw drone data, the values of the statistical indices for the quadratically adjusted drone data are shown in bold, these values indicate that there are improvements, especially in correlation, when the adjustment of the data is applied. The same pattern is repeated for all variables, WRF generates better values in terms of correlation, but the error metrics indicate that the magnitude of the drone values is the most accurate with all values below one in all variables. AIRS data show in general a better representation than MODIS, except in temperature, where the metrics show big differences compared to the reference data.

4.2. Upper Troposphere Analysis

Section 4.1 detailed the findings obtained by analyzing the different profiles up to 300 magl., which was done mainly to evaluate the performance of the measurements in the lower part of the troposphere and to validate the drone data. However, lower atmosphere studies require analysis at higher altitudes, so that the hPBL and the atmospheric behavior up to the tropopause can be determined. Therefore, in this section, the results of the profiles up to 15 km determined with radiosonde (reference), the MOD07_L2 product of the MODIS sensor, AIRS satellite data, and WRF are presented. Figure 6 shows the vertical profiles of the temperature up to 15 km above ground level for 10 measurement events. In contrast to the profiles observed up to 300 magl., it can be observed that the temperature representation up to 15,000 magl. is more accurate in terms of vertical behavior and magnitude. It can be clearly observed that in most of the analyzed dates, MODIS profiles are the ones that present more differences in relation to the reference profile (black), while AIRS and WRF follow a close path to the reference data. In addition, as detailed in the previous section, there are three dates on which MODIS data could not be obtained due to cloud coverage, and four data gaps in AIRS data.

Figure 7 shows the TD profiles for the 10 dates analyzed in this article. The profiles show a strong vertical variability. It is observed that MODIS, AIRS, and WRF are not able to capture the largest variations as in the case of 10 December 2021, or 21 January 2022. However, WRF is the best adapted to these variations because there are certain similarities in relation to the reference data, although it does not manage to represent very well the strongest changes in magnitude. MODIS, due to its low vertical resolution, underrepresents the vertical variability of the radio sounding measurements. AIRS profiles show a better performance than those of MODIS, however, the sharp variations of dewpoint temperature in the vertical are not well represented because of the low vertical resolution of the data. The profiles for Q and PT variables are presented in Appendix A.3 and Appendix A.4, respectively. They show that WRF and AIRS are the best to represent the variables in the atmosphere, MODIS show stronger variations, being the least accurate for all variables.

A statistical comparison of the effectiveness of the two profiling methods is presented in Figure 8 with Taylor plots for all the variables under study. In general, it is observed that the representation of all the variables up to 15 km above the surface is highly accurate since most of the points are close to the reference point, which indicates that the profiles obtained from MODIS, AIRS, and WRF are effective overall. As in the analysis up to 300 magl., the most accurate variable is T (a), in this case with a correlation higher than 0.95 and standard deviations lower than ±0.75. Regarding TD (b), the measurement effectiveness is also high, and it has low variability, the correlations for this variable are higher than 0.9, which indicates that the profiles are also very accurate. For both variables (T and TD), it is not possible to determine exactly if MODIS, AIRS, or WRF profiles have better results, because in general, all of them are equally good. When analyzing Q and PT (c and d), there are also high levels of accuracy with correlations higher than 0.85 and 0.9, respectively, but it is observed that the MODIS profiles (squares) generally show higher standard deviation, so it is inferred that the AIRS and WRF profiles have greater accuracy when compared to the radiosonde data.

As detailed above, Taylor plots determine the effectiveness of the profiling methods. However, there are cases, such as the one in Figure 8a, where all values are close to the reference point and the efficiency cannot be accurately quantified.

Table 4. Average values of the statistical metrics of drone, MODIS, and WRF profiles for the four variables analyzed up to 15,000 magl.

MODIS			WRF			AIRS
RMSE	MAE	KENDALL	RMSE	MAE	KENDALL	RMSE	MAE	KENDALL
Temperature (°C)
4.07	3.48	0.98	1.62	1.20	0.98	2.39	1.51	0.98
Dewpoint Temperature (°C)
7.87	6.63	0.94	8.36	6.15	0.91	7.20	5.20	0.93
Mixing Ratio (g/kg)
1.48	0.81	0.89	2.98	1.53	0.88	0.77	0.41	0.89
Potential Temperature (°C)
25.37	22.45	0.94	2.66	1.94	0.99	4.83	2.80	0.99

shows the values of the evaluation metrics selected for this study for each variable up to 15,000 magl. It can be determined that the variables that are best represented in the profiles are T and Q. The profiles modeled and extracted from the WRF are generally accurate and show high correlation values and relatively low errors, except in the TD profiles, where the errors are significant, this is due to the large variations in the observed profiles in Figure 7. The MODIS profiles have a poor representation, especially in PT the errors that reach values higher than 20 °C, so it is evident that the profiles are not accurate to be used as an alternative instead of radiosonde data in the study area. The AIRS profiles have a good representation of T and Q, showing a better performance than MODIS profiles, and are comparable to WRF values.

4.3. hPBL Analysis

The hPBL is a fundamental parameter for air quality and meteorology studies in the lower troposphere; therefore, its estimation was considered as part of the variables compared in this study. Figure 9 shows the methods of potential temperature (PT), relative humidity (RH), and mixing ratio (Q) gradients for the calculation of the hPBL using data of 10 globe radio soundings. It is generally observed that the values are consistent between them, with some minor differences. All values are in a range between 3603 and 6858 masl., considering that the height of the IZO station is 3048 masl, the values are up to 3810 magl. The radiosondes were performed in a time range between 10:20 am and 12:48 pm, local time (GMT-5).

Table 5. Results synthesis of the hPBL calculation at the IZO station.

IZOBAMBA STATION (3048 masl)
No	Date and Local Time	CAPE [J/kg]	PT Gradient	RH Gradient	Q Gradient	hPBL
1	6 August 2021—11:50	87	2055	2430	2430	2430
2	17 September 2021—12:48	47	2025	2025	2025	2025
3	22 October 2021—10:53	740	3810	3735	3735	3735
4	29 October 2021—10:58	380	1665	1590	1590	1590
5	12 November 2021—10:48	592	2775	2775	2775	2775
6	19 November 2021—10:54	232	2760	2760	2895	2760
7	26 November 2021—10:20	33	2580	2160	1665	2160
8	10 December 2021—10:45	323	555	555	555	555
9	17 December 2021—11:11	386	3645	3315	3315	3315
10	21 January 2022—11:14	44	2145	2145	2145	2145

summarized the hPBL values for all dates with three methods tested and presents the reference hPBL value to be taken for the evaluation of the tested methods. Table 5 also contains the value of the Convective Available Potential Energy (CAPE), which is a measure of the amount of energy available for convection in the analyzed place, this parameter indicates the degree of instability of the atmosphere.

The hPBL values shown in Table 5 are primarily congruent between methods. The most recurrent value was selected as the reference value, or in the case of not having similar values, the average value among the results was taken. Among the hPBL values shown in Table 5, particularly stand out 3735 magl. occurred on 22 October 2021, which may be explained by a very low dew point depression and high humidity, unstable conditions that may help to develop such a high value of PBL (Figure 10). Another high hPBL value correspond to 3315 magl. occurred 17 December 2021, in this case the dew point depression is also very low producing high humidity and consequently high atmospheric instability. Although there may be computational uncertainties in the calculation of these values, the use of three methods for calculating the hPBL added to the evaluation of the vertical profiles gives a more robust estimation.

To verify that these results are correct, validation was performed with the temperature and dewpoint temperature profiles as a function of the pressure values, and they were compared with the profiles obtained from WRF, AIRS and MODIS, as shown in Figure 10, Figure 11 and Figure 12, respectively. It is observed that in most cases, the hPBL height can be obtained from the plots. In cases where an adequate representation is not achieved, a recalculation of the value was performed with reference to the radiosonde profiles (6 August 2021 and 26 November 2021 profiles). Comparing these results with the profiles extracted from the WRF (Figure 10), it can be noted that the representation is good and mostly match with the variations of the reference values. On the other hand, the values from AIRS (Figure 11) have some cases where the profiles are close to the reference, but the vertical resolution does not allow a proper hPBL determination. Additionally, there are just 6 of the 10 dates with data, and the profiles from 12 October are completely out of line with the reference values. The values extracted from MODIS (Figure 12) also fail to capture the variations and have inconsistencies in terms of magnitude and variability of profiles, so the use of AIRS and MODIS values for the calculation of the hPBL was discarded and only the comparison of the corrected radiosonde values and the values directly extracted from the WRF was performed.

The validated hPBL values, corrected when needed, were compared with the hPBL values extracted directly from the WRF modeling, as it is a resulting parameter (M1), and the hPBL values calculated with the same methods as the ones used for radiosondes (M2), which results are shown in Appendix A.5.

Table 6. Comparison of hPBL values at the IZ station between Radiosonde and WRF.

	IZOBAMBA STATION (3048 masl)
N°	Date and Local Time	Radiosonde hPBL	WRF hPBL M1	WRF hPBL M2	Difference
1	6 August 2021—11:50	1960	1144	2360	−816/+400
2	17 September 2021—12:48	2025	1831	3200	−194/+1175
3	22 October 2021—10:53	3735	1581	1800	−2164/−1935
4	29 October 2021—10:58	1590	1940	2080	+350/+490
5	12 November 2021—10:48	2775	1498	1520	−1277/−1255
6	19 November 2021—10:54	2760	1027	2080	−1733/−680
7	26 November 2021—10:20	1665	1029	1240	−636/−425
8	10 December 2021—10:45	555	1184	2640	+629/+2085
9	17 December 2021—11:11	3315	647	3480	−2668/+165
10	21 January 2022—11:14	2145	1331	2640	−814/+495
Mean difference					−931/+52
RMSE					1367/1112

shows the values compared for each date and the existing difference. To obtain the WRF hPBL value, the half-hour closest to the time of the sounding was taken, since the modeled data has a temporal resolution of half an hour. In contrast to the results obtained in the profiles analyzed in Section 4.1 and Section 4.2, the WRF hPBL value extracted from the wrfout file is not accurate with respect to the reference values. In general, it is observed that the WRF underestimates the hPBL values with differences of up to 2668 m below the reference value. The mean difference between the modeled data and the reference values is –932 m, and the RMSE is up to 1367. With the second method (M2), which used the same hPBL determination criteria, the results are closer to the reference, however, they still have a RMSE up to 1112, and differences up to 2085 m. The wide differences found reveal that possibly the PBL parametrization used in this research (YSU scheme) might not be appropriate for the study area. The YSU scheme has a nonlocal mixing approach which was analyzed by Parra (2017) in the DMQ. This scheme produced less favorable results than other schemes early in the morning due to its nonlocal characteristics; however, it was selected for this study because it considers different layers with the overlap of large and small eddies. Thus, a complete sensitivity analysis should be performed to determine the scheme that best represents the hPBL in the DMQ.

5. Concluding Remarks

The characterization of the vertical atmospheric behavior is of great importance in several areas of knowledge, the main ones being air quality and meteorology studies. Despite technological advances, the most efficient technique for vertical data collection is still radiosondes, which use high-precision sensors, a surface data station, and a weather balloon that ascends with the sonde up to ~35 km. The great limitation of these measurements, especially in developing countries, is their high cost, which makes it difficult to conduct these measurements frequently, although they provide valuable data to meteorological centers. For this reason, the search for alternative methods with a similar accuracy has been constant in the field of research. To test the effectiveness of new methods, it is essential to make a tough analysis of the results in relation to radiosondes, through evaluation metrics and statistical graphs that allow quantifying and qualifying the accuracy of these new methods. In this context, this study was proposed to evaluate different atmospheric profiling methods at the Izobamba station in the capital of Ecuador, Quito.

The constant development of technology has allowed more efficient measurement methods to be used more frequently. Satellites are becoming more accurate and have worldwide coverage, and despite the limitations presented, such as data loss in case of high cloud coverage or data gaps, they generate usable profiles although the vertical scale still represents a drawback to expand their use, as observed in Quito, where under the conditions of the study MODIS profiles accuracy was not comparable to other methods, and AIRS profiles, which, although significantly better than those of MODIS, have a vertical and horizontal resolution that does not allow the development of fine studies at local level. To avoid these problems, modeling can be performed with programs of increasing efficiency, resolution, and complexity, such as WRF, which allows both horizontal and vertical modeling of the atmosphere. Simulation models have more parameters and updates that offer a wide range of possibilities; however, at the same time, this makes their usage more complex and require in-depth research to find the best configurations depending on the geographical area in which their use is required. Finally, there are also direct measurement techniques, such as radiosondes, but more economical as unmanned aerial vehicles and/or surface measurement sensors. These methods are the most recently studied because over time they have been improved to meet the specific needs of this field of knowledge. However, their use also requires the evaluation of the effect of the equipment components.

In this study, four vertical profiling methods were tested compared with a campaign of 10 radiosondes performed at the Izobamba station near Quito, Ecuador. It was determined through vertical profiles, Taylor diagrams, and evaluation metrics that the most effective data in terms of magnitude are the profiles captured using the drone equipped with the same sensor of the radiosondes. Although in magnitude it is the one that best represents the profile, there are limitations, such as the low correlation of the data, that can be attributed to several factors that should be further evaluated, however, the sharp variations in temperature data can be simply corrected with a quadratic adjustment that significantly increases the correlation of the profile, for the other variables, a more complex adjustment should be further analyzed. The profiles extracted from the WRF modeling also represent a viable alternative for atmospheric profiling; these data have a high correlation with the radiosonde profiles, though they present errors in magnitude greater than those obtained with the drone. The parametrizations used generated a favorable scenario in the study area although the hPBL values were not accurate. The time of launching was a limiting factor to determine the values of the hPBL, due to its low height in the early morning hours, and its fast increase, especially when there are clear sky conditions. AIRS constitutes a viable source of information, especially with its better accuracy with variables influenced by humidity, such as Q, but its resolution and data gaps are still a limitation. Finally, the profiles extracted from MODIS product MOD07_L2 are the poorest performer, even though in specific variables, such as temperature, have a fair representation with respect to the radiosondes. The limitations regarding cloudiness and the low vertical resolution of the profiles can be further analyzed with the option of interpolation. However, in this study, it is not considered a viable alternative to radiosondes.

The limitations of the analyzed techniques make their use challenging, and it is important to perform a specific analysis for each of them to be considered viable alternatives to radiosondes. Referring to drones, it is of utmost importance to characterize the effects of factors, such as their ascent speed, the maximum height they reach, the effect of turbulence during ascent, the interference they may cause in the connection between the sonde and the station, the batteries ‘heating’, the blades influence on the measurement, among others. A complete study should be conducted to evaluate these effects and test drones that are specific for vertical measurements at high altitudes. Regarding the WRF, a test of parametrizations and a sensitivity analysis to changes in all the profiles generated by using different microphysics schemes and changing parametrizations should be performed. This would allow improving the model’s ability to represent variations in the profile or specific variables such as the hPBL. Additionally, the modeling approach for Quito’s meteorological variables requires further studies on the influence of other surface-layer and land-surface schemes and other physics options. About satellite data, the interpolation in the vertical profiles, and the interpolation of missing data should be further analyzed to determine the effectiveness and the best techniques for a more robust result.

In addition to future considerations for specific studies with each technique, it is important to note and emphasize that this study was conducted between August 2021 and January 2022, so the results are very specific for that period of the year and the behavior of the profiles in the rest of the year should be also evaluated to see if the results are congruent throughout the year, or if seasonality has an impact on their effectiveness.