Measurement of turbulent water vapor fluxes using a lightweight unmanned aerial vehicle system

We present here the first application of a lightweight unmanned aerial vehicle (UAV) system designed to measure turbulent properties and vertical latent heat fluxes (λE). Such measurements are crucial to improve our understanding of linkages between surface moisture supply and boundary layer clouds and phenomena such as atmospheric rivers. The application of UAVs allows for measurements on spatial scales complimentary to satellite, aircraft, and tower derived fluxes. Key system components are: a turbulent gust probe; a fast response water vapor sensor; an inertial navigation system (INS) coupled to global positioning system (GPS); and a 100 Hz data logging system. We present measurements made in the continental boundary layer at the National Aeronautics and Space Administration (NASA) Dryden Research Flight Facility located in the Mojave Desert. Two flights consisting of several horizontal straight flux run legs up to ten kilometers in length and between 330 and 930 m above ground level (m a.g.l.) are compared to measurement from a surface tower. Surface measured λE ranged from −53 W m−2 to 41 W m−2, and the application of a Butterworth High Pass Filter (HPF) to the datasets improved agreement to within +/ −12 W m−2 for 86 % of flux runs, by removing improperly sampled low frequency flux contributions. This result, along with power and co-spectral comparisons and consideration of the differing spatial scales indicates the system is able to resolve vertical fluxes for the measurement conditions encountered. Challenges remain, and the outcome of these measurements will be used to inform future sampling strategies and further system development.


Introduction
Vertical water vapor transport in the planetary boundary layer is an important component of the Earth's energy systems, particularly in the marine environment.This moisture supply is a key contributor to boundary layer cloud liquid water content (LWC), which in turn impacts the cloud albedo, lifetime and radiative effects -a large source of climatic uncertainty.Recent work has shown the assumption of a constant LWC with aerosol perturbations is not apparent in cloud systems inherent to the boundary layer and its dynamics, e.g.(Xue and Feingold, 2006;Roberts et al., 2008;Sandu et al., 2009).Stevens and Feingold (2009) describe a number of linkages between water content, cloud morphology and microphysics, aerosol properties, radiative forcing and boundary layer dynamics, which act to buffer cloud responses to perturbations of these elements.Observations of water vapor fluxes (as well as cloud, aerosol, and boundary layer properties) on local spatial scales with a high temporal resolution will help constrain the water budget and aid understanding of such mechanisms, facilitating their inclusion into more complete cloud-resolving models.
Studies of atmospheric Rivers (AR, ribbon-like structures extending thousands of kilometers contained within the lowest 3 km of the troposphere) also benefit from water vapor flux observations.ARs are a critical pathway for meridional moisture transport (Zhu and Newell, 1998) and play a key role in Californian flooding events (Ralph et al., 2006).Regular water vapor flux observations on local scales in AR development regions would improve understanding of their formation, maturation and ultimately help to improve forecasting algorithms.
Published by Copernicus Publications on behalf of the European Geosciences Union.

R. M. Thomas et al.: Measurement of turbulent water vapor fluxes
The major advantage of aircraft measurements is the potential to provide flux measurements which are more spatially representative than tower point sources and with greater temporal resolution and accuracy than those derived from satellite observations (e.g.Smith et al., 2011).They can also cover the vertical extent of the boundary layer and produce reliable data in relatively short time periods.Slower moving, smaller aircraft offer similar temporal and spatial resolutions, but with much less disturbance than manned aircraft (Crawford and Dobosy, 1992).Unmanned Aerial Vehicles (UAVs) offer the potential to sample small-scale dynamic turbulent structures within the lower troposphere (e.g.across a 20 m thick entrainment layer and through small clouds); they are cheaper to own and maintain and they can fly in multiple formations (Ramanathan et al., 2007) and in areas where it is difficult or dangerous for manned aircraft to fly.
Eddy covariance has become a well established method for the direct measurement of the vertical exchange of gases and/or particles in the atmosphere, suitable for use in a variety of environments (Nemitz et al., 2008;Famulari et al., 2010;Spirig et al., 2005;Martensson et al., 2006).The turbulent transfer flux (F ) through a horizontal plane at the measurement height is given by the covariance between the instantaneous deviations (denoted by prime) of vertical wind velocity (w) from the averaging period (T m ) mean (denoted by overbar) and those of the tracer of interest (in this case water vapor, q), such that: (1) q = q + q (2) Latent heat (λE) is derived by multiplying the water vapor flux by the enthalpy of vaporization (λ).Ground based measurements are typically made from a tower within the so-called "constant flux" layer usually present in the lower 50 m or so of the boundary layer.Here, the flux is considered independent of height, and kinetic energy is conserved and cascaded from larger to smaller eddies with a −5/3 power law (Kolmogoroff, 1941).A measurement frequency of 10 Hz and T m ≈ 30 min are generally considered acceptable for tower based instruments to capture the frequency bandwidth of eddy sizes contributing to the flux in the surface layer, without introducing errors due to mesoscale influences (Vesala et al., 2008).
The key to successful aircraft flux measurements lies in the translation of accurately measured, aircraft-orientated, wind vectors to earth-referenced orthogonal wind vectors (Lenschow, 1986).This requires accurate measurement of the aircraft velocities and attitude with respect to the ground, which has been achieved with increasing accuracy over the years, particularly with the advent of GPS and differential GPS (DGPS) technology coupled with INS systems (Inertial Navigational Systems) e.g. through Karman filtering (Leach and Macpherson, 1991).
Scaling such systems to lightweight UAVs poses further size, mass and power challenges when developing flux instrumentation.For turbulence measurements, recent progress has been made by the development of the meteorological mini UAV (M 2 AV), which has shown promising measurements of the wind vector (Van den Kroonenberg et al., 2008).However, to acquire measurements complimentary to turbulence, work would be required to miniaturize existing scalar flux instrumentation to satisfy the 6 kg gross takeoff weight.
Payloads for lightweight UAVs have been developed previously by C 4 to measure aerosol, radiation, cloud, and meteorological properties.These measurements, when coupled with the UAV's versatility, have allowed investigation of the atmospheric heating rates of black carbon using stacked UAVs (Ramanathan et al., 2007), developed links between cloud microphysics and albedo (Roberts et al., 2008), and established insights into the long range transport of aerosols and their influence on solar absorption (Ramana et al., 2010).
Ideally, for direct comparison with surface tower fluxes, flying at low altitude over long flight legs over a uniformly homogeneous surface is desirable.The low altitude reduces vertical divergence, long legs enable capture of the low fluxcontributing eddy frequencies, and the surface homogeneity simplifies horizontal flux interpretation.If these conditions are met, aircraft flux systems will sample a turbulent wind field broadly equivalent to that advected past a tower, but on a much shorter averaging time (in the form of straight and level horizontal runs) due to the rapid motion of the aircraft through the assumed "frozen" turbulent wind field (Taylor, 1938).Airborne systems therefore require higher frequency response instrumentation than their stationary counterparts, in order to capture the smallest eddies contributing to the flux.In reality, such conditions are rare, and research is progressing to reconcile (λ)E horizontal flux variability with surface inhomogeneities (Kiemle et al., 2011;Samuelsson and Tjernstrom, 1999;Mahrt et al., 2001;Desjardins et al., 1992).
At altitudes within a convective boundary layer (CBL) between the constant flux surface layer and capping inversion height (z i ) a series of slow moving or stationary convective cells tend to form, with dimensions and movement dependent upon z i , stability, topography and wind velocity.Contributions to variances and fluxes are not only limited to those scales associated with turbulent fluxes e.g.k ≈ 10 −2 cycles km −1 , but continue out towards 0.1 cycles km −1 (Lenschow and Sun, 2007).For aircraft, extended sampling paths of such structures is estimated to be required to properly capture the spatial variability, ideally with run legs on the order of 100 km, and repeated sampling is required to drive down the random variability, especially when sampling towards the upper portion of the boundary layer (Mahrt et al., 2001).For towers, this pseudo-organization produces regions where stationary instruments could produce long periods of positive or negative w if situated within a consistently down-or up-draught region, respectively (Mahrt, 1998).Inclusion of low frequency flux contributions to the total derived fluxes is possible through the use of linear detrending, and more robustly at flux tower sites where a horizontal plane can be defined over the long-term (Wilczak et al., 2001).In summary, there is a measurement balance to be met between long flux legs, which will sample lower frequencies, but lack stationarity, and shorter repeated runs to reduce the random error; all limited by the available flight time.High pass filtering can be used to exclude poorly sampled turbulent scales and aid comparison with surface measurements (Desjardins et al., 1992).
The UAV system presented in this paper is designed to measure turbulence and water vapor fluxes for the reasons mentioned above, and to integrate with existing C 4 systems to improve their ability to address cloud/atmospheric dynamics/aerosol interactions.First we describe the new system and then its application in a continental boundary layer experiment comparing UAV fluxes with surface tower measurements.

System description
In this section we describe the UAV platform utilized, the flux system, calibrations, and some results from surface based tests undertaken during system development.

UAV
The platform used is BAE's Manta UAV offering a compact, durable and aerodynamic platform with extended flight endurance (Fig. 1).The aircraft has a wingspan of 2.6 m, a fuselage length of 1.9 m and a maximum take-off weight of 27.7 kg, of which up to 5 kg is designated for scientific payloads.With only slight modifications, it can carry both internal and external instrumentation and sensors.
Scientific missions using these aircraft have seen flight times of up to 5-h, typically cruising at a groundspeed of around 110 kph, and with flight ceilings of up to 4 km (Ramana et al., 2010).The C 4 Manta aircraft are equipped with DGPS capability and can perform automated takeoff and landing when requested.The DGPS gives the aircraft the ability to control its flight path to within less than 1 m and permits tight coordination in time and space.Iridium satellite communication is used for beyond the horizon missions.

Flux payload
Fast measurements of air velocities with respect to the aircraft are obtained from measurements of attack (α) and sideslip (β) angles (relative to the UAV) made using a precision-engineered 5-hole differential pressure aeroprobe gust probe extending beyond the shockwave point in front of the UAV.It has a diameter of 6.35 mm, length 201.66 mm, and weighs <10 g.Six pressure ports sharing a common manifold for static pressure measurement are located 30 mm back from the tip.Four of the tip holes are arranged in a cruciform pattern and connected to two differential pressure transducers (DPT) providing α and β.A centrally located port corresponds to the nominal stagnation point of the measured airflow and is coupled with the static port to a third DPT to resolve true airspeed (TAS).An absolute pressure transducer (range 0-1034 mb) measures static pressure.Connecting tubing length for all ports is kept to a minimum (76 mm of 0.794 mm ID tubing) to enable a near constant response for frequencies <100 Hz.DPTs are manufactured by All Sensors and have a range of +/− 12.5 mb, low hysteresis (0.5 %) and a measured variability of 0.007 mb.Calibrations are performed before and after field measurements using a precision manometer.Aircraft attitude and groundspeed are measured using a C-Migits-III tactical sensor which offers up to 100 Hz outputs of aircraft attitude and velocity data, and outputs Kaman filtered data from its internal GPS.This capability can be extended to use the DGPS system for greater accuracy (e.g.Khelif et al., 1999).A Campbell KH 2 O open path sensor is used to measure water vapor fluctuations of up to 100 Hz by directly measuring UV light absorption by water molecules at 123.58 and 116.49nm wavelengths, emitted by a krypton UV lamp.There is no reported pressure sensitivity of this instrument and it is widely used in surface based measurements as well as on aircraft (Khelif and Friehe, 2008).It is situated on top of the UAV fuselage and the UV path is in a location observed during windtunnel tests to be well outside of the aircraft's boundary layer, even at extreme aircraft pitch and roll angles.Absolute concentrations are compromised due to signal attenuation by scaling of the sensor windows; q absorption and is recalibrated using an accurate RTD/RH probe (Honeywell HIH-4602-C) with a 1 s response time mounted on the underside of the airframe.Gust probe pressure transducer voltages and KH 2 O data are subject to a Sallen-Key, 3-pole Low-Pass Butterworth filter with a cutoff of 20 Hz.To reduce aircraft noise prior to logging analogue data are logged on one of 16 balanced channels on a 9205 module as part of a National Instruments (NI) CRIO data logger comprising of a Field-Programmable-Gate-Array (FPGA) hardware with a real-time controller.NI Labview 8.6 logging code is directly compiled onto the FPGA, improving reliability and allowing for high determinism (25 ns of timing accuracy) with <0.5 ms of jitter.Digital GPS/INS data are simultaneously logged at 100 Hz using the NI 9870 RS232 CRio Module.
A Laser Technology Incorporated Trupulse ™ 200 rangefinder was modified and incorporated onto the UAV to provide additional height information to the aircraft avionics in addition to the GPS/DGPS instruments.It is accurate to within 0.3 m above high quality targets and has a maximum range of 655 m.A resistance thermometer is included in the instrument payload to monitor the payload bay temperature.
The system can run for 6 h powered by two 6000Ah Lipoly batteries.The total weight of the flux payload is 4.3 kg.

Calibrations
A precision look-up table of gust probe pressures with respect to probe orientation in wind fields within the operational air speeds of the UAV was provided by the gust probe manufacturer and can be used to derive calibrated α, β and total and static pressure required to characterize the local flow relative to the aircraft.See Telionis et al. (2009) for an in-depth discussion of current probe technology.However, the manufacturer's calibrations are performed with the standalone probe; additional calibrations were performed by mounting the probe on the UAV fuselage (minus wings) on two aerodynamic pivots in the 0.91 × 1.2 m low speed wind tunnel at the Aerospace Engineering department at San Diego State University.In broad accord with the methods presented in Garman et al. (2006), the UAV was stepped through α and β from +/−15 and +/−11 respectively at windspeeds spanning the operational flight speed.The results indicated a near-constant offset in β meas − β actual and a wind speed dependence on α meas − α actual which were interpolated across the measured windspeed and incorporated into the gust probe calibrations (Fig. 2).
At least one new RH/RTD probe is used during each campaign, and the manufacturer's sensor-specific calibration information is checked in the laboratory.The 9205 analogue input module voltages are self-calibrated by an internal routine which corrects for temperature differences between those measured when the manufacturer's external calibration was last performed.The KH 2 O probe uses the manufacturer's calibration for the UV laser absorption and the absolute water vapor concentration is linearly adjusted using humidity probe data.A standard correction procedure is applied for oxygen absorption following data collection.Respective linear calibrations are applied to measured analogue voltages from the KH 2 O, temperature, and humidity probes, and the pressure transducers prior to flux calculation.

Preliminary tests
Tests were performed on two occasions with the UAV mounted on a frame on a motor vehicle to check system logging capability, positional information from the INS/GPS, and gust probe performance.A sonic anemometer was mounted on the frame, with a horizontal displacement of approximately 0.75 m.The vehicle velocity was maintained close to 60 mph (28 m s −1 ) along the I5 freeway in San Diego.Positional information detailed the vehicle location and orientation perfectly throughout the tests with the exception of some incorrect altitudinal information during periods where there was a limited view of the sky and/or horizon, a typical issue with GPS systems.We could not be sure that the anemometer and the gust probe were out of the truck's boundary layer, but the resulting wind velocities and λE agreed well (Fig. 3a) with expectations.Despite the many imperfections with the setup, the UAV demonstrated values more or less consistent (at least in sign) with expectations along vehicular path (Fig. 3b) with positive (upward flux) in the vicinity of the Coronado bridge over San Diego Bay and over the San Diego River.

Measurements of convective boundary layer fluxes
Following successful ground based tests this short experiment was designed to collect turbulence and water vapor flux data in the continental boundary layer and compare with existing surface techniques, paving the way for use alongside the existing C 4 instrumented UAV fleet.

Experiment description
Two test flights were conducted on 27 May 2010 at the NASA Dryden Flight Research Center (NDFRC) located within Edwards Air Force Base (EAFB) in the western Mojave Desert, California.EAFB has designated a 112 km 2 UAV test airspace area (within the FAA restricted airspace zone R-2515) above the smooth surface of Rogers lake bed, which thus offers a multi-directional runway (Fig. 1).
The UAV data collection attempted to attain the best possible balance between the ideal conditions noted above and those specific to the UAVs and the lakebed setting.For example, the UAV work area airspace allows a maximum horizontal run length of 8.7 km, which along with run repeatability, is ultimately determined by the flight duration.However, the relative surface homogeneity and expected flux outcome combined with the relatively easy access to Dryden's controlled UAV airspace and facilities makes this setting ideal for this continental boundary layer experiment.
Risks from piloting a UAV at near surface altitudes from several km away were deemed worth taking only when the system is proven capable of the desired measurements; thus during this developmental experiment flux run altitudes were kept above 300 m.
To provide additional meteorological information and to offer an established eddy covariance flux measurement technique for comparison with UAV measurements, a λE flux measurement system (a GILL WindMaster Pro Sonic anemometer, a Licor 7500 open path Infra-red gas analyzer, and a Vaisala HMP235 temperature and RH sensor) was installed on a 10 m mast.These sensors were provided by NOAA's Earth System Research Laboratory and similar to ones used for measuring marine fluxes from ships (Bradley and Fairall, 2006).Working on the lakebed involves due care and attention to ensure the safety of the many users who depend upon a uniform runway surface.With this in mind, the tower was set up with minimum ground intrusion and close to the launch site at 34.954 • N, 117.857 • W (Figs. 1 and 4a, b).Two flights were scheduled on the 27 May 2010.The first flight, FT A , departed at 09:13 a.m. and lasted 2 h and 41 min and the second flight, FT B , departed at 12:48 p.m. and landed 1 h and 24 min later.FT B was of reduced duration due to safety concerns brought about by increasing wind speed and gustiness, and meant the airborne wind calibrations could not be completed.
Both flights followed a similar racetrack pattern; and aimed to collect turbulence data during straight and level flux runs averaged at 330, 520, 720, 930 m a.g.l.within the workspace area (Fig. 4a and b).Two to three patterns were conducted at each altitude and a summary of SW run information is presented in Table 1.The system was not configured to transmit real-time scientific data back to base over the communication link, therefore altitudes were selected prior to flight based on analysis of available meteorological data.

Flux data conditioning and analysis
For the UAV data, the 100 Hz analogue and GPS/INS UAV data are first screened for spikes (typically points >3.5 σ from a 5000 point mean) and data replaced using an interpolative replacement method similar to Højstrup (1993).Signal data is then smoothed to 50 Hz to reduce noise, and time lags between channels are investigated by locating the maximum correlation attained between staggered time series data and corrected as necessary.Geo-referenced u, v, and w wind components are then calculated from the well adopted equations of Lenschow (1986) using gust probe α and β, and average GPS/IMU measured pitch, roll and yaw angles and surface velocities.Suitable T m is determined by ogive inspection (i.e. the cumulative co-spectum).The real part of the co-spectrum between the vertical wind and the scalar tracer of interest (χ ) indicates the flux transported by turbulent eddies of that characteristic frequency, and is given by: where f denotes the frequency, W indicates a Fourier windowing function (in this case Hanning) and the asterix denotes the complex conjugate.
When analyzing and comparing aircraft data to ground measurements, it is usual to convert f to a wavenumber, k, by normalizing to the aircraft's speed, u, using: (5) Data here are not detrended, but subjected to high-pass filtering (0.04 Hz for the UAV and 0.01 Hz for the tower) to remove insufficiently sampled large eddies and facilitate comparison between the UAV and surface measurements.Integrated across the frequency domain, the ogive ideally displays low and high frequency asymptotes bounding a frequency bandwidth denoting the region in which the majority of the flux is transported.The high frequency asymptote corresponds to the inertial subrange, where turbulent energy is cascaded down towards smaller scales with a −5/3 power law.The reciprocal of the lower bandwidth limit represents the minimum T m at which the vast majority of the low frequency eddies contributing to the turbulent portion of the flux are included in their calculation.For aircraft measurements, this is related to the minimum run length required for a given altitude via the relationship given in Eq. ( 5).Ogives are calculated and inspected to allow selection of a suitable T m with which to parse the straight and level flux run data into segments for flux calculation.Figure 5 indicates a maximum averaging time of ∼1000 s (0.001 Hz) and >289 s (0.00346 Hz) is required to capture the majority of low frequency eddies contributing to the turbulent λE flux at this location, and also demonstrates the anticipated bandwidth narrowing of UAV Fig. 5. Integrated Co < w q > plots (ogives) from unfiltered data normalized to total covariance (< w q >) for 520 m UAV run at 10:40 a.m. and tower data at 10:17 a.m.Indicating T m of 289 (1/0.00346Hz) and 1000 s is suitable for UAV tower flux calculations, respectively.Tower measurements have a greater influence of the higher frequency turbulent data, and take longer to reach an asymptote at the low frequency region; the narrower flux bandwidth for the UAVs is due to aircraft movement through the turbulent field.measured flux frequencies relative to tower measurements.
To investigate λE changes due to surface morphology, fluxes are calculated over T m for the entire run as a moving average with a δt in start time of typically 10 or fewer seconds.Such fluxes are not used in ensemble averages unless the turbulent integral length (L m ) is used as the interval between averages therefore maintaining statistical independence (Buzorius et al., 2006).
Flux analysis of tower data was performed using slightly modified aircraft algorithms to allow calculation of 1 s sliding average windows.Wind speed and direction, T and RH measurements were also available from a NASA tower situated ∼0.5 km upwind of the launch location.

Error analysis
There are a wide range of corrections and calculations one can apply to aircraft derived flux data to assess the limitations of this technique (Mahrt, 2010).The most important measurement is w -the lynchpin of scalar flux measurements.Using the methods of Garman et al. (2006) we can derive a minimum resolvable w of 0.17 m s −1 .For estimation closer to our application in the continental boundary layer, we here use the methods of Lenschow and Sun (2007), by first estimating the typical signal level required under the encountered experimental conditions from: TAS is the true airspeed, m s −1 .We estimate the peak signal magnitude, σ wm , and wavenumber, k m , from power spectra (Fig. 9) and derive a minimum requirements for w rate measurement of ∂w/∂t <0.052 m s −2 .To calculate measurement error from dominant sources in the system we adopt: and, where w p is the aircraft vertical velocity, and θ is the aircraft's pitch angle relative to the local earth plane (+ve for nose up).The first term is dominated by drift in the differential pressure transducer, the second term is a combination of INS/GPS pitch accuracy and drift in the measured attack angles.Error in TAS is assumed dominated by the 0.31 mb dynamic pressure error, and was generally <6 • .We use the manufacturer's stated pitch accuracy and a measured TAS of 28 m s −1 to compute the pitch accuracy to be <2.7 × 10 −6 m s −2 .Measured attack angle drift is difficult to quantify, but can be estimated by solving Lenschow and Sun's (2007) Eq. ( 10) also using σ wm of 0.1 m s −1 (from Fig. 6), and peak K m measured during this campaign.The absolute pressure sensitivity is calculated to be well within the required drift rate of 0.071 mb s −2 during flux legs, based on a static pressure transducer accuracy estimate of 0.036 mb.Systematic flux errors specific to each UAV flight are estimated from equations given in Mann and Lenschow (1994) by: where L is the length of the fight leg, L wq is the integral turbulent length scale inherent to each level of flight and derived via the relationship of L wq ≈ (L w L q ) 0.5 (Mann and Lenschow, 1994) where L w and L q are integral lengths of w and q respectively and calculated from autocorrelation functions (Lenschow and Stankov, 1986) to give L ws on the order of 100 and 195 m for the FT A and FT B respectively.Systematic flux error is then calculated to be on the order of 5 % for the measurements presented here.

System performance
Except for a DGPS failure due to an electrical fault, the flight and flux systems performed well on the 27th, enabling the successful collection of high-frequency data in both the morning and afternoon flights.Roll and pitch angles during flight runs mostly varied between 3 • to 5 • and 0 • to 10 • respectively (e.g.Fig. 6.Composite smoothed power spectra, kS(k), of unfiltered wind components (u, v, w) and water vapor, q, from flux legs at 520 m a.g.l.compared to surface spectra between 09:30-10:00 a.m., a period when ζ < 1.Similar spectral structures appear at both altitudes across the bandwidth with magnitudinal departures at the low and high frequency ranges.
According to the equations of Lenschow et al. (1986) corrected wind components (dependent on pitch, roll, α and β) are subtracted from the aircraft's vertical velocity to derive w. Figure 7b displays these components at the two altitudes samples in FT B ; the autopilot was able to maintain level runs resulting in domination of the derived w by the gust probe signal.The differences in scales of vertical motions is apparent -smallest scales dominate at the surface level, thermal events are seen at mid altitudes, whilst larger undulations were recorded at elevations at the highest measured altitude.

Regional and local meteorology
An eastwardly moving and deepening low-pressure system of 1008 mb (at 03:00 a.m.PST, 27 May 2010) situated 330 km NE of DFRC coupled with a strengthening oceanic high of 1016 mb, resulted in marine air being drawn in over land and up over the San Gabriel mountains.Winds increased in intensity over the course of the day.Hysplit back trajectory (Fig. 4a) shows the source of moisture.The day initially started off cloud free with cumulus cloud cover becoming increasingly congested towards the afternoon.
Temperature, q, wind speed and directional data measured by the NOAA, NASA and EAFB towers, and by the UAV are presented in Fig. 8.Here we can see an increase in temperature from 12.9 to 18.3 • C from 08:00 a.m. until 02:00 p.m., and wind speeds ranging from 3.9-5.5 m s −1 in the morning increasing to 20.2 m s −1 over the same period.Absolute water vapor concentration was slightly lower on the lakebed than at EAFB meteorological station, and displayed variability between the NOAA and NASA towers towards the end of the measurement period.The sonic anemometer on the NOAA tower allows the calculation of friction velocity, u * = (u w 2 + v w 2 ) 0.25 , Turbulent kinetic energy, TKE = 0.5[u 2 + v 2 + w 2 ]), sensible heat flux, H , and the Monin-Obukhov stability parameter, L = u 3 * /[κ(g/θ v )H ], where κ is the von-Kármán constant, g is acceleration due to gravity, and θ v is the virtual potential temperature.These parameters were calculated for a one second sliding window of averaging period 15 min from the NOAA tower data (Fig. 8b).H increased throughout the measurement period from 50 to reach a 200 W m −2 plateau by noon.This declined sharply in the afternoon to <100 W m −2 for a period between 02:30-03:00 p.m., possibly influenced by the increased cloud cover during this period.Friction velocity and TKE generally increase which corresponds with the increasing wind speeds.The dimensionless stability parameter, (defined here as ζ = −z/L, where z is the tower height of 10 m a.g.l. in this case), indicates the surface layer to be unstable throughout the day (ζ > 0) with periods of free convection (ζ > 1) occurring periodically.
Periods tending towards neutrality (ζ → 0) are apparent in the morning, and most notably during the second half of FT B , coinciding with a reduced H , and u * of >0.6 m s −1 .
Vertical profiles of potential temperature calculated from flight ascent data (Fig. 9a) indicate an inversion at z i = 505 m a.g.l. in the morning rising to 705 m a.g.l. in the afternoon.A super adiabatic layer present above the lakebed surface displays an increasing lapse rate in θ from 19.6 • C km −1 at 08:30 a.m. to 36 • C km −1 just before 01:00 p.m., above this layer the potential temperature indicates a near-neutral mixed layer present to just below the capping inversion.Absolute humidity (Fig. 9b) displays a decreasing gradient with altitude over the course of the measurement period of −1.7 g kg −1 km −1 .UAV measured wind speeds are reasonably consistent with altitude, and agree with the tower measured mean of 7.9 m s −1 in the morning, but approximately 1 m s −1 greater at altitude than the tower in the afternoon (Fig. 9c).Wind directions veer with altitude compared to the NOAA tower measurements before backing above z i (Fig. 9d).Also apparent in the profiles is the vertical broadening of the super-adiabatic layer over the course of the day which manifests as a near-neutral to slightly stable surface layer shown in the descent profile of the UAV following completion of FT B , at 02:10 p.m.

Flux Data
With allowance for aircraft turns, straight and level flux runs covered a mean distance of 7.7 km within the approximately 10 km maximum possible path in the designated UAV airspace.The mean time for runs in a north easterly direction was 161 s, compared with 414 s for SW runs travelling into the wind.Runs of less than 290 s in length were excluded from flux analysis, leaving the majority of the SW runs only.Of these runs, 2-3 were conducted at each altitude shown in Table 1, the maximum height reached corresponded to a maximum of 922 m a.g.l.
Averaged turbulent power spectra calculated using equivalent periods of tower and UAV data for the wind components u, v, w and q agreed more in both magnitude and form during periods where ζ < 1. Figure 6a-d displays averaged power spectra for the period 09:45-10:15 a.m., the period calculated when the airmass sampled in the 2nd and 3rd 520 m a.g.l.runs was advected past the tower; we consider this a reasonable comparative approach given the largely consistent vertical wind speed profile.These spectra adhere to the theoretical −2/3 slope expected for f S(f ) plot in the inertial subrange and also indicate significant contributions to the overall measured variance from lower frequency components; an expected result given the progressively larger energy containing turbulent eddies with altitude.This difference in measured scales ultimately hinders airborne/surface measurement comparison for the reasons mentioned above (i.e.insufficient UAV run length and the stationarity limits on tower averaging times), thus high pass filtering (HPF) is applied to the data to limit the influence of larger eddies, whose sampling was not statistically valid over the short flux run paths.The HPF used for UAV data is a 0.04 Hz Butterworth HPF and an equivalent (0.01 Hz) filter was also applied to the tower data.
We display the results of unfiltered and HPF UAV and surface data calculated for the UAV using T m = 290, a ten second sliding window interval in Fig. 10.In considering this plot it is important to bear in mind differences can be expected due to horizontal and vertical variability which are discussed later.In general, the unfiltered UAV data span a large range of values, particularly evident towards the afternoon, with a maximum range on the order of 150 W m −2 noted during the second 330 m a.g.l.run during FT B ; such a range is not uncommon in aircraft studies (Mahrt et al., 2001;Samuelsson and Tjernstrom, 1999;Song and Wesely, 2003).This λE variability is reflected in the surface data, particularly in the afternoon, which varies between −51 to 87 W m −2 towards the end of FT B .HPF data results in a much closer agreement with surface measurements, with maximum variations from the surface fluxes during runs A5 and B4 of 17.6 and 60 W m −2 , respectively and aside from these runs a mean maximum deviation of 5.1 W m −2 .The surface values themselves undergo sign changes following removal of the lower frequency contributions, resulting in near-zero or consistently positive upward fluxes more expected of the turbulent scale range.For the 930 m a.g.l.UAV unfiltered and HPF analysis, much lower variability is seen, particularly during FT A (σ 2 = 1.1 W m −2 for HPF data).This is consistent with decreased turbulence expected above the boundary layer inversion, and the increased afternoon variance (σ 2 = 12.2 W m −2 ) at this altitude is most probably due to increased penetration of this capping inversion by increasing thermal convection.At lower altitudes, and within the boundary layer, agreement with the surface data is also apparent, but with more variability.The simplest ideal of vertical water vapor flux divergence is an assumption of a monotonic decrease with height; fluxes produced at the surface layer decrease towards zero at the top of the boundary layer where typically drier air is entrained.In reality, aircraft and ground based profiling studies have found λE profiles to display no divergence, (Gioli et al., 2004) and monotonic decrease (Samuelsson and Tjernstrom, 1999), and also much vertical variability ascribed to entrainment rate changes and cloud effects (Giez et al., 1999;Linne et al., 2007;Kiemle et al., 2007).Average measured vertical λE profiles for FT A and FT B (Fig. 11), using HPF data and a time interval equivalent to L wq values given above, are likewise consistent with the boundary layer structure.The level of random noise is indicated by 1 σ .At the lowest level, 320 m a.g.l. the mean morning λE of 2.4 W m −2 corresponds with the surface flux of 0.9 W m −2 , whilst they are seen to decrease to around zero at the highest level in both flight profiles.The positive λE during FT A centered around 520 m a.g.l.displays a large influence due to the inclusion of the run A5 (Table 1), removal this run would result in a largely uniform profile although it is worth noting that at the time of measurement 520 m a.g.l. is both the capping inversion height and the calculated lifting condensation level (LCL), situated in the region of rapidly decreasing q with altitude (Fig. 9b).Although difficult to dissect using these measurements alone, this suggests the UAV location within the entrainment zone itself may be the reason for this result.The same thing however, cannot be said for the low level outliers during FT B , B3, which at an altitude of 330 m a.g.l.
(inversion and LCL were >705 m a.g.l.) were not influenced by such features.Closer towards the surface one can expect a greater contribution from surface elements as one moves below the mixing height, such features are also easier to investigate by using HPF data as it is the higher frequency events which are dominantly produced from surface elements (large eddies form in the region above the mixing height).Flux runs into the wind in the afternoon were at a slightly different angle and as such were able to squeeze a little more distance out of the UAV work area.When each run altitude λE for morning and afternoon is averaged into 100 m horizontal intervals and plotted against distance from SW to NE along the flight track (Fig. 12), these longer runs display a clear increase in λE with distance of ≈40 W m −2 for the 330 m a.g.l.runs, and reflected also in the 930 m a.g.l.run, providing a surface influenced reason for the increased variability discussed in reference to the flux time series.Although too short to really discern a pattern, FT A runs hint at this increase close to 5.6 km.The obvious candidate to explain this change is the edge of the lakebed, making the assumption that even a dry lakebed surface contains more moisture than the surrounding desert area.Extrapolating this step change to the  surface indicates an occurrence at around 4.8 km along the track, which is close to the central part of the lake bed where we could expect more water near the surface, further investigation by flux source modeling would aid this interpretation, but is beyond the scope of this study.

Co-spectral analysis
Normalized cospectral plots of HPF data indicate contributions to the flux from different frequency regions.These are presented in Fig. 13 for each altitude sampled during FT A for < q w > (calculated using an averaging time of 290 across all legs at the same altitude for UAV data, and using T A = 900 s for tower measurements, averaging across the period of runs A8-A10).These spectra demonstrate the general decrease with altitude of peak k from 23 m at the surface to close to 60 m at the highest altitude.This shift is expected due to surface limitations on the eddy size reaching the tower (Rissmann and Tetzlaff, 1994).The general form of (particularly the lowest two) spectra indicate a bimodal distribution in the fluxes, with this hump becoming less clear with altitude and the peaks becoming sharper.Conceptually, these cospectra are consistent with large scale motions of a convective boundary layer (Fig. 14, adapted from Shao (2008) supported by the presented measurements.The observed super adiabatic lapse rate is indicative of dry soils coupled with strong surface heating.Such thermal instability typically produces small plumes which merge into large thermals which rise to the top of the boundary layer, before sinking through the mixed layer, with a length scale limited by z i .These small plumes are seen in the cospectral plots, and also in the w time series.Similarly, although the larger eddies have been removed by HPF of the data, they are also discerned in the w time series plot (Fig. 7b).Occasionally, plumes may break through the inversion, but generally cause undulations.This lack of penetration is seen in the lack of variability along average flux track in the uppermost FT A leg, but strong surface heating around noon of the lake bed increases the ability of these plumes to break through the inversion, leading to the increase in variability in the afternoon flight seen in the plot of λE vs. distance over the lakebed (Fig. 12).The decrease in q with altitude in the vertical profiles (Fig. 9) reflects entrainment of drier air from above the inversion.The tower measurements are predominantly dominated by the smaller scale turbulent convective motions which rise up to form large eddies.Occasionally, when the neutral layer extends down to the surface, the tower sampling frequencies are more akin to those of the overlying boundary layer, resulting in power spectral agreement such as that in Fig. 6.
The lower frequency portions of the UAV co-spectra demonstrate both positive and negative periods, which can be expected with a limited horizontal leg extent; longer legs will reduce this variability by evenly sampling the large-scale atmospheric motions.In gaps between clouds one can expect down-draughts to dominate, with up-draughts dominating beneath clouds (Smith and Jonas, 1995).Based on the low-frequency portion of the < q w > co-spectra sampled in the neutral layer below the inversion we estimate the clouds to have an average spacing on the order of 440 m.Future experiments could investigate cloud effects on climate by the concurrent deployment of existing cloud microphysics and/or aerosol/radiation payloads.

Conclusions and further work
We have developed and demonstrated, for the first time a UAV platform based turbulent water vapor flux measurement system and have presented a snapshot of the recent development phase.Measurements of λE made in the lower troposphere, both below and above the capping inversion are in broad agreement with the meteorological situation of a convective boundary layer.High pass filtering of the data allow the removal of the majority of insufficiently sampled large eddies and, alongside spatial and temporal considerations, improves comparison with surface tower measurements to within 12 W m −2 for 86 % of flux runs.This consistency implies that the aircraft turbulent data collection system and data analysis are sufficient to measure λE for the conditions encountered here.We aim to continue the development of the data collection and particularly processing aspect of this system to continue learning from the extremely informative and detailed work undertaken by the manned aircraft flux community (Khelif et al., 1999).For example, airborne calibrations which were not undertaken due to the shortened nature of FT B would further enhance the probe characterization.Similarly, a Kalman filter system designed to implement the DGPS data in w derivation may improve accuracy.
Another important factor for the use of this system is gaining the flexibility to perform long flux runs which would more adequately sample large eddy flux contributions.Sampling strategies are a challenge when attempting to make multiple runs with sufficient length within the available flight time.Sending real-time data such as q and θ over the communication link would allow for in-flight changes by providing situation dependent data.Steps have been taken to achieve this, and it is also within the capability of the CRIO system to calculate real-time averaged cospectra which can be used to assess flux scales across the depth of the boundary layer in real time.We also strive to add a fast response temperature probe to the system to enable derivation of the sensible heat flux and stability information.

Fig. 1 .
Fig. 1.Image of the manta UAV on Roger's lake bed in preparation for afternoon take off on the 27th May 2010.Also shown are the water vapor flux instruments on the Manta (inset left) and the surface tower (inset right).

Fig. 3 .
Fig. 3. Vehicular test results indicating (a) time series plots of sonic anemometer and gust probe vertical wind, and (b) calculated fluxes along truck measurement track.

Fig. 4 .
Fig. 4. (a) Regional view looking west of the NASA Dryden Flight Research Center (NDFRC) situation displaying the tower location (Triangle), flight positional data from FT A , and HYPLIT Modelled back trajectories ending at 1055, 1445 and 1645 m a.g.l., 12:00 p.m. PST.The distance from HYSPLIT site end points to the coast is approximately 115 km.(b) Plan view of FT B flight path, colored according to water vapor concentration.Both figures use Google earth imagery.
Fig.7a).The GPS/INS demonstrated excellent agreement with the flight computer system; greater variability in the vertical velocity of the GPS/INS was expected, indicative of its higher sensitivity.The GPS/INS maintained an average of eight satellites during flights.
Fig. 7. (a) Time series comparison between on-board flight computer attitudes and ground velocities with those measured by the GPS/INS device during 20 min of FT A resampled to 1 Hz.(b) Eight minutes vertical velocities at two altitudes during FT B , showing 50 Hz Airplane INS and gust probe data used to derive w.Tower data are also shown.High frequency fluctuations dominate the surface data.Thermal events are clearly seen at 330 m a.g.l., and lower frequency signals are present above the boundary layer at 930 m a.g.l.Gust probe data are inverted for clarity.
Fig. 8. (a) Time series of surface and UAV meteorological data.(b) From top to bottom: H , u * , TKE and ζ derived from the sonic anemometer data on the surface tower.For a confirmation calculation techniques, H is shown calculated with the standard NOAA and SCRIPPS UAV algorithms.

Fig. 9 .
Fig. 9. UAV ascent profile plots averaged over 20 m intervals during FT A (open symbols, yellow, 09:13 a.m.-11:54 a.m.)/FT B (closed symbols, red, 12:48 p.m.-02:12 p.m.) for (a) potential temperature θ (also shown is the descent profile of FT B ), and (b) water vapor mixing ratio q; shaded areas denote 1-σ .Plots (c) and (d) display flux run averages of wind speed and direction.Ground station measurements corresponding to the launch periods are included in each plot as well as EAFB 03:00 a.m.(1100 z) radiosonde data.Inset plots are of radiosonde (a) derived θ and (b) q for up to 6 km altitude.

Fig. 10 .Fig. 11 .
Fig.10.Time series of high pass filtered and unfiltered UAV and tower λE data.UAV λE are calculated with T m = 290 s and a 10 s sliding period, colored according to altitude.Tower data are calculated with T m = 1000 s also for 10 s intervals.High Pass Filtering (HPF) data (0.04 Hz for the UAV and 0.01 Hz for the tower), removes poorly sampled low frequency contributions to the measured flux and results in closer agreement between the two measurement techniques and horizontal spatial scale differences contribute to the residual variability.

Fig. 12 .
Fig. 12. Mean of λE of flight runs in a SW direction at each altitude for FT A and FT B .North east is towards 0 m and the position of the tower along the track is shown.The longer afternoon runs indicate an increase in λE at approximately 5 km, hinted at in FT A data, and most likely reflecting the change in surface properties from the surrounding drier desert to the lakebed surface.

Fig. 13 .
Fig.13.Averaged Cospectral plots of < w q > normalized to the total variance from UAV FT A and surface data.Tower data are for a 30-min period from 11:17.

Fig. 14 .
Fig. 14.Structure of the entrainment zone capping the convective atmospheric boundary layer.Horizontal scale is on the order of 10 km.Adapted from Shao (2008).

Table 1 .
Properties of straight and level southwesterly runs during Fts A&B with t > 300.Showing start time (t start ), duration t, Direction, surface velocity U g , horizontal length x, altitude z, and corresponding windspeed, U tow , measured during the run at the NOAA tower.