Interactive comment on “ Application of Gauss ’ s Theorem to quantify localized surface emissions from airborne measurements of wind and trace gases

The paper by Conley et al. presents and validates a technique to infer point-source emission rates from in-situ aircraft observations of the atmospheric concentrations in a cylindrical volume around the source. The techniques is approached theoretically through LES modelling as well as experimentally through the analysis of actual aircraft observations. The observation part is, however, rather short deserving more detailled discussions.


Introduction
Accurate national inventories of greenhouse gas emissions (primarily carbon dioxide -CO 2 ; methane -CH 4 ; and nitrous oxide -N 2 O) is of paramount importance in developing strategies to understand global emissions.The multitude of sources, however, are so often highly variable in area, emission magnitude, height above ground, and duration that rigorous verification is exceedingly difficult.Nevertheless, measurement techniques have improved markedly in the past decade, and these are being employed to an unprecedented extent in an effort to evaluate and refine emission inventories (Nisbet and Weiss, 2010).Most so-called "bottom-S.Conley et al.: Application of Gauss's theorem to quantify localized surface emissions up" inventories are developed by aggregating statistical correlates of individual process emissions to such mapping variables as population density, energy consumption, head of cattle, etc., extrapolating to total emissions using a relatively small number of direct measurements.On the other hand, atmospheric scientists have long striven to use measurements from global surface networks, aircraft campaigns, and satellites to try to determine emissions based on the amounts and build-up rates of observed trace gases.Aircraft and satellites, the "top-down" approach, conveniently integrates the multitude of sources, but is heavily reliant on a detailed knowledge of atmospheric transport.Top-down methods also suffer from difficulties attributing sources and generalizing measurements made over a relatively short time period.Attempts to reconcile these two distinct methods on global (Muhle et al., 2010) and continental scales (Gerbig et al., 2003;Miller et al., 2013) have often indicated an apparent underestimation by the bottom-up methods of a factor 1.5 or more.
In principle, the aircraft top-down measurements can be conducted at all the atmospheric scales to better understand and identify the emissions at comparable scales.For longlived greenhouse gases, which readily disperse throughout the atmosphere, the global scale is very instructive.The seminal experiment began with Keeling's acclaimed CO 2 curve (1960), and has continued through more contemporary techniques by Hirsch et al. (2006) and Neef et al. (2010) for CH 4 and N 2 O, respectively.At progressively smaller scales more details of the source strengths and apportionment can be made: from synoptic or continental scales which can help constrain national inventories (Bergamaschi et al., 2005) or specific biogeographic regions (Gallagher et al., 1994), to mesoscale investigations that estimate emissions from urban areas (Mays et al., 2009;Turnbull et al., 2011;Wecht et al., 2014) or specific oil-and gas-producing fields (Karion et al., 2013;Petron et al., 2014) and even down to individual point/area sources on the order of 10-100 m size (Denmead et al., 1998;Lavoie et al., 2015;Roscioli et al., 2015).
Aircraft in situ measurements are particularly useful for top-down methods at the sub-mesoscale because they can be used to measure the air both upwind and downwind of a source region.However, deployments tend to be costly and thus sporadic.As far as we know, the aircraft methods used so far can be categorized into three types.First, there is the eddy covariance technique that is carried out at low altitudes wherein the vertical fluxes of gases carried by the turbulent wind are measured by tracking rapid fluctuations of both concentrations and vertical wind (Hiller et al., 2014;Ritter et al., 1994;Yuan et al., 2015).This method is generally thought to be the most direct, but it is limited to small footprint regions which must be repeatedly sampled for sufficient statistical confidence, requires a sophisticated vertical wind measurement, and can be subject to errors due to flux divergence between the surface and the lowest flight altitude and acceleration sensitivity of the gas sensor.The second and by far the most common approach is what chemists usually refer to as "mass balance" and what is known in the turbulence community as a "scalar budget" technique.Many different sets of assumptions and sampling strategies are employed, but the overall goal is to sample the main dispersion routes of the surface emissions as they make their way into the overlying atmosphere after first accumulating near the surface.The scales that can be addressed by this method are from a few kilometers (Alfieri et al., 2010;Hacker et al., 2016;Hiller et al., 2014;Tratt et al., 2014) to tens of kilometers (Caulton et al., 2014;Karion et al., 2013;Wratt et al., 2001) to even potentially hundreds of kilometers (Beswick et al., 1998;Chang et al., 2014), and this approach has been the focus of recent measurements in natural gas production basins.These basins present a source apportionment challenge in that emissions from multiple sources (agriculture, oil and gas wells, geologic seepage, etc.) commingle as the air mass travels across the basin.The third method of source quantification is to reference measurements of the unknown trace gas to a reference trace gas with a metered release (tracer) or otherwise-known emission rate and assume that the tracer and the scalar of interest have the same diffusion characteristics.Typically this tracer release technique is applied to small scales of tens to hundreds of meters (Czepiel et al., 1996;Lamb et al., 1995;Roscioli et al., 2015), but the principle has been attempted at the basin (Peischl et al., 2013) and continental (Miller et al., 2012) scales using a reference trace gas with a suitable known emission rate such as CO 2 or CO.
Here we describe a new airborne method borne out of a necessity to identify and quantify source emissions to within 20 % accuracy in a large heterogeneous field of potential sources.The novel technique applies an aircraft flight pattern that circumscribes a virtual cylinder around an emission source and, using only observed horizontal wind and trace gas concentrations, applies Gauss's theorem to estimate the flux divergence through that cylinder.By integrating the outward horizontal fluxes at each point along the circular flight path, the flux contributions from enclosed sources can be accounted for.Making an accurate estimate, however, requires the selection of an appropriate circling radius based on the micrometeorological conditions inferred in flight from measurements onboard the aircraft.The pattern must be far enough downstream for the plume to mix sufficiently in the vertical, yet not so far that the trace gas plume enhancements do not stand out sufficiently from the background concentration.
In this study we first present the general analytical method used to derive emission estimates using airborne measurements.Next, we investigate the structure of a generalized dispersing plume using large eddy simulation (LES) to better understand the optimal sampling strategies for quantifying near-surface gas sources.Because the wind fields of turbulent flows cannot be predicted in detail, we do not attempt to compare specific features of our observations with specific LES results, but rather we use the numerical experiments to guide the development of the observational methodology.For example, by investigating the LES flux divergence profiles in the layer below the lowest flight altitude, we are able to estimate the contribution of this unmeasured component to the overall source strength.We then evaluate the accuracy of the approach using coordinated planned release experiments and by applying the method to CO 2 emitted from several power plant plumes to compare with reported emissions.
2 Data collection

Airborne instrumentation
The airborne detection system is flown on a fixed wing single-engine Mooney aircraft, extensively modified for research as described in Conley et al. (2014).Ambient air is collected through ∼ 5 m of tubing (Kynar, Teflon, and stainless steel) that protrudes out of backward-facing aluminum inlets mounted below the right wing.In situ CH 4 , CO 2 , and water vapor are measured with a Picarro 2301f cavity ring down spectrometer as described by Crosson (2008), which is operated in its precision mode at 1 Hz.In situ ethane (C 2 H 6 ) is measured with an Aerodyne methane/ethane tunable diode infrared laser direct absorption spectrometer (Yacovitch et al., 2014).There is a 5-10 s time lag in both analyzers that depends on the flow rate and tubing diameter.We use a 1/8 in.OD (3.175 mm) stainless line for the Picarro (∼ 0.2 slpm flow rate), and a 6.3 mm (1/4 in.)Teflon line for the CH 4 /C 2 H 6 spectrometer (∼ 4 slpm flow rate).This results in lag times of ∼ 5 s for the Aerodyne and ∼ 10 s for the Picarro.The lag time for the Picarro is calculated using a "breath test", whereby we exhale into the air inlet and measure the time required for the CO 2 measurement to peak.The ethane lag time is adjusted to maximize the correlation between the ethane and Picarro methane time series in plumes where both gases are emitted.Both lag times are slightly dependent on pressure, i.e., with a typical altitude change of ∼ 1 km, the change in lag time is less than 10 %, and is inconsequential when applying this method within a few hundred meters from the surface.The horizontal wind speed and direction, sampled at 1 Hz, is based on a standard aircraft pitot-static pressure airspeed measurement and a dual GPS compass that determines aircraft heading and ground speed.The accuracy of the horizontal wind measurement is about 0.2 m s −1 (Conley et al., 2014).The horizontal wind is cal-ibrated periodically by flying ∼ 5 km L-shaped patterns in the free troposphere; a heading rotation and airspeed adjustment is made to the wind calculation to minimize the dependence of the wind on aircraft heading.These adjustments typically amount to less than 2 • rotation and 3 % adjustment of the airspeed.In flying the tight circle patterns described below, the pilot does not adjust the rudder trim to use the same calibration coefficients in the wind measurement calculation throughout the flight.

Large eddy simulations
In order to study the plume behavior of surface emissions as it relates to sampling in the stacked circles, we use the LES module of WRF V3.6.1.WRF-LES explicitly resolves the largest turbulent eddies by filtering the Navier-Stokes scalar conservation equations at some scale in the inertial subrange, and allowing the smaller motions beyond the cutoff to be modeled using a sub-grid (also called a sub-filter) scale turbulence parameterization that is based on properties of the larger-scale, resolved flow.Because the aircraft data is typically sampled at 1 Hz and the true airspeed is around 70 m s −1 , we use an LES horizontal grid size roughly half (40 and 50 m) the distance between aircraft data samples.Because periodic lateral boundary conditions are imposed on the WRF-LES variables, care must be taken to ensure that the effluent does not reach the lateral boundaries of the simulation domain.On the other hand, WRF-LES does not allow for parallelized computation, making the simulations quite expensive in terms of computation time.We therefore struck a balance between a large enough domain in horizontal extent (6 and 8 km) such that the effluent would not reach the downwind boundary before the end of our simulation, while maintaining a grid size small enough to resolve scales of the aircraft observations.The vertical domain needs to be large enough to encompass a developing convective boundary layer (CBL), while at the same time containing substantial free tropospheric flow above to serve as a reservoir that can feed momentum and free-tropospheric scalars to the CBL.Moreover, the stable region (potential temperature lapse rate dθ/dz = 5 C km −1 ) between the CBL inversion base and the top of the domain had to be large enough to damp any wave activity before it could reflect off the upper boundary and create spurious motions throughout the domain.
The standard WRF-LES module is not set up to allow for effluent release, so we implemented a modified version of the WRF source code (S.-H.Chen, personal communication, 2015) that includes a surface effluent release with a specified position and release rate.Three different convective simulations were run with varying resultant mean wind speeds in the boundary layer, and each was allowed 4-5 h spin-up dynamically before the effluent was released at a rate of between 2.9 and 3.5 kg h −1 .The exact release time was selected to give reasonably stationary CBL depths and turbulent ki-  netic energy.The conditions for the three simulations are listed in Table 1, and based on the different wind speeds they range from moderate to strongly convective boundary layers (-z i /L MO from ∼ 50 to ∼ 200, where L MO is the Monin-Obukhov length and z i is the CBL depth.) 3 Methods

Theory of measurement using Gauss's theorem
We use an integrated form of the scalar budget equation for a passive, conservative scalar in a turbulent fluid to estimate the emission of a gas of interest within a cylindrical volume V .The volume is circumscribed by a series of closed aircraft flight paths (typically circular) flown around the emission source over a range of altitudes.The altitudes encompass the lowest safe flight level (usually 75-150 m a.g.l.) up to an altitude where no discernable change in the trace gas mixing ratio, χ, is observed around the flight loop, z max .The scalar in our case is the mass concentration (i.e., density of a chemically unreactive species in a turbulent flow field, u = ui + vj + wk; its Reynolds decomposition is c = C + c , where C is the mean concentration around each loop and c is the departure from the loop mean.Figure 1 shows an actual example of the effluent sampled by the aircraft in a sequence of stacked paths l that circumscribe an area, A, enclosing the source in a volume, V .The effluent is carried downwind as it mixes upward in the CBL.A virtual surface circumscribed by the circular flight tracks is assumed enclosing the source and extending above the vertical extent of the plume so that there is no vertical transport above that level.To estimate the source strength, we start with the integral form of the continuity equation: where denotes an average over the volume V , Q c is the sum of the internal sources and sinks of c within V , and m is the total mass of c within the volume V .At this point, we recognize that the flux divergence is composed of two terms (2) In Sect.3.2 we perform a scale analysis of the terms on the right-hand side (rhs) of Eq. ( 2) and show that the second term, which is proportional to the horizontal wind divergence, may be neglected under our normal flight protocol.This is fortunate because of the difficulty in accurately estimating the horizontal wind divergence from aircraft measurements (Lenschow et al., 2007).The vertical flux across the top of the flight cylinder is assumed to be zero, and the flux from the bottom (ground) is the surface source we are measuring.This leaves us with only the horizontal flux, i.e., cu h , where u h (= ui + vj ).In order to minimize the contribution from the horizontal wind divergence term, we remove the loop mean concentration, C, which does not alter the first term on the rhs because ∇C = 0, so that Eq. ( 2) becomes Next, we use Gauss's theorem to relate the volume integral to a surface integral around the volume that is sampled by the aircraft flight loops: where S is the surface enclosing V and n is an outward pointing unit vector normal to the surface.
The surface integral can be broken into three elements: a cylinder extending from the ground up to a level above significant modification by the emission, the ground surface circumscribed by a low-level (virtual) circular flight path (z = 0), and a nominally horizontal surface circumscribed by a flight path above the level modified by the source (z = z max ).We assume there is no significant flux (other than the source of interest) into or out of the ground.Next, the surface integral is estimated solely from a sequence of closed path integrals measured by the aircraft at multiple flight levels to estimate the right side of Eq. ( 5) (blue dashed lines in Fig. 1), where l is the flight path.Combining Eqs. ( 4) and ( 5) leads to the result that is the basis for this measurement technique where a series of horizontal loops at different altitudes are flown around a source region: Along each path the instantaneous outward flux is computed and summed over the loop to yield the mean flux divergence via Gauss's theorem.A temporal trend of the total mass within the volume ( ∂m ∂t ) can be estimated from the flight data and added to the flux divergence integral to obtain the emission rate.

Divergence uncertainty
In order to estimate the relative error in the horizontal divergence term that we are eliminating, we perform a scale analysis of the relative size of the two terms in Eq. ( 2), using some typical values of the CBL parameters (convective velocity scale w * = 1 m s −1 , boundary layer depth, z i = 1000 m), and sampling geometry (flying at a radius of 1 km around a point source.)The Taylor (1922) statistical theory of dispersion in a homogenous and stationary turbulent fluid predicts that the root mean square lateral (σ y ) and vertical (σ z ) dispersion parameters increase linearly with time, or equivalently advection distance, downwind in the near field.Weil (1988) shows several examples of the growth of both of these parameters downwind to be ∼ 0.5w * , which we use here for a rough estimate of a conical plume spreading to quantify the dilution of the source's emission as it travels downwind to be intercepted by the aircraft.We use a large background mixing ratio characteristic of global CH 4 (∼ 1.9 ppmv), estimate the mean gradient by the plume concentration divided by the distance downwind, and assume a conservatively large horizontal wind divergence of 10 −5 s −1 , which may in fact be typical for our small sampling region (Stull, 1988).The results are shown in Fig. 2 and, for all but the smallest sources of a few kilograms per hour and wind speeds below 1 m s −1 , the divergence term is at least an order of magnitude smaller than the gradient term.

Applying the theory to the LES results
We calculated a comparable estimate of Q c in the LES domain from the air density, concentration, and wind along circular flight paths as a virtual aircraft would fly.Willis and Deardorff (1976) generalized results of their convection tank experiments to downwind dispersion in the convective boundary layer (CBL) in terms of a dimensionless length scale X, the ratio of the horizontal advection time to the large eddy turnover time: where x is the downwind distance and U is the vertically averaged mean wind speed.3), as a function of wind speed and source magnitude for methane, using a typical global background of 1.9 ppm and divergence of 10 −5 s −1 .
Figure 3 shows the crosswind-integrated concentration profile for the plume release in the UCD50B WRF-LES run as function of X, and normalized height, Z = z/z i .Because of the time limitation due to the periodic boundary conditions, the plume is averaged for only ∼ 15 min of simulation time which is just under a large eddy turnover time for the conditions of the run.The results displayed in Fig. 3 are in good qualitative agreement with the results of Willis and Deardorff (1976) and Weil et al. (2012), save for the release being at the surface in our LES study, and at Z = 0.067 for the above studies (see Figs. 1 and 2 of Weil et al., 2012).Figure 3 shows the maximum concentration being lofted near X ∼ 0.2 and leveling off near Z ∼ 0.8 around X ∼ 0.6; beyond X > 1.5 the plume is fairly well mixed throughout the extent of the boundary layer.

The upwind-directed turbulent flux
Horizontal turbulent fluxes are generally ignored in boundary layer budget studies due to the fact that while they are often sizeable in magnitude they do not change significantly over horizontal length scales under consideration (the horizontal homogeneity assumption).In the vicinity of a point source, however, this is not likely.The method outlined here estimates source emissions using a measured horizontal flux that incorporates wind and scalar measurements at 1 Hz sample rate, resolving scales of ∼ 70 m (Conley et al., 2014), which should include nearly all of the turbulent contributions to the horizontal flux.Here we consider the nature of this turbulent flux and the error in emission estimates if only the mean transport were considered.We start with the budget equation for a horizontal scalar flux in a horizontally homogeneous turbulent flow where the molecular diffusive/viscous term has been neglected (Wyngaard, 2010), where ρ is density and p is the pressure fluctuation.We then assume stationarity and integrate across the source from a point just upwind to a point within the plume and obtain We further assume that although the scalar field is not homogeneous the flow field is, and the background horizontal c-flux upwind is much smaller than the flux induced by the point source.This results in an equation for the in-plume flux The first three terms on the rhs of Eq. ( 10) are negative within the plume with their largest magnitudes on the upwind side and diminishing downwind.On the largest, boundarylayer-filling eddy scales, the mean concentration of C downwind of a source is greater than in the upwind region, (C x − C bckg ) > 0, and therefore the first term is negative, but decreases in magnitude with distance downwind.However, this term is also positive on smaller scales within the plume where the mean gradient is directed upwind towards the source, and is most likely responsible for the specious intuitive impression that the horizontal turbulent flux should transport the plume downwind from the source along with the mean wind advection.Moreover, the second and third terms on the rhs are negative because the momentum flux, u w , and mean vertical gradient, ∂C/∂z, are negative while the vertical turbulent flux, c w , and wind shear, ∂U/∂z , are positive.
Based on the vertical concentration profiles shown in Weil et al. (2012;their Figs. 3 and 4) it can be inferred that the vertical concentration gradient, ∂C/∂z, changes from negative to positive near X ∼ 1 and becomes negligible for X > 2-3.Similarly, in the third term, the vertical flux, c w , decreases with fetch.Thus the counter-directed flux (c u < 0) will fade with distance downwind.Wyngaard et al. (1971) have shown that the third-moment, turbulent transport terms (4 and 5 on rhs of Eq. 10) in the horizontal heat flux equation are small in the surface layer compared to the source terms, so we assume the same holds for this scalar flux.Finally, the remaining pressure covariance term is believed to be the main sink in the budget equation working to decorrelate the wind and the scalar as was shown in the surface layer measurements of Wilczak and Bedard (2004).Therefore, the dominant production terms for negative c u (terms 1-3 on the rhs of 10) must be balanced by the pressure-correlation term leading to an upwind-directed horizontal turbulent flux within the plume that decreases in magnitude in the downwind direction.This conclusion is supported by several previous studies.For example, in a wind-tunnel study of flux-gradient relationships, Raupach and Legg (1984) reported that the mean stream-wise horizontal heat flux calculated by multiplying the mean wind by the mean temperature overestimates the total heat flux by approximately 10%, which suggests that the turbulent component of the horizontal heat flux is negative; that is, the turbulent flux is upwind, directed counter to the mean flow.Other researchers have reported an even larger disparity.Field experiments by Leuning et al. (1985) indicate that the horizontal turbulent flux of a trace gas is ∼ 15 % of the mean flux, while Wilson and Shum (1992) suggest it may be 20 %.A recent LES study of particle dispersion over a plant canopy by Pan et al. (2014) indicates magnitudes of 20 % or more for the negative turbulent component of scalar fluxes in the vicinity of the source and decreasing with downwind fetch.We therefore conclude that when sampling a near-surface point source at X of order unity or less, if only the mean concentration difference is measured, a significant overestimate of the scalar source is likely to occur.4: there is an obvious peak at the mean flight loop frequency (usually ∼ 100 s period) followed by a smaller negative dip at higher frequencies within the meandering effluent plume.We believe this to be good evidence that our method captures this important component of the overall flux away from the source, which cannot be obtained with a traditional mean wind and an integrated concentration enhancement measurement that is so often employed in airborne source estimates (Ryerson et al., 2001;White et al., 1976).

Choosing the downwind sampling distance
Determining the optimal sampling distance from the targeted point source is a balance of several factors.First, not surprisingly, the largest plume signal occurs closest to the source (Fig. 3).Second, a high degree of confidence in the results is contingent upon sampling the majority of the plume at and above the lowest flight altitude, which only occurs downwind after a sufficient time has elapsed to loft the initially nearsurface plume.And third, an attempt is made to sample the plume before it reaches the top of the boundary layer so that the vertical turbulent entrainment flux does not become appreciable violating the assumption of negligible flux through the top of the volume V as discussed in Eq. ( 2).Finally, close to the source, the fluctuations in concentration will be very large, intermittent, at small scales, and highly variable.To gain further insight into the second feature of the dispersing plume, Fig. 5 shows the average horizontal flux divergence profiles derived from the three WRF-LES runs.
Here we discuss a dimensionless R, which is identical to X, to emphasize that this scaled downwind distance from the source is a radius of a flight loop.The flux divergence values are made dimensionless by the boundary layer height, z i , and the source emission rate, Q. Very close to the source, before the plume has had a chance to loft, the flux divergence profile exhibits a strong gradient below the minimum safe flight altitude, making that term difficult to measure directly, as shown in Fig. 5. Farther from the source, the signal becomes weaker with increasing altitude and eventually becomes increasingly influenced by entrainment fluxes.We therefore seek a sampling distance that is far enough to allow sufficient vertical lofting yet close enough so that plume crossings are easily observable against the background variability and instrument noise, and are not yet influenced by entrainment mixing.
Based on the simulation results presented in Fig. 5, we see the gradient below the lowest flight safe altitude typically becomes very small for R > 0.4, and therefore we attempt to target that distance to minimize the extrapolation error from the flight data to the surface.We do not currently measure all the necessary parameters to estimate R in flight (primarily the surface heat flux (w θ v ) 0 which is required to estimate w * ).Instead, we estimate w * based on the observed boundary layer height, standard deviation of wind speed, and a parameterization for w * = σ u /0.6 (Caughey and Palmer, 1979).By 15 laps, the emissions estimate (blue line) has stabilized to 2.5 kg h −1 compared to the actual leak rate (red line) of 2.9 kg h −1 .Dimensionless distance R = 0.25, 50 realizations.Gray area represents 1 SD (standard deviation).

Minimum number of passes
The atmospheric boundary layer is a turbulent medium, meaning that two passes across a plume at the same altitude and distance downwind will likely make very different measurements of the trace gases of interest.A natural question arises as to how many passes are required to develop a statistically sound estimate of the emission rate.We investigate the number of passes required to obtain a statistically robust estimate using the WRF-LES results and a controlled release experiment.By calculating the horizontal flux divergences with a virtual airplane flying through the simulated tracer field, and then randomly sampling the flux divergences from each of the legs and plotting the resultant estimated emission rate as a function of the number of samples used, we obtain the results presented in Fig. 6.The gray region around the red line mean represents the standard deviation of estimates based on a random set of loops.Figure 7 shows results from an analysis of actual flight data from the controlled ethane release test near Denver, Colorado, on 19 November 2014.It is evident from both the simulation data and the field data that a statistically stable estimate seems to be achieved somewhere between 20 and 25 loops around the source.

Discretization and altitude binning the flux divergence data
Measurements of the relevant scalars (e.g., CH 4 ) and meteorological variables are sampled at discrete time intervals.
For our analyses, we interpolate all measurements including GPS (3 Hz), methane (1 Hz), and temperature (1 Hz), and computed variables including horizontal wind (1 Hz), and air density (1 Hz) onto a synchronous 1 Hz time series.Next, we estimate the path integral for each individual loop of the flux normal to the flight path by summing up the flux contri- butions times the sample length around each loop and then summing over the height intervals, where ρ is the scalar air density, u n is the wind speed normal to the flight path, s is the distance covered during the 1 s time interval of each measurement, and L is the distance covered in one complete circuit.The outer summation sums each of the discrete vertical laps from the bottom (z = 0) to the highest lap (z = z t ).If all laps were sampled at equidistant altitudes, the total divergence could be calculated as the average divergence of all laps multiplied by the top altitude.However, because there is greater horizontal transport and variability at lower altitudes, as demonstrated by the widening standard deviations approaching the surface in the theoretical flux divergence profiles shown in Fig. 4, more sampling laps at lower altitudes increase the statistical validity of the largest horizontal transport values.To ensure that all altitudes are nearly equally weighted, we divide the vertical range into six equally spaced bins, save for the lowest bin which is extended to the surface, and then average the measurements from the laps within each bin.The total emission is the sum of the flux times the path length in each bin multiplied by the bin width.We also performed six flights where we sampled equally at all altitudes to allow for comparison of the direct average vs. the binned results, and in all of these flights the values derived by the two methods agreed to within 5 %.

Error analysis
Our method assumes a stationary emission source.The legto-leg variability is primarily driven by the stochastic nature of turbulence (e.g., we may sample the plume on one lap and miss it on another).By aggregating the laps into vertical bins, we can use the standard deviation of the horizontal fluxes within each bin as an estimate of the uncertainty within that bin.Then the total uncertainty in the estimate of the flux divergence is simply estimated by adding up the individual bin uncertainties in quadrature.The first term on the rhs of Eq. ( 6) is the time rate of change of the scalar mass within the cylindrical flight volume.This storage term is estimated by performing a least squares fit of the methane density with time and altitude.The resulting uncertainty in the time rate of change is then combined (summed in quadrature) with the uncertainty from the altitude bins to achieve a total uncertainty in the measurement.

Results and discussion
We use measurements from three sets of flights to characterize the accuracy of this estimation method.We flew 2 days measuring an controlled ethane release provided by Aerodyne Research, Inc., 4 days measuring a controlled natural gas release provided by the Pacific Gas & Electric Company (PG & E), and six power plant flights where our estimates are compared with reported hourly power plant CO 2 emissions.

Controlled ethane releases
Two experiments with known/controlled ethane releases were performed in collaboration with the Aerodyne Mobile Laboratory team.Pure ethane was released and measured with a flowmeter by the Aerodyne ground crew.The Colorado site (November 2014) was in a remote area approximately 170 km NE of Denver.This site was chosen because of the flat terrain and lack of other nearby ethane sources that could pollute the controlled release plume.The flux profiles for both releases are shown in Figure 8 and indicate that, in both cases, the aircraft successfully flew above the ethane plume (measurements tend toward zero with increasing altitude).An example of an individual lap is shown in Fig. 9 and indicates a clear plume signal downwind of the release.
As the aircraft climbs, eventually the signal disappears, as shown in the figure.Agreement was excellent, with the estimated emission just 2 % over the actual controlled release rate.The second Aerodyne controlled release in Arkansas on 3 October 2015 was performed at a site surrounded by nearby emission sources and an elevation change (∼ 70 m) within the aircraft flight path.The aircraft-derived ethane emission estimate was 25 % higher than the controlled release rate and the calculated uncertainty was significantly higher than on other sites (Table 2).A significant upwind ethane source was observed during the Arkansas experiment.This source was evident on roughly half of the upwind passes, suggesting that techniques which rely on a limited number of upwind passes to characterize the background could have a large random error and thus erroneously estimate the upwind source strength.A similar problem would affect those techniques that employ a downwind transect, using the edges of that transect lying outside the plume to estimate the background concentration.These observations demonstrate the complication (and bias) that can arise from nearby sources.Since this method integrates all the emission sources in the area within the flight circle and a small distance upwind of the circle depending on the vertical mixing, estimates from Gauss's method may be biased high if there are sources within that area.The average error of the two ethane releases is 13 %.

Controlled natural gas releases
In conjunction with PG & E, we performed two sets of 2day ground-level controlled release experiments from existing PG & E facilities, exactly 1 year apart.The first set was performed southeast of Sacramento near the town of Rio Vista, CA at the Rio Vista "Y" station and the second set near Bakersfield, CA.For the Rio Vista test, the release rate was not calibrated with a flow meter but, based on the size of the orifice and the upstream pressure, the release rate was estimated at 15.2 ± 1.5 kg h −1 .This release rate is an estimate of the total gas being released which is a combination of primarily CH 4 and C 2 H 6 .We use the regression fit of ethane to methane (averaging 0.085 by mass) to estimate the actual release rate of each scalar.
In comparison with the controlled C 2 H 6 release, controlled CH 4 releases suffer from the effect of small enhancements relative to the background concentration.During the Rio Vista release, the largest enhancement that we measured was 100 ppb, with 30-40 ppb being typical.Using a typical background level of 1.9 ppm, a 40 ppb enhancement represents 2 % of the background.In contrast, for ethane the enhancements are as large or larger than the background.The results of the controlled methane release tests are shown in Table 3 and indicate aircraft-derived estimates within 17 % of the controlled release rate.This large background results in increased uncertainty in the emission calculation.The average difference between the estimated emission and the calculated flow rate is 13 %.

Power plant flights
Power plants in the US are required to report CO 2 emissions to United States Environmental Protection Agency (EPA; https://ampd.epa.gov/ampd) on an hourly basis.The accuracy of the reported CO 2 emissions has been determined to be ±10.8-11.0% based on reported US average differences between Energy Information Administration (EIA) fuel-based estimates and EPA continuous emission monitoring-based estimates (Ackerman and Sundquist, 2008;Peischl et al., 2010;Quick, 2014).Also, Peischl et al. (2010) determined an accuracy of power plants reporting CO 2 emissions in Texas of ±14.0 % based on differences between observed downwind SO 2 / CO 2 and NO x / CO 2 emission ratios and those reported via EPA continuous emission monitoring (Peischl et al., 2010).Here, we use the slightly larger uncertainty from Peischl et al. (2010).Power plant emissions are "hot" gases and very buoyant, in contrast to a surface emission source that is typically not buoyant.An additional uncertainty arises from temporal emission variability (hourly averaged reported CO 2 emissions vs. < 1 h power plant flights that may cover parts of two reported consecutive hourly values).We estimate the total reported uncer-tainty by summing in quadrature the Peischl estimate and the relative difference between two reported consecutive hourly CO 2 emission values closest to the time of the power plant sampling.The aircraft frequently encountered power plants during oil and gas monitoring campaigns, but usually did not have the flight time to perform a full emissions characterization of the power plant.Here we limit our comparison to days when the aircraft performed a minimum of 10 laps around the plant, thus excluding the quick fly-bys where uncertainties would be unacceptably large.The results are presented in Table 4 and indicate very good agreement between Gauss's method and the reported CO 2 emissions with the averaged difference being 10.6 %.A comparison plot of the reported vs. measured CO 2 emissions is shown in Fig. 10.The average difference between the reported and measured emissions for the five power plants is 11 %.

Conclusion
This technique was developed out of the necessity to identify and quantify individual well pads in an extensive oil and gas production field.Consequently the frequent tracking of the upwind and downwind side of the source provides a very ac- curate determination of the location and magnitude of a given emission site.The main uncertainty arises from the effluent below the lowest flight altitude, but this is minimized by targeting a downwind distance determined by LES studies to provide very little change in the plume flux divergence from the lowest loop to the ground.In addition to the controlled release experiments, hundreds of sites have been measured using this technique with varying levels of success.Ideal conditions include flat terrain, ample sunlight to promote vertical mixing, consistent winds, and no nearby competing sources.
Under optimal conditions we have demonstrated that measurement uncertainties are quite low, often better than 10 %.
As the conditions deteriorate from the ideal to situations involving complex terrain, variable winds or nearby upwind sources, measured uncertainties can increase to be as large or larger than the emission estimates themselves.In the worst case of stably stratified conditions (winter or nighttime), for instance, the lack of vertical mixing may preclude the trace gases emitted at the surface from reaching the minimum safe flight altitude.Complex terrain provides a challenge to the method because the aircraft is unable to maintain a constant altitude above the ground.A possible future refinement of this technique to be applied in complex terrain would be to fit the measurements of both wind and mixing ratio to a uniform 3-D surface surrounding the source, where the grid passes through the terrain and then integrate the flux normal to this irregular virtual flight path.This would not assume level loop flight legs and would, in principle, account for individual loops being flown at differing altitudes and thus more closely track mass continuity near the terrain elevation.

Figure 1 .
Figure 1.Map of the airplane flight pattern sampling a methane plume emanating from an underground storage facility.Wind direction is indicated by the white arrow and the methane mixing ratio is given by the color bar to the right.This flight was conducted on 28 June 2016 and took place between 12:46 and 13:52 LT at altitudes ranging from 91 to 560 m with a loop diameter of approximately 3 km.The measured methane emission rate was 763 ± 127 kg h −1 .

Figure 2 .
Figure 2. Graphical representation of the relative magnitude (%) of the contribution of the horizontal wind divergence to the horizontal advective terms in Eq. (3), as a function of wind speed and source magnitude for methane, using a typical global background of 1.9 ppm and divergence of 10 −5 s −1 .

Figure 3 .
Figure 3. Relative cross wind integrated concentrations of an effluent plume released at the surface in the UCD50B simulation.The data are averaged over 15 min of simulation time and normalized by the maximum concentration.

Figure 4 .
Figure 4. Average cospectrum of the outward-directed component of the observed wind and the methane concentration from 70 laps around a point source near San Antonio, Texas.The peak at 10 −2 Hz corresponds to the period of the circle.

Figure 5 .
Figure 5. Dimensionless flux divergence profiles generated from averaging over 3 different WRF-LES runs using 30 time steps for each one.The horizontal flux per unit altitude (d = F / z) is normalized by the boundary layer height, z i , and source strength, Q.The colored profiles are averages at various dimensionless distances, R = 0.1, 0.2, 0.3, and 0.4, and the gray areas represent 1 standard deviation about the mean.The horizontal dashed lines are the approximate lowest safe flight altitude.

Figure 6 .
Figure 6.Rate of convergence toward the final leak rate estimation as a function of the number of loops for LES CASE UCD508.By 15 laps, the emissions estimate (blue line) has stabilized to 2.5 kg h −1 compared to the actual leak rate (red line) of 2.9 kg h −1 .Dimensionless distance R = 0.25, 50 realizations.Gray area represents 1 SD (standard deviation).

Figure 7 .
Figure7.Averaged LES estimates for the Aerodyne case.This leak shows a slightly higher number of laps before convergence (∼ 25 laps).This simulation was performed using the conditions for the Aerodyne controlled release near Denver, Colorado, on 19 November 2014.

Figure 8 .
Figure 8. Ethane horizontal transport profiles for the Aerodyne controlled releases near Denver, Colorado, on 19 November 2014 ((a)) and in Bee Branch, Arkansas, on 3 October 2015 ((b)).Blue dots represent individual flight loop measurements and the red circles represent the bin average values for altitude intervals represented by the red bars.

Figure 9 .
Figure 9. (a) Time series of methane (blue) and ethane (red) along with (b) the geographic distribution of ethane (color bar) and instantaneous winds (arrows) from a single flight loop during the second controlled ethane release.

Figure 10 .
Figure 10.Comparison of aircraft vs. reported power plant emissions.

Table 1 .
Domain and micrometeorological parameters for the three WRF-LES experiments in this study.L represents the Monin-Obukhov length.

Table 3 .
Controlled natural gas release.

Table 4 .
Power plant estimates.The mid-point of the measurements (Hour UTC) is indicated in the third column (Hour).The reported emissions from the hour before to the hour after that time were averaged to derive the "Reported" emissions in column 5. Emissions are reported in units of metric tons (t) per hour.