Upscaling point measurements from micrometeorological towers is a challenging task that is important for a variety of applications, for example, in process studies of convection initiation, carbon and energy budget studies, and the improvement of model parameterizations. In the present study, a technique was developed to determine the horizontal variability in sensible heat flux H surrounding micrometeorological towers. The technique was evaluated using 15-min flux observations, as well as measurements of land surface temperature and air temperature obtained from small unmanned aircraft systems (sUAS) conducted during a one-day measurement campaign. The computed H was found to be comparable to the micrometeorological measurements to within 5–10 W m−2. Furthermore, when comparing H computed using this technique with H determined using large-eddy simulations (LES), differences of <10 W m−2 were typically found. Thus, implementing this technique using observations from sUAS will help determine sensible heat flux variability at horizontal spatial scales larger than can be provided from flux tower measurements alone.
Eddy covariance observations from micrometeorological towers have long been used to determine exchanges of heat, mass, and momentum between the land surface and the overlying atmosphere (e.g., Baldocchi et al. 1988; Foken and Wichura 1996; Aubinet et al. 2000; Baldocchi et al. 2001). Oftentimes, questions arise as to the representativeness of these point measurements and how best to upscale these measurements to larger areas (e.g., Baldocchi 2014; Xu et al. 2017). Acquiring this information is necessary so that these measurements can be more reliably used in, for example, process studies of convection initiation (e.g., Trier et al. 2004; Kang and Bryan 2011), carbon budget studies (e.g., Baldocchi et al. 2001), and the improvement of model parameterizations (e.g., Gryanik and Hartmann 2002; LeMone et al. 2008).
To help determine the horizontal variability in fluxes surrounding micrometeorological towers, small unmanned aircraft systems (sUAS) may be used. Over the past 10–20 years, sUAS have emerged as an important tool for atmospheric science research (e.g., Holland et al. 2001; Spiess et al. 2007; Houston et al. 2012), as they have been used to make quasi-continuous vertical and horizontal profiles of, for example, temperature (e.g., van den Kroonenberg et al. 2012), wind speed and direction (e.g., Bonin et al. 2013; Palomaki et al. 2017), and aerosol concentrations (e.g., Corrigan et al. 2008). More recently, multispectral and thermal cameras on board sUAS have been used to estimate sensible and latent heat fluxes, generally at scales of up to a few hundred meters (e.g., Hoffmann et al. 2016a,b; Ortega-Farías et al. 2016). However, few studies have used sUAS for determining fluxes surrounding micrometeorological towers. In the present study, we developed a new approach to estimate fluxes using sUAS and measurements from a surface flux station. We evaluated our approach using observations obtained from a site in eastern Tennessee coupled with large-eddy simulations (LES) for this case.
2. Datasets and models
a. Surface meteorological observations
Micrometeorological observations were obtained from a site (36.1259°N, 83.7944°W) that was located 23 km northeast of Knoxville and 3 km south of Corryton. The study site was characterized by a flat grassy field encompassing an area of approximately 0.6 × 0.6 km2. Two 2-m tripods were installed at the study site. Tripod 2 (328 m MSL) was approximately 300 m northeast of tripod 1 (330 m MSL). Both of these tripods were outfitted with an R. M. Young sonic anemometer, a downward-pointing Apogee infrared temperature sensor, and a platinum resistance thermometer (PRT) that was enclosed within an aspirated shield to minimize radiative errors. Sonic anemometer measurements were sampled at 10 Hz and were used to compute sensible heat flux H at 15-min intervals using the eddy covariance method. Samples at 1 Hz from the infrared sensor and PRT were used to compute 15-min means. We chose 15 min as the averaging time scale for the fluxes and for the meteorological measurements in order to coincide with the average length of the sUAS flights, discussed in the next section.
b. sUAS platforms
Six sUAS flights (Table 1) were conducted on 15 December 2016. During the experiment, we used two different sUAS: a DJI S-1000 and a Microdrone MD4-1000. The DJI S-1000 is an eight-rotor sUAS with a wingspan of approximately 1 m. The platform can carry a payload of up to 4.5 kg and has a mean flight time of 15 min (see Dumas et al. 2016 for more details). The Microdrone MD4-1000 is a four-rotor sUAS capable of carrying a 1.2-kg payload for up to 25 min. Internal diagnostic data, including the sUAS position, velocity, and height, were recorded at 192 Hz on the DJI S-1000. These same sUAS diagnostic data were recorded on the MD4-1000 at 100 Hz. On board each sUAS were two iMet-XQ sensors, manufactured by International Met Systems Inc., that were mounted on top of each sUAS to sample temperature, pressure, and relative humidity at a frequency of 1 Hz. The iMet-XQ sensors have a manufacturer-stated accuracy of ±0.3°C, ±5%, and ±1.5 mb for temperature, relative humidity, and pressure, respectively. Air temperatures from the two iMet-XQ sensors during the six flights agreed to within 0.17° ± 0.35°C (r = 0.99, p < 0.01) of each other (not shown), and the mean of the temperatures from the two sensors was used to compute air temperature from the sUAS.
On the underside of each sUAS was a downward-pointing FLIR Tau 2 infrared camera. The camera has a 7.5-mm lens and 336 × 256 pixel resolution and recorded surface temperature at 7.5 Hz (Dumas et al. 2016, 2017). Because of differences in the sampling frequencies of the sUAS internal diagnostics and the FLIR Tau 2 camera, data from the sUAS and infrared camera were subsampled to 1 Hz to coincide with the temporal resolution of the iMet-XQ sensors and to calculate over the area where the sUAS was flown.
Vertical profiles were conducted using the DJI S-1000 starting approximately 10 m downwind of each 2-m tripod to minimize the inflow effects generated by the sUAS rotors on the tripods’ measurements. Vertical profiles were performed between the surface and 365 m AGL, which is the maximum altitude allowed by NOAA/ARL/Atmospheric Turbulence and Diffusion Division (ATDD)’s agreement with the Federal Aviation Administration. In addition to these vertical profiles with the DJI S-1000, the MD4-1000 was flown between the two tripods at the following heights above the surface: 25, 75, 125, 175, 225, and 275 m AGL. These heights were chosen to coincide with the vertical model levels in the LES, discussed in the next section.
c. Numerical modeling
We performed LES using the Collaborative Model for Multiscale Atmospheric Simulation (COMMAS; e.g., Wicker and Wilhelmson 1995; Coniglio et al. 2006; Buban et al. 2012) to provide an independent evaluation of our proposed technique of deriving from sUAS, which is described in section 3. COMMAS is a cloud-resolving and nonhydrostatic model and includes a weakly diffusive fifth-order horizontal advection scheme (Wicker and Skamarock 2002), a 1.5-order parameterization for turbulent kinetic energy, and a modified force–restore land surface–atmosphere exchange scheme (Deardorff 1978; Peckham et al. 2004; Buban et al. 2012). The simulations had a horizontal and vertical grid spacing of 100 and 50 m, respectively, with the lowest model level 25 m AGL. The domain size was 36 km × 36 km × 6 km for the x, y, and z dimensions, respectively, and periodic lateral boundary conditions were applied.
We initialized the LES using 15-min 2-m surface and air temperature from tripod 1 and rawinsonde observations from a Graw DFM-09 rawinsonde that was launched from the site at 1100 LST 15 December. A second rawinsonde was launched at 1300 LST and was used to evaluate the LES. Although the planetary boundary layer (PBL) depth was overestimated in the LES by about 200 m, LES output generally showed good agreement with the rawinsonde observations (Fig. 1). Potential temperature, water vapor mixing ratio, wind speed, and wind direction agreed to within ±0.34 K, ±0.015 g kg−1, ±1.06 m s−1, and ±4.1°, respectively, over the lowest 3000 m of the profile. The good comparison between the observations and LES provides us with confidence in the use of our LES to help evaluate our proposed technique of deriving heat fluxes from sUAS, which we discuss in the next section.
3. Technique for determining sensible heat fluxes from sUAS
We determined using the conditional sampling technique (Businger and Oncley 1990). Although the conditional sampling technique was developed and has been used to determine fluxes of passive scalars (e.g., Cobos et al. 2002; Meyers et al. 2006), the technique can be expressed for heat flux using the following equation:
In Eq. (1), is an empirical parameter obtained from flux tower measurements, is the standard deviation of vertical velocity, is the surface temperature, and is the air temperature at a given altitude. Term σw is measured from a flux tower, and is determined by rearranging Eq. (1) as
Assuming that calculated over the tower is representative of the entire flight path, at any location () over which the sUAS has flown can be determined by solving Eq. (2) for H and using 1) and obtained from a flux tower (i.e., and , respectively); 2) obtained from a downward-pointing infrared camera on board the sUAS [i.e., ] and obtained from the sUAS while the sUAS is flown near the flux tower [i.e., ]; 3) over the flight path of the sUAS [i.e., ]; and 4) and for each location over which the sUAS is flown [i.e., and , respectively]:
Furthermore, if we 1) neglect horizontal variations in and assume that is constant over the sUAS flight path, and 2) neglect vertical variations in , which is reasonable given the small range of heights over which the sUAS can be flown (cf. section 2b), then . Thus, cancels from Eq. (3), and thus the final expression for has no dependence on . Therefore, we introduced a new empirical parameter, β, that does not depend on . We derived β while the sUAS was flown over a flux tower and sampling air temperature and surface temperature at a frequency of 1 Hz. In the present paper, we evaluated how β varies both vertically and horizontally while the sUAS was flown, which is expressed as
After quantifying the horizontal and vertical variability in β, Eq. (3) can be rewritten using β to determine the heat flux at any point along the sUAS flight path. Once β is known for a given height, the sUAS is flown at this height so that does not depend on , and can be computed as
4. Results and discussion
a. Surface observations
We found good agreement between measurements of (r = 0.82, p < 0.01), (r = 0.80, p < 0.01), (r = 0.97, p < 0.01), and (r = 0.99, p < 0.01) from the two tripods (Fig. 2). The mean difference in between the two tripods was 9.8 ± 27.1 W m−2, and the mean difference in between the two tripods was 0.008 ± 0.02 m s−1. The negligible difference in between the tripods provided us with confidence that can be assumed to be constant in the area over which the sUAS was flown and justified the assumption to neglect when calculating ’s horizontal variability [cf. Eqs. (3)–(5)].
b. Vertical variability in
To determine the vertical profile of from the sUAS flights, we selected only the portion of the flight when the sUAS was ascending at >1 m s−1 to eliminate any possible effects of the sUAS inflow generated by the propellers on , which was measured by the iMet-XQ sensors on board the sUAS. The vertical profiles indicated the largest values of were near the surface, exceeding 10 W m−2 K−1 during each flight (Fig. 3). In all flights, decreased with height, and the largest changes in occurred below 30 m AGL. Above approximately 100 m AGL, only marginal decreases in with height were observed. During the vertical profiles, varied from 5 W m−2 K−1 during flight 6 to 10 W m−2 K−1 during flight 1 above 100 m AGL.
Vertical profiles of obtained from the LES, computed using a domain average, also showed that the largest changes in occurred near the surface. However, there was a positive offset between the LES profile and the profile observed from the sUAS. This offset varied from about 8 W m−2 K−1 at the lowest model level (25 m AGL) to about 3 W m−2 K−1 near the top of the sUAS profile. The larger values for in the LES than in the sUAS observations were caused by the larger values of in the LES.
c. Horizontal variability in
The stacked horizontal flights between the two tripods that were performed with the MD4-1000 sUAS provided evidence was relatively constant when the sUAS was flown at different heights parallel to the land surface. The flight duration for each of the six heights ranged from 50 to 70 s. During flight 3, showed little variation in height; the mean of was between 13 and 16 W m−2 K−1 for the different heights at which the sUAS was flown (Fig. 4a). Additionally, standard deviations in were <1.5 W m−2 K−1, or about 5%–10% of . Although the magnitude of was lower during flight 6 than during flight 3 (Fig. 4b), with mean values for between 6 and 8 W m−2 K−1 for the different heights, standard deviations in were proportionally smaller and were typically <10% of . Thus, even though the magnitude of varied between these two flights, the small standard deviations in observed at all heights during both flights indicate 1) the largest changes in occur with height, rather than at the horizontal spatial scales at which the sUAS is flown; and 2) once is known for a given time and height, this value of can be used to infer H in the area over which the sUAS is flown.
d. Evaluation of technique
We evaluated our technique of deriving H from sUAS by comparing measurements from each of the sUAS flights with the meteorological and flux measurements made at tripod 2 during each of these flights. We calculated the means over each entire flight in order to correspond with the averaging time scales for the fluxes from tripod 2 (cf. section 2a). We found good agreement between obtained from the sUAS and obtained from the Apogee infrared temperature sensor at tripod 2 and averaged during the sUAS flight (r = 0.83, p = 0.04), with a mean difference between from tripod 2 and from the sUAS of −1.72° ± 0.87°C (Fig. 5a). There was also good agreement between H calculated from tripod 2 and H calculated from the sUAS (r = 0.93, p < 0.01); the mean difference between the two values of H was −8.14 ± 19.0 W m−2 (Fig. 5b). The small differences and high correlations between the sUAS measurements and the tripod 2 measurements provided us with confidence in our technique and in its ability to determine the horizontal variability in , which we do in the next section.
e. Horizontal variability in
We used measurements from flight 3 as an example to help us illustrate the horizontal variability in obtained using the technique we developed in the present manuscript. To this end, we used obtained from 75 m AGL during flight 3 when the sUAS was flown between tripod 1 and tripod 2. Over the field of view of the infrared camera, which was 61 m × 46 m when the sUAS was 75 m AGL, varied between 5° and 11°C (Fig. 6a), and H varied between 120 and 170 W m−2 (Fig. 6b).
Over the LES domain, modeled was comparable to the observed during flight 3 and very close to the measured at both tripods at 1300 LST. However, the modeled varied by only 0.5°C over the LES domain, which was about an order of magnitude smaller than the observed variability. These differences between the LES and observations arose because of the homogenous conditions that were used to initialize the LES. Furthermore, H, which we calculated by implementing the technique described in section 3 on the LES output and using a determined from the second model level (75 m AGL) for comparison with the results shown in Fig. 4b, was about 50–70 W m−2 larger (Fig. 7) in the LES than in the observations. The larger H in the LES is likely caused by the advection of thin high cirrus clouds over the study region between late morning and early afternoon. Since the LES was initialized using a cloud-free initial sounding and using periodic boundary conditions, the LES was unable to advect cloud cover across the boundaries and thus did not suppress the H that was seen in the observations. This is evident by the deeper modeled PBL (cf. Fig. 1) and lower dewpoints in the LES than in the observations between 500 and 700 mb (not shown).
However, when comparing H computed using our technique with H output from the LES (Figs. 7a and 7b, respectively), there was good agreement between these two approaches for determining H. Differences over a 5 × 5 km2 subset of the LES domain were as large as 30 W m−2 over a few isolated locations, but they were generally ±10 W m−2 (Fig. 7c). Also, the differences between the calculated and simulated H values were randomly distributed and not tied either to structures in (Fig. 7d) or to simulated H. This lack of bias provides us with additional confidence in the robustness of our technique.
5. Conclusions and outlook
In the present study, we developed a new technique to compute from sUAS. The technique used surface flux measurements and and obtained from the sUAS over an instrumented tripod to compute . This was then used with additional and sUAS measurements to derive H. Comparisons between computed using this technique with independent observations from a second instrumented tripod agreed to within 5–10 W m−2. Furthermore, computed by implementing our technique into an LES compared well with LES output of , as differences between these fluxes were generally <10 W m−2 over the LES domain.
Important to note, though, is that we demonstrated the use of this technique using measurements from a one-day field campaign. Additional studies are needed to further evaluate the technique developed in the present paper and may include, for example, using additional independent flux tower measurements, evaluating the technique over different land surface types and in different climatic regimes, etc. Nonetheless, the present technique shows promise for the use of sUAS instrumentation for determining the horizontal variability in fluxes surrounding micrometeorological towers and for deriving heat flux variability at spatial scales that are relevant to many applications, most notably for improving surface and PBL parameterization schemes where knowledge of the horizontal variability in heat fluxes is required.
We gratefully acknowledge Mark Heuer, Randall White, Tom Wood, and Joshua Blackwell at NOAA/ARL/ATDD for assisting with instrument deployment during the 15 December 2016 experiment. We also thank John Kochendorfer and Tilden Meyers at NOAA/ARL/ATDD and Paul Kelley of NOAA/ARL, whose insightful comments and suggestions helped us to improve an earlier draft of this manuscript. Finally, we acknowledge and thank the two anonymous reviewers, whose comments and suggestions helped us improve the manuscript.