Ultra-Light Airborne Measurement System for Investigation of Urban Boundary Layer Dynamics

Winter smog episodes are a severe problem in many cities around the world. The following two mechanisms are responsible for influencing the level of pollutant concentrations: emission of pollutants from different sources and associated processes leading to formation of secondary aerosols in the atmosphere and meteorology, including advection, which is stimulated by horizontal wind, and convection, which depends on vertical air mass movements associated with boundary layer stability that are determined by vertical temperature and humidity gradients. The aim of the present study was to evaluate the performance of an unmanned aerial vehicle (UAV)-based measurement system developed for investigation of urban boundary layer dynamics. The evaluation was done by comparing the results of temperature, relative humidity, wind speed and particulate matter fraction with aerodynamic diameter below 10 μm (PM10) concentration vertical profiles obtained using this system with two reference meteorological stations: Jagiellonian University Campus (JUC) and radio transmission tower (RTCN), located in the urban area of Krakow city, Southern Poland. The secondary aim of the study was to optimize data processing algorithms improving the response time of UAV sensor measurements during the ascent and descent parts of the flight mission.


Introduction
Air quality in urban and suburban areas is a major long-term public health problem. Recent reports of air quality in Europe [1] indicate that many Polish cities do not meet the required standards. Many actions have been undertaken to improve the air quality in Poland. However, this is a complex problem that is connected to both the distribution and efficiency of pollution emission sources and local meteorological conditions. Particulate matter (PM) concentration is a parameter that is used to determine air quality. A high PM concentration affects human health [2,3], decreases visibility [4] and affects the global climate.
The concentration of individual compounds in ambient air depends on many factors, such as the location of emission sources, land forms, degree of urbanization and particular local environmental conditions. Studies of severe haze episodes in Beijing indicated that they are largely driven by meteorological conditions. High aerosol concentrations were recorded with weak wind conditions and air mass stagnation in the city [4]. The next step to better understand the dispersion of pollutants in urban areas is to investigate the boundary layer (BL) dynamics. A key element in such studies is to determine the vertical atmospheric profiles. Studies of vertical atmospheric profiles can be performed using devices such as balloon sounding [5], light detection and ranging (LIDAR) [6,7] and sound detection article presents configuration of the system, calibration of the sensors and full sequence of a data postprocessing procedure consisting of data smoothing and correction of sensors' response time. Finally, the performance of the system is demonstrated in the example of three atmospheric profiling measurement campaigns by comparing obtained results with other available static measurement points available in the vicinity of flight location. According to our literature studies, this is one of the few articles presenting a complete and comprehensive description of the full data flow algorithm validated by the analysis of real cases [21,30,31]. The system has been developed by the Environmental Physics Group at the AGH University of Science and Technology.
The second section presents the architecture of the system, sensor calibration and numerical procedures aimed at improving data quality consisting of time response error correction and smoothing procedure. Section Results presents the performance of the system demonstrated by the example of three measurement campaigns performed in Krakow, Southern Poland. These examples were aimed at the characterization of the transition phase of the boundary layer between stable and convective after sunrise and convective and stable after sunset. The presented results allowed us to determine the impact of the formation of convective and stable nocturnal boundary layers during the transition phase on vertical pollutant distribution in urban areas. To date, only a few studies presenting the formation of convective and stable nocturnal boundary layers using a drone measurement system in the city have been reported, e.g., [32].

AirDust Measurement System
Figure 1a presents an unmanned aerial vehicle (UAV) equipped with air quality and meteorological component sensors. The ultra-light measurement system dedicated for temperature, relative humidity, pressure and particulate matter (PM 10 fraction) measurements is presented in Figure 1b. It was built based on the Arduino MKR ZERO microcontroller, which is responsible for communication with the sensors, storing the measurements on the memory card and sending information in real time to the ground station using a 433-MHz radio line. The system was powered by a drone battery, and for this purpose, a step-down converter reducing the input voltage from ca. 22 V to 5 V at 2 A was used. Temperature, relative humidity and atmospheric pressure were measured using a digital humidity, pressure and temperature sensor (BME280, Bosch Sensortec GmbH, Reutlingen, Germany). It communicated with the microcontroller via an I2C interface. Additionally, the AirDust ultra-light measurement system (b) for temperature, relative humidity, atmospheric pressure, particulate matter and GPS position equipped with a data logger. The system was powered by a drone battery, and for this purpose, a step-down converter reducing the input voltage from ca. 22 V to 5 V at 2 A was used. Temperature, relative humidity and atmospheric pressure were measured using a digital humidity, pressure and temperature sensor (BME280, Bosch Sensortec GmbH, Reutlingen, Germany). It communicated with the microcontroller via an I2C interface. Additionally, the measurement system was equipped with a thermocouple type T to monitor rapid temperature changes. For PM 10 measurements, a PMS7003 low-cost sensor (Plantower, Beijing, China) based on the light scattering method was used. To reduce the impact of water vapor condensation at high humidity and low air temperature on PM measurement, sensor air inlet was heated. Additionally, to reduce the air turbulences at the air inflow to the sensor, the air inlet was equipped with a cap located in the upper part of the drone providing a laminar flow during the drone's climb. The measurement system was supplemented with the NEO-7 GNSS module (u-blox AG, Thalwil, Switzerland) allowing monitoring of the horizontal and vertical sensor positions. Additionally, the system is equipped with an electrochemical methane (CH 4 ) sensor for detection of high concentration natural gas plumes, but this part was not used in the presented study. A list of components and their technical parameters provided by a producer is included in Figure 1 and Table 1. The whole system is housed in 11 × 9 × 5 cm 235 g of weight box fixed to the drone platform.

TriSonica Wind Sensor
The second element of the presented airborne measurement system is an acoustic wind sensor using a TriSonica TM Mini Wind and Weather Sensor (Anemoment LLC, Longmont, CA, USA). This module is also equipped with a data logging module based on the Arduino MKRZERO microcontroller and the NEO-7 GNSS module (u-blox AG, Thalwil, Switzerland) allowing monitoring of the horizontal and vertical sensor positions. Wind sensor is a 50 g miniature 3D sonic anemometer supplemented with temperature, pressure and humidity sensors (located at the bottom side of the TriSonica sensor), magnetometer and 3D accelerometer ( Figure 2).

TriSonica Wind Sensor
The second element of the presented airborne measurement system wind sensor using a TriSonica TM Mini Wind and Weather Sensor (Anemome mont, CA, USA). This module is also equipped with a data logging modul Arduino MKRZERO microcontroller and the NEO-7 GNSS module (u-blo Switzerland) allowing monitoring of the horizontal and vertical sensor p sensor is a 50 g miniature 3D sonic anemometer supplemented with temper and humidity sensors (located at the bottom side of the TriSonica sensor), and 3D accelerometer ( Figure 2).  The anemometer measures wind speed (0-50 m·s −1 ) and wind direction (0-360 • ) along with the actual position of the sensor with a 10 Hz data rate. To avoid measurement disturbance due to turbulences generated by the drone's propellers, a wind sensor is fixed on a small boom ca. 30 cm above propellers' level (see Figure 2). The idea of wind direction correction using an onboard magnetometer during rotation of the drone during the flight was tested. The results show that magnetic field disturbance generated by the currents supplying the drone's motors makes it impossible to use this sensor for a correction procedure. Reliable wind direction data can be acquired only by avoiding drone rotation during the flight. A list of components and their technical parameters is included in Table 2. Test of air temperature, air humidity and atmospheric pressure sensors in the TriSonica sonic anemometer showed that their inertia is more significant than sensors used in the AirDust system (not shown in the article). Sensors are located inside the anemometer shield which causes their poor ventilation and longer time response. Due to this fact, only measurements of wind components and sonic temperature were used in the analysis.

Sensors Calibration
Calibration of meteorological and wind sensors was performed by comparison with the stationary meteorological station at the Faculty of Physics and Applied Computer Science, AGH University of Science and Technology (AGH UST). The meteorological station was located on the roof of the faculty building. Air temperature, relative humidity and atmospheric pressure, at the level of 20 m a.g.l., were measured using a VAISALA WXT520 platform. Additionally, particulate matter fraction measurements (PM 10 ) were conducted at the same point by a TSI 8530 optical aerosol monitor calibrated using the gravimetric method. Figure 3 presents the comparison of tested sensors with reference instruments for three days of tests with one-minute resolution (VAISALA WXT520 platform saves data once per minute). The calibration equation and correlation coefficient R 2 were obtained using linear regression. In the plots in Figure 3, there are presented linear regression equations for individually tested parameters. For calibration of the BOSH atmospheric pressure and PM 10 concentration sensor, the intersection points of the straight lines matching the graph were set to 0. Tests of the Plantower sensor showed the correct measurement for a concentration close to 0 µg·m −3 . Correlation coefficient R 2 calculated for all sensors except the BOSH air temperature sensor was higher than 0.9. The highest correlation coefficient value was obtained for the PM 10 concentration sensor and the BOSH air humidity sensor (above 0.99). Calibration of sensor Plantower indicated that concentration of PM 10 was overestimated, therefore the slope of a regression line is equal to 0.85. Accuracy of all sensors used in the study calculated based on the accuracy provided by the manufacturer and calibration were equal for: temperature sensors from 0.6 • C (thermocouple) and 0.7 • C (BOSH temperature sensor) to 1.2 • C (TriSonica TM -sonic temperature) in the range from 0 to 20 • C; for BOSH air pressure sensor, it was equal to 0.6 hPa (in the range of 300-1100 hPa, 0-65 • C); 2% for BOSH air humidity sensor (in the range of 20-90%); 0.1 m·s −1 for TriSonica TM wind speed sensor (range from 0 to 10 m s −1 ); 6 µg·m −3 for PM10 sensor error in the range from 0 to 100 µg·m −3 .
Formulas of linear regression were used to correct the systematic error of measurements presented in this article.
Sensors 2021, 21, 2920 7 of 32 was overestimated, therefore the slope of a regression line is equal to 0.85. Accuracy of all sensors used in the study calculated based on the accuracy provided by the manufacturer and calibration were equal for: temperature sensors from 0.6 °C (thermocouple) and 0.7 °C (BOSH temperature sensor) to 1.2 °C (TriSonica TM -sonic temperature) in the range from 0 to 20 °C; for BOSH air pressure sensor, it was equal to 0.6 hPa (in the range of 300-1100 hPa, 0-65 °C); 2% for BOSH air humidity sensor (in the range of 20-90%); 0.1 m⋅s −1 for TriSonica TM wind speed sensor (range from 0 to 10 m s −1 ); 6 μg⋅m −3 for PM10 sensor error in the range from 0 to 100 μg⋅m −3 . Formulas of linear regression were used to correct the systematic error of measurements presented in this article.

Smoothing Procedure
To remove the noise related to the sensors, the local polynomial regression fitting method was used. Locally estimated scatterplot smoothing (LOESS) is one of the most common representations of the above procedure that combines multiple regression models in a k-nearest-neighbor-based meta-model [33]. At each point within the data set range, a low-degree polynomial was fitted to a data subset. The polynomial was fitted using weighted least squares, giving more weight to the points near the area where the response was being estimated and less weight to points that were further away. The value of the regression function for the point was then obtained by evaluating the local polynomial using the explanatory variable values for that data point. The LOESS fit was complete after the regression function values had been computed for each from the n data points. The method was implemented using R software, and it was called the LOESS function [34]. The vertical profile of air temperature without smoothing (raw measured data) and the results of the LOESS method are presented in Figure 4. The degree of smoothing was set to 0.1, which means that 10% of the data were used in the calculations for each point. Using higher values caused data smoothing that was too strong, and therefore, too much valuable information about the atmospheric vertical profiles was lost. The degree of the polynomial used in the calculation is a default value that was equal to 2. The smoothing method was used to remove noise from the vertical profiles of the air temperature, air humidity, air pressure and PM 10 concentration profile. set to 0.1, which means that 10% of the data were used in the calculations for each point. Using higher values caused data smoothing that was too strong, and therefore, too much valuable information about the atmospheric vertical profiles was lost. The degree of the polynomial used in the calculation is a default value that was equal to 2. The smoothing method was used to remove noise from the vertical profiles of the air temperature, air humidity, air pressure and PM10 concentration profile.

Correction of Sensors' Response Time
Calibration of sensors presented in Figure 2 was made for slow-changing meteorological parameters, therefore inertia of the tested sensor was invisible. During the evening or early morning, measurements of vertical profiles of atmosphere differences between ground measurements and measurements at a few hundred meters can be significant. Therefore, measurements of meteorological components [35] and gas concentrations [36] in the atmosphere are subject to several sensor-related sources of error. One of the most

Correction of Sensors' Response Time
Calibration of sensors presented in Figure 2 was made for slow-changing meteorological parameters, therefore inertia of the tested sensor was invisible. During the evening or early morning, measurements of vertical profiles of atmosphere differences between ground measurements and measurements at a few hundred meters can be significant. Therefore, measurements of meteorological components [35] and gas concentrations [36] in the atmosphere are subject to several sensor-related sources of error. One of the most common in atmospheric measurements is time-lag error, which occurs in radiosonde profiles [37] or mobile measurements [38]. It is caused by the fact that sensors need a certain time in order to reach equilibrium with their ambient environment. To correct this error for humidity and temperature sensor, the procedure based on Equation (1) [39] time-dependent response presented below was developed.
U a and U m are ambient and measured parameter values, respectively; the coefficient k is sensor time constant dependent on atmospheric conditions. "The sensor time constant is not the same for all atmospheric conditions. For temperature sensors, the sensor time constant is found to be primarily dependent on the atmospheric pressure (air density) and the ventilation of the sensor. For relative humidity sensors, the sensor time constant is additionally a function of absolute temperature and water vapor content. A decrease in temperature and atmospheric pressure leads to increased sensor response times, and so does a lack of ventilation and low water vapor content" [39]. Time-lag error in vertical profiling depends also on the vertical velocity of sensors.
Solving Equation (1) for U a provides a simple way to compute a true ambient parameter from measurements [36]. Before calculation of a derivative, original datasets had to be smoothed. Earlier studies of this issue confirmed the need to smooth the data before correction [37][38][39]. To compute derivative dU m /dt, a 4th-order Savitzky-Golay filter was used [40,41] with a length of 11 samples, without losing detailed information Sensors 2021, 21, 2920 9 of 32 of the vertical structure of the atmosphere. Calculated derivative dU m /dt was smoothed to omit sudden changes (peak values). It should be noted that to compute a derivative, homogenous-in-time data with known time resolution are needed.
To remove the noise related to the sensors, the local polynomial regression method was used, and detailed information about this method is presented in Section 2.3.
The coefficient k was determined for each flight by an optimization that minimizes the difference between the corrected humidity and temperature profiles for ascent and descent profiles.
In order to validate the calculations, tests of time response of temperature and humidity sensors for big amplitudes of temperature and humidity in indoor-outdoor tests in the cold season were conducted. For this case, the sensors were moved inside and outside the window assuming that the time course of temperature and humidity should be a quadratic wave. By comparing results with the theoretical time course, k coefficient was calculated. Tests indoor-outdoor lasted several minutes until sensors stabilized. Obtained time curves of temperature and humidity were smoothed using 10% of the data. Calculation of derivative dU m /dt was obtained by using 11 points. Calculated derivatives' curves were also smoothed by using 10% of data. Figure 5 presents an example plot of air humidity and temperature time series before and after correction. At the moment of putting outdoor or indoor humidity sensor, a sudden peak of humidity was observed (caused by condensation and drying of the plate on which the sensor is located)-Figure 5e,f. In order to omit enhancement of the sudden peak of humidity observed during the change of air conditions, some data were removed, and results are presented in Figure 5c,d. Optimal k coefficient was found to equal 20, 120 and 220 for the thermocouple, BOSH temperature sensor and BOSH humidity sensor, respectively.
Tests of temperature change rate at indoor-outdoor conditions pointed out that time response and temperature change rate of the thermocouple was significantly faster than the BOSH sensor. After applying a correction method to the thermocouple and BOSH time course of temperature, differences between both sensors' responses were smaller.
Detailed information about the time response of the sensor, length of sensor stabilization and temperature/humidity amplitude are included in Appendix A Table A1. Analysis of results presented in Figure 5 pointed out that the time response of temperature sensors on condition change was faster even by 60 s; before correction of the time course of the thermocouple, there was a visible sharp change of temperature in the first few seconds. Due to the fact that the optimal k coefficient for the thermocouple is relatively small compared with the coefficient for the BOSH temperature sensor, the rate change of temperature in 1 s is higher from the reference value by 15%-20%. Significant differences between reference time course and time course after correction are visible at air temperature and humidity courses from the BOSH sensor (Figure 5a-d). Length of the period during temperature or humidity change in the range from 10% to 90% of maximum amplitude is two or even three times shorter than the reference value.
However, differences of the rate of parameter change in time after correction are visible between heating and cooling of the sensors, from 15-20% lower for the thermocouple and air humidity BOSH sensor to 140% lower for temperature BOSH sensor for heating.
Presented procedures were applied to vertical profile correction. During the flight operation, vertical drone speed can differ slightly, among others during landing and takeoff or caused by the disturbances (wind gusts), and due to this fact density of points in the vertical profile slightly fluctuates. Because of uneven distribution of measuring points as a function of height, smoothing original data and calculated derivative were made to the altitude, not time. Smoothing procedure was applied to the original data and calculated derivative. Depending on the amount of measurement points, the length of the sample used for smoothing varied from 20% to 30% of all data in the profile. temperature or humidity change in the range from 10% to 90% of maximum amplitude is two or even three times shorter than the reference value.
However, differences of the rate of parameter change in time after correction are visible between heating and cooling of the sensors, from 15%-20% lower for the thermocouple and air humidity BOSH sensor to 140% lower for temperature BOSH sensor for heating. Presented procedures were applied to vertical profile correction. During the flight operation, vertical drone speed can differ slightly, among others during landing and takeoff or caused by the disturbances (wind gusts), and due to this fact density of points in the vertical profile slightly fluctuates. Because of uneven distribution of measuring points as a function of height, smoothing original data and calculated derivative were made to the altitude, not time. Smoothing procedure was applied to the original data and calculated  The flight corresponding to the ascent and descent of the drone was treated as two separate sets of data; for one data set containing both phases, the method did not give promising results.
In order to determine the convergence of the vertical profile measured during the ascent and descent parts of the profile, the mean square error (RMSE) value was determined. Calculations of RMSE were conducted for the pair of points with an altitude difference of less than 4 m.
Because during the last phase of landing at a low altitude disturbances generated by drone are significant, measurements at the lowest layer up to altitude 10-15 m above the ground level were omitted. Figure 6 presents the examples of temperature and humidity vertical profiles before and after correction.

Study Area
Krakow is the second-largest city in Poland (area: 327 km 2 , population: 771,000 inhabitants), located in the Lesser Poland region, in a large valley of the Vistula River. The city's area includes three different regions and geological structures, i.e., the Polish Uplands, the Western Carpathians and the basins of the Carpathian Foredeep [42]. The dominant part of the city is located in the third region, with a diversified natural environment. The hilltops bordering the city to the north and the south reach about 100 m above the river valley floor, similar to the hilltops in the western part of the valley which means that the city is located in a semi-concave land form (open only to the east) and sheltered from the prevailing western winds (Figure 7). This location is responsible for frequent air temperature inversions and the prevailing weak western winds, which contribute to the poor ventilation in the city. A clear reduction in the wind speed in the city is reflected by the frequent incidence of atmospheric calm. Analysis of data from 2001-2010 indicated that there was atmospheric calm in Krakow during an average of 16.8% of the year [43]. Mobile measurements of air temperature during calm and cloudless nights indicate that most probably katabatic flows occur on the slopes of the hills which surround the city, but the resulting cold air pool is formed only in rural areas, i.e., those flows do not enter the ur- Coefficient k for the thermocouple and temperature BOSH sensor was set to 24 and 44, respectively. Optimal coefficient for the BOSH air humidity sensor was found to be 50. After the correction method, we could observe that differences between ascent and descent vertical profiles were significantly smaller. Differences between the BOSH sensor and thermocouple at an altitude between 0 and 200 a.g.l. were reduced by the presented procedure. Air humidity profiles during descent and ascent of the drone after correction procedure still differed at the lowest part of the atmosphere, but a significant improvement was observed. Relative humidity profiles presented in Figure 6b show that this method has still some limitations.

Study Area
Krakow is the second-largest city in Poland (area: 327 km 2 , population: 771,000 inhabitants), located in the Lesser Poland region, in a large valley of the Vistula River. The city's area includes three different regions and geological structures, i.e., the Polish Uplands, the Western Carpathians and the basins of the Carpathian Foredeep [42]. The dominant part of the city is located in the third region, with a diversified natural environment. The hilltops bordering the city to the north and the south reach about 100 m above the river valley floor, similar to the hilltops in the western part of the valley which means that the city is located in a semi-concave land form (open only to the east) and sheltered from the prevailing western winds (Figure 7). This location is responsible for frequent air temperature inversions and the prevailing weak western winds, which contribute to the poor ventilation in the city. A clear reduction in the wind speed in the city is reflected by the frequent incidence of atmospheric calm. Analysis of data from 2001-2010 indicated that there was atmospheric calm in Krakow during an average of 16.8% of the year [43]. Mobile measurements of air temperature during calm and cloudless nights indicate that most probably katabatic flows occur on the slopes of the hills which surround the city, but the resulting cold air pool is formed only in rural areas, i.e., those flows do not enter the urbanized area and have no impact on the urban canopy layer inside the city [44].
ences between measurement results that originated from flights and the meteorological station may be caused by the influence of microclimatic conditions in both locations. The flight point is located in open space, while the station is located between three-story University buildings.
Vertical profile measurements of air temperature and relative humidity corresponding to altitudes of 50 m a.g.l. and 100 m a.g.l. were compared with the measurements from the RTCN tower, located 2.7 km north from the flight point. The location of the RTCN tower is shown in Figure 7. Temporal resolution of meteorological parameters from the JU Campus and RTCN tower stations was equal to 10 min; for the AGH UST station, temporal resolution was equal to 1 min; measurements from the synoptic station of Balice had hourly resolution. Location of measurement points used in the study is included in Table  3.   Explanations: numbers as in Table 3. Topographic data used in Figure 7 come from Shuttle Radar Topography Mission database provided by the National Aeronautics and Space Administration (https://www2.jpl.nasa.gov/srtm/ (accessed on 18 April 2021)). Vertical profile measurement data were collected on the Jagiellonian University (JU) Campus (UAV place; Figure 7). The flight point was located in the western part of the city; the terrain was covered with grass, and altitude differences of the surrounding area were up to 15 m (average value, 209 m a.s.l.). The flights were performed up to a height of 500 m a.g.l., and the maximum flight height was determined by the limitations of the aircraft. The location of the vertical measurements using a UAV was chosen based on safety requirements (e.g., distance from buildings and people), and it was a location that was close to the city center that was representative for incoming air masses.
Vertical profile measurements representing the lower part of the profile (altitude up to 5 m a.g.l.) were compared with the observations from a meteorological station located at a distance of 300 m east from the flight point on the JU Campus ( Figure 7). Small differences between measurement results that originated from flights and the meteorological station may be caused by the influence of microclimatic conditions in both locations. The flight point is located in open space, while the station is located between three-story University buildings.
Vertical profile measurements of air temperature and relative humidity corresponding to altitudes of 50 m a.g.l. and 100 m a.g.l. were compared with the measurements from the RTCN tower, located 2.7 km north from the flight point. The location of the RTCN tower is shown in Figure 7. Temporal resolution of meteorological parameters from the JU Campus and RTCN tower stations was equal to 10 min; for the AGH UST station, temporal resolution was equal to 1 min; measurements from the synoptic station of Balice had hourly resolution. Location of measurement points used in the study is included in Table 3. Explanations: numbers as in Table 3. Topographic data used in Figure 7 come from Shuttle Radar Topography Mission database provided by the National Aeronautics and Space Administration (https://www2.jpl.nasa.gov/srtm/ (accessed on 18 April 2021)).

Measurement Periods
From December 2017 to November 2018, more than 60 vertical profile test measurements of the atmospheric parameters were performed in Krakow. In the present study, measurements taken on the following 2 days (12 flights) were analyzed: The aim of the measurements was to study the dynamics of the atmospheric boundary layer during transition phases between the convective boundary layer (CBL) and stable boundary layer (SBL) after sunset and between SBL and CBL after sunrise. The measurements were made during weather conditions controlled by a high-pressure system, with clear sky and weak wind at the ground level. High-pressure weather in the summer season is favorable for fog formation and a strong thermal inversion at night and strong convective mixing during the day.
In 2019, the drone system was expanded by a TriSonica sonic anemometer sensor. An example of wind profiles associated with other measured parameters in profiles is presented for the flight campaign performed in March 2019.
Observations during the measurements, synoptic comments (available at http:// www.meteo.pl/ (accessed on 18 April 2021)) and surface pressure charts (available at https://danepubliczne.imgw.pl (accessed on 18 April 2021)) were used to describe the meteorological conditions during the measurement campaigns.

Vertical Profile Measurements
Each measurement campaign consisted of several vertical flights performed over the same location approximately every hour. The duration of each flight lasted from 5 to 15 min depending on the maximum altitude and the climb velocity, which ranged from 1 m s −1 to 2 m s −1 . The maximum drone climb velocity (2 m s −1 ) was determined based on the assumed flight altitude and battery operating time. It was tested that, for velocities below 2 m s −1 , the obtained results after applying all postprocessing algorithms described in Section 2 reflected actual variability in the vertical profile. Because of strong turbulence generated by the multirotor propellers below the aircraft, the sensors were located on the upper part of the drone, and only the ascending part of the flight was analyzed.

Results
In this section, vertical profile measurements from three campaign days are analyzed. Figures presented in Sections 3.1-3.3 show vertical profiles of the air potential temperature, relative humidity and PM 10 concentration for all campaigns. Additionally, for the third campaign wind speed profiles are presented in Section 3.3. In subsections below are presented the time series of air temperature, humidity and wind speed recorded at the surface station of the JU Campus. The red points added to charts represent measurements conducted by the UAV platform at the start of each flight. Thermal stratification of the atmosphere in Krakow during the measurement days are presented in the last Figures at Sections 3.1-3.3, using data from two different levels at the RTCN tower (50 and 100 m a.g.l.). Points presented in those Figures represent measurements of air temperature made by a drone during each flight at the altitudes that were closest to the tower sensors' elevations.
Due to the fact that the measurement campaign was conducted during nighttime and daytime conditions, temporal variation of vertical PBL structure was observed. During each campaign period, ground thermal inversion with and ground layer with increased air humidity were observed. Depending on temporary conditions, the depth of both layers could be different.
For this purpose, the vertical gradient of temperature/humidity was calculated. A location where significant gradient change is observed in the vertical profile of potential temperature is assigned as estimated ground inversion layer height (see Figure 8). Studies of the vertical air humidity profile and its vertical stratification will be continued in further work in order to better understand the phenomenon. Atmospheric stability was estimated based on the vertical gradient of potential temperature [45]. Based on the analysis of all vertical profiles from the measurement campaign of potential air temperature and climatological analysis of the analyzed region (not shown) for neutral and stable conditions, potential temperature gradient interval was set to ±0.7 • C/100 m and >0.7 • C/100 m, respectively. Cases of unstable atmospheric conditions with the potential gradient lower than −0.7 • C/100 m were not observed. General information presenting characteristics of relative air humidity and potential air temperature vertical profiles are included in Appendix A, Tables A2-A4. of the vertical air humidity profile and its vertical stratification will be continued in further work in order to better understand the phenomenon. Atmospheric stability was estimated based on the vertical gradient of potential temperature [45]. Based on the analysis of all vertical profiles from the measurement campaign of potential air temperature and climatological analysis of the analyzed region (not shown) for neutral and stable conditions, potential temperature gradient interval was set to ± 0.7 °C/100 m and >0.7 °C/100 m, respectively. Cases of unstable atmospheric conditions with the potential gradient lower than −0.7 °C/100 m were not observed. General information presenting characteristics of relative air humidity and potential air temperature vertical profiles are included in Appendix A, Tables A2-A4.

18 September 2018-Evening Campaign
Between 17 and 21 September 2018, Poland was under the influence of a high-pressure system that was moving slowly from the west to the east of the European continent. Polar, warm maritime air masses were passing over Poland. During the measurement period of 18 September 2018, the sky was almost cloudless, and the wind in the afternoon was weak (below 2 m s −1 ). At night, local radiation fog was observed close to the ground level at flight point. Sunset for the city of Krakow was at 16:50 UTC, and flights were

18 September 2018-Evening Campaign
Between 17 and 21 September 2018, Poland was under the influence of a high-pressure system that was moving slowly from the west to the east of the European continent. Polar, Sensors 2021, 21, 2920 15 of 32 warm maritime air masses were passing over Poland. During the measurement period of 18 September 2018, the sky was almost cloudless, and the wind in the afternoon was weak (below 2 m s −1 ). At night, local radiation fog was observed close to the ground level at flight point. Sunset for the city of Krakow was at 16:50 UTC, and flights were conducted from 15:00 UTC until 21:00 UTC. Figure 9 presents vertical profiles obtained during this campaign. The maximum altitude of the measurements varied between individual flights from 200 m a.g.l. up to 500 m a.g.l.  Air temperature during the first four flights measured by a drone differed from the measurements conducted at the Campus JU station, and the differences ranged from 1.8 ± 0.6 First measurements were made 2 h 40 min before the sunset. The vertical profile of the potential temperature and relative humidity at the first flight was almost constant, and maximum differences in the potential temperature were approximately 0.8 • C. Starting from the third flight (16:32 UTC), formation of ground thermal inversion was observed. The process started before the sunset, because of strong cooling of the ground at clear sky conditions.
The observed mean gradient of the potential temperature from the ground level up to 300 m a.g.l. at the sixth flight was equal to 3.3 • C/100 m. The largest change in potential air temperature was measured during the third flight. Differences in the potential air temperature between 7 m a.g.l. and 40 m a.g.l. were 6.1 • C (the potential temperature gradient was equal to 18 • C/100 m). Similar phenomena were presented in Banta et al. 2007 [46] defined as a very stable boundary layer and as a shallow boundary layer with a thickness of 10-30 m with weak intermittent turbulence within the strong surface-based radiation inversion. Because of strong surface cooling at night, we observed radiation fog, which was also visible in the vertical profiles of air humidity after 18:00 UTC. Air humidity measurements suggested that the upper limit of the fog was 60 m a.g.l. Subsequent vertical air humidity profiles above 290 m a.g.l. did not differ significantly for flights 2 and 6. Measurements of vertical profiles of potential temperature and relative humidity indicated that in the valley, a cold air pool developed below valley depth. Significant decrease of air humidity and potential temperature in the layer between 60 and 300 m a.g.l. indicates that at the upper layers above the valley, advection of drier and warmer air masses occurred.
The PM 10 concentration in the vertical profile measured at 15:09 UTC also confirmed that the lower part of the planetary boundary layer (PBL) was homogeneous, and values from the ground up to 400 m a.g.l. were less than 20 µg m −3 . PM 10  Possible explanation of this phenomenon is occurrence of pollutant transport from distant areas or influence of residual turbulence on pollutants trapped in the ground inversion layer and residual layer located above.
Second possible scenario is occurrence of a nocturnal low-level jet above valley depth, which can modify the dispersion of pollutants. The long-term studies for the period 1994-1999 with the use of SODAR for this region showed the presence of a significant wind direction change during the day in the cold season and in the warm season during the night. The largest wind direction change was observed in the layer of 100-150 m a.g.l. Additionally, attention should be paid to the frequent occurrence of a significant wind speed increase in the layer 50-170 m a.g.l. during the night in the warm season (nocturnal low-level jet) [29]. Figure 10 presents the air temperature, humidity and wind conditions on 18 September 2018. During the drone measurements, the mean wind speed was below 2 m s −1 , with a maximum speed reaching 3 m s −1 ; between 14:00 UTC and 16:30 UTC, the wind direction was south, but between 16:30 UTC and 17:30 UTC, the direction changed to west-south-west. Maximum air temperature measured on the JU Campus was 26.4 • C at 13:40 UTC. Air temperature during the first four flights measured by a drone differed from the measurements conducted at the Campus JU station, and the differences ranged from 1.8 ± 0.6 • C to 3.6 ± 0.6 • C. The differences during the afternoon hours may be caused by well-developed thermal turbulence, which generates significant microclimatic variability in the air temperature and relative humidity. The turbulence is the result of the differences in air temperature over various active surfaces which can be found around the measurement site and the campus. Each active surface shows a different structure of the heat balance [47]. Air parcels of various properties are generated over each type of the active surface, and then they are moved freely so the sensors measure the properties of accidental air parcels. During the night, when the thermal turbulence is much weaker, the air temperature measurements from both sources were consistent, and the maximum difference in air temperature in flights 5 and 6 was 0.7 ± 0.6 • C. Measurements of air humidity during almost all the flights did not differ significantly from measurements made at the station, except for measurements during flight 4 when the difference was 15%.
Sensors 2021, 21, x FOR PEER REVIEW 18 of 33 and then they are moved freely so the sensors measure the properties of accidental air parcels. During the night, when the thermal turbulence is much weaker, the air temperature measurements from both sources were consistent, and the maximum difference in air temperature in flights 5 and 6 was 0.7 ± 0.6 °C. Measurements of air humidity during almost all the flights did not differ significantly from measurements made at the station, except for measurements during flight 4 when the difference was 15%.
(a) (b) (c) (d) Figure 10. Wind direction (a), wind speed (b), air temperature (c) and relative humidity (d) measured at the JU Campus meteorological station on 18 September 2018. Red dots show measurements from the drone that were made at the start of each flight. Figure 11 presents measurements of air humidity and temperature on 18 September 2018 that were measured at the RTCN tower and using a drone at the following two levels: 50 m a.g.l. and 100 m a.g.l. Observations of air temperature and humidity that were conducted using a drone during the day (first flight) were slightly different from the tower site for both levels. Temperatures measured using a drone were approximately 2 °C higher than the tower, and relative humidity was approximately 8% lower, but the observed gradients between 50 and 100 m a.g.l. were the same for both sites. For the later flights, the temperature that was measured using the drone was lower, and the relative humidity was higher than the tower results, especially for 50 m a.g.l. The differences in meteorological conditions between the two places were possibly caused by different surface radiation budgets because of differences in land cover that was adjacent to both sites. The second possibility of significant differences between layer 50 and 100 m a.g.l. is topographic relief; in the west from the RTCN tower are hilltops with altitudes that reach from 300 to 400 m a.s.l. (Figure 7, RTCN tower-point no. 3), while at the nearest region to the flight place, the altitude varies between 200 and 250 m a.s.l. Due to this fact, the tower station may be under the influence of topographic shading at selected atmospheric conditions (i.e., western advection). On the other hand, inconsistences between relative humidity measurements can be partially caused by insufficient response time correction of the sensor on the drone (see Figure 6b).  Figure 11 presents measurements of air humidity and temperature on 18 September 2018 that were measured at the RTCN tower and using a drone at the following two levels: 50 m a.g.l. and 100 m a.g.l. Observations of air temperature and humidity that were conducted using a drone during the day (first flight) were slightly different from the tower site for both levels. Temperatures measured using a drone were approximately 2 • C higher than the tower, and relative humidity was approximately 8% lower, but the observed gradients between 50 and 100 m a.g.l. were the same for both sites. For the later flights, the temperature that was measured using the drone was lower, and the relative humidity was higher than the tower results, especially for 50 m a.g.l. The differences in meteorological conditions between the two places were possibly caused by different surface radiation budgets because of differences in land cover that was adjacent to both sites. The second possibility of significant differences between layer 50 and 100 m a.g.l. is topographic relief; in the west from the RTCN tower are hilltops with altitudes that reach from 300 to 400 m a.s.l. (Figure 7, RTCN tower-point no. 3), while at the nearest region to the flight place, the altitude varies between 200 and 250 m a.s.l. Due to this fact, the tower station may be under the influence of topographic shading at selected atmospheric conditions (i.e., western advection). On the other hand, inconsistences between relative humidity measurements can be partially caused by insufficient response time correction of the sensor on the drone (see Figure 6b).

21 September 2018-Morning Campaign
On 21 September 2018, Poland was under the influence of a high-pressure system that was located on the border between Ukraine and Moldova. In the morning, the sky was clear, and fog was observed from the ground level to a height of several dozen meters. At the synoptic station of Balice, visibility until 7 UTC did not exceed 10 km; in the valley, atmospheric silence occurred (wind speed below 1 m s −1 ). Figure 12 presents the results of vertical profiles obtained during the flights from September the 21st. The first flight started before sunrise, and profiling was continued until forenoon. Sunrise for that day for Krakow was at 4:22 UTC, and measurements were conducted from 4:00 UTC to 9:30 UTC. Maximum altitude of the measurement ranged from 200 m a.g.l. to 295 m a.g.l.
During the early morning hours, before sunrise, atmospheric conditions were stable (cold air pool developed in the valley), which was confirmed by the vertical profiles of potential air temperature from the first and second flights (the depth of the ground inversion layer at the first and second flights was equal to 230 and 190 m. a.g.l., respectively) (Appendix A, Table A3). The gradient of potential air temperature for both flights was 5.2 °C/100 m (between 0 m a.g.l. and 295 m a.g.l.) and 6 °C/100 m (between 0 m a.g.l. and 195 m a.g.l.), respectively. During subsequent flights, the stable boundary layer gradually disappeared as a result of surface heating. Vertical profiles from flights 3 and 4 indicated the disappearance of the near-ground inversion (neutral conditions in the lowest layer with a depth of ca. 100m a.g.l.) and the appearance of the upper inversion layer. Vertical potential temperature gradient in the layer between 200 and 100 m a.g.l was equal to 5.4 °C/100 m and 3.9 °C/100m during the third and fourth flights, respectively. After sunrise, the air layer close to the ground quickly heated up, which triggered convection processes. Vertical profiles of potential air temperature from flights 5 and 6 (after 8:00 UTC) presented a developed mixed layer with constant potential temperature throughout the profile. Significant differences of the vertical profile of air humidity between the first and second flights at a layer between 70 and 200 m a.g.l. indicate that advection of dry warmer air masses caused destruction of the cold air pool from upper layers.
During the last flight, we observed a thermal plume between 190 and 250 m a.g.l. The maximum difference between the potential temperature in the plume and the surrounding air was 2.8 °C, while the maximum value was measured at 240 m a.g.l.
Accurate analysis of the subsequent pollutant concentration profiles provided valuable information. First flight in the early morning indicated that PM10 concentration at the ground layer up to 150 m a.gl. slightly increased with the altitude; at the layer above were visible significant fluctuations of PM10 concentration. After the sunrise, vertical profiles of PM10 concentration changed because of convective movements.

21 September 2018-Morning Campaign
On 21 September 2018, Poland was under the influence of a high-pressure system that was located on the border between Ukraine and Moldova. In the morning, the sky was clear, and fog was observed from the ground level to a height of several dozen meters. At the synoptic station of Balice, visibility until 7 UTC did not exceed 10 km; in the valley, atmospheric silence occurred (wind speed below 1 m s −1 ). Figure 12 presents the results of vertical profiles obtained during the flights from September the 21st. The first flight started before sunrise, and profiling was continued until forenoon. Sunrise for that day for Krakow was at 4:22 UTC, and measurements were conducted from 4:00 UTC to 9:30 UTC. Maximum altitude of the measurement ranged from 200 m a.g.l. to 295 m a.g.l.
During the early morning hours, before sunrise, atmospheric conditions were stable (cold air pool developed in the valley), which was confirmed by the vertical profiles of potential air temperature from the first and second flights (the depth of the ground inversion layer at the first and second flights was equal to 230 and 190 m. a.g.l., respectively) (Appendix A, Table A3). The gradient of potential air temperature for both flights was 5.2 • C/100 m (between 0 m a.g.l. and 295 m a.g.l.) and 6 • C/100 m (between 0 m a.g.l. and 195 m a.g.l.), respectively. During subsequent flights, the stable boundary layer gradually disappeared as a result of surface heating. Vertical profiles from flights 3 and 4 indicated the disappearance of the near-ground inversion (neutral conditions in the lowest layer with a depth of ca. 100m a.g.l.) and the appearance of the upper inversion layer. Vertical potential temperature gradient in the layer between 200 and 100 m a.g.l was equal to 5.4 • C/100 m and 3.9 • C/100m during the third and fourth flights, respectively. After sunrise, the air layer close to the ground quickly heated up, which triggered convection processes. Vertical profiles of potential air temperature from flights 5 and 6 (after 8:00 UTC) presented a developed mixed layer with constant potential temperature throughout the profile. Significant differences of the vertical profile of air humidity between the first and second flights at a layer between 70 and 200 m a.g.l. indicate that advection of dry warmer air masses caused destruction of the cold air pool from upper layers.
During the last flight, we observed a thermal plume between 190 and 250 m a.g.l. The maximum difference between the potential temperature in the plume and the surrounding air was 2.8 • C, while the maximum value was measured at 240 m a.g.l.
Accurate analysis of the subsequent pollutant concentration profiles provided valuable information. First flight in the early morning indicated that PM 10 concentration at the ground layer up to 150 m a.gl. slightly increased with the altitude; at the layer above were visible significant fluctuations of PM 10 concentration. After the sunrise, vertical profiles of PM 10 concentration changed because of convective movements.
A significant change in the PM 10 concentration on the first flight was measured at 150 m a.g.l. Due to the fact that the PM 10 sensor inlet was heated, the error caused by significant air humidity change was reduced (provided by laboratory experiments not shown in the article). Sudden change of PM 10 concentration close to the transition layer observed at the vertical profile of air humidity was probably caused by new particle formation (NPF), which was also observed at the top of the atmospheric boundary layer (ABL) in different studies [32]. Research by Platis et al. (2016) [48] indicated that the thermodynamic conditions such as high thermodynamic fluctuations, strong gradients of temperature and humidity near the top of the ABL were favorable for the NPF. Vertical profiles of PM 10 from the second and third flights showed that at an altitude of ca. 200 m a.g.l. transport of pollutants occurred; in subsequent flights, increased PM 10 concentration was not visible at this altitude, which is probably caused by a change of meteorological conditions (wind vertical profile and upward convection movements). The PM 10 concentration in the vertical profile up to 250 m a.gl. was significantly reduced after 9:00 UTC, and the values did not exceed 20 µg m −3 , except local maximum at 100 m a.g.l. (24 ± 6 µg m −3 ). A significant change in the PM10 concentration on the first flight was measured at 150 m a.g.l. Due to the fact that the PM10 sensor inlet was heated, the error caused by significant air humidity change was reduced (provided by laboratory experiments not shown in the article). Sudden change of PM10 concentration close to the transition layer observed at the vertical profile of air humidity was probably caused by new particle formation (NPF), which was also observed at the top of the atmospheric boundary layer (ABL) in different studies [32]. Research by Platis et al. (2016) [48] indicated that the thermodynamic conditions such as high thermodynamic fluctuations, strong gradients of temperature and humidity near the top of the ABL were favorable for the NPF. Vertical profiles of PM10 from the second and third flights showed that at an altitude of ca. 200 m a.g.l. transport of pollutants occurred; in subsequent flights, increased PM10 concentration was not visible at this altitude, which is probably caused by a change of meteorological conditions (wind vertical profile and upward convection movements). The PM10 concentration in the vertical profile up to 250 m a.gl. was significantly reduced after 9:00 UTC, and the values did not exceed 20 μg m −3 , except local maximum at 100 m a.g.l. (24 ± 6 μg m −3 ).   Figure 13 shows meteorological conditions on 21 September 2018. The wind direction during the morning was changing from north to west-southwest at 8:30 UTC, mean wind speed during the night (until 8:00 UTC) was less than 1 m s −1 , and maximum wind speed during the measurements was 8 m s −1 at 10:00 UTC. Measurements of air humidity and temperature that were obtained using the drone were consistent with measurements from the meteorological station.
Sensors 2021, 21, x FOR PEER REVIEW 21 of 33 Figure 13 shows meteorological conditions on 21 September 2018. The wind direction during the morning was changing from north to west-southwest at 8:30 UTC, mean wind speed during the night (until 8:00 UTC) was less than 1 m s −1 , and maximum wind speed during the measurements was 8 m s −1 at 10:00 UTC. Measurements of air humidity and temperature that were obtained using the drone were consistent with measurements from the meteorological station.  Figure 14 presents measurements of air humidity and temperature on 21 September 2018 at the RTCN tower and using a drone at two levels 50 m a.g.l. and 100 m a.g.l. Measurements at both levels differed significantly in the first two flights when air temperature from the drone was lower than the observations from the RTCN tower. Vertical profiles of air temperature and air humidity from the tower station indicated that at the western part of the Vistula Valley, during the night, ground inversion occurred. Differences between layer 50 and 100 m a.g.l. in air humidity points out that during the night there was advection of warmer, drier air close to the valley depth. After sunrise, observations from both places did not differ significantly. Lower air humidity at the flight location in comparison with measurements from the tower station after the third may have been caused by ventilation condition differences of both places.  Figure 14 presents measurements of air humidity and temperature on 21 September 2018 at the RTCN tower and using a drone at two levels 50 m a.g.l. and 100 m a.g.l. Measurements at both levels differed significantly in the first two flights when air temperature from the drone was lower than the observations from the RTCN tower. Vertical profiles of air temperature and air humidity from the tower station indicated that at the western part of the Vistula Valley, during the night, ground inversion occurred. Differences between layer 50 and 100 m a.g.l. in air humidity points out that during the night there was advection of warmer, drier air close to the valley depth. After sunrise, observations from both places did not differ significantly. Lower air humidity at the flight location in comparison with measurements from the tower station after the third may have been caused by ventilation condition differences of both places.

7 March 2019-Whole-Day Campaign
On 7 March 2019, air advection to Poland from the southwest was caused by the influence of the high-pressure system above Russia and the low-pressure system over the British Isles. The movement of both pressure systems to the east resulted in the inflow of a cold front on that day in the afternoon hours. Measurements from the synoptic station of Balice indicated that high cloudiness during the measurement day maintained up to 14 UTC (6-7 octants), and between 14 and 17 UTC, a sudden decrease in cloud cover up to 2 octants and subsequent increase in cloud cover occurred. Sunrise on that day occurred at 5:12 UTC, until 9 UTC there was a weak wind below 2 m s −1 , and afterward, there was an increase in wind speed in the Vistula River Valley.
The last presented campaign was done on March the 7th 2019 between 6:00 and 15:00 UTC. The aim was to observe the development of foehn wind occurring in the Tatra Mountains (ca. 100km south from Krakow). The first flight ( Figure 15) indicated that the thermal inversion reached up to 150 m a.g.l. (potential temperature gradient equal to 5 °C/100 m); during subsequent flights (2-4), its height decreased up to 100 m a.g.l., and the vertical gradient of relative humidity and potential temperature decreased as a result of soil heating and advection of warm air masses associated with the occurrence of the foehn wind. Flights 4-11 indicated that at strong vertical mixing (well-mixed air), vertical gradient of potential temperature and relative humidity during these flights were equal on average to −0.1 °C/100 m and 1.9%/100 m, respectively (Figures 15 and 16). Air humidity in the vertical profile during the measurement campaign did not vary significantly; however, in the first three vertical profiles up to 100 m a.g.l., relative air humidity varied between 40% and 50%. The vertical profiles of PM10 for the first two flights (Figure 17a) showed moderate values (up to 50 μg m −3 ) with fine structure including local concentration peaks at 80 m a.g.l. which was probably connected to the change of vertical gradient of relative air humidity from −29%/100 m in the ground layer up to 100 m a.g.l. at first flight and −33%/100 m in the ground layer up to 80 m a.g.l. to almost 0%/100 m at the upper layer. Later flights until 10:00 UTC showed low constant values, which was an effect of foehn wind intrusion into the valley and exchange of air masses (Figure 17a,b). After 11:00 UTC, the PM10 concentration started to rise (up to 150 μg m −3 ), and profiles become irregular (Figure 17b,c). One possible explanation is long-range transport of particulate matter by strong foehn wind. The second possibility is that ascending airflow lifted particles from the ground, and therefore at the 10th and 11th flights, high concentration peaks in the layer up to 150 m a.g.l were observed. Final interpretation requires deeper analysis to be performed in another study. The wind speed profile showed the evolution of strong wind gradually penetrating deeper parts of the valley where the city is located. The first five vertical profiles (Figure 18a,b), with the exception of flight 2, indicated that horizontal wind speed in the layer up to 80-100 m a.g.l. did not exceed 2 m s −1 , Figure 14. Air temperature (a) and relative humidity (b) at 50 and 100 m a.g.l. measured at the RTCN tower and by UAV on 21 September 2018.

7 March 2019-Whole-Day Campaign
On 7 March 2019, air advection to Poland from the southwest was caused by the influence of the high-pressure system above Russia and the low-pressure system over the British Isles. The movement of both pressure systems to the east resulted in the inflow of a cold front on that day in the afternoon hours. Measurements from the synoptic station of Balice indicated that high cloudiness during the measurement day maintained up to 14 UTC (6-7 octants), and between 14 and 17 UTC, a sudden decrease in cloud cover up to 2 octants and subsequent increase in cloud cover occurred. Sunrise on that day occurred at 5:12 UTC, until 9 UTC there was a weak wind below 2 m s −1 , and afterward, there was an increase in wind speed in the Vistula River Valley.
The last presented campaign was done on March the 7th 2019 between 6:00 and 15:00 UTC. The aim was to observe the development of foehn wind occurring in the Tatra Mountains (ca. 100km south from Krakow). The first flight ( Figure 15) indicated that the thermal inversion reached up to 150 m a.g.l. (potential temperature gradient equal to 5 • C/100 m); during subsequent flights (2-4), its height decreased up to 100 m a.g.l., and the vertical gradient of relative humidity and potential temperature decreased as a result of soil heating and advection of warm air masses associated with the occurrence of the foehn wind. Flights 4-11 indicated that at strong vertical mixing (well-mixed air), vertical gradient of potential temperature and relative humidity during these flights were equal on average to −0.1 • C/100 m and 1.9%/100 m, respectively (Figures 15 and 16). Air humidity in the vertical profile during the measurement campaign did not vary significantly; however, in the first three vertical profiles up to 100 m a.g.l., relative air humidity varied between 40% and 50%. The vertical profiles of PM 10 for the first two flights (Figure 17a) showed moderate values (up to 50 µg m −3 ) with fine structure including local concentration peaks at 80 m a.g.l. which was probably connected to the change of vertical gradient of relative air humidity from −29%/100 m in the ground layer up to 100 m a.g.l. at first flight and −33%/100 m in the ground layer up to 80 m a.g.l. to almost 0%/100 m at the upper layer. Later flights until 10:00 UTC showed low constant values, which was an effect of foehn wind intrusion into the valley and exchange of air masses (Figure 17a,b). After 11:00 UTC, the PM 10 concentration started to rise (up to 150 µg m −3 ), and profiles become irregular (Figure 17b,c). One possible explanation is long-range transport of particulate matter by strong foehn wind. The second possibility is that ascending airflow lifted particles from the ground, and therefore at the 10th and 11th flights, high concentration peaks in the layer up to 150 m a.g.l were observed. Final interpretation requires deeper analysis to be performed in another study. The wind speed profile showed the evolution of strong wind gradually penetrating deeper parts of the valley where the city is located. The first five vertical profiles (Figure 18a,b), with the exception of flight 2, indicated that horizontal wind speed in the layer up to 80-100 m a.g.l. did not exceed 2 m s −1 , which was partially caused by topographic relief (Vistula River Valley depth). Significant increase of wind speed from the ground level up to 150 m a.g.l. (local wind maximum) at second flight was caused by a sudden intrusion of foehn wind from the East into the valley (eastern wind at 7 UTC at the JU Campus station). Sudden temporal intrusion of foehn wind is characteristic of this phenomenon which was confirmed by other studies [49,50].
Sensors 2021, 21, x FOR PEER REVIEW 23 of 33 which was partially caused by topographic relief (Vistula River Valley depth). Significant increase of wind speed from the ground level up to 150 m a.g.l. (local wind maximum) at second flight was caused by a sudden intrusion of foehn wind from the East into the valley (eastern wind at 7 UTC at the JU Campus station). Sudden temporal intrusion of foehn wind is characteristic of this phenomenon which was confirmed by other studies [49,50]. During the last flight, a strong wind speed exceeding 10 m s −1 was observed in the whole profile (Figure 18c).   which was partially caused by topographic relief (Vistula River Valley depth). Significant increase of wind speed from the ground level up to 150 m a.g.l. (local wind maximum) at second flight was caused by a sudden intrusion of foehn wind from the East into the valley (eastern wind at 7 UTC at the JU Campus station). Sudden temporal intrusion of foehn wind is characteristic of this phenomenon which was confirmed by other studies [49,50]. During the last flight, a strong wind speed exceeding 10 m s −1 was observed in the whole profile (Figure 18c).   Figure 19 shows meteorological conditions at the JU Campus station on 7 March 2019. The wind direction during the morning was changing from east to west at 9:00 UTC, wind speed during the night (until 9:00 UTC) was less than 1.5 m s −1 , and maximum wind speed during the measurements reached 10 m s −1 at 14:30 0 UTC.
Increase of wind speed being the effect of foehn wind was associated with temperature increase and significant drop of relative humidity between 6 and 9 UTC at the ground station. Measurements of air humidity and temperature that were obtained using the drone were consistent with measurements from the meteorological station. Due to the fact that temporal resolution of the wind component data from the ground meteorological station is the mean values of the 10-min period, Figure 19b presents mean and maximum wind speed measurements. Wind speed data from the sensor located in the UAV have 1  Figure 19 shows meteorological conditions at the JU Campus station on 7 March 2019. The wind direction during the morning was changing from east to west at 9:00 UTC, wind speed during the night (until 9:00 UTC) was less than 1.5 m s −1 , and maximum wind speed during the measurements reached 10 m s −1 at 14:30 0 UTC.
Increase of wind speed being the effect of foehn wind was associated with temperature increase and significant drop of relative humidity between 6 and 9 UTC at the ground station. Measurements of air humidity and temperature that were obtained using the drone were consistent with measurements from the meteorological station. Due to the fact that temporal resolution of the wind component data from the ground meteorological station is the mean values of the 10-min period, Figure 19b presents mean and maximum wind speed measurements. Wind speed data from the sensor located in the UAV have 1 During the last flight, a strong wind speed exceeding 10 m s −1 was observed in the whole profile (Figure 18c). Figure 19 shows meteorological conditions at the JU Campus station on 7 March 2019. The wind direction during the morning was changing from east to west at 9:00 UTC, wind speed during the night (until 9:00 UTC) was less than 1.5 m s −1 , and maximum wind speed during the measurements reached 10 m s −1 at 14:30 0 UTC. s resolution, therefore during windy conditions, single measurement points can fluctuate significantly. Measurements from the drone conducted at 10 m a.g.l. at the first five flights during weak wind conditions were consistent with observations from the ground station. During the subsequent flights, temporal measurements of wind speed from the drone varied between mean and maximum wind speed observed at the JU ground station for most of the cases, except the sixth and seventh flights.  Figure 20 presents measurements of air humidity and temperature on 7 March 2019 at the RTCN tower and using a drone at two levels, 50 m a.g.l. and 100 m a.g.l. In the first part of the record, a clear vertical gradient of temperature and relative humidity was observed. The same effect was confirmed by the data from flights. In the later part, after increase of wind speed, a strong turbulence in the atmosphere caused the vertical gradients to disappear. Temperature data from the RTCN and drone were consistent while the relative humidity recorded by the drone during a higher wind speed period was lower by 10% compared to tower observations at both elevations (50 and 100 m a.g.l.). Increase of wind speed being the effect of foehn wind was associated with temperature increase and significant drop of relative humidity between 6 and 9 UTC at the ground station. Measurements of air humidity and temperature that were obtained using the drone were consistent with measurements from the meteorological station. Due to the fact that temporal resolution of the wind component data from the ground meteorological station is the mean values of the 10-min period, Figure 19b presents mean and maximum wind speed measurements. Wind speed data from the sensor located in the UAV have 1 s resolution, therefore during windy conditions, single measurement points can fluctuate significantly. Measurements from the drone conducted at 10 m a.g.l. at the first five flights during weak wind conditions were consistent with observations from the ground station. During the subsequent flights, temporal measurements of wind speed from the drone varied between mean and maximum wind speed observed at the JU ground station for most of the cases, except the sixth and seventh flights. Figure 20 presents measurements of air humidity and temperature on 7 March 2019 at the RTCN tower and using a drone at two levels, 50 m a.g.l. and 100 m a.g.l. In the first part of the record, a clear vertical gradient of temperature and relative humidity was observed. The same effect was confirmed by the data from flights. In the later part, after increase of wind speed, a strong turbulence in the atmosphere caused the vertical gradients to disappear. Temperature data from the RTCN and drone were consistent while the relative humidity recorded by the drone during a higher wind speed period was lower by 10% compared to tower observations at both elevations (50 and 100 m a.g.l.).

Discussion
The results obtained during the presented campaigns using an airborne multisensor measurement platform enabled detailed analysis of temporal variation of vertical PBL structure at a local scale during three different synoptic conditions. Results from morning campaigns (21 September 2018 and 7 March 2019) pointed out that in the layer with strong gradients of temperature and humidity near the top of the ground inversion layer, peaks of PM10 concentration were observed. During stable conditions, PM10 concentration in the ground layer up to 100 m a.g.l. (valley depth) did not vary significantly, while at the upper layers, strong fluctuations of PM10 concentration were observed. Possible explanation is weakening of mechanical turbulence in the valley due to topographic shading and occurrence of nocturnal low-level jet above the valley which can induce mechanical turbulence and modify pollutant dispersion.
Other research presenting the occurrence of a layer with increased concentration of air pollution at a certain height points out that these conditions were connected to vertical diffusion-increased concertation of pollutants at a layer a few dozens of meters above the ground [51]-and diffusion of plumes from elevated sources-upper layers at an altitude of a few hundred meters above the ground [52] or downward flows of regional pollutants from the upper layer [53].
Studies of vertical profiles of potential air temperature and relative humidity pointed out that an important factor that affects cold air pool break is advection of warmer air above the valley during morning hours. At the first campaign (September 18, 2018), advection of warmer air strengthened the ground inversion layer during the night. Observed potential temperature gradient up 100 m a.g.l. at sixth flight (after sunset) was equal to 6.3 °C/100 m, while at first flight on 21 September 2018 (before sunrise) was equal to 5.4 °C/100 m.
The temperature data obtained during the first and second campaigns were used for comparison of atmospheric processes occurring during convective and stable nocturnal boundary layer formation. In particular, a potential temperature gradient and the altitude associated with its rapid change were analyzed for subsequent flights. While the first two profiles obtained during the campaign representing the formation of stable nocturnal boundary (15:09 UTC and 16:56 UTC) were almost constant (temperature gradient below 0.3 °C/100 m), next profiles showed temporal increase of nocturnal boundary layer (NBL) depth (defined as the altitude where potential temperature is stabilizing) starting from ca. 100 m a.g.l observed during flight no. 3 and reaching the end value of ca. 300 m a.g.l. for flight 6 (Figure 9b). In the first phase of NBL formation (flight 3), a very steep gradient was observed close to the ground level (ca. + 18 °C/100 m) which gradually disappeared during the next flights when the NBL depth was increasing. The gradient value in the

Discussion
The results obtained during the presented campaigns using an airborne multisensor measurement platform enabled detailed analysis of temporal variation of vertical PBL structure at a local scale during three different synoptic conditions. Results from morning campaigns (21 September 2018 and 7 March 2019) pointed out that in the layer with strong gradients of temperature and humidity near the top of the ground inversion layer, peaks of PM 10 concentration were observed. During stable conditions, PM 10 concentration in the ground layer up to 100 m a.g.l. (valley depth) did not vary significantly, while at the upper layers, strong fluctuations of PM 10 concentration were observed. Possible explanation is weakening of mechanical turbulence in the valley due to topographic shading and occurrence of nocturnal low-level jet above the valley which can induce mechanical turbulence and modify pollutant dispersion.
Other research presenting the occurrence of a layer with increased concentration of air pollution at a certain height points out that these conditions were connected to vertical diffusion-increased concertation of pollutants at a layer a few dozens of meters above the ground [51]-and diffusion of plumes from elevated sources-upper layers at an altitude of a few hundred meters above the ground [52] or downward flows of regional pollutants from the upper layer [53].
Studies of vertical profiles of potential air temperature and relative humidity pointed out that an important factor that affects cold air pool break is advection of warmer air above the valley during morning hours. At the first campaign (18 September 2018), advection of warmer air strengthened the ground inversion layer during the night. Observed potential temperature gradient up 100 m a.g.l. at sixth flight (after sunset) was equal to 6.3 • C/100 m, while at first flight on 21 September 2018 (before sunrise) was equal to 5.4 • C/100 m.
The temperature data obtained during the first and second campaigns were used for comparison of atmospheric processes occurring during convective and stable nocturnal boundary layer formation. In particular, a potential temperature gradient and the altitude associated with its rapid change were analyzed for subsequent flights. While the first two profiles obtained during the campaign representing the formation of stable nocturnal boundary (15:09 UTC and 16:56 UTC) were almost constant (temperature gradient below 0.3 • C/100 m), next profiles showed temporal increase of nocturnal boundary layer (NBL) depth (defined as the altitude where potential temperature is stabilizing) starting from ca. 100 m a.g.l observed during flight no. 3 and reaching the end value of ca. 300 m a.g.l. for flight 6 (Figure 9b). In the first phase of NBL formation (flight 3), a very steep gradient was observed close to the ground level (ca. + 18 • C/100 m) which gradually disappeared during the next flights when the NBL depth was increasing. The gradient value in the upper part of NBL was almost constant for flights 4-6 (ca. +3.3 • C/100 m). During the second campaign representing the formation of the convective boundary layer, potential temperature profiles looked different. In the initial state before sunrise (first flight at 4:14 UTC), a clear thermal inversion layer reaching the height of ca. 300 m a.g.l. was observed. The potential temperature gradient within the NBL layer was equal to ca. 5 • C/100 m. Starting from flight 3 (6:23 UTC), a formation of the convective layer was observed close to the ground. The surface temperature was increasing for subsequent flights and remained constant in the first part of the profile. The depth of the convective layer (defined as the maximum altitude where the temperature was close to the surface) was increasing with time starting from ca. 80 m a.g.l for flight 3 and ending up with ca. 200 m a.g.l. for flight 6. Above the observed convective layer, a residual layer was still present, and the potential temperature gradient within this layer remained constant.
The comparison of the results of the wind speed from the drone with observations from the JU Campus meteorological station showed good agreement for weak wind conditions (below 2 m s −1 ), while during stronger wind speed occurrence, drone measurements were in between average and maximum wind speed recorded by the JU station, being in most cases closer to the maximum values at a certain moment. Such an effect can be explained by local reduction of wind speed at the JU station location due to the influence of the surrounding campus buildings while the flights were made in the open field area.
Analysis of the period with occurrence of foehn wind in the valley allowed us to identify two opposite effects associated with elevated wind speed. On the one hand, an improvement of the air quality induced by the intense dispersion of pollution stimulated by strong turbulences associated with high wind speed in the first phase of foehn wind (wind speed below 5 m s −1 ) was observed (Figure 17a,b, flights 3-6), but on the other hand, deterioration of air quality by resuspension of deposited particles from the ground was observed in the second phase when wind speed exceeded the value of 5 m s −1 (Figure 17b,c, flights 7-11). A similar effect was reported by Kishcha et al. [54] for the Dead Sea Valley region.
An attempt was made to calculate the turbulent kinetic energy (TKE) parameter which can be very useful to identify thermal winds occurrence in the valley; however, drone position fluctuations during the flights in turbulent conditions masked recorded wind components fluctuations required for its calculation, and the obtained values of TKE were extremely noisy and did not show any pattern along the vertical profiles. Future studies are required to test if the analysis of drone fluctuations itself can be useful for determination of the TKE parameter.

Conclusions
The availability of high-quality atmospheric measurements by using a UAV provides unquestionable value to meteorological studies. However, attention should be paid to the time resolution of the equipment used to measure meteorological (air temperature, wind components, air pressure and air humidity) and air quality parameters (e.g., PM 10 concertation, amount of particulate matter at selected size). The recommended time resolution of meteorological parameters should be 1s or less [30]. Recent studies indicated that the location of the measuring equipment on the drone is very significant; it can cause two significant problems: bias measurement related to exposure to external radiation sources (e.g., insufficient solar shielding) and heat (e.g., engines or electronics) [30,55]. Position of the sensors on the drone should provide proper ventilation undisturbed by drone propeller movement (quick sensor response to changing environmental conditions).
Several other studies that investigated the effects of a multicopter on air sampling pointed out similar conclusions. The results of computational fluid dynamics (CFD) simulation for the quadrotor with added sampler presented in the article by McKinney et al. 2019 [56] and field studies of the flow of a four-rotor drone using colored pyrotechnical smoke cartridges [57] confirm that the largest airflow disturbances exist at the propellers (especially below the propellers) and the smallest in the middle of the drone above the propeller position. Hence the optimal position of the sensors is in the middle part of a UAV and elevated above the propeller height.
The smaller size of a UAV makes operation much simpler and easier but may reduce the resistance capability of strong winds and increase the instability of measurements significantly [32]. It was confirmed that fluctuations of drone position influence wind speed variability making reliable calculation of, e.g., the TKE parameter, impossible. Additionally, tests of low-cost air pollution sensors indicated influence of air humidity on the measurement of PM 10 concentration (e.g., sensor Mie pDR-1500 [32]). It was demonstrated that equipping the sensor with a heated inlet solves this problem.
Atmospheric measurements by using a UAV can be applied in high-resolution research, verification of modeling results [58] and as an additional source for data assimilation in the lower troposphere [59]. Recent study indicates that UAVs that are able to accurately measure three-dimensional wind might be used as a cost-effective and flexible addition to measurement masts and LIDAR scans [60].
The atmospheric vertical profile results obtained during three campaigns using an airborne measurement system demonstrate the usefulness of this tool in investigating the dynamics of the lower part of the PBL. The results present three different pictures of the transition phase between the nocturnal boundary layer and the convective layer and influence of the foehn wind on the PBL structure in the city located inside the river valley. Campaign 1 showed the transition from convective to stable stratified atmosphere. Campaign 2 showed vertical atmosphere dynamics during the opposite process. The last campaign was dedicated to studying the influence of foehn wind formation on urban PBL structure.
Concentrations of the PM 10 fraction in air pollution that are associated with temperature and relative humidity measurements along the vertical profile help to better understand meteorological processes influencing the formation of high PM 10 concentration. Results of experiments during the night indicate that development of a stable nocturnal boundary layer affected higher PM 10 concentrations at the lowest part of the boundary layer. After sunrise, the PM 10 concentration above the ground decreases because of convection that dilutes pollutants in the larger air volume and transports them to the higher levels. The last case showed that, during strong wind periods, a long-distance pollutant transport can significantly influence the air quality.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.