Estimation of Sensible Heat Flux and Atmospheric Boundary Layer Height Using an Unmanned Aerial Vehicle

Abstract

In this work, sensible heat flux estimated using a bulk transfer method was validated with a three-dimensional ultrasonic anemometer or surface layer scintillometer at various sites. Results indicate that it remains challenging to obtain temperature and wind speed at an appropriate reference height. To overcome this, alternative observations using an unmanned aerial vehicle (UAV) were considered. UAV-based wind speed and sensible heat flux were indirectly estimated and atmospheric boundary layer (ABL) height was then derived using the sensible heat flux data. UAV-observed air temperature was measured by attaching a temperature sensor 40 cm above the rotary-wing of the UAV, and UAV-based wind speed was estimated using attitude data (pitch, roll, and yaw angles) recorded using the UAV’s inertial measurement unit. UAV-based wind speed was close to the automatic weather system-observed wind speed, within an error range of approximately 10%. UAV-based sensible heat flux estimated from the bulk transfer method corresponded with sensible heat flux determined using the eddy correlation method, within an error of approximately 20%. A linear relationship was observed between the normalized UAV-based sensible heat flux and radiosonde-based normalized ABL height.

1. Introduction

Indirect methods of estimating sensible heat flux from measurable mean wind speed and temperature in the surface layer or the whole atmospheric boundary layer (ABL) are based on appropriate flux–profile relationships. The simplest and most widely used method is the bulk aerodynamic approach [1]. However, it is not easy to obtain meteorological variables at an assumed standard height of 10 m. The ABL height is usually determined using vertical profiles of specific humidity, potential temperature, and CO₂ mixing ratios from atmospheric soundings [2,3]. New data sources (e.g., radiometers, wind profilers, lidars, sodars, tethered balloons, and aircraft) have resulted in the development of additional methods to understand ABL structure and processes [4,5]. While these data sources can acquire stable vertical profiles of temperature and wind speed, there are challenges regarding their widespread application. For example, the use of expensive instrumentation is not always economically viable, requiring skilled handling and labor-intensive operations.

Since the 1970s, interest in unmanned aerial vehicles (UAV) for atmospheric research into convective processes [6], ABL observations [7], surveillance, and weather forecasting [8] has grown, especially in situations that are too dangerous or remote for manned aircraft. In a previous study, a UAV equipped with shortwave and longwave radiometers was used to determine radiometric flux [9]. More recently, UAVs have been equipped with various weather observation sensors to obtain quasi-continuous vertical and horizontal data, including temperature, humidity, wind speed, wind direction, and aerosol concentration. Although atmospheric science studies mostly use fixed-wing UAVs [10,11,12,13], Shimura et al. [14] and Palomaki et al. [15] have confirmed the possibility of measuring vertical profiles of temperature, humidity, and wind vectors within 1 km using an ultrasonic anemometer installed on a small hexacopter UAV. UAVs are more flexible in areas they can cover, unlike meteorological masts or towers.

UAV-based wind vectors have been studied since early 2000 [13,16,17,18,19]. For example, Van den Kroonenberg et al. [13] used a five-hole wind probe on a small UAV to measure mean wind vectors; however, their approach was based on aerodynamics that are not always available and require an additional airspeed sensor. Palomaki et al. [15] estimated atmospheric winds from a UAV-mounted anemometer with a direct approach and an indirect approach based on inertial measurement unit (IMU) data. They reported that both approaches had a root mean-squared error (RMSE) of about 0.5 m s⁻¹ when compared with independent wind measurements in calm to moderate wind conditions (0–5 m s⁻¹). However, Brosy et al. [20] have reported that a UAV-mounted anemometer cannot capture the full range of wind vectors because the volume of the rotary-wing UAV is relatively larger than the ultrasonic anemometer. Neumann and Bartholomai [21] estimated wind vectors by means of a wind triangle (navigation equation) using attitude data (pitch, roll, and yaw angles) stored in an IMU sensor in a rotary-wing UAV and GPS data (ground speed and ground direction). Brosy et al. [20] have suggested a method similar to that of Neumann and Bartholomai [21] but without wind tunnel tests. UAV-based wind estimations are also an alternative to ultra-high frequency wind profiler radar that are unable to measure wind vectors within 100–200 m.

Sensible and latent heat flux were estimated using data from the UAV small multifunction research and teaching sonde (SMARTSonde) [22] and showed that UAV-based meteorological data can be used to examine how the early evening transition affects the thermodynamic structure of the lower ABL. As suggested by Reineman et al. [23], UAVs are useful for direct measurements of sensible, latent, and momentum flux within the ABL while simultaneously detecting surface topography. Atmospheric turbulent heat fluxes were estimated over Terra Nova Bay in September 2009 using UAV-observed temperature, wind speed, and relative humidity by the integral method [24]. They highlighted sensitivities in UAV observations that were important to the methodology used to estimate the turbulent flux. Båserud et al. [25] used the micro-remotely piloted aircraft system (RPAS) small unmanned meteorological observer (SUMO) to investigate turbulent kinetic energy during the boundary-layer late afternoon and sunset turbulence (BLLAST) field campaign in 2011. They reported limitations in capturing the largest turbulent scales and turbulent production scales.

Indirect methods of estimating sensible heat flux from UAV-observed temperature and UAV-based wind speed have advantages. We suggest an indirect approach to estimate wind speed, sensible heat flux, and ABL height by applying UAV-based IMU data observed from flight, UAV-observed temperature, and surface measurement, which is more economic and less labor-intensive compared to methods of earlier studies. It has the advantage of being able to determine the horizontal variability of sensible heat flux at a regional scale and low-altitude flights at no risk to human pilots and can capture the evolution of ABL over transition zones such as air-sea boundaries. This paper is organized as follows. Section 2 describes the methods for estimating sensible heat flux using a bulk transfer method. Section 3 evaluates UAV-based estimates of wind speed using UAV attitude data and UAV-based estimates of sensible heat flux and considers the relationship between the latter and radiosonde-observed ABL height. The conclusions made from these evaluations are given in Section 4.

2. Methods

2.1. Bulk Method

When using the sensible heat flux H (W m⁻²) in practical applications, it can be measured directly with data from eddy covariance equipment and Equation (1a). It is often related to the mean atmospheric variables via so-called bulk transfer parameterization, i.e.,

$$\begin{array}{ccc}\hfill \mathrm{H}={\mathsf{\rho }\mathrm{c}}_{\mathrm{P}}\overline{{\mathsf{\theta }}^{\prime }{\mathrm{w}}^{\prime }}& \simeq & {\mathsf{\rho }\mathrm{c}}_{\mathrm{P}}{\mathrm{C}}_{\mathrm{H}}{\mathrm{u}}_{\mathrm{r}}\left({\mathsf{\theta }}_0-{\mathsf{\theta }}_{\mathrm{r}}\right)\hfill \\ \hfill \left(\mathsf{a}\right)& & \left(\mathsf{b}\right)\hfill \end{array}$$

(1)

where $\mathsf{\rho }$ is the density of air (1.25 kg m⁻³), ${\mathrm{c}}_{\mathrm{P}}$ (1004 J K⁻¹ kg⁻¹) is the specific heat at a constant pressure, ${\mathrm{u}}_{\mathrm{r}}$ is the wind speed (m s⁻¹) measured at a reference height (10 m), $\mathsf{\theta }$ is the potential temperature, ${\mathsf{\theta }}_0-{\mathsf{\theta }}_{\mathrm{r}}$ is the difference in potential temperature between the surface and the air, and C_H is the dimensionless exchange coefficient for sensible heat. The latter is often given as a function of wind velocity [26] and particularly on a surface, such as an ocean or a lake. However, atmospheric stability has an important influence on land surface. The heat transfer coefficient C_H needs to be determined before the surface sensible heat flux can be estimated using Equation (1b). The well-known Monin–Obukhov similarity theory is mostly used to close the bulk transfer parameterization [1,27,28,29], i.e.,

$$C_H=k^2{\left[ln\left(\frac{z_r-d_0}{z_{0m}}\right)-{\mathsf{\psi }}_m\left(\frac{z_r-d_0}{L}\right)\right]}^{-1}{\left[ln\left(\frac{z_r-d_0}{z_{0m}}\right)-{\mathsf{\psi }}_h\left(\frac{z_r-d_0}{L}\right)\right]}^{-1}$$

(2)

where $\mathrm{k}$ is the Von Karman constant (0.4) and ${\mathrm{d}}_0$ is the zero-plane displacement as a function of vegetation height (h₀). Zero-plane displacement is determined by Equation (3) [30], ${\mathrm{z}}_{\mathrm{r}}$ is the reference measurement height, and ${\mathrm{z}}_{0\mathrm{m}}$ is the roughness length, which is calculated with Equation (4) [31] using the heat transfer coefficient C_H proposed by Verkaik and Holtslag [32].

$$ln\left(d_0\right)=0.98\times ln{\left(h_0\right)}^2-0.15$$

(3)

$$z_{0m}=\frac{1}{n}{\displaystyle \sum }_{i=1}^{n}z exp{\left[\frac{0.4\times u_r}{u_{*}}+{\mathsf{\psi }}_m\left(\frac{z_r-d_0}{L}\right)\right]}^{-1}$$

(4)

${\mathsf{\psi }}_m$ and ${\mathsf{\psi }}_h$ are different similarity functions related to ${\mathsf{\varphi }}_m$ and ${\mathsf{\varphi }}_h$, which are the dimensionless functions of atmospheric stability parameters for momentum and sensible heat flux, respectively, and can be expressed as, depending on the atmospheric stability [33],

$$\frac{z_r}{L}\ge 0,{\mathsf{\psi }}_m={\mathsf{\psi }}_h=-\frac{5\left(z_r-d_0\right)}{L}$$

(5)

$$\frac{z_r}{L}<0,{\mathsf{\psi }}_m=ln\left[\left(\frac{1+x^2}{2}\right){\left(\frac{1+x}{2}\right)}^2\right]-2ta{n}^{-1}x+\frac{\pi }{2}$$

(6)

$$\frac{z_r}{L}<0,{\mathsf{\psi }}_h=2ln$$

where x is ${\left(1-15\left(\frac{z_r-d_0}{L}\right)\right)}^{0.25}$.

$$u_{*}=\frac{k\left(u_r-u_0\right)}{ln\left(\frac{z_r-d_0}{z_{0m}}\right)-{\mathsf{\psi }}_m\left(\frac{z_r-d_0}{L}\right)}$$

(7)

The Monin–Obukhov atmospheric stability parameter $z_r/L$ indicates the relative significance of buoyancy versus shear effects. The Obukhov length (L) can be expressed as Equation (8), i.e.,

$$L=-\frac{\mathsf{\rho }{c}_p{u}_{*}^3\overline{\mathsf{\theta }}}{kgH}$$

(8)

where g is the gravitational acceleration (m s⁻²), $\overline{\mathsf{\theta }}$ is the mean potential temperature (K) between two heights, and ${\mathrm{u}}_{*}$ is the friction velocity (m s⁻¹). The friction velocity and sensible heat flux must be estimated to determine L and vice versa. Therefore, L was calculated using an iterative method. ${\mathsf{\psi }}_{\mathrm{m}}$ and ${\mathsf{\psi }}_{\mathrm{h}}$ will equal 0 when the initial value of L is assumed to be the neutral atmospheric stability (=10¹⁵). Based on these variables, the friction velocity, sensible heat transfer coefficient, and sensible heat flux can be calculated in an iterative manner using Equations (1) to (7) and be substituted into Equation (8) to obtain a new value of L under the condition that the ratio $\left({\mathrm{L}}_{\mathrm{n}}-{\mathrm{L}}_{\mathrm{n}-1}\right)/{\mathrm{L}}_{\mathrm{n}-1}$ is less than 1% [34].

The sensible heat flux from the bulk transfer method was validated with direct flux measurement at various sites (

Table 1. Observation periods and characteristics for each site.

Site	Period	Surface	d0 (m)	z0m (m)	Synoptic Condition	Main Instruments
GH1	2006.08.08–2006.08.11	Reed	2.1	2.8 × 10−1	North pacific anticyclone	USA
GH2	2005.02.23–2005.02.25	Reed	1.7	1.0 × 10−1	Siberian anticyclone	USA
GR1	2012.10.08–2012.10.11	Rice	0.6	8.5 × 10−2	Migratory anticyclone	USA
GR2	2012.10.14–2012.10.16	Cut paddy	0.01	5.6 × 10−2	Migratory anticyclone	USA
KHR	2012.10.01–2012.10.03	Water	0	1.0 × 10−3	Migratory anticyclone	USA, SLS20
SYR	2013.03.14–2013.03.16	Water	0	1.0 × 10−3	Migratory anticyclone	USA, SLS20
CHR	2016.07.09–2016.07.10	Asphalt	0	3.0 × 10−2	North pacific anticyclone	RS, USA, UAV
WON	2016.07.11–2016.07.14	Stone	0	0.2 × 10−2	North pacific anticyclone	RS, USA, UAV
SOK	2017.02.09–2017.03.15	Stone	0	2.0 × 10−2	Migratory anticyclone	RS, USA, UAV
BOS	2018.09.10–2018.09.12	Grass	0.3	1.2 × 10−2	Migratory anticyclone	RS, USA, UAV
YHI	2018.11.27–2018.12.05	Trees	2.0	1.7 × 10−1	Siberian anticyclone	RS, USA, UAV

Legend: RS, radiosonde; USA: 3D ultrasonic anemometer; UAV, unmanned aerial vehicle.

). The results were divided into three categories according to differences in the observation data.

The air temperature, surface skin temperature, and wind speed were based on data collected at six observation sites using an automatic weather system (AWS) (Table 1). Sites GH1 and GH2 refer to the same reed field (36°36′45″ N, 127°13′05″ E) in Goheung Bay, Jeollanam-do, Korea, and are classified as summer and winter, respectively. GH1 consists of green reeds, approximately 2.5 m tall, and is characterized by a reed cover density higher than that of GH2. Site GH1 extends approximately 10 km north to south and 3 km east to west, with no other terrain or buildings in the surrounding area. At GH2, the reeds were dry, brown, sparse, and few in number compared with GH1. Measurements were collected from a 10 m meteorological tower equipped with a three-dimensional ultrasonic anemometer (CSAT3A, Campbell Sci., Logan, UT, USA), an aerovane (05103, R. M. Young, Traverse City, MI, USA), and a thermo-hygrometer (HMP45C, Campbell Sci., Logan, UT, USA). GR1 and GR2 refer to the same agricultural field (35°49′50.77″ N, 128°27′33.72″ E) at Goryeong-gun, Gyeongsangbuk-do, Korea, and were classified as before and after harvest, respectively. GR1 is a rice field, with crop heights of approximately 70 cm, while GR2 is a rough field that yields rice harvests. A 10 m meteorological tower operated by the Korea Meteorological Administration (KMA) near the field site is equipped with a three-dimensional ultrasonic anemometer (CSAT3A, Campbell Sci., Logan, UT, USA), an aerovane (05103, R. M. Young, Traverse City, MI, USA) and a thermo-hygrometer (HMP45C, Campbell Sci., Logan, UT, USA). Sensible heat fluxes were measured simultaneously at both sites at a height of approximately 2 m using a three-dimensional ultrasonic anemometer (SATI-3K, Applied Technologies, Inc., Longmont, CO, USA) and a surface layer scintillometer (SLS20, Scintec, Rottenburg am Neckar, Germany) to facilitate comparisons between the different levels. The heat fluxes were in good agreement, with a mean bias of 2.87 W m⁻² and RMSE of 24.3 W m⁻². The sensible heat fluxes based on the anemometer and surface layer scintillometer were compared and yielded reasonable agreement (i.e., a mean bias of 19.8 W m⁻² and RMSE of 34.9 W m⁻²). This confirmed the demonstrated reliability of the SLS20 scintillometer when performing measurements of various surfaces [35,36,37,38].

Site KHR (130 m wide; located at 35°51′00″ N, 128°28′13″ E) is part of the Kumho River, which flows southwest of Daegu, Korea, and is approximately 2.5 km from sites GR1 and GR2. Site SYR (90 m wide; located at 35°11′26″ N, 129°06′52″ E) is part of the Suyeong River, which flows south into the city of Busan. Air temperature and wind speed were measured at 10 m using a weather transmitter sensor (WXT520, Vaisala, Vantaa, Finland). Water temperature was measured using a water temperature sensor (LTC Levelogger, Solinst, Georgetown, Canada) and the sensible heat flux across the width of the water surface was measured using a SLS20 scintillometer, which calculates turbulence parameters using variations in the refractive index of light based on changes in temperature, air density, momentum, and moisture between the receiver and transmitter.

We then estimated the sensible heat flux using the bulk transfer method with UAV-observed air temperature and AWS-observed wind speed and surface skin temperature data at three KMA observation sites: Cherwon (CHR; 38°08′53″ N, 127°18′16″ E) is surrounded by farmland, Wonju (WON; 37°20′15″ N, 127°56′48″ E) is a complex area including residential buildings, and Sokcho (SOK; 38°15′02″ N, 128°33′52″ E) is located in the East sea, approximately 30 m southwest off the coast of Korea. Anticyclonic conditions were dominant during the UAV-based observation periods (Table 1).

Finally, we estimated the UAV-based wind speed and evaluated it with ground observations. We estimated the sensible heat flux using the UAV-based wind speed and UAV-observed air temperature and observed surface skin temperature at two sites: Boseong (BOS; 34°45′48″ N, 127°12′53″ E) is a standard weather observatory operated by the KMA that conforms to international standards set by the World Meteorological Organization; site BOS was designated as a leading observation station in January 2012, with the site consisting of a horizontal homogeneous agricultural area (154,500 m²) which is in contact with the sea to the south.

Site YHI (37°15′50″ N, 126°29′44″ E) is located along the southern coast of Youngheung Island, near the city of Incheon (Table 1). Site YHI has an area of 23.46 km² and a 42.2 km coastline. A meteorological tower was installed and equipped with a weather transmitter (WXT520, Vaisala, Vantaa, Finland), surface temperature sensor, and three-dimensional ultrasonic anemometer.

2.2. Unmanned-Aerial-Vehicle-Based Wind Speed

The UAV used in this study (Figure 1a) has four rotors and is capable of stable flight for up to 30 min with a fully charged battery. The UAV frame (Tarot X4, Wenzhou Aviation Technology Co., Wenzhou, China) is 960 mm × 960 mm × 320 mm with a weight of 1.6 kg. It was equipped with a Pixhawk 2 (3DR, USA) autopilot consisting of an external GPS (HERE + RTK GNSS, 3DR, Berkeley, CA, USA), a three-axis accelerometer, a three-axis gyroscope, a three-axis magnetometer, a temperature sensor, and a barometer to ensure stable position control. The IMU, which consisted of a three-axis accelerometer, and gyroscope automatically measured the UAV’s position data during the flight. Kalman filter calibrations were carried out in flight to decrease errors. The IMU and flight measurement unit (FMU) systems were separated to reduce sensor interference and to filter out high-frequency vibrations, reducing errors in the IMU data. The attitude data were recorded at 100 Hz for flight safety and GPS data (ground speed (${\mathrm{s}}_{\mathrm{g}}$) and ground direction (${\mathsf{\theta }}_{\mathrm{g}}$)) were collected at 50 Hz. Flight stability was achieved using open source Mission Planner software with a horizontal hovering accuracy of ±0.01 m. UAV maximum elevation, horizontal movement, wind speed, and horizonal and vertical flight speeds were 2,000 m, 1,000 m, 10 m s⁻¹, 6 m s⁻¹, and 5 m s⁻¹, respectively. To measure the air temperature, a thermo-hygrometer (HOBO H8 Pro, ONSET, Wareham, MA, USA) was installed on the rotary-wing of the UAV at 40 cm above the rotors. The rotary-wing UAV moves through the air by setting a tilt angle that is roughly proportional to speed. This tilt angle varies in real time to maintain stable flight. Therefore, the UAV-based wind vector can be indirectly estimated using the pitch, roll, and yaw angles measured by onboard sensors (IMU) for UAV’s attitude control without additional sensors for wind vector measurements [13,14,15,16,17,18,19,20,21].

The UAV-based wind vector was able to be calculated using the wind triangle (Figure 1b), which is widely used in navigation and estimates the wind vector using a flight and ground vector. The ground vector is measured directly from the UAV GPS receiver while the flight vector can be measured using an optical flow sensor or pitot tube [13,17,19]; however, these methods are not suitable for rotary-wing UAVs due to their flight speed (${\mathrm{s}}_{\mathrm{f}}$) (which is lower than fixed-wing UAV), inconsistent flight direction (${\mathsf{\theta }}_{\mathrm{f}}$), and a wide range of tilt angle changes except for an examination [39].

Neuman and Bartholomai [21] have found a relationship between the tilt angle and flight speed ($s_f$) in wind tunnel experiments while Brosy et al. [20] have estimated flight speed assuming ground speed should be equal to flight speed under calm wind conditions where the wind speed is negligible (less than 1 m s⁻¹) and have developed a relationship between the tilt angle and flight speed. We developed a relationship between the tilt angle and flight speed using the method given in Brosy et al. [20]. This is explained in detail in Section 3.2.

$${\overrightarrow{e}}_{pitch}=\left(\begin{array}{c}cos\left(pitch\right)\\ 0\\ -sin\left(pitch\right)\end{array}\right),{\overrightarrow{e}}_{roll}=\left(\begin{array}{c}0\\ \cos \left(roll\right)\\ \sin \left(roll\right)\end{array}\right)$$

(9)

The tilt angle ($\mathsf{\alpha }$) can be estimated from the inverse scalar product of the normal unit vector ${\overrightarrow{\mathrm{n}}}_{\mathrm{XY}}=\left(0,0,1\right)$ to the XY-plane which is parallel to the ground and the cross product of the unit vectors as shown in Equation (10).

$$\mathsf{\alpha }=co{s}^{-1}\left(\frac{{\overrightarrow{n}}_{XY}·\left({\overrightarrow{e}}_{pitch}\times {\overrightarrow{e}}_{roll}\right)}{\left|{\overrightarrow{n}}_{XY}\right|·\left|{\overrightarrow{e}}_{pitch}\times {\overrightarrow{e}}_{roll}\right|}\right)$$

(10)

The angle $\mathsf{\beta }$ between the projection of the vector ${\overrightarrow{\mathrm{e}}}_{\mathrm{pitch}}\times {\overrightarrow{\mathrm{e}}}_{\mathrm{roll}}$ onto the XY-plane, and the viewing direction of the UAV, defined as the negative normal vector $-{\overrightarrow{\mathrm{n}}}_{\mathrm{YZ}}=\left(-1,0,0\right)$ is calculated using Equation (11).

$$\mathsf{\beta }=co{s}^{-1}\left(\frac{-{\overrightarrow{n}}_{Yz}·{\left({\overrightarrow{e}}_{pitch}\times {\overrightarrow{e}}_{roll}\right)}_{XY}}{-\left|{\overrightarrow{n}}_{YZ}\right|·\left|{\left({\overrightarrow{e}}_{pitch}\times {\overrightarrow{e}}_{roll}\right)}_{XY}\right|}\right)$$

(11)

To calculate the flight direction (${\mathsf{\theta }}_{\mathrm{f}}$) whether the orthogonal vector ${\overrightarrow{\mathrm{e}}}_{\mathrm{pitch}}\times {\overrightarrow{\mathrm{e}}}_{\mathrm{roll}}$ is located in the left or right direction of the UAV’s viewing direction $-{\overrightarrow{\mathrm{n}}}_{\mathrm{YZ}}=\left(-1,0,0\right)$, Equation (12) is used. The flight direction can be estimated using $\mathsf{\beta }$ and the compass angle (${\mathsf{\delta }}_{\mathrm{c}}$) of the UAV viewing direction using Equation (13).

$${\overline{n}}_{XZ}·{\left({\overrightarrow{e}}_{pitch}\times {\overrightarrow{e}}_{roll}\right)}_{XY}=\{\begin{array}{cc}<0,\hfill & \mathrm{if}{\overrightarrow{e}}_{\mathrm{pitch}}\times {\overrightarrow{e}}_{\mathrm{roll}}\mathrm{is}\mathrm{left}\hfill \\ >0,\hfill & \mathrm{if}{\overrightarrow{e}}_{\mathrm{pitch}}\times {\overrightarrow{e}}_{\mathrm{roll}}\mathrm{is}\mathrm{right}\hfill \\ =0,\hfill & \mathrm{otherwise}\end{array}$$

(12)

$${\mathsf{\theta }}_f=\{\begin{array}{c}360°-\mathsf{\beta }+{\mathsf{\delta }}_c,\mathrm{if}\mathrm{Equation}\left(11\right)<0\\ \mathsf{\beta }+{\mathsf{\delta }}_c,\mathrm{otherwise}\end{array}$$

(13)

Finally, the UAV-based wind speed (${\mathrm{s}}_{\mathrm{w}}$) is calculated using the wind triangle (Equation (14)) with the estimated flight vector (${\mathrm{s}}_{\mathrm{f}}$, ${\mathsf{\theta }}_{\mathrm{f}}$) and measured ground vector (${\mathrm{s}}_{\mathrm{g}}$, ${\mathsf{\theta }}_{\mathrm{g}}$).

$$s_w=\sqrt{s_f^2+s_g^2-2s_f·s_g·cos\gamma }$$

(14)

The drift angle, $\mathsf{\gamma }$, is equal to the difference between the ground direction (${\mathsf{\theta }}_{\mathrm{g}}$), which is measured by the GPS receiver and flight direction (${\mathsf{\theta }}_{\mathrm{f}}$). For detailed explanations of the wind direction (${\mathsf{\theta }}_{\mathrm{w}}$) calculation, refer to Neuman and Bartholomai [21].

3. Results and Discussion

3.1. Sensible Heat Flux

Weather conditions at the study sites were dominated by anticyclones during the observation period (Table 1). In vegetated areas (GH1, GH2, GR1, and GR2), the mean wind speed was 2.0 to 2.8 m s⁻¹ during the day, which weakened overnight. On rivers (KHR and SYR), weak winds (<2 m s⁻¹) were observed during the day and night. For all sites, the heat flux and heat transfer coefficients were calculated indirectly using the bulk transfer method and validated using anemometer and SLS20 scintillometer measurements. As vertical fluctuations of horizontal wind and the thermal instability due to buoyancy vary with time, heat transfer coefficients clearly vary in response to temporal changes in atmospheric stability (Figure 2a). In vegetated areas, the heat transfer coefficient varied depending on atmospheric stability and roughness length z_0m (m) from 6.0 × 10⁻³–3.0 × 10⁻² during the day to 2.0 × 10⁻³–10.0 × 10⁻³ overnight, which is consistent with values summarized according to surface conditions by Stull [40]. The heat transfer coefficient was larger in summer compared to winter. This is because the temperature difference between two heights is higher in summer than in winter. When the atmospheric stability was unstable ((z_r−d₀)/L < 0), the differences between heat transfer coefficients was larger than when the atmospheric stability was stable. However, the atmospheric stability was always unstable at river sites, because the water temperature was higher than the air temperature, even overnight. Consequently, heat transfer coefficients did not show significant differences between day and night.

Sensible heat flux was estimated using a heat transfer coefficient depending on the atmospheric stability ((z_r−d₀)/L), zero-plane displacement d₀, and roughness length z_0m, and compared with measurements from the 3D ultrasonic anemometer and SLS20. Figure 2b shows good agreement between directly and indirectly derived values, with a correlation coefficient of 0.94, mean bias of −1.26 W m⁻², and RMSE of 19.9 W m⁻², which is less than approximately 10% of the maximum measured value. We confirm that the bulk transfer method is useful to estimate UAV-based sensible heat flux.

3.2. Temperature and Wind Speed from Rotary-Wing UAV

The measurements from the UAV-mounted sensors (i.e., the anemometer and thermo-hygrometer) required validation because they were acquired from a moving body and were therefore liable to be affected by rotor wash. Palomaki et al. [15] have examined the effect of rotors on the wind vector from an anemometer located 30 cm above the rotary-wing UAV and found that the average bias due to the rotors was 0.5 m s⁻¹ in indoor tests; however, rotor disturbance was found to be negligible at a height of 40–47 cm from the top of the rotary-wing UAV [14,41].

For the first part of the validation process, the profiles of temperature and humidity measured at 40 cm above the rotary-wing UAV were compared with those from the 10 m high meteorological tower that was equipped with an aerovane, a thermo-hygrometer, and surface temperature sensor (KWT1002, Wellbian system, Seoul, Korea) at sites CHR, WON, and SOK.

Figure 3 compares the air temperature and relative humidity between UAV-observed and AWS-observed data at a height of 10 m while hovering at a distance of 15 m from the AWS. UAV-observed temperatures are seen to be in good agreement with the AWS-observed temperatures, with a correlation coefficient of 0.98, mean bias of −1.1 °C, and RMSE of 1.79 °C. UAV-observed relative humidity measurements had lower agreement, with a correlation coefficient of 0.9, mean bias of 6.66%, and RMSE of 9.86%. AWS-observed relative humidity tended to be higher than the UAV-observed relative humidity during the day and lower during the night. The difference in agreement could be attributed to the accuracy of the temperature sensor (±0.2 °C) being higher than the accuracy of the relative humidity sensor (±2.5%). The UAV-observed atmospheric data was measured based on the volume of the air around the UAV rather than observations at a single point. Increasing the response time of the sensor could overcome this discrepancy [20].

Radiosondes (DFM-09, Graw, Nürnberg, Germany; M10, Meteomodem, Ury, France) and a rotary-wing UAV were deployed simultaneously. Radiosondes provide fine resolution (one second) profiles of barometric pressure, temperature, wind, and relative humidity. Figure 4a–c shows vertical profiles of radiosonde-observed air temperature and vertical profiles of UAV-observed air temperature at CHR, WON, and SOK, respectively. Although radiosonde data loss and error were confirmed near the ground, there were small difference patterns between the profiles at nighttime. However, there were large difference patterns between the profiles during the daytime. These large difference patterns could be due to the time difference between the profiles. The radiosonde-observed air temperature had a temporal resolution of 1 s with an ascent rate of more than 6 m s⁻¹, but UAV-observed temperatures were recorded while hovering at the same height for 1 min. When the surface is heated by solar radiation, the linearity of the air temperature gradient is not maintained. Martin et al. [11] have reported that it is possible to measure temperature, humidity, wind direction, wind speed, and even turbulence flux at a 10 cm vertical resolution using fixed-wing UAVs equipped with a very fast response speed sensor (30 Hz). These measurements were compared with measurements from the weather tower, wind profiler, and radiometer for verification. Because these observations were difficult to summarize from existing standard observation sensors, they reported only the usefulness and utility of the rotary-wing UAV. We examined the estimated surface skin temperature using linear extrapolation of the profile of UAV-observed air temperature and compared it with the observed surface temperature sensor (KWT1002, Wellbian system, Seoul, Korea) at CHR and WON. The correlation coefficient, mean bias, and RMSE were approximately 0.9, 0.2 °C, and 3.3 °C, respectively (Figure 4d); the latter shows that the largest difference between the temperatures was 10.5 °C at 10:00 local solar time (LST) on July 10, 2016. This is likely because of the specific heat difference between air and land surface and the land characteristics. After sunrise, land surface temperature usually begins to increase faster than air temperature due to the impact of solar radiation. The difference between surface temperatures at CHR is greater than those at WON because the land surface at the latter comprises stone and sand but that at the former comprises asphalt. Katz and Zhu [28] state that the surface skin temperature should be that defined at the height of thermal roughness length and that it may change from time to time even at a fixed location. A slight change in the height of the sensors can result in a large bias in the skin temperature because there is a sharp temperature gradient close to the surface. Thus, a proper extrapolation method to estimate surface temperature needs to be developed and could be the focus of further studies. In this study, an observed surface temperature was used for estimating UAV-based sensible heat flux from the bulk transfer method.

UAV-based wind speed was estimated by applying the flight vector (${\mathrm{s}}_{\mathrm{f}},{\mathsf{\theta }}_{\mathrm{f}}$) and the ground vector (${\mathrm{s}}_{\mathrm{g}},{\mathsf{\theta }}_{\mathrm{g}}$) to the wind triangle (Equation (14)) for the rotary-wing UAV at BOS and YHI. Notably, BOS has the second highest meteorological tower (307 m) in Asia (Figure 5). The accuracy of UAV-based wind speed depends on the accuracies of the flight and ground vectors. The ground vector is measured directly by the UAV’s GPS receiver. Therefore, a precise GPS receiver is very important for estimating UAV-based wind speed and flight performance. In this study, GPS equipped with real time kinematic (RTK) technology was installed to reduce error between the travel route and actual travel route measured by the GPS receiver to 3 cm from approximately 30 cm. For time synchronization between the GPS data (ground speed and ground direction) measured at 50 Hz and attitude data measured at 100 Hz, attitude data closest to the time of the GPS data were selected. The tilt angle was calculated using synchronized pitch and roll angles (50 Hz).

The ground speeds (50 Hz) recorded in the RTK GPS receiver were replaced with the flight speeds (${\mathrm{s}}_{\mathrm{f}})$. The tilt angle and flight speed were each averaged for one minute and then classified into 0.5 m s⁻¹ intervals in this study. The tilt angle was also classified according to flight speed class and the number of flight speed classes was different for each case. Flight speeds from the highest class to the third highest class were averaged to determine the representative value of a case. Tilt angles were averaged in a similar manner.

Figure 6 shows scatter plots of the calculated tilt angles using Equation (10) and flight speeds that were replaced by ground speeds recorded in the RTK GPS receiver for 73 cases of calm winds (<3 knots) at BOS and YHI. The relationship between the tilt angle and the flight speed was estimated using the least squares method (Equation (15)) and is depicted by the solid black line in Figure 6. The difference in the relationships is possible due to hover characteristics. These are dependent on the diameter and weight of the UAV frame, tuning parameters that differ between UAV designs, and the quality of the GPS signal. Our relationship considers low flight speeds (<2 m s⁻¹) and tilt angles (<4°).

$$flightspeed\left(s_f\right)=0.0372{\mathsf{\alpha }}^{2.856}$$

(15)

For another 54 cases with wind (≥3 knots) at BOS and YHI, the tilt angles were calculated using pitch and roll angles measured during hovering at a distance of 15 m (for YHI this was 10 m) from the meteorological tower using Equation (10), and the flight speeds were calculated using the relationship we present in Equation (15). The calculated flight speeds and ground speeds were applied to Equation (14) to estimate the UAV-based wind speeds, which were then compared to those from the AWS at heights of 10, 20, and 40 m at BOS and 10 m at YHI.

Figure 7 shows the UAV-based wind speeds at BOS and YHI and the five-minute moving average AWS-observed wind speed. The correlation, mean bias, and RMSE between the UAV-based and AWS-observed wind speeds were 0.66, −0.15 m s⁻¹, and 0.4 m s⁻¹, respectively, at BOS, and 0.88, 0.34 m s⁻¹, and 0.72 m s⁻¹, respectively, at YHI. The relationship between the tilt angle and the flight speed was found to yield a reasonable flight speed for estimating UAV-based wind speed when compared to the AWS-observed wind speed. Examining the largest mean bias, the RMSEs between the UAV-based and AWS-observed wind speeds for each time period were −0.04 m s⁻¹ and 0.7 m s⁻¹, respectively, occurring at approximately 11:00–13:00 LST. The mean bias and RMSE (–0.7 m s⁻¹ and 0.9 m s⁻¹, respectively) between the UAV-based and AWS-observed wind speeds at 40 m were larger compared to at 10 m (0.1 m s⁻¹ and 0.4 m s⁻¹, respectively). These errors could have been caused by many factors, including atmospheric conditions such as wind gusts (i.e., the wind speed cannot be negligible as assumed), sharp wind variation over the distance between the rotary-wing UAV and anemometer, and time discrepancies between the GPS and IMU due to delays in GPS data. Shimura et al. [14] noted an error when turbulence was strong due to increased wind shear. Neumann and Bartholomai [21] reported an RMSE of ±0.6 m s⁻¹ when hovering 5 m away from an ultrasonic anemometer and an RMSE of ±0.3 m s⁻¹ at 2 m. Palomaki et al. [15] reported an RMSE of 0.5 m s⁻¹ when the wind speed was 0–5 m s⁻¹ while Brosy et al. [20] found an RMSE of 0.7 m s⁻¹ and a maximum standard deviation of 2 m s⁻¹. Palomaki et al. [15] have stated that stable hovering is important for accurate wind speed estimation. For this, it is necessary to optimize weight, diameter, quality of GPS signal, and various tuning parameters of the UAV.

3.3. Heat Flux Based on Rotary-Wing UAV Hovering

Estimation of the sensible heat flux based on rotary-wing UAV hovering was tested at two wind speeds. The AWS-observed wind speed was tested at SOK during international collaborative experiments for the Pyeongchang 2018 Olympic Winter Games (https://weather.msfc.nasa.gov/sport/icepop2018/) in February and March 2018 while UAV-based wind speeds were tested at BOS in September 2018 and at YHI in November and December 2018. The BOS site is the largest observatory in Korea with a 307 m meteorological tower (Figure 5) equipped with a thermo-hygrometer (5628, Fluke; HMP155, Vaisala) and two-dimensional ultrasonic anemometer (UA-2D, Thies Clima, Göttingend, Germany; 05103, R. M. Young, U.S.A.) with 11 heights that have projecting rods oriented in three directions. A three-dimensional ultrasonic anemometer, infrared gas analyzer, and net radiometer are also installed at the surface (2.5 m), as well as at 60, 140, and 300 m. The heat transfer coefficient and sensible heat flux were estimated using the bulk transfer method with AWS-observed surface temperature, UAV-observed air temperature, AWS-observed wind speed, or UAV-based wind speed.

UAV-based heat fluxes (H_{b_AWSUAV}) were calculated by applying the heat transfer coefficient, AWS-observed wind speed, AWS-observed surface temperature, and UAV-observed air temperature to Equation (1b) and comparing the result with the sensible heat flux H measured directly with eddy covariance equipment and calculated according to Equation (1a) (Figure 8a). The raw data measured at 10 Hz by the anemometer were subjected to peak inspection and direction correction according to the quality control method of Vickers and Mahrt [42]. The eddy covariance method was applied to calculate heat flux following double coordinate conversion. The correlation coefficient between H and H_{b_AWSUAV} (0.77), mean bias (0.51 W m⁻²), and RMSE (30.98 W m⁻²) were within 20% of the maximum H.

UAV-based heat fluxes (H_{b_UAVUAV}) were calculated by applying the heat transfer coefficient, UAV-based wind speed, AWS-observed surface temperature, and UAV-observed air temperature to Equation (1b) and comparing the results with H. Figure 8b presents the five-minute moving averaged H_{b_UAVUAV} and Husing different symbols for time. The correlation coefficient (0.95), mean bias (−0.75 W m⁻²), and RMSE (10 W m⁻²) between H and H_{b_UAVUAV} were within 5% of the maximum H. The largest mean bias and RMSE between the H_{b_UAVUAV} and H were −3.6 W m⁻² and 18.5 W m⁻², respectively, occurring between 11:00–13:00 LST. This error can be explained by Equation (1b). We confirmed that the mean bias and RMSE between UAV-based and AWS-observed wind speeds were relatively the largest at that time, and those were −0.04 ms⁻¹ and 0.7 ms⁻¹. The mean AWS-observed wind speed was 3 ms⁻¹ between 11:00–13:00 LST, which was the strongest compared to the other times. The mean bias and RMSE were −0.9 °C and 1.2 °C between UAV-observed and AWS-observed temperatures, respectively, and were the largest compared to the other times. Knuth and Cassano [24] calculated the sensible heat flux by applying UAV-observed pressure, temperature, and wind speed to the integral method and examined the sensitivity of each UAV-observed measurement on the calculated heat flux using random perturbations for 500 repetitions. They found that the standard deviation for most flight days was within 5–15 W m⁻². Lee et al. [43] estimated the horizontal variability in heat flux surrounding meteorological towers by applying the rotary-wing UAV-observed pressure, temperature, and relative humidity to the conditional sampling technique. Comparisons between UAV-based with independent eddy covariance sensible heat flux agreed to within 5–10 W m⁻². When comparing UAV-based and large eddy simulation-modeled heat flux, differences over a few isolated locations were as large as 50–70 W m⁻². Lee et al. [43] have also reported that these discrepancies are likely caused by the advection of high thin cirrus clouds, because large eddy simulations were unable to advect the cloud cover across the sensible heat flux that was seen in the observations.

3.4. Atmospheric Boundary Layer Height Based on Rotary-Wing UAV Hovering

Reineman et al. [23] have demonstrated measuring UAV-based turbulent flux and the development of the ABL. Furthermore, Zhang et al. [44] have reported that low soil moisture, high surface sensible heat flux, high Bowen ratio, and low latent heat flux generally facilitate ABL height growth in the afternoon. Low soil moisture tends to favor ABL growth by partitioning more solar radiation to sensible heat flux [45]. Seidel et al. [46] calculated ABL height using seven common methods and found large uncertainties with respect to these methods. Turbulence is generally responsible for the mixing processes in the ABL, and the turbulent motions affect the vertical redistribution of momentum and heat. Regardless of the synoptic condition, the ABL height is mostly determined by the activation of atmospheric turbulence by the supply of thermal energy from the surface to the atmosphere in areas where the heat flux is strong. Therefore, accurate sensible heat flux measurements are important for understanding ABL structural processes. The most widely used and simplest indirect method of estimating heat flux is via the bulk aerodynamic approach, which is based on the bulk transfer formula (Equation (1b)). In radiosonde data, the height of the maximum gradient of vertical change in potential temperature is defined as the ABL height. The relationship between UAV-based heat flux and radiosonde-based ABL height was found for nine cases at SOK, BOS, and YHI.

Figure 9a shows radiosonde-based ABL heights and UAV-based sensible heat flux over time. The ABL height increased with sensible heat flux due to surface heating following sunrise. The UAV-based sensible heat fluxes were between 20–30 W m⁻² and radiosonde-observed ABL heights were between 20–250 m between 8:00 and 9:00 LST. However, UAV-based heat flux was 30.5 W m⁻² and radiosonde-observed ABL height was 400 m at 9:00 LST on November 30, 2018. It has been shown that ABL heights can grow quickly in the morning and reach around 1200 m because of the cloud-free convective conditions, with surface skin temperatures reaching 30 °C. UAV-based heat flux increased to more than 100 Wm⁻² between 11:00 and 12:00 LST. As a result, the ABL height of 20–400 m rapidly increased to 750–1,050 m. When the UAV-based heat flux was 190 W m⁻² at 14:00 LST on September 11, 2018 in BOS, the ABL grew further to 1250 m. These results are consistent with a previous study where a heat flux of 321 W m⁻² increased the ABL height from 457 m to 625 m, i.e., an increase of 39% [25]. The increase in ABL height was due to the increase in heat flux, which was caused by an increase in the temperature difference between the two profiles [24]. By contrast, the radiosonde-based ABL height was 350 m at 15:00 LST on December 15, 2018 at YHI. This occurred because the surface sensible heat flux was low due to clouds and clouds are generally associated with lower ABL heights.

When the surface skin cools due to nocturnal radiation processes, the sensible heat flux becomes negative instead of positive. However, while the UAV-based heat flux was 17 W m⁻², the ABL grew to 1100 m at 18:00 LST on September 11, 2018 in BOS, caused by strong winds blowing from the south coast. Strong wind shear within the ABL facilitated dynamic growth, turbulent kinetic energy was generated in the surface layer, and heat flux increased with height. Many previous studies have stated that cloudiness, heating of the ground surface, and variability in surface characteristics can affect the evolution of the ABL. Furthermore, daytime ABL height is positively correlated with surface skin temperature [47]. Conversely, turbulent flux was measured from aircraft following approximately linear gradient profiles [48,49,50].

The linear relationships between normalized heat flux (dimensionless) and normalized ABL height are shown in Figure 9b. Heat flux is assumed to decrease with height and reach zero at the top of the ABL. As shown by the gray lines, the slopes of the straight lines formed by normalized heat flux and ABL heights are the rate of change in the heat flux according to height [40]. By averaging all these slopes, the representative heat flux rate of change was approximately 1.6, as shown by the dashed black line. Therefore, it is possible to estimate ABL height using the UAV-based sensible heat flux.

Our results confirm that the ABL height varies linearly with the surface heat flux and indirect estimation of the ABL height using a UAV-based heat flux is possible; however, it is necessary to analyze the relationship between sensible heat flux and ABL height using more UAV observation cases in future research.

4. Conclusions

Accurate measurement of vertical profiles of temperature, moisture, and wind in the lower ABL at high spatial and temporal resolutions are important for understanding ABL development mechanisms. In this study, we show that these profiles, and estimates of wind speed, heat flux, and ABL height can be achieved using a rotary-wing UAV equipped with a temperature/humidity sensor approximately 40 cm above the UAV. We confirmed that UAV-observed temperatures and relative humidity were in good agreement with the AWS-observed temperatures and relative humidity, with correlation coefficients of 0.98 and 0.9, and a mean bias of −1.05 °C and 6.66%, respectively. Our own relationship (Equation (15)) between flight speed and tilt angle was estimated under calm wind conditions and compared with previous studies. Flight speed was calculated for wind conditions using our relationship and the UAV-based wind speed was estimated using the flight and ground speeds with the wind triangle. The UAV-based wind speeds agreed with the AWS-based wind speeds within an error of 10%. Accurate wind measurements provided excellent UAV-based heat flux with a mean bias of −0.75 W m⁻² and RMSE less than 10 W m⁻². Radiosonde-based ABL height increased with the UAV-based heat flux by heating the surface after sunrise and vice versa. The relationship between UAV-based heat flux and radiosonde-based ABL height can be used to estimate ABL height. We confirmed the usefulness of estimating UAV-based wind speed, heat flux, and ABL height; however, additional studies are needed to further evaluate our newly proposed approach. The well-known Monin–Obukhov similarity theory is mostly used to close bulk transfer parameterization. However, it is not easy to obtain an accurate surface skin temperature to determine surface flux parameterization based on the Monin–Obukhov similarity and bulk transfer method. The estimation of surface skin temperature was examined using linear extrapolation of the profile of the UAV-observed temperature and evaluated with the observed surface temperature sensor. The correlation coefficient and RMSE showed agreement (0.9 and 3.5 °C, respectively), but the largest difference between the temperatures was 10.5 °C at 10:00 LST on July 10, 2016. This could be due to specific heat difference between the air and land surface and land characteristics. This is a limitation of our newly proposed approach based on the bulk transfer method, which requires surface temperature. Thus, future studies are needed to estimate surface temperatures using a UAV. Utilizing a thermal camera (Duo pro R 640, FLIR, Wilsonville, OR, USA) to observe surface temperature directly might be a viable solution. Other solutions include using the gradient method [33] or the profile method [51], which do not require surface temperature and roughness length. Utilizing UAVs to estimate wind speed, heat flux, and ABL height has advantages over existing standard observation methods. They have the advantage of low-altitude flights at no risk to human pilots and can capture the evolution of ABL over transition zones such as air-sea boundaries. UAV can also help determine the horizontal variability of heat flux at a regional scale. Consequently, UAV can play an important role in acquiring ABL height and near surface atmospheric flux, especially in extreme conditions where satellite-derived flux products are unreliable. UAV-based wind vectors supplement the data of the area that cannot be observed by the UHF wind profiler.