Towards vertical wind and turbulent flux estimation with multicopter uncrewed aircraft systems

Vertical wind velocity and its fluctuations are essential parameters in the atmospheric boundary layer (ABL) to determine turbulent fluxes and scaling parameters for ABL processes. The typical instrument to measure fluxes of momentum and heat in the surface layer are sonic anemometers. Without the infrastructure of meteorological masts and above the typical heights of these masts, in situ point measurements of the three-dimensional wind vector are hardly available. We present a method to obtain the three-dimensional wind vector from avionic data of small multicopter uncrewed aircraft systems (UAS). To achieve a good accuracy in both average and fluctuating parts of the wind components, calibrated motor thrusts and measured accelerations by the UAS are used. In a validation campaign, in comparison to sonic anemometers on a 99 m mast, accuracies below 0.2 m s⁻¹ are achieved for the mean wind components and below 0.2 m² s⁻² for their variances. The spectra of variances and covariances show good agreement with the sonic anemometer up to 1 Hz temporal resolution. A case study of continuous measurements in a morning transition of a convective boundary layer with five UAS illustrates the potential of such measurements for ABL research.

1 Introduction

The nature of turbulence in the atmospheric boundary layer (ABL) is multifaceted and depends on internal and external parameters such as thermal stratification and shear, surface heat fluxes, and surface roughness. Turbulent fluxes of momentum and heat are the main driver for mixing in the ABL and thus determine its diurnal cycle. Determining turbulent fluxes requires measuring the three-dimensional wind vector and temperature with a temporal resolution that covers all turbulent scales that contribute significantly to the process. The eddy covariance (EC) method is a way to determine the momentum fluxes from the covariances of the wind components $\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}$, $\stackrel{\mathrm{‾}}{v^{\prime }{w}^{\prime }}$, $\stackrel{\mathrm{‾}}{u^{\prime }{v}^{\prime }}$ and the heat flux from the covariance between potential temperature and vertical wind $\stackrel{\mathrm{‾}}{{\mathit{\theta }}^{\prime }{w}^{\prime }}$ (Aubinet et al., 2012). From the fluxes, scaling parameters for the ABL such as shear stress u^*, convective velocity scale w^*, Obukhov length L, and Richardson number Ri can be derived (Stull, 1988), which are essential for a profound understanding of boundary layer processes. It is thus evident that the measurement of vertical wind velocity and its fluctuations is crucial.

In practice, the most prominent instrument to measure turbulent fluxes in the ABL is the three-dimensional sonic anemometer (Kaimal and Businger, 1963) which has undergone extensive development and commercialization since the 1990s and is nowadays the essential part of every surface energy balance station (Mauder and Zeeman, 2018). Operating sonic anemometers is mostly limited to the surface layer (Baldocchi et al., 2001) or requires meteorological masts to reach heights up to 300 m above ground, which involves significant costs and infrastructure (Wolfe and Lataitis, 2018). To obtain representative quantities for a certain area with single point measurements, assumptions of homogeneity and ergodicity need to be made. It has recently been shown that coherent structures such as secondary circulations can violate the ergodicity assumption and lead to errors in the EC method (Mauder et al., 2020; Morrison et al., 2022). Only few experiments exist where multiple sonic anemometers have been installed over a wide area and heights above the surface layer, due to the large logistical effort. Some examples are the Great Belt Coherence experiment over water (Mann et al., 1991) and the Perdigão 2017 experiment in complex terrain (Fernando et al., 2019).

One possibility to obtain area-averaged fluxes are airborne measurements with research aircraft incorporating flow probes and fast temperature sensors. Such measurements can cover many square kilometers in a relatively short amount of time but also require assumptions of stationarity throughout the flight which can lead to sampling errors (Mahrt, 1998). To overcome limitations of crewed aircraft when flying at low altitudes and at lower costs, fixed-wing uncrewed aircraft systems (UAS) have been deployed in the past to measure turbulent fluxes (van den Kroonenberg et al., 2008; Wildmann et al., 2015).

To combine the benefits of the flexible deployment of UAS and allow point measurements instead of area-averaged measurements, multicopter UAS are a possible solution. Wetz et al. (2021) showed that such systems can be deployed in small fleets and obtain distributed horizontal wind measurements with a good accuracy. They become increasingly popular for vertical profiling of the ABL when temperature and humidity sensors are carried (Koch et al., 2018; Segales et al., 2020). The measurement of vertical velocity with such systems is, however, challenging due to the obvious distortion of the flow by the thrust of the rotors. Thielicke et al. (2021) showed that with a sonic anemometer installed on a multicopter in sufficient distance to the rotors, measurements of vertical wind are feasible. However, such a design increases the complexity of the flight system, its weight, and its cost. In order to operate a fleet of small and lightweight multicopter UAS, we propose a method to obtain the three-dimensional wind vector with avionic data alone, including calibrated motor thrust measurement. This is the first method that allows vertical wind and turbulent flux measurement with multicopter UAS without an external sensor.

In Sect. 2, we describe the UAS and the methodology to filter sensor data, calculate the forces that act on the UAS and calibrate the vertical velocity measurements. Section 3 describes the setting of the calibration and validation campaign at the boundary layer field site (Grenzschichtmessfeld, GM) Falkenberg of the German Meteorological Service (DWD). Section 4 contains the results of three-dimensional wind measurements, second-order statistics, and the derivation of momentum fluxes in comparison to measurements of sonic anemometers on a 99 m mast. A case study with measurements during a morning transition of a convective boundary layer illustrates the potential of the method. We conclude and give an outlook on future research in Sect. 5.

2 Methods

2.1 System description

We operate a fleet of 35 UAS of type Holybro QAV250 within the SWUF-3D project (Simultaneous Wind measurements with Unmanned Flight Systems in 3D). The main characteristics are described in Wetz et al. (2021). With a weight below 1 kg and a frame size of 0.25 m, the UAS are small enough to be operated in the open category of the European Regulation 2019/947. Operation in a fleet with fewer pilots than the UAS requires special permission from the flight authorities. A Pixhawk® 4 Mini autopilot is used as the flight controller and allows access to all avionic data, which are stored at high sampling rates on an internal SD card. The sampling rate depends on the type of data and is presented in more detail in Sect. 2.5. In the current configuration, flight times up to 17 min can be achieved. The UAS are operated in a weather vane mode, which means that a controller will minimize the roll angle so that the front side of the UAS will always face the wind in the main wind direction. Figure 1 shows two UAS in flight.

Figure 1 Two Holybro QAV250 UAS in flight. In the background the 99 m mast with a sonic anemometer on the boom at 50 m can be seen.

2.2 Motor thrust calibration

In order to determine the forces that act on the UAS, it is important to know the actuator forces of the four rotors. There are no direct sensors for rotational speed, torque, or even thrust in the system. However, the motor controller PPM (pulse-position modulation) signal and the battery voltage and current are measured by the autopilot. The PPM signal is directly related to the voltage that is applied to the motor and is in theory directly proportional to the revolution speed. The voltage at each rotor is

$$\begin{array}{}\text{(1)}& U_i=U_{\mathrm{b}}{\displaystyle \frac{s_i-s_{\min }}{s_{\max }-s_{\min }}},\end{array}$$

with U_b the battery voltage and s_i the servo PPM command (s_min for zero and s_max for full throttle). It is important to include the battery voltage here as it decreases throughout the flight. Assuming a linear relationship for the electric servo motors, the revolution speed of each motor can be determined as

$$\begin{array}{}\text{(2)}& {\mathrm{\Omega }}_i=k_V{U}_i.\end{array}$$

The motor thrust T_i can be calculated from the rotational speed Ω_i, the rotor radius R, the rotor swept area A=πR², air density ρ, and the thrust coefficient C_t (Leishman, 2016, Eq. 2.31):

$$\begin{array}{}\text{(3)}& T_i=C_{\mathrm{t}}\mathit{\rho }A{\mathrm{\Omega }}_i^{\mathrm{2}}{R}^{\mathrm{2}}.\end{array}$$

In order to determine this relationship experimentally, in addition to the thrust coefficient of the rotor, which relates thrust to revolution speed, a hardware-in-the-loop test stand for a single motor–rotor setup is used. A schematic of the setup is shown in Fig. 2. We can thus determine the motor constant k_V and the thrust coefficient C_t experimentally.

Figure 2 Schematic of the motor test stand. The rotor is producing thrust in the free stream and the counteracting force of the equal length lever arm is measured with a scale on a table.

2.3 Quadrotor dynamics and rotor aerodynamics

To describe the motion of a quadrotor in three dimensions in its body-fixed frame of reference (x to the front, y to the right, z downwards), the rigid body model as described in Wetz et al. (2021) can be used. The force balance equations for translational motion within this model are as follows:

$$\begin{array}{}\text{(4)}& {\displaystyle }m(\ddot{x}+q\dot{z}-r\dot{y})& {\displaystyle }=-mg\left[\sin \right(\mathit{\theta }\left)\right]+F_x,\text{(5)}& {\displaystyle }m(\ddot{y}+p\dot{z}-r\dot{x})& {\displaystyle }=-mg\left[\cos \right(\mathit{\theta }\left)\sin \right(\mathit{\phi }\left)\right]+F_y,\text{(6)}& {\displaystyle }m(\ddot{z}+p\dot{y}-q\dot{x})& {\displaystyle }=mg\left[\cos \right(\mathit{\theta }\left)\cos \right(\mathit{\phi }\left)\right]+F_z-T,\end{array}$$

with m the mass of the UAS; p, q, and r the gyroscopic terms; g the acceleration due to gravity; and θ and φ the pitch and roll angle. The external forces $\mathbf{F}=(F_x,F_y,F_z)$ are caused by wind. Equation (6) contains the total actuator thrust T, which counteracts all the external forces due to wind and gravity. The concept of operation for the UAS fleet is to hover at fixed positions so that motions due to maneuvering can be neglected. The variance of gyroscopic terms during hover periods is 2 orders of magnitude smaller than the variance of the acceleration terms. Gyroscopic terms of the force balance equation are thus neglected in the following. In Wetz et al. (2021), only the forces in the x direction in Eq. (4) were considered since the UAS yaws into the main wind direction, and the wind speed magnitude can thus be approximated by this one-dimensional simplification. In Wetz and Wildmann (2022), the two-dimensional wind vector was determined, including Eq. (5) in the algorithm and improving the dynamic response via the inclusion of the acceleration terms. In order to measure the full three-dimensional wind vector, Eq. (6) needs to be included. The single-motor thrusts T_i can be obtained from calibration as described in Sect. 2.2 and can be summed to obtain the total acting force $T=\sum T_i$. With increasing wind speed, the quadrotor lift not only generated rotor thrust, but aerodynamic effects will also start to occur. In theory, these effects are well described if the forward flight speed is much larger than the velocity that is induced by the rotor v_i (V_∞>>v_i). Glauert's high-speed approximation (Leishman, 2016) gives an equation for the total resulting thrust for a single-rotor helicopter:

$$\begin{array}{}\text{(7)}& T=\mathrm{2}\mathit{\rho }A{v}_i{V}_{\mathrm{\infty }}.\end{array}$$

For the quadrotors used in this study, the induced velocity is on the order of 10 m s⁻¹. This means that Glauert's approximation can not be used for hover flight and wind speeds of the same order of magnitude. Nevertheless, lift effects are observed by reduced power consumption and less required thrust in dependence of wind speed. We thus include a lift term F_L in Eq. (6) to account for this effect but need to determine the relationship experimentally.

Table 1 Sensor data overview and initial raw sampling rates.

With the assumption of hover state and the additional term for aerodynamic lift, we can then solve Eqs. (4)–(6) for the external wind forces, so that

$$\begin{array}{}\text{(8)}& {\displaystyle }{F}_x& {\displaystyle }=mg\sin \left(\mathit{\theta }\right)+m\ddot{x},\text{(9)}& {\displaystyle }{F}_y& {\displaystyle }=mg\cos \left(\mathit{\theta }\right)\sin \left(\mathit{\phi }\right)+m\ddot{y},\text{(10)}& {\displaystyle }{F}_z& {\displaystyle }=-mg\cos \left(\mathit{\theta }\right)\cos \left(\mathit{\phi }\right)+m\ddot{z}+T+F_L\left(F_x\right).\end{array}$$

Figure 3 Map of the GM Falkenberg site with the location of measurement points for the different flight patterns (stars) and the meteorological mast (black triangle). Background map data are taken from ©OpenStreetMap contributors, SRTM (map style: ©OpenTopoMap CC BY-SA).

2.4 Determination of the three-dimensional wind vector

To translate the forces that act on the quadrotor to wind speeds, a straightforward approach, as applied in Wetz et al. (2021), is the Rayleigh drag equation. It is, however, evident that the complex aerodynamic system of a multicopter is not perfectly represented by the physical relationships for solid bodies. Wetz and Wildmann (2022) have therefore allowed a calibration curve that is a power law with a free exponent. With this approach, higher accuracies are achievable. We apply this method here for all three dimensions, so that the wind components in the body frame of the quadrotor are

$$\begin{array}{}\text{(11)}& {\displaystyle }\left(\begin{array}{c}{u}_{\mathrm{b}}\\ v_{\mathrm{b}}\\ w_{\mathrm{b}}\end{array}\right){\displaystyle }=\left(\begin{array}{c}{c}_x{F}_x^{b_x}\\ c_y{F}_y^{b_y}\\ c_z{F}_z^{b_z}\end{array}\right),\end{array}$$

where the parameters b and c need to be calibrated for the respective direction. For the vertical direction, it is conceivable that upward and downward motion of the quadrotor relative to the air is subject to significantly different drag. For this reason, a distinction of cases is done for the z`component so that

$$\begin{array}{}\text{(12)}& {\displaystyle }{w}_{\mathrm{b}}{\displaystyle }=\left\{\begin{array}{ll}{c}_{z\uparrow }{F}_z^{b_{z\uparrow }}& F_z\ge \mathrm{0}\\ c_{z\downarrow }{F}_z^{b_{z\downarrow }}& F_z<\mathrm{0}\end{array}\right..\end{array}$$

The final step to obtain the meteorological wind vector $\mathit{U}=(u,v,w)$ is a rotation of the wind speeds from the body frame into the geodetic coordinates with the Euler angles of the quadrotor and a subtraction of the linear velocities (v_x, v_y, v_z), which are measured by the autopilot on the basis of GNSS (global navigation satellite system) information:

$$\begin{array}{}\text{(13)}& {\displaystyle }\left(\begin{array}{c}u\\ v\\ w\end{array}\right){\displaystyle }=\mathbf{R}(\mathit{\phi },\mathit{\theta },\mathit{\psi })\left(\begin{array}{c}{u}_{\mathrm{b}}\\ v_{\mathrm{b}}\\ w_{\mathrm{b}}\end{array}\right)-\left(\begin{array}{c}{v}_x\\ v_y\\ v_z\end{array}\right).\end{array}$$

2.5 Data preparation and filtering

The raw sensor data of the UAS are stored on the SD card with different sampling rates depending on the data type but with a synchronized timestamp. Table 1 presents the raw data types and variables that are used for the wind calculation and their initial sampling rate. Reasons for the different sampling rates are the capabilities of the sensors and the requirements for the flight controller. Except for the battery data, all data are sampled with a rate of 10 Hz or higher. The raw data can be subject to significant noise in the high-frequency range due to sensor noise or vibrations. Preconditioning steps are taken before further processing of the data to obtain clean and synchronized data. Accelerometer data are filtered with a Butterworth filter to prevent anti-aliasing effects from frequencies above 10 Hz. All data are then interpolated and synchronized to a common time step with Δt=0.1 s. After this step, noise can still be recognized in accelerometer and actuator data, so that for the wind calculation, frequencies above 2 Hz are cut off with a fast Fourier transform (FFT) filter in the frequency domain.

Figure 4 Calibration results of the motor test stand. (a) The calibration curve of revolution speed versus applied voltage and (b) the thrust coefficient curve. The fitted calibration curves are shown as red lines.

2.6 Calculation of turbulence parameters

From the three-dimensional wind data, we can calculate momentum fluxes, friction velocity, and turbulence kinetic energy (TKE). Momentum fluxes are calculated as the covariances between the wind components, i.e., $\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}$, $\stackrel{\mathrm{‾}}{v^{\prime }{w}^{\prime }}$, and $\stackrel{\mathrm{‾}}{u^{\prime }{v}^{\prime }}$. The friction velocity u^* is calculated from the momentum fluxes according to Stull (1988):

$$\begin{array}{}\text{(14)}& u^{*}={\left({\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{v^{\prime }{w}^{\prime }}}^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{4}}.\end{array}$$

TKE is defined as half of the sum of the variances of the three wind components:

$$\begin{array}{}\text{(15)}& \mathrm{TKE}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\stackrel{\mathrm{‾}}{u^{\prime }}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{v^{\prime }}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{w^{\prime }}}^{\mathrm{2}}\right).\end{array}$$

3 Field experiment

A calibration and validation campaign for the SWUF-3D fleet was carried out in the framework of the FESSTVaL (Field Experiment on submesoscale spatiotemporal variability in Lindenberg) campaign in June/July 2021 at the GM Falkenberg site of the DWD. Figure 3 gives a map of the location with an indication of the measurement positions, including the location of the 99 m meteorological mast at the site. Two sonic anemometers (Metek UAS-1) are installed at 50 and 90 m height and serve as the reference for calibration and validation. Their measurements are distorted by the tower for wind directions between 345 and 50^∘ via north and are thus disregarded for such conditions.

Calibration flights were done with five UAS at 50 and 90 m respectively, in close proximity (Δx=20 m) to the mast (see black stars in Fig. 3 and Wetz and Wildmann, 2022). During a morning transition of a convective boundary layer on 28 June 2021, continuous measurements were done with two sets of five UAS stacked in five heights (10, 50, 90, 150, 200 m) at the location of the blue star in Fig. 3. The two sets were measuring alternately in close separation of 5 m and with a short overlap of 1 min to guarantee continuous observations. Spatially distributed measurements (red stars) were performed to calculate coherence of the flow but are not analyzed in this study.

4 Results

4.1 Motor test stand calibration

With Eq. (1), the voltage at the motor can be estimated from the actuator outputs. Results from the hardware-in-the-loop experiment (Fig. 4a) reveal a nearly linear relationship between voltage and revolution speed. The coefficient that relates voltage to revolution speed is $k_V\approx \mathrm{30.7}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{V}}^{-\mathrm{1}}$.

From Fig. 4b we find that the thrust curve T_i as a function of revolution speed Ω_i is best approximated with a polynomial of the second order. Such behavior is quite common and has been described before (Bangura et al., 2016). We determine the calibration curve as follows:

$$\begin{array}{}\text{(16)}& T_i=\mathrm{4.6}\cdot {\mathrm{10}}^{-\mathrm{5}}{\mathrm{\Omega }}_i^{\mathrm{2}}-\mathrm{1.68}\times {\mathrm{10}}^{-\mathrm{3}}{\mathrm{\Omega }}_i.\end{array}$$

4.2 In-flight wind vector calibration

For calibration of the dynamic equations as described in Sect. 2.3, we use UAS#12 because it was operating in most of the calibration flights (i.e., 14 flights, with a wind speed range of 1–8 m s⁻¹). The results of the calibration of the horizontal wind components are presented in Wetz and Wildmann (2022). In this study, we focus on the calibration of the vertical velocity component. This calibration is done in two steps. First, the lift force, which reduces the necessary thrust in hover flight, is determined as a function of the wind speed in the x direction. Second, thrust minus lift is calibrated against the vertical wind velocity as measured by the sonic anemometers at the 99 m mast.

4.2.1 Lift force

To estimate the lift force F_L, we analyze the measured thrust force during all calibration flights in dependency of the horizontal wind. Horizontal wind is not known a priori, so we use the drag force in x-direction in the body frame F_x as a proxy which depends on the wind speed as described by Eqs. (8) and (11). It is not possible in a field experiment to perfectly decouple the lift effects from the wind in x-direction from the drag and lift in z-direction. For calibration purposes, we assume that the x and z-component of the wind are uncorrelated and that vertical wind is normally distributed around zero. To reduce the scatter, a 10 s moving average is applied to the data. Figure 5 shows the point cloud of lift versus drag force in x-direction for all calibration flights of UAS#12. The red line gives a polynomial fit. It can be seen that significant lift is generated up to F_x≈1 N lift, before it abruptly decreases, as it is typical for stall conditions. This force in x-direction corresponds to a horizontal wind speed of approximately 8 m s⁻¹. Beyond that wind speed, only a few data are available, which yields a high uncertainty in the calibration at higher wind speed conditions.

Figure 5 Calibration curve of thrust generated by lift F_L versus drag force in x direction F_x in the body frame. The greyscale represents the probability of occurrence.

The determined polynomial for the lift correction is

$$\begin{array}{}\text{(17)}& F_{\mathrm{L}}=T_{\mathrm{0}}+\mathrm{0.9}{F}_x-\mathrm{0.27}{F}_x^{\mathrm{5}},\end{array}$$

where T₀ is a thrust offset value that can depend on changes in UAS mass, changes in motor thrust coefficients, or unmodeled aerodynamic effects.

Figure 6 Scatterplot of vertical wind of the sonic anemometer versus the acting force on the UAS in the body frame. The greyscale represents the probability of occurrence in the dataset. The solid red line shows the estimated calibration curve. The dashed red line shows a linear fit for comparison.

4.2.2 Thrust offset

We observe that the absolute value of thrust T₀ that is needed to maintain altitude in the hovering state varies between single flights and especially between individual UAS. The reasons can be differences in UAS mass, motor performance, or aerodynamic effects. At the present state of development and hardware robustness, we find that T₀ needs to be determined in each single flight as the average $\stackrel{\mathrm{‾}}{T}$ over the whole flight for best performance. This is justified in flat terrain where zero average vertical velocity can be assumed in most cases. For cases with non-zero vertical velocity, determining T₀ from the average over a single flight can yield errors in the absolute value of w but will have little effect on the variance, which is most important for the calculation of TKE and fluxes.

4.2.3 Vertical velocity

In the second step of the calibration, vertical wind is calibrated against thrust (including aerodynamic lift, see above) with data from the FESSTVaL field campaign. Figure 6 shows the scatterplot of all calibration flights of UAS#12. Since measurements of UAS and sonic anemometer are comparatively far apart (≈20 m) and absolute values of vertical wind are rather small, it is not surprising that significant scatter is found in the comparison. It is nevertheless possible to fit a curve according to the theory of Eq. (18). The fitted equation is

$$\begin{array}{}\text{(18)}& {\displaystyle }{w}_{\mathrm{b}}{\displaystyle }=\left\{\begin{array}{ll}\mathrm{3.3}|F_z{|}^{\mathrm{0.85}}& F_z\le \mathrm{0}\\ -\mathrm{1.6}|F_z{|}^{\mathrm{0.6}}& F_z>\mathrm{0}\end{array}\right..\end{array}$$

Figure 7 Example of time series of u, v, and w for a flight on 2 July 2021 by UAS#13. UAS data (red) are compared to the sonic anemometer at the same height (black).

In comparison to a linear fit, the exponential fit appears to follow the measurements increasingly better for higher negative values of w. However, high absolute values of w are not measured with a high probability.

4.3 Wind vector validation

Figure 8 Scatterplots of flight-averaged UAS measurements versus sonic anemometer for u, v, and w (a–c), as well as the corresponding variances (d–f) and covariances (g–i).

The parameters that are determined for UAS#12 are now applied to the whole fleet. We want to compare the wind components u, v, and w in the geodetic coordinate system according to Eq. (13). We rotate the vector horizontally into the main wind direction for both the sonic anemometer and UAS. The coordinate transformation allows us to separate motion around the pitch axis (streamwise) from motion around the roll axis (lateral), which are treated with different calibration coefficients. The wind direction as measured by the UAS is used for the rotation of both UAS and sonic anemometers. The time series in Fig. 7 demonstrates the quality of the wind retrieval in all three dimensions and at a high resolution for a flight in a turbulent ABL with horizontal wind speeds ranging from u=2.5–7.5 m s⁻¹. Figure 8 gives scatterplots of the comparison of flight-averaged wind components; their variances; and the three covariances $\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}$, $\stackrel{\mathrm{‾}}{v^{\prime }{w}^{\prime }}$, and $\stackrel{\mathrm{‾}}{u^{\prime }{v}^{\prime }}$. Since the wind vector is rotated into the main wind direction during the flight time, the lateral wind component v is zero by definition (Fig. 8b). Flight-averaged vertical velocities are smaller than the uncertainty of the measurement so that a linear regression does not yield a good correlation (p value for test of null-hypothesis P>0.1). The root-mean-square (rms) deviation is ${\mathit{ϵ}}_w=\mathrm{0.25}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$. Better parameters to evaluate the quality of the wind retrieval algorithm are the variances of the wind components in lateral and vertical direction. For those, the linear regressions show a very good correlation between UAS and sonic anemometer results (R²>0.9, Fig. 8e, f). The covariance measurements are a good indicator as to whether an independent measurement of horizontal and vertical wind can be guaranteed with the proposed method. The high correlation proves the validity of the method and the rms deviation of approximately $\mathit{\epsilon }=\mathrm{0.1}{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{2}}$ for all directions can be considered the prevailing uncertainty of the UAS measurements (Fig. 10g, h). The lower R² value for $\stackrel{\mathrm{‾}}{v^{\prime }{w}^{\prime }}$ is likely due to the imperfect rotation into wind direction for both UAS and sonic anemometer measurements, which introduces larger relative errors for the overall small values. In Fig. A1 we additionally present the validation results in a meteorological wind coordinate system (u is eastward, v is northward, w is up).

4.4 Validation of resolved turbulence scales

A closer evaluation of the resolved turbulence scales by the UAS can be achieved with the analyses of power spectra of the three wind components in comparison to the sonic anemometer. We choose two representative flights of UAS#12, one in low-wind and low-turbulence conditions in the morning of the 26 June (blue), the other in a turbulent convective boundary layer on 2 July (red). The spectra of the three wind components in Fig. 9 show that an inertial subrange of locally isotropic turbulence with a slope of $-\mathrm{5}/\mathrm{3}$ (dashed grey line) can be observed in both cases. For the low-turbulence case, the spectra as measured by the UAS deviates earlier from the curve as measured by the sonic anemometer than in the high-turbulence case. It can be read from the plots that frequencies above 1 Hz and spectral power below ${\mathrm{10}}^{-\mathrm{2}}{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{2}}{\mathrm{Hz}}^{-\mathrm{1}}$ can barely be resolved by the UAS but are dominated by noise. This effect is also observed for the cospectra of the wind components (Fig. 10). For conditions with low turbulence, the noise contributions are limiting the resolution of variance and covariance measurements with the UAS. From the noise level, we estimate this limit in resolution for the variances and covariances to $\approx \mathrm{0.1}{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{2}}$, which is almost equal to the root-mean-square error (RMSE) that was determined in Sect. 4.3.

Figure 9 Comparison of power spectra of all three wind components between sonic anemometers and UAS.

Figure 10 Comparison of power cospectra between sonic anemometers and UAS.

Figure 11 Time series of horizontal wind speed U, wind direction Ψ, temperature T, friction velocity u^*, and TKE at 50 m above ground level during the morning transition of 28 June 2021 in comparison between UAS and sonic anemometers (left). The continuous lines show 10 min moving averages, whereas the dots give the 30 min mean for each quantity. On the right, Hovmöller diagrams for the observed height range from 10–200 m are shown. The black line indicates the height of the mixed layer below the lowest temperature inversion as observed from vertical profiles of potential temperature θ. For wind speed U, wind direction Ψ, and θ, continuous ascent vertical profiles from a single drone are shown in the background of the measurements from the hovering UAS at 10, 50, 90, 150, and 200 m.

4.5 Momentum flux during a morning transition

A good case to study the validity and capability of the turbulence and flux measurements from the UAS is a transition phase of a convective boundary layer (CBL), which features some significant variability and an increase in turbulence at constant heights. In the morning of 28 June 2021, conditions were present that allowed for the observation of a morning transition of a CBL that was mainly driven by surface heating. In the cloudless night, a stable nighttime inversion had developed, and no clouds were present during the period of observation from 04:00 to 08:00 UTC. A steady increase in shortwave downward radiation from ≈100 to $\approx \mathrm{700}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ was observed.

Continuous measurements during this period with five stacked UAS at the heights 10, 50, 90, 150, and 200 m were achieved by launching two separate sets of five UAS in turns, with a short overlap of 60 s. The UAS were located close to the meteorological mast as indicated by the blue star in Fig. 3. In total, 23 flights with five UAS each were performed in this pattern. Additionally, a sixth UAS was operating in close proximity (≈5 m apart), profiling the same height range with continuous ascents and descents. The ascents of these flights allow us to obtain profiles of wind speed, wind direction, temperature, and humidity with a high vertical resolution but are not used for the turbulence retrieval. Figure 11 shows the results of the measurements during the morning transition. On the left, the time series of the UAS at 50 m are shown in comparison to the sonic anemometer at the mast at the same level. The data are post-processed with a 10 min moving average for better visualization. Between 07:30 and 07:45 UTC, data are missing because one UAS did not take off as planned. The dots in the plots show the 30 min average with error bars representing the uncertainties that were estimated in Sect. 4.3. We find a good agreement between UAS and sonic anemometer for all variables. The wind direction bias is higher than the determined RMSE from the calibration, especially during those times when the wind speed is lowest. In such conditions, the weather vane controller does not align the UAS very well with the mean wind direction. Around 05:30 UTC, the mixed CBL is growing through the measurement height and wind speed is very weak, and thus the estimation of wind direction deviates by almost 20 ^∘ between the UAS and sonic anemometer. In addition, temperature measurements differ most during this period. TKE only picks up significantly after 06:30 UTC, with a large peak just before 08:00 UTC that can be seen in both the UAS and sonic anemometer results. The friction velocity u^* is quite variable, also increasing after 06:30 UTC.

In the right column of Fig. 11, Hovmöller diagrams show the transition period in dependency of time and height. From the profiles of potential temperature, an estimate of the boundary layer height as the height of the mixed layer below the lowest inversion is derived and shown as a black line. A low-level jet (LLJ) was present at 100 m above ground level (a.g.l.) that disappears as the mixed layer grows and is the reason for the observed higher wind speeds in the early morning at 50 m and above. It can be clearly seen how low turbulence and high wind speeds above the inversion are separated from low wind speeds and higher turbulence in the mixed layer. Despite the large number of individual UAS that were used to take these measurements (i.e., 12 UAS), no significant biases or discontinuities can be found in the data, which is a good validation of the accuracy of the calibration.

5 Conclusions

In this work, we present a method to obtain the three-dimensional wind vector with small multicopter UAS. To obtain the vertical wind velocity, motor thrust is included in the equations of motion, and the resulting forces that act on the UAS are calibrated against sonic anemometers to obtain the wind velocity with an rms deviation below 0.2 m s⁻¹ in all three dimensions. We showed that the variance of the wind components can be obtained particularly well with a remaining rms deviation below 0.2 m² s⁻² in all dimensions. The range and limits of the resolvable scales of turbulence are presented with spectra and cospectra of the wind components. It is evident that there is a noise limit in the temporal resolution of approximately 1 Hz and a limit for low-turbulence conditions with spectral power density below ${\mathrm{10}}^{-\mathrm{2}}{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{2}}{\mathrm{Hz}}^{-\mathrm{1}}$. The fact that we could exercise the calibration with a single UAS and apply it to a fleet of more than 20 UAS without a big loss of accuracy shows the robustness of the method and the applicability in practice. From these findings, it can be concluded that UAS of the presented size and weight are capable of resolving eddies of only a few meters. Bigger and heavier UAS will necessarily disturb the flow more, meaning that a small and lightweight design is essential for a high resolution of wind measurements.

In the presented dataset, the range of horizontal wind velocities is limited to $\approx \mathrm{8}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$. We see that above such wind speeds, aerodynamic effects change, and more studies will be necessary in future to extend the calibration to higher wind speeds. We remove a thrust bias for every individual flight in this study because the assumption of zero average vertical wind speed can be made in the flat terrain in which the experiment was carried out, and this approach will lead to the best estimates of vertical velocity variance. For conditions with a significant average vertical velocity component, different strategies will need to be developed in the future, which may incorporate more robust hardware configurations, an integrated weight estimator in the Kalman Filter, or individual calibration of motor thrusts for UAS. For both constraints, the obtained dataset from the FESSTVaL campaign cannot be used to improve the calibration. In general, calibration with field data is prone to large uncertainties due to the complex ABL flow. The calibration dataset is thus subject to significant scatter, and the polynomials that are used for the relationship between forces and wind speed are not strictly based on physical principles, as simplifying assumptions are made. In the future, we plan to conduct wind tunnel experiments to study the aerodynamic effects in a broader range of horizontal and vertical wind velocities and in more detail. Hardware modifications that allow for hover flights indoors are currently being prepared for this purpose. We are also preparing temperature sensors with a faster response time to be able to measure sensible heat flux in a broad range of atmospheric conditions in the future.

We have shown the potential of multicopter UAS measurements for fundamental ABL research with a case study during a morning transition of a CBL. In this case study, the calibration was applied to 12 individual UAS and a seamless time series over 4.5 h could nevertheless be realized at six measurement heights. The evolution of wind speed and temperature, as well as TKE and shear stress, is well captured in comparison to the reference from sonic anemometers throughout the transition period. This example demonstrates the advantage of operating a fleet of multiple UAS to obtain spatiotemporal in situ data at flexible measurement points throughout the ABL. The possibilities for deploying such a system for ABL studies in complex flow environments such as wind farms or mountainous terrain are numerous, and an observational gap can potentially be closed with this approach.

Appendix A

Figure A1 Scatterplots of flight-averaged UAS measurements versus sonic anemometer for u, v, and w (a–c), as well as the corresponding variances (d–f) and covariances (g–i) in the meteorological frame of reference.