Testing of a Lagrangian model of dispersion in the surface layer with cattle methane emissions

Agricultural and Forest Meteorology 150 (2010) 1428–1442

Contents lists available at ScienceDirect

Agricultural and Forest Meteorology journal homepage: www.elsevier.com/locate/agrformet

Testing of a Lagrangian model of dispersion in the surface layer with cattle methane emissions Johannes Laubach ∗ Manaaki Whenua – Landcare Research, P.O. Box 40, Lincoln 7640, New Zealand

a r t i c l e

i n f o

Article history: Received 12 April 2010 Received in revised form 15 July 2010 Accepted 19 July 2010 Keywords: Backward-Lagrangian stochastic model Gas dispersion Cattle CH4 emissions Atmospheric surface layer Horizontal flux profiles Turbulent Schmidt number

a b s t r a c t Results from an experiment measuring methane emissions from a herd of cattle are used to investigate the performance of a backward-Lagrangian stochastic model (distributed under the name WindTrax). The availability of simultaneous mass-budget measurements of the emission rate, together with a unique setup geometry, allow to compare modelled and measured normalised concentration profiles and horizontal flux profiles with five sensor heights, z, and for four horizontal source–sensor distances, x. Simulated emission rates differ typically by 10–20% to those obtained from the mass-budget measurements, which is in agreement with previous tests of the accuracy of WindTrax. Thus, the idealisation of a herd of animals as a homogeneous area source at ground level does not seriously affect the model’s applicability to infer emission rates. The profile comparison suggests that WindTrax may overestimate the speed of vertical dispersion. As a consequence, for this experiment an ideal z/x ratio exists where the modelled emission rate is unbiased. Its value is about 0.080 in unstable and 0.067 in stable stratification. Using concentration measurements taken above or below this z/x threshold leads to emission rates that are slightly under- or overestimated, respectively. Simultaneous measurements with an open-path methane laser are compatible with this finding. Possible causes of the apparent overestimate of vertical dispersion rates are discussed, leading to the cautious suggestion that it may stem from the choices for the Kolmogorov constant and/or the normalised dissipation rate in the model, which reflects gaps in our understanding of the atmospheric surface layer. It is argued that this notion does not contradict the earlier results from a number of controlled tracer-release experiments that had been designed to test WindTrax. © 2010 Elsevier B.V. All rights reserved.

1. Introduction “Lagrangian stochastic” (LS) models are the “natural and most powerful means” (Wilson and Sawford, 1996) to describe atmospheric dispersion processes. Such models can be run backwards in time in order to efficiently infer emission rates of gases from confined sources, using measurements of concentrations downwind of the sources as inputs. One particular implementation of such a backward-Lagrangian stochastic model (BLS) for the atmospheric surface layer is available under the name WindTrax (Thunder Beach Scientific, Nanaimo, Canada). WindTrax has recently become popular to estimate gas emissions from farm operations (see references later in this section). There is, therefore, a need to assess the performance of this model not only in ideal test setups, such as controlled-release experiments, but also in real-world farm situations. This study attempts such an assessment.

∗ Tel.: +64 3 321 9999; fax: +64 3 321 9998. E-mail address: [email protected]. 0168-1923/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.agrformet.2010.07.006

The mechanics of WindTrax have been described in detail by Flesch et al. (1995, 2004). At the core is Thomson’s (1987) wellmixed three-dimensional LS model for Gaussian inhomogeneous turbulence. Wilson and Sawford (1996) point out that for the vertical dispersion in this type of turbulence, the Thomson model is the unique LS model. It contains one empirical parameter, which can be expressed either as a Lagrangian timescale or the Kolmogorov constant, as discussed by Wilson et al. (2001, 2009). In the horizontal dimensions, the Thomson model is not unique, and it is well-known that any model that uses only surface-layer parameters is prone to error because streamwise and lateral wind fluctuations contain larger-scale components originating outside the surface layer. For this reason, Flesch et al. (2004) considered line-averaged concentration measurements preferable to one-point measurements: provided that the line average includes the entire width of the gas plume, the BLS results will be insensitive to modelling error of horizontal dispersion. Flesch et al. (1995, 2004) showed that WindTrax derives emission rates from measured concentrations with about 20% accuracy for a wide range of conditions in the atmospheric surface layer. This makes it arguably the best available model for this purpose

J. Laubach / Agricultural and Forest Meteorology 150 (2010) 1428–1442

in obstacle-free flow over flat terrain. Subsequently, WindTrax has been used to determine the trace gas emissions from livestock: of methane (CH4 ) by Laubach and Kelliher (2005a,b), McGinn et al. (2006), Laubach et al. (2008), of ammonia (NH3 ) by Denmead et al. (2004), Sommer et al. (2005), Flesch et al. (2007), McGinn et al. (2007), Harper et al. (2009), and of both gases simultaneously by Loh et al. (2008) and van Haarlem et al. (2008). Laubach et al. (2008), henceforth L08, compared the BLS results to those of two other micrometeorological methods, using the same measured concentrations, wind speeds and turbulence parameters, and to an animal-scale tracer ratio method which determined the daily CH4 emissions from each animal. The two other micrometeorological methods were an integrated horizontal flux (IHF) method (also known as mass-budget) and a flux-gradient (FG) method (with footprint correction accounting for the limited source area extent), both described in detail by Laubach and Kelliher (2004). Considering the animal-scale method as an absolute reference, L08 concluded that FG and IHF “provide unbiased means if the minimum distance between the animals and the instruments is of order 20 m”, while BLS was “the most consistent method over time, but a bias of order 20% needs to be corrected for”. This latter result is in qualitative agreement with two other experiments on CH4 emissions from grazing cattle (Laubach and Kelliher, 2005a,b). However, it is at odds with the results of tracer-release experiments which showed no significant bias of the mean gas recovery rate (Flesch et al., 2004; Gao et al., 2009a). A somewhat more complicated picture emerges from the results of two other release experiments. First, McBain and Desjardins (2005) observed a dependence of gas recovery rates on measurement height, with emissions overestimated when using concentrations at the lowest height, and underestimated when using the uppermost height, respectively. In contrast to this, Laubach and Kelliher (2005a) had previously found the calculated emission rates to increase with the height of the concentration measurement, in three experiments with cattle herds. Second, Gao et al. (2009b) reported a stability dependence of the gas recovery rate, underestimating the true emission rate in unstable and overestimating it in stable stratification. The present paper aims to understand whether such discrepancies can be attributed to the performance of the BLS model itself, i.e. the Thomson (1987) model, or whether they are due to operational differences between the tracer-release experiments and the animal-herd experiments. Possible factors are: differences in source geometries, differences in sensor geometries, and differences in the implementation of WindTrax. These are briefly described in the following. First, there are some principal differences in the nature and distribution of the sources between the two groups of experiments. In the tracer-release experiments a true area source at ground level is created as well as possible, by a pipe grid releasing the trace gas uniformly in the horizontal dimensions, and constant in time. In the animal-herd experiments, the assumption of an area source is an approximation for a large number of point sources, the cattle, which can move freely within the confines of the paddock. Observations of the cattle behaviour show that they are generally evenly spread, but nevertheless the source distribution is patchier and less homogeneous in space and time than in the tracer-release experiment. In addition to the variability of emission locations due to animal movement, there is generic variability of emission rates between animals, and there are temporal variations related to the digestion processes. Another idealisation of the animal-herd experiments is that the emissions are assumed to be at ground level. While this is true for a major fraction of the emissions (whenever an animal exhales while its mouth is near the ground, either because the animal grazes or because it is lying down to sleep or ruminate), it is not true for all of the emissions, since the mouth of a standing cattle would typically be about 1 m above ground. A third difference

1429

between the two groups of experiments is the size of the source area, typically a few meter in each horizontal dimension for the tracer releases, compared with tens to hundreds of meters for the animal-herd experiments. Second, the quoted tracer-release experiments and the animal-herd experiments differ principally in the distribution of concentration measurements. The latter employed vertical profiles of point measurements, realised by intake tubes leading, via ballast volumes, to the same gas analyser, which received air from these intakes sequentially in a switching cycle. By contrast, the tracer-release experiments all employed arrays of line concentration sensors, and, except in the study of McBain and Desjardins (2005), these sensors were all deployed at the same height. (Though Gao et al. (2009a) made measurements at six heights, they reported recovery rates obtained with the BLS method only for one.) Third, WindTrax offers the user some choices how the wind and turbulence parameters are implemented. While the relationships between wind profile and wind statistics in the model are always forced to obey Monin–Obukhov similarity theory (MOST), there are different ways how these parameters can be entered. In McGinn et al. (2006, 2007), Flesch et al. (2007), and Harper et al. (2009), all required turbulence parameters are taken from a single sonic anemometer and WindTrax uses these to derive a vertical wind profile such that the roughness length, z0 , is fitted in each run for consistency between measured mean wind speed u, friction velocity u* and Obukhov length L. As a consequence of inevitable measurement error as well as inhomogeneity and instationarity effects, z0 can vary considerably and unrealistically from run to run. Where implausible values of z0 occur, they are an indication that the shape of the wind profile near the ground may be inaccurate, which is a concern because all the emitted gas passes through the near-ground air layer on its way towards the measurement location, and consequently the speed of its horizontal advection may be in error. Alternatively, Laubach and Kelliher (2005a,b) and L08 first determine realistic z0 values from the whole dataset of a fivepoint cup anemometer profile, corroborated by u* and L from a sonic anemometer and visual assessment of site characteristics. These z0 values do not vary at short time scale, only in response to observed variations of vegetation height (and where necessary, they vary with wind direction). This z0 together with the wind speed from the uppermost cup anemometer and L from the sonic is entered into WindTrax. The model then computes u* and ignores the value measured by the sonic anemometer, which is typically subject to 10% random error (and sometimes biased due to imperfect mounting). The different ways of entering the flow field parameters will lead to some differences in the calculated emission rates. The previously quoted studies all assess the quality of the BLS results by comparing the modelled emission rate to either a known release rate or an alternative estimate of the unknown emission rate obtained by a different method. Here, a second aspect is added to the model assessment. We will investigate the field of the ratio of concentration, C (minus background, Cb ) to emission rate, Qc , in the alongwind-vertical plane (x,z) that contains the measurement mast (y = 0): B(x, z) =

C(x, z) − Cb Qc

(1)

The normalised concentration field B(x,z) is unambiguously constrained in WindTrax for a given source geometry and a given set of turbulence parameters. The specific experimental geometry of L08 provides the – so far unique – opportunity to compare the simulated B for five values of z and four values of x with reference B values obtained from measured concentrations and independently derived emission rates. Hence, the data of L08 will be re-analysed in the following to suit this purpose.

1430

J. Laubach / Agricultural and Forest Meteorology 150 (2010) 1428–1442

Fig. 1. Schematic of the experimental setup, reproduced from Laubach et al. (2008). Four grass strips were grazed on consecutive days. Each strip was subdivided into N (uppermost) and S half, containing 29 steers each. Wind direction was predominantly SW to NW. Each strip is labelled by rmin , the minimum distance between its fenceline and the profile mast (measured along the dashed line).

2. Properties of the BLS model

as (Flesch et al., 1995):

WindTrax is described comprehensively by Flesch et al. (2004). It determines B(x,z) as follows. First, the flow is assumed stationary and horizontally homogeneous, over a flat surface with uniform z0 . Its characteristics are prescribed by a set of measurable parameters, detailed below. Using these, the model simulates the flight paths of air parcels from the detection point backwards. The velocity of each air parcel is varied for each time step, to simulate how the air parcel responds to the forces exerted by other air parcels that it encounters as it changes position in the flow. This process is described by the Langevin equation:

B(x, z) =

uk = (˛k + ˇ Pk )t

(2)

where t is the discrete time step, uk (u1 , u2 , u3 = u, v, w) the associated longitudinal, lateral, and vertical velocity changes, Pk a Gaussian random number, from a population with a mean of zero and variance of one, and ˛k and ˇ are model coefficients. For each time step, the position of the air parcel, xk (x1 , x2 , x3 = x, y, z), changes by the amount: xk = uk t

(3)

and the model process is iterated by re-applying (2). Over many time steps, a complicated trajectory develops. The simulation is repeated to represent a large number of air parcels, n, each yielding a different trajectory because of the random term in (2). When an air parcel touches the surface, it is reflected upwards, i.e. its momentary vertical speed at the touchdown, w0 , changes sign. If a touchdown occurs within the source area, the air parcel contributes to the surface flux of the emitted gas, at a rate proportional to w0 . Touchdowns outside the source area do not contribute to the flux. Hence the normalised concentration B at a given location (x,z) can be constructed as the average of flux-contributing and noncontributing air parcels passing through this location. B then results



1 n

2 . |w0 |

(4)

within emitting area The measurable parameters required to describe the flow are: mean wind speed u at one point, friction velocity u* , an estimate of stability, e.g. L, the standard deviations of the three wind components,  u ,  v ,  w , and the dissipation rate of turbulence kinetic energy (TKE), ε (Flesch et al., 2004). The first three together determine the shape of the mean wind profile unambiguously. In practice, the model allows replacement of one of these three by z0 . The approach chosen here is to prescribe u, L, and z0 , and have u* computed by WindTrax. The other four parameters, together with u* , describe the turbulent properties of the flow. The standard deviations of the wind are entered to the model in normalised form as  u /u* ,  v /u* ,  w /u* , or, if these are not available from measurements, then WindTrax uses a default parameterisation with similarity relationships (Flesch et al., 2004, Appendix). The coefficients ˛k in (2) are functions of  u ,  v ,  w , u* , and also ˇ, which is a function of ε (see Eq. (5)). The TKE dissipation rate cannot be entered by the user of the model. It is parameterised as ε = (u3∗ /z)ε , where ␧ is a dimensionless function of z/L and  the von-Kármán constant. It determines the random coefficient ˇ, as: ˇ = (C0 ε)1/2

(5)

where the Kolomogorov constant, C0 , needs to be set on an empirical basis. Wilson et al. (2001) show that the choice of C0 is equivalent to prescribing a Lagrangian timescale, and also sets the turbulent Schmidt number, Sc (the ratio of momentum diffusivity to gas diffusivity). In WindTrax, the default value is C0 = 4.41, implying Sc = 0.64 at neutral stability (Flesch et al., 2004), thus the turbulent vertical “diffusion” of gas is assumed 1.56 times more rapid than that of momentum. This choice provides consistency with the data analysis of Wilson et al. (1981) for the Project Prairie Grass Experiment (Wilson et al., 2001).

J. Laubach / Agricultural and Forest Meteorology 150 (2010) 1428–1442

3. Experimental The data were collected over a 3-week period in austral spring 2006 on a research farm 20 km inland from the West coast of the North Island of New Zealand (40◦ 20 S, 175◦ 28 E, 13 m a.s.l.). The dominant wind direction (in general and in this experiment) is from W, including frequent afternoon sea-breezes when synoptic wind patterns are weak. The surrounding terrain consisted of flat paddocks, free of flow obstacles to W, S and E for at least 500 m. To the N, there was a shelterbelt at ca. 200 m distance. 3.1. Grazing area and setup geometry Fifty-eight steers of average weight 325 (±20) kg grazed highquality ryegrass pasture ad libidum. Each day, the cattle moved to a fresh rectangular strip, measuring 35 m (W–E) by 100 m (N–S). The strip was subdivided in the centre of the longer axis, with each half containing 29 steers, to partially prevent clustering of the animals (Fig. 1). In order to estimate the CH4 emissions from the cattle herd, the grazing strip is considered as an area source with spatially uniform emission rate (varying in time between measurement runs, each of 20 min duration). The CH4 concentrations as effected by the steers’ emissions were measured with two separate instrumentations: an open-path laser system, which will only briefly be considered below, and a vertical profile system, which provides the main dataset. Every fifth day, all instruments were moved to a new area providing a fresh set of four grazing strips as sketched in Fig. 1. Four cycles associated with fixed instrument locations were completed. Suitable wind directions prevailed on 13 days, with some nighttime gaps. 3.2. Vertical profile system and turbulence measurements The profile included five measurement heights (zj = 0.7, 1.3, 2.3, 4.2 and 7.3 m for j = 1–5, respectively). Air from these heights was fed via a datalogger-controlled valve-switching system to a flame ionisation detector (inside a gas chromatograph, HP 6890, Wilmington, Delaware), which measured CH4 concentration. In addition, air was sampled from two locations, NW and S of the herd, via 300 m long intake tubes. Either of these provided the local background concentration (upwind of the herd), depending on wind direction. Two pumps were employed: one to draw air from all seven intake lines continuously, in parallel, into separate ballast volumes (one per intake line), the other to then draw air from the selected intake to the CH4 sensor. That way, even though the intake lines were sampled sequentially, each sample represents a time period about five times longer than the actual sampling time, due to the premixing in the ballast volumes. Of every 20-min run, 2 min were reserved for a two-point calibration check of the CH4 sensor with known concentrations from gas cylinders. These checks were used to remove signal drifts in postprocessing. Wind speed was measured by cup anemometers (Model A101 M, Vector Instruments, Rhyl, Clwyd, UK) with matched calibrations. They were installed at the same heights zj as the intakes for the CH4 concentration samples. Wind direction, velocity statistics and the Obukhov length were provided by a sonic anemometer (RM Young 81000 V, Traverse City, Michigan) on a separate mast, east of the grazing strips and 3.86 m above ground. Leaving the instruments in the same positions for four consecutive grazing days had the consequence that every day the horizontal distance between the grazing cattle and the profile mast changed, while the shape of the grazing area stayed virtually identical, as shown in Fig. 1. This makes the layout physically equivalent to one in which the profile mast would have been relocated horizontally while the area source remained fixed; or, if one selected for runs of similar wind vector and turbulence parameters from

1431

different days, to a layout in which profile masts in different horizontal locations would have operated simultaneously. It is this translational symmetry that will be exploited here to investigate the two-dimensional field B(x,z). For the purpose of run classification, the horizontal distance is labelled by the minimum distance between the cattle and the mast, rmin , as in Fig. 1. In the easternmost strip, rmin = 5 m, then successively to the W, rmin increased to 22, 57 and 92 m, respectively. The actual x for each run is not a single value but an interval xmin < x < xmax where xmin is the distance from the measurement mast to the intersection of the alongwind axis (in upwind direction) with the relevant near fenceline, and xmax the distance from the mast to the intersection with the relevant far fenceline. For a given run, obviously xmin ≥ rmin , dependent on wind direction. The ranges of x for neighbouring grazing strips can partially overlap, but wind direction statistics and run selection procedures ensure that average x systematically increases with increasing rmin , roughly in proportion, which suffices for the present analysis. 3.3. Path-averaged concentration measurements The open-path laser system (GasFinder MC, Boreal Laser Inc., Spruce Grove, Alberta, Canada) measured line-integrated CH4 concentrations along each of four paths located outside the grazed area’s fencelines to the N, W, S and E, respectively (Fig. 1). Here only data from the E path, mounted 1.1 m above ground, will be considered. This path was 97 m long during the first four and 66 m during the last three measurement days (in between, the laser system was offline for repair). The minimum distances between this path and the animals were xmin = 6.5, 41, 77 and 112 m for the four grazing strips labelled in Fig. 1 with rmin = 5, 22, 57 and 92 m, respectively. 4. Data collection and processing 4.1. Emission rates and C/Qc ratios for BLS model and reference L08 compared altogether five methods to obtain the CH4 emission rates from this experiment. Four of these were micrometeorological methods, and three of these based on identical data from the vertical profile mast: the integrated horizontal flux (IHF) method (referred to as “mass-budget (MB)” in L08), the fluxgradient (FG) method, and the BLS method with concentrations from all five heights as simultaneous inputs, which leads to an overdetermined system of equations that the WindTrax software solves with a linear optimisation algorithm. The other two methods were BLS using concentration data from a CH4 laser system (with much smaller data yield due to instrument failure), and an animal-scale tracer ratio method which delivered daily averages of CH4 emissions from each animal. Here, we are interested only in the BLS method, and specifically in the question if it models the field B(x,z) correctly. To investigate this, the execution of WindTrax was repeated for each zj with a single concentration input from that height, Cj , predicting an estimate of the emission rate that is denoted Qcj . Combining model input and output using (1) therefore allows to retrieve: B(x, zj ) =

Cj − Cb Qcj

(6)

for each run. If the Cj were error-free and the parameterisation of B(x,z) exactly correct, then all five Qcj would be equal (and accurately represent true Qc ). In practice, neither condition is fulfilled, so the Qcj differ from each other, with the differences being caused by the combined measurement errors of concentrations and flow parameters plus the model’s error (which includes inaccuracies and approximations in the source geometry). Simultaneously, the

1432

J. Laubach / Agricultural and Forest Meteorology 150 (2010) 1428–1442

model can provide predictions of the concentrations1 at the other heights, Cjk , with k = / j, which analogously allows one to obtain: B(x, zk ) =

Cjk − Cb Qcj

(7)

The five repetitions of WindTrax thus provide five predictions of the concentration profile, with the value for k = j equal to the actually measured value. Any two of these profiles differ from each other linearly with a gain factor equal to the ratio of the respective Qcj and an offset representing the error in (Cj − Cb ). Because of the compounded errors, neither the Qcj nor the Cjk are easily used to assess the accuracy of the BLS model. Two combinations of these variables suggest themselves for this purpose. The first is obtained from dividing (7) by (6), which gives the ratio of a simulated profile (for any j) to the measured profile (Cjk − Cb )/(Cj − Cb ). Ideally this ratio should be identical to 1 for all k, in which case the shape of the measured concentration profile would be fully consistent with the given flow field. In practice, deviations from 1 can be interpreted either as a measure of disagreement between measured and simulated flow field, or as a measure of inconsistency between the Cj from different heights. Since Qcj drops out by division, no assumption is required on what the true emission rate is. However, because it is unknown which of the five concentration measurements comes closest to the truth, the profile ratio cannot measure the absolute error of the simulation, for any height. The other useful variables for model error assessment are the B(x,zj ). These are identical for the different repetitions of WindTrax because their parameterisation depends only on the set of flow parameters, which are identical for all repetitions. Together, for fixed x, a set of B(x,zj ) represents the flux-normalised concentration profile at that horizontal position. The error of this normalised profile depends only on the measurement errors of the flow parameters and on the model’s error, but not on concentration errors. In order to assess the accuracy of the B(x,zj ) profile, it is necessary to compare it to a reference which is assumed to represent the true profile as closely as possible. This requires a reference concentration profile and a reference emission rate Qc for each run. The former is the measured profile; for the latter, there are a number of possible choices, all with their inherent uncertainties (since the true emission rate is not known on a run-to-run basis). Here, the emission rates obtained by L08 with the IHF method are chosen as the reference, for the following reasons. First, the analysis of error propagation and sensitivity tests have shown that the IHF method has smaller measurement error than the FG method (Laubach and Kelliher, 2004) and the BLS method (Laubach and Kelliher, 2005a). Second, the IHF method uses the background concentration in exactly the same way as the BLS method, as a constant additive term, hence the error of Cb affects both methods equally. Third, the IHF method makes direct use of the complete concentration profile, thereby providing internal consistency between the (Cj − Cb ) and the Qc for each run. Fourth, the major terms contributing to the IHF estimate of Qc do not depend on any flow parameters

other than the mean wind profile – only the correction term for turbulent backflow does, which is of order 10% of Qc (see Laubach and Kelliher, 2004 for details). For each individual run the field B(x,z) describes the essence of the model’s output, in principle sufficient to assess model performance. However, when pooling a number of runs into a class, as done below, a number of parameters cause variations in (C − Cb ) even when Qc remains constant. A class-averaged B(x,z) will therefore have an additional error component that is due to uneven weighting of these parameters between runs. The first such parameter is wind speed: the larger u, the smaller (C − Cb ) in proportion, so runs of low wind speed are overrepresented in the average B. The second parameter is the alongwind contact distance of the mean flow with the gas-emitting area, x = xmax − xmin (Laubach and Kelliher, 2004). Due to the setup geometry, x varies with wind direction, also causing variations of B between runs. The weighting effect of u and x can be avoided if horizontal flux profiles H(z) are considered instead of B profiles, as shown in the next section. A third parameter influencing (C − Cb ) is stability, which will be explicitly accounted for in the class selection. Run-to-run variability in the  u /u* ,  v /u* ,  w /u* input values has a minor effect on average B and will be neglected. 4.2. Normalised horizontal fluxes The horizontal flux, Fj , at height zj , without turbulent backflow correction, is given by: Fj (x, zj ) = uj (Cj − Cb )

  Fj zj x 6

QcIHF =

It was found that the concentration outputs of the WindTrax version used by L08 (2.0.6.6) and the version which is now available as Freeware (2.0.8.1) did not agree. Further tests showed that changes of Cb by a constant offset changed the Cjk in Version 2.0.8.1 by exactly the same offset, as it should be, but in Version 2.0.6.6 by slightly different values. The observed differences were small relative to Cjk , but often significant relative to (Cjk − Cb ). It is concluded that Version 2.0.6.6 contains a processing error that occurs for situations with non-zero background concentration, and that this error has been corrected in more recent versions of WindTrax. In the present dataset, the effect of this error on calculated emission rates is negligible (<10−3 Qcj ), hence the results of L08 are still valid. Interestingly, the earlier Version 1.6.0 of WindTrax, used by Laubach and Kelliher (2005a,b), exactly reproduced the results of Version 2.0.8.1 and thus appears free of this processing error.

(9)

j=1

where  is the turbulent backflow correction factor. Normalising (9) with QcIHF and using (6) yields: x = 

6 

uj Bj zj

(10)

j=1

x depends on the setup geometry (Fig. 1) and wind direction, with xmin ≥ rmin , and as the wind direction fluctuates within a run, so does x. Details of how to compute x under the assumption of a Gaussian distribution of the crosswind fluctuations are given by Laubach and Kelliher (2004). Here, we wish to remove the geometrical dependence in order to obtain classes of similar runs that can be averaged to assess model performance. To do so, we arbitrarily define a standard contact distance as the West-East extension of the paddock rectangles (Fig. 1), i.e. x0 = 35 m, and introduce: b=

x x0

(11)

which can be considered a “fetch normalisation factor”. Using that: Hj = uj Bj b−1

1

(8)

where uj is the mean wind speed at zj , and the x dependence of Cj has been dropped for the sake of shorter notation. In this particular realisation of the IHF method:

(12)

is the normalised horizontal flux at zj , which can be compared between runs as well as between the IHF and BLS method. Vertical integration of Hj retrieves x0 /. 4.3. Meteorological parameters and run selection L08 obtained 383 valid 20-min runs with the following two selection criteria. First, wind speed had to be fast enough (by demanding that u* > 0.12 m/s) and steady enough (mean wind speed, at sonic anemometer height, exceeding 5 u* ). This is necessary because neither the IHF nor the BLS method are expected to

J. Laubach / Agricultural and Forest Meteorology 150 (2010) 1428–1442

Fig. 2. Frequency distribution of the inverse Obukhov length for the 307 selected runs, in classes of width 0.005 m−1 .

work in very weak or variable winds. Second, in order to capture a sufficient quantity of CH4 gas from the animals, data were accepted if the grazed strip accounted for more than 10% of the mast’s footprint, i.e. when the source area weight function, f, exceeded 0.1. The calculation of f combined the model of Schuepp et al. (1990) in the alongwind direction, including a stability correction, with a Gaussian plume description of the lateral dispersion and the actual fenceline geometry (Laubach and Kelliher, 2004). For this study, as a further quality requirement it is demanded that the emission rates obtained with the IHF method and with BLS using the concentration at the lowest height as input do not differ by more than 40% of the larger of the two, i.e. 3/5 < Qc1 /QcIHF < 5/3. There are a number of reasons that may cause Qc ratios outside these bounds: very small emission rates causing large error in the measured (Cj − Cb ), disturbed concentration profiles producing biased Qc , large error in one or more of the turbulent parameters affecting Qc1 but not QcIHF , and possibly others. By excluding such dubious runs the dataset is restricted to normal operating conditions for which the BLS model has been designed. 307 runs passed this additional quality test. These runs were sorted into four classes according to the minimum source–sensor distance, rmin , as shown in Fig. 1. In order to average only runs with similar profile shapes, these were further subdivided according to stability and roughness, as follows. A histogram of the inverse stability parameter, 1/L, shows that the stability distribution is heavily dominated by near-neutral conditions (Fig. 2). The total range spans only −0.103 < 1/L < 0.069, so (after the necessary exclusion of runs of low wind speed) no instances of extreme (in-)stability remain, and only a small number of moderately stable or unstable runs. The dataset is therefore

1433

divided into only two stability classes, “stable” and “unstable” according to the sign of L. It makes sense to assess these two classes separately because WindTrax, like any MOST-based model, uses different empirical parameterisations for stable and unstable conditions. For rmin = 92 m, no stable runs were found to satisfy the quality criteria, so in effect three stable classes and four unstable classes were obtained. The roughness length, z0 , was determined as a single value for each grazing strip (grazing day), as follows. First, the measurement heights were reduced by the zero-plane displacement d, assumed as d = 0.13 m (constant throughout the experiment), which was 2/3 of the average grass height. Due to the uniform flatness of the site, no directional dependence of d or z0 was apparent. Run-by-run estimates of z0 from the cup anemometer profiles were selected for near-neutral stability and averaged on a daily basis. The resulting values agreed well with an alternative estimate calculated from u and u* of the sonic anemometer. They ranged from 0.016 to 0.025 m. During the course of the experiment, z0 tended to increase as the grass on the ungrazed parts of the paddock grew. z0 also tended to be larger (≥0.022 m) for rmin = 5 m than for the other classes, indicating that the presence of the cattle near the measurement mast increased the roughness somewhat. Data from different days in one class were averaged together if z0 differed by no more than 0.002 m between days. Otherwise, the class was split into two. This split occurred for rmin = 22 and 57 m, both stable and unstable. Consequently, a total of 11 classes resulted (Table 1). The average z0 of these classes all assumed one of the following values: 0.016, 0.017, 0.023 or 0.025 m. The first two of these will be referred to as “smooth”, the last two as “rough”, however, it should be kept in mind that the difference is only a gradual one. 5. Results The first point of the analysis is to compare the shape of the concentration profile simulated by the BLS model to that of the observed concentration profile. This is done by displaying ratios of simulated and observed concentrations (Fig. 3), for each class. For deeper analysis of the observations, three types of profiles are then shown in parallel: profiles of wind speed, normalised concentration B, and normalised horizontal flux (Fig. 4 for unstable and Fig. 5 for stable stratification). Along with the discussion of these figures, the significance of each type of profile is explained in a brief section. Implications for the calculated emission rates and results from the path-averaging concentration sensor are discussed in subsequent sections. 5.1. Concentration ratio profiles Simulated and observed concentration profiles are compared as follows. Instead of prescribing an emission rate, the concentration at the lowest height, C1 , is given and used to predict the emission

Table 1 Classes of 20-min runs distinguished by minimum mast-cattle distance, rmin , sign of the stability parameter, 1/L, and roughness length, z0 . Columns give class-means of these parameters, number of runs, N, and figure panels where results are shown. Description

rmin (m)

1/L

z0 (m)

N

Profiles shown in Figs.

Nearest – unstable – rough Nearest – stable – rough 2nd-nearest – unstable – smooth 2nd-nearest – stable – smooth 2nd-nearest – unstable – rough 2nd-nearest – stable – rough 2nd-farthest – unstable – rough 2nd-farthest – stable – rough 2nd-farthest – unstable – smooth 2nd-farthest – stable – smooth Farthest – unstable – smooth

5 5 22 22 22 22 57 57 57 57 92

−0.0187 0.0105 −0.0086 0.0123 −0.0108 0.0192 −0.0019 0.0052 −0.0099 0.0109 −0.0061

0.023 0.023 0.017 0.017 0.025 0.025 0.025 0.025 0.016 0.016 0.016

46 61 64 34 22 12 12 6 21 9 20

3, 4 (top panels) 3, 5 (top panels) 3, 4 (top panels) 3, 5 (top panels) 3, 4 (centre panels) 3, 5 (centre panels) 3, 4 (centre panels) 3, 5 (centre panels) 3, 4 (bottom panels) 3, 5 (bottom panels) 3, 4 (bottom panels)

1434

J. Laubach / Agricultural and Forest Meteorology 150 (2010) 1428–1442

The concentration ratios in Fig. 3 tend to be systematically >1 (except for the lowest height where the ratio is 1 by definition). For the smallest source–sensor distance (top panel), the largest Csim /Cobs ratio occurs at z ≈ 4 m, for the medium distances the largest ratios occur towards the top of the observed height range, and for the largest distances, the ratio drops back towards one. Fig. 3 thus shows that simulated concentrations deviate in a small, yet measureable and systematic, way from the observed concentrations, and they do so in the absence of any assumption on what the true emission rates were. From this it follows that the observed deviations do not originate from any computational details in determining B or H as presented in the following sections. 5.2. Wind profiles For each run, WindTrax constructs a complete, horizontally homogeneous wind profile from the given flow parameters, using MOST. Here, the given parameters are u5 (mean wind speed at z5 − d = 7.14 m), L and z0 . On the left-hand side of each panel of Figs. 4 and 5, the modelled wind profiles are compared to those measured with cup anemometers, again for each of the 11 classes defined earlier. The modelled profiles are anchored to the measured ones at z5 − d. Below that, any differences between measured and modelled profiles can be due to a combination of cup anemometer error (of order 1%), mismatch of modelled and true u* (due to error of d, z0 or L), and profile disturbances caused by obstacles (the latter considered negligible here). The figures show very good agreement between measured and simulated wind speeds. In Fig. 4 (unstable stratification) the modelled u at z2 , z3 and z4 are up to 3% larger than the measured ones, probably due to an overestimate of u* . In Fig. 5 (stable stratification) the measured and modelled wind speeds agree within 1% (indistinguishable in the graphs). It is concluded that the horizontal (advective) gas transport has been modelled correctly. It also follows that the simulated u* is consistent with the observed wind shear (apart from the slight overestimate in unstable conditions). 5.3. Normalised concentration profiles

Fig. 3. Concentrations (minus background), simulated for the upper four heights when given at the lowest height, divided by observed concentration. The top two panels show mean profiles of four classes each and the bottom panel of three classes, from Table 1, with rmin increasing from top to bottom. Note the logarithmic scale on the concentration ratio axis.

rate, Qc1 , for each run. Simultaneously, the concentration values at the upper four heights, C1k , are simulated. The result is a concentration profile that is self-consistent for the given flow parameters and the simulated emission rate. For each of the 11 classes defined previously (dependent on rmin , z0 , and stability), an average profile is computed, after subtracting the background concentration from the C1k . The observed concentration profiles are reduced by the same Cb values and class-averaged, too. Simulated and observed profiles are then compared by forming their ratio (Fig. 3). The choice of the bottom height for the input concentration is arbitrary: choosing another height would linearly transform the resulting profile, as explained below Eq. (7), and therefore change the coordinates on the concentration ratio axis, but not the profile shape.

In the centre of each panel of Figs. 4 and 5, the class-averaged normalised concentrations according to Eq. (6) are shown as observed and as simulated by WindTrax. In the top panels of both figures, i.e. for the smaller values of rmin , the simulated B profiles agree well with the observed profiles. In the centre panel of Fig. 4, for the medium range of rmin combined with large z0 values, systematic differences emerge: below ca. 4 m height, the measured B exceeds the simulated value, while at larger heights, the reverse is true. The same is observed for stable stratification (centre panel of Fig. 5), perhaps with arguable significance at rmin = 22 m, but clearly recognisable at rmin = 57 m. In the bottom panel of Fig. 4 (large rmin combined with small z0 ), the measured B exceeds the simulated value for all of the visible profile, i.e. up to at least 7 m. In Fig. 5 (same for stable stratification), measured B exceeds simulated for heights below ca. 5.5 m, and the reverse is true aloft. It therefore appears that the model consistently overestimates the rate of vertical mixing, with the result that as downwind distance increases (represented by the proxy rmin ), according to the model a larger fraction of the released gas reaches the upper levels than in reality. The random sampling error of Bj can be estimated as the standard error of the mean (SE), for each class and height. For z1 to z4 , the relative error (SE/mean) ranges from 2% to 20%, with a median of 6%. For z5 , the relative error tends to be larger, but the absolute error (SE itself) is smaller than at the lower heights. These random errors are not indicated in Figs. 4 and 5 so as not to distract from the finding that the observed and simulated systematically crossover, in the majority of classes.

J. Laubach / Agricultural and Forest Meteorology 150 (2010) 1428–1442

1435

Fig. 4. Profiles, for six classes of unstable stratification, as observed and as simulated by WindTrax for wind speed u (left), normalised concentration B (centre), and normalised horizontal flux u B b−1 with b from Eq. (11). Each panel shows two classes from Table 1. rmin increases from top to bottom. In the top panel, z0 = 0.023 and 0.017 m for the classes with rmin = 5 and 22 m, respectively. Horizontal dashed lines indicate layer depths for the IHF method, Eq. (9).

5.4. Normalised horizontal flux profiles The normalised horizontal flux profiles H = u B b−1 are shown on the right-hand side of each panel in Figs. 4 and 5. Simulated wind speeds were used for both the measured and modelled profiles. Thus, for each run the factor u b−1 is a simple scaling factor, equal for measured and modelled profile, and the only cause of

non-equivalence to the B profiles in the middle of the panels is the different relative weight of each run within a class. Since the height integrals of the H profiles are equal to x0 /, they are the same for all classes except for small variations of the backflow correction factor , see Eqs. (10)–(12). Average  was found, by the method of Laubach and Kelliher (2004), to be 0.90 (± 0.03 standard deviation). Consequently, the H profiles can be interpreted as a direct visuali-

1436

J. Laubach / Agricultural and Forest Meteorology 150 (2010) 1428–1442

Fig. 5. Same as Fig. 4 but for five classes of stable stratification. No stable runs with sufficiently large footprint proportion were available for rmin = 92 m (bottom panel).

sation of how the normalised emission rate is redistributed in the vertical as the gas plume is swept along by the mean wind profile. Comparing them for increasing rmin then shows how this vertical spread changes with downwind distance from the source area. The random relative error of Hj is generally similar to or slightly less than the relative error of Bj , typically 5% (not shown). The normalised horizontal flux profiles reiterate the main findings from the normalised concentration profiles. At the bottom of the profiles the measured values exceed the modelled values (except for the rmin = 5 m classes, both unstable and stable), and vice

versa at the top of the profiles. The crossover height, where measured and modelled H are equal, increases with increasing rmin . It is somewhat lower for the stable classes than for the corresponding unstable classes. Again, a consistent explanation of this behaviour would be that the model somehow overestimates the speed of the vertical dispersion. The maximum of the H profile, which can be interpreted as the centre of the gas plume, is usually located below the crossover height, hence the model tends to underestimate the concentration at the plume centre, by up to 30% for rmin = 92 m (bottom panel of

J. Laubach / Agricultural and Forest Meteorology 150 (2010) 1428–1442

1437

Fig. 4). Whether the model over- or underestimates the height of the plume centre cannot be unambiguously resolved with the current 5-point profiles. The suspected mechanism of somehow overestimated speed of vertical dispersion would suggest that the model overestimated the height, but since the plume shape is not static, the same mechanism may simultaneously lower the position of the maximum by moving gas above the maximum upwards at a faster rate than it is replaced from below. 5.5. Location of the crossover point Figs. 4 and 5 indicate that the field B(x,z) is simulated by WindTrax with a systematic error that depends on position. Let us assume that the effective average source–sensor distance can be approximated as x = rmin + 17 m (half the width of the grazing strip added), giving x = 22, 39, 74 and 109 m for the four values of rmin . (For rmin = 5 m, this assumption is clearly incorrect for W winds, because then xmax = 22 m, but since for these classes the acceptable wind directions range from 180 to 360◦ , they include runs with xmax up to 50 m also.) For each class, the height of the crossover point, zcr , is then read off the right-hand panels of Figs. 4 and 5 (the H profiles), with an error range defined as the height range in which simulated and observed H differ by less than 5%. (The crossover heights of the B profiles are very similar to those of the H profiles, and always within the given error ranges.) Crossover heights are thus obtained for 10 of the 11 classes, since the first class in Fig. 4 (nearest–unstable–rough) does not exhibit crossover, and are displayed as a function of x in Fig. 6. Fig. 6 suggests that the relationship between x and zcr is approximately linear, at least within the considered x range. It makes therefore sense to define the crossover location by the zcr /x ratio. For unstable stratification (open symbols), this ratio tends to be somewhat larger than for stable stratification (full symbols). Considering these two groups separately, zcr /x is estimated as the slope of linear regression in Fig. 6, giving 0.080 for unstable and 0.067 for stable stratification. The present data therefore imply that if the location of a single concentration input C is chosen such that z/x < zcr /x, then WindTrax will overestimate the emission rate, since B is underestimated and, from (6), Qc = (C − Cb )/B. Conversely, if z/x > zcr /x then Qc will be underestimated. Whether this statement can be generalised for other experiments remains to be investigated. 5.6. Emission rate calculations from multiple concentration inputs When WindTrax is used with several simultaneous concentration inputs Cj , as e.g. by L08, then the model computes a linear best-fit solution for Qc from an overdetermined system of equations. The resulting bias of Qc is then a weighted average of the biases of the individual Bj , with the weight in proportion to the relative weight that each Bj has in the optimal solution. In the case of L08, this plays out as follows. At rmin = 5 m, all simulated Bj are slightly larger than their observed counterparts, hence Qc from the linear optimisation should be less than Qc from the IHF method.

Fig. 6. Height of the crossover point between simulated and observed profiles of the normalised horizontal flux, H, as a function of average source–sensor distance. The error bars indicate the height range where simulated and observed H differ by <5%. The solid and dashed line are linear regressions for stable (R2 = 0.59) and unstable (R2 = 0.87) classes, respectively.

At rmin = 22 m, two measurement heights are located well below the crossover point, one near the crossover point, and two well above, so the distribution appears well-balanced with respect to bias, and Qc from the linear optimisation should agree with Qc from the IHF method within statistical error. At rmin = 57 m and 92 m, the majority of measurement heights are located below the crossover point, consequently Qc from the linear optimisation should exceed Qc from the IHF method. Table 2 generally confirms these predictions, showing a clear trend of QcBLS /QcIHF to increase from 0.91 for rmin = 5 m to 1.20 for rmin = 92 m. (Note that the QcBLS and QcIHF values are slightly different from those in L08 because of the stricter data selection, using about 20% fewer runs, but the trends with rmin are not changed by this.) A little surprising is that already for rmin = 22 m the BLS method overestimates the emission rate by 13% on average, despite the “well-balanced” distribution of concentration measurements. Perhaps this is an indication of WindTrax giving stronger weight to the bottom heights in the linear optimisation algorithm. 5.7. Comparison to line sensor measurements Additional points of the B(x,z) field can be obtained from the path-averaged concentration CL provided by the downwind (E) laser path. With this and background Cb (matched to Cb for the IHF method) as the only concentration inputs, WindTrax was executed BLS (additional subscript “L” to compute emission rate estimates QcL for “laser” or “line”) which were then compared to QcIHF , assumed again to be the true value. From (1) it follows that the true and

Table 2 CH4 emission rates obtained with the IHF method, and with the BLS method using concentrations measured at five heights as simultaneous model inputs for linear optimisation. Four classes are distinguished by the minimum horizontal source–sensor distance rmin , see Fig. 1. For each class, means of emission rates are shown (with standard error of the mean in parentheses), as well as the ratio of BLS mean to IHF mean (with propagated error). An emission rate of 200 g CH4 d−1 head−1 is equivalent to a surface flux of 38 ␮g CH4 m−2 s−1 (based on 58 cattle and a paddock area of 3500 m2 ). rmin (m)

No. of runs

QcIHF (g CH4 d−1 head−1 )

QcBLS (g CH4 d−1 head−1 )

QcBLS /QcIHF

5 22 57 92

107 132 48 20

250 (14) 194 (8) 168 (9) 179 (7)

228 (13) 219 (10) 193 (10) 214 (7)

0.91 (0.07) 1.13 (0.07) 1.14 (0.09) 1.20 (0.06)

1438

J. Laubach / Agricultural and Forest Meteorology 150 (2010) 1428–1442

Table 3 Ratios of the emission rates obtained with the BLS method, using concentrations measured by a line concentration sensor, to those obtained with the IHF method (assumed to be the “true” reference). Data are given as averages of four classes, labelled by the ratio of sensor height to average source–sensor distance z/x. In the third column, emission rates are first averaged for each class, then divided; in the fourth column, they are first divided then averaged. z/x

No. of runs

BLS QcL /QcIHF

BLS (QcL /QcIHF )

0.041 0.017 0.010 0.0075

49 37 3 5

1.31 1.85 1.76 1.29

0.88 1.63 1.22 0.79

in practice it is not advisable to choose very large x because then C − Cb will become very small and may not be resolved by the concentration sensor. BLS /Q IHF values in Perhaps more reliable than the individual QcL c Table 3 is the trend with z/x, which is the same for the mean of the ratios and the ratio of the means. In both cases it is indicated that as z/x decreases to somewhere between 0.017 and 0.010, the emission rate from the BLS method with a single concentration input becomes increasingly positively biased, while beyond this point, the bias decreases again. 6. Discussion

the simulated B are inversely proportional to the corresponding emission rates: Btrue Bsim

=

BLS QcL

QcIHF

(13)

A crucial difference to the vertical profile data is that there is no implicit consistency between CL and QcIHF , because the latter was obtained completely independent of the former. Consequently, the errors of both combine in the ratio (13), and it is of concern that BLS /Q . any potential bias of either (CL − Cb ) or QcIHF will bias QcL c BLS To avoid runs of dubious quality, it is demanded that QcL and QcIHF do not differ by more than 70% of the larger of the BLS /Q IHF < 10/3. This test is passed by 94 runs. two, i.e. 3/10 < QcL c Because of the small sample size, classes are formed solely according to xmin , ignoring differences in roughness or stability. Assuming x = xmin + 17 m gives x = 23.5, 58, 94 and 129 m for the four classes. The corresponding z/x ratios (using zL − d = 0.97 m) are given in BLS /Q IHF computed in two ways, as the ratio Table 3, along with QcL c of the mean emission rates and as the mean of the run-wise ratios. There are considerable differences between these two estimates, which are a consequence of large run-to-run variability of BLS /Q IHF . QcL c For all classes, z/x indicates that the measurement is taken below the crossover points observed in Figs. 4 and 5. Consequently, one BLS to overestimate the true Q in all cases. This is indeed expects QcL c found in the third column of Table 3 (ratio of the means), so the laser results are compatible with the results from the profile measurements. The fourth column (mean of the ratios) confirms the BLS for the two middle classes, while for the largest positive bias of QcL and smallest z/x, it indicates an underestimate of Qc rather than an overestimate. For the largest z/x (grazing strip nearest to the laser path) the likely explanation of the ambiguous result is that the actual z/x range spans from 0.023 to 0.15, so the nearer parts of the grazing strip represent z/x above the crossover point, while the farther parts represent z/x below the crossover point. Indeed, if one used the mean of the extreme z/x (instead of z divided by the mean of the extreme x) to define a representative value, this value would be 0.086, slightly above the crossover point, so the simulated B should be subject to very little bias. Curiously, the data of this class show a BLS /Q IHF <(>) 1 for u <(>) 5 m s−1 , trend with wind speed: usually QcL 5 c but the significance of this observation is unclear. For the smallest z/x (grazing strip farthest from the laser path) the ambiguous result in Table 3 may well be due to random error, given the small number of runs (N = 5). However, there is another possibility. Intuitively, one would expect that for very large x the z/x dependence disappeared, as both the modelled and the true B profile would smear out vertically to become an almost height-constant addition to the background concentration, BLS /Q IHF should approach 1. (Fig. 4 shows, for so theoretically, QcL c rmin = 92 m, that modelled and measured H at the bottom two heights are much closer together than at z3 and z4 , indicating the beginning of this approach process.) Though it may be tempting,

The mean emission rates obtained with BLS and IHF methods differ by 20% or less for each of the four groups distinguished by source–sensor distance in Table 2. This level of accuracy is of the same magnitude as found in other, stricter, tests of WindTrax (Flesch et al., 1995, 2004; McBain and Desjardins, 2005; Gao et al., 2009a,b). It seems, therefore, that the idealisation of a herd of animals as a homogeneous area source at ground level does not seriously affect the model’s applicability to infer emission rates. The same level of accuracy was not achieved with the laser data in Table 3, but this is partly explained by larger uncertainty of C–Cb with this instrument,2 and partly by the small number of runs. Comparison of the B and H profiles obtained with BLS and IHF provides a more detailed picture. These profiles differ in a systematic way that has been shown to introduce a small positiondependent bias to the inferred emission rates. The following subsections discuss possible causes of this bias. The first two of these are concerned with the effect of the crude idealisation of source locations and shapes applied here. The next two deal with implementation and accuracy of the flow parameters. The fifth discusses inner mechanics of the model, and the last reviews the literature for indications whether the presented findings are of generality. 6.1. Elevation of sources One assumption underlying this particular application of the BLS model is that the gas sources are assumed to be at ground level. Yet the real sources (cattle mouths) can be located anywhere between 0 and about 1 m height. The effect of this approximation is that some fraction of the gas plume, by virtue of originating from elevated sources, should in reality have been transported higher up than the BLS model simulates. Consequently, one would expect to find the centre of the measured gas plume somewhat higher than that of the modelled plume. This is the opposite of what is observed in this experiment. Further, the effect of an elevated release point is likely to diminish with increasing travelling distance, x. However, the observed mismatch between measured and simulated B and H profiles increases with x. Elevation of sources can therefore not explain this mismatch. Source elevation is most likely to have an observable effect near the sources, i.e. here for the classes with rmin = 5 m. It is indeed only for these classes (unstable and stable, top panels of Figs. 4 and 5) that the simulated B and H at the lowest two measurement heights

2 For both the IHF and BLS method, the resolution of emission rates is determined by the resolution of the downwind-upwind concentration differences, Cj – Cb . For the HP flame ionisation detector, the resolution was estimated by Laubach and Kelliher (2004) as about 3 ppbV. Observed Cj − Cb were one to three (typically two) magnitudes larger, so the typical relative error of Cj − Cb is 1–2%. For the Boreal Laser GasFinder MC, precision of the path-integrated concentration is 2 ppmV m, so for a 66 m long path, the uncertainty of a single concentration measurement is 30 ppbV. Hence, CL − Cb has a resolution of ca. 40 ppbV (from standard error propagation), and a typical relative error of order 20%.

J. Laubach / Agricultural and Forest Meteorology 150 (2010) 1428–1442

1439

(0.7 and 1.3 m) are slightly larger than their measured counterparts, with the consequence that there is no crossover point where measured and modelled H would agree. Elevation of some of the emission sources may therefore be the reason why simulated H near the ground exceeds measured H for rmin = 5 m. By contrast, at the upper three measurement heights, the explanation for overestimated H is likely to be the same as for the other rmin classes, discussed below. 6.2. Idealisation of point sources as an area source One may ask whether the idealisation of representing a collection of moving point sources (animals) within a fenced area by a uniform area source can introduce bias to the model results. This idealisation makes two assumptions: uniformity and equivalence of point and area sources. As to the first, the animals were generally observed to spread evenly across the paddock. They all were of similar age and size, and consequently their CH4 emission rates were quite uniform: on each of nine days of application of the SF6 tracer method, the standard deviation of the emission per head was between 10.2 and 14.9% of the herd mean (Knight et al., 2006). Hence, the uniformity assumption appears reasonable in the present experiment. For any given measurement situation, it should be carefully assessed if that is the case. Supposing even spread, each steer was responsible for the CH4 emissions from an area of 60 m2 , which implies a typical point source separation of about 8 m. This represents a rather lumpy source, compared to the model prescription of a uniform area source. If this mismatch between true and modelled source configuration mattered, then one would expect the largest observable effect near the sources, and that differences disappear at large x. Since in the present experiment the mismatch between measured and modelled H profiles increases with x, it seems unlikely that nonequivalence of point and area sources is the cause of this mismatch. Rather, the overall acceptable agreement of QcBLS and QcIHF suggests that the BLS model is a legitimate tool to estimate emissions from animal herds on a per-area basis. 6.3. Implementation of flow parameters Flesch et al. (2004) tested in a controlled-release experiment if QcBLS /Qc depended on how the flow field was parameterised: using wind profile data only, using sonic anemometer data only, or a combination of both. They found that if a profile of cup anemometer speeds was used and the turbulent flow parameters inferred from that, the BLS model tended to underestimate emission rates that were computed from concentration measurements at small x, and overestimated those from large x. They argued this was due to internal overestimation of u* and  w based on the usage of averaged wind speed rather than averaged alongwind vector component, an interpretation that is in agreement with the findings of the following section. Alternatively, they used a full set of the means, variances and covariances of the wind vector plus temperature flux as input parameters (obtained with the sonic anemometer), from which the model computed u* , L and z0 internally. With this approach, the dependence of QcBLS /Qc on fetch was strongly reduced. The same was true when a combination of sonic and wind profile data was used as input parameters, and results from the latter two approaches differed little. For the present data, the same was tested, by replacing the combination approach of L08 (in which wind speed from the cup anemometer at z5 together with L and a fixed z0 were prescribed) with the approach recommended by Flesch et al. (2004), using only sonic anemometer data. The latter does not lead to any significant changes in the shapes of the simulated profiles of B and H

Fig. 7. Ratio of friction velocity as simulated by WindTrax (from u5 , L, z0 ) to that measured with the sonic anemometer, dependent on wind speed.

(not shown). The simulated H profiles still systematically exhibit crossover points with the observed ones, at virtually the same heights as in Figs. 4 and 5, so these crossover points are not artefacts of a particular implementation method for u* and the wind profile. The sensitivity of QcBLS to errors of the parameters L, z0 ,  u /u* ,  v /u* , and  w /u* was tested by Laubach and Kelliher (2005a) and found to be minor, except for  w /u* , to which it was roughly proportional. This is addressed next. 6.4. Sensitivity to error of  w /u* The simplest possible explanation of an overestimated rate of vertical dispersion in the model is an overestimate of  w from the  w /u* input parameter. This could come about either from  w itself being overestimated by the sonic anemometer, or from a mismatch between the u* measured by the sonic anemometer, which is not used in the present implementation of flow parameters, and the u* computed from u5 , L and z0 , which is used instead and here, for brevity, called “simulated” u* (although strictly speaking, it is not simulated but enforced by MOST from the given parameters). If simulated u* overestimates true u* , then multiplication with measured  w /u* would result in too large  w . To investigate this possibility, simulated u* was compared to “measured” u* (i.e. computed from the wind covariances measured by the sonic anemometer, following Stull, 1988). The mean ratio of simulated to measured u* was 1.05, and the comparison revealed a slight trend of this ratio to increase with wind speed (Fig. 7). This supports the suspicion noted earlier that overestimated u* causes u2 , u3 and u4 in unstable stratification to be slightly overestimated. The scatter observed in Fig. 7 is compatible with the common statistical error of u* measurements with a sonic anemometer, of 10 to 15%. Errors of u* affect a number of parameters in the BLS model by error propagation, but the test described in the previous section shows that such errors cannot explain the systematic appearance of the crossover points between simulated and measured H profiles. Since u* is overestimated most commonly at high wind speeds,  w is then most likely to be overestimated, too. To investigate the effect of this on the profiles in Fig. 4, a sensitivity study was conducted with 12 runs in which 10 m s−1 < u5 < 12 m s−1 . These runs form the class with rmin = 57 m, z0 = 0.025 m and unstable strati-

1440

J. Laubach / Agricultural and Forest Meteorology 150 (2010) 1428–1442

Fig. 8. Same as Fig. 4 but for one class with high wind speeds only. The observed and the original simulated profiles are repeated from Fig. 4, centre panel, in identical colours and symbols, with the axis scaling optimised for this class. In addition, the results of a simulation with reduced  w /u* are shown (thin dashed lines).

fication. For this class, significant differences between measured and simulated H profiles had been observed (light colours in centre panel of Fig. 4). To compensate for biased u* , the original input parameters  u /u* ,  v /u* ,  w /u* were divided by the ratio of simulated to measured u* for these runs, and the BLS computations were repeated with these changed parameters. The resulting B and H profiles are shown in Fig. 8 (thin dashed lines). They are between the original simulation results and the measured profiles, indicating that the correction of  w /u* works in the right direction and that the profiles shapes are indeed sensitive to the details of the vertical dispersion parameterisation. However, the likely overestimate of  w /u* in the original model runs accounts only for about 40% of the gap between measured and modelled H profile in this class, and by implication for less than 40% in the other classes, where the u* bias was smaller. The remainder of the H mismatch requires a different explanation. It should also be noted that the correction of  w /u* has very little effect on the location of the crossover point between measured and modelled H profile. 6.5. Potential error sources in the BLS model Since the previous four sections, dealing with idealisations and implementations of this particular experiment, have not identified a likely cause for the apparent overestimate of vertical dispersion speed, it is briefly discussed here which assumptions of the BLS model itself may cause such an overestimate. The correctness of the model depends on the correctness of (1) the MOST parameterisations used, which set the height and stability dependence of the flow statistics  u ,  v ,  w and ε, (2) the specification of the coefficients ˛k as functions of these flow statistics, u* and ˇ, and (3) the value of the universal constant C0 . 6.5.1. Coefficients in the Thomson model equations Addressing (2) first, Thomson (1987) showed that there is only one LS model for vertical dispersion in a homogeneous surface layer that satisfies the well-mixed condition. From this it follows that the only freedom in the choice of ˛3 concerns the value of ˇ, given by Eq. (5) and discussed below. By contrast, the horizontal acceleration coefficients, ˛1 and ˛2 , are not uniquely constrained, and Thomson’s expressions, used in WindTrax, may not be the best possible representation of reality. Yet, horizontal motion in the surface layer is controlled by processes of the outer boundary layer. In most experimental situations it is impossible to know these precisely, and therefore futile to attempt formulation of a “better” horizontal dispersion model. Model error with respect to horizontal disper-

sion would affect the vertical profiles of the present experiment, because too much or too little CH4 may have been carried in the x or y direction prior to arrival of the simulated air parcel at the point of measurement. Error of ˛2 would not affect the results of experiments where a line concentration measurement integrates across the whole width of the gas plume, which is why Flesch et al. (2004) recommend line sensors over point sensors. However, this advantage must in practice be balanced against the accuracy of either type of instrumentation to measure the upwind-downwind concentration differences (in the present experiment, the CH4 point sensor was clearly more accurate than the line sensor). In order to explain the result of a systematic crossover point between simulated and observed H profiles, though, the horizontal dispersion model error would need to have a strong height dependence, which could then cause underestimated concentrations at small z simultaneously with overestimated concentrations at larger z. The author believes this to be unlikely. 6.5.2. Parameterisations using similarity theory The MOST expressions for the dimensionless parameters  u /u* ,  v /u* ,  w /u* and ␧ consist of their neutral-stability values and modifying factors that are functions of z/L. The first three are direct inputs to the model, at the location of the sonic anemometer. These inputs supersede the neutral default values of WindTrax and are correct within measurement uncertainty at this location. The extrapolation to other locations assumes no height dependence for  u /u* and  v /u* . Given that the QcBLS outputs were found insensitive to variation of the  u /u* and  v /u* inputs (Laubach and Kelliher, 2005a), it appears extremely unlikely that the spatial extrapolation of these input values could cause noticeable model error. The height dependence of  w /u* is empirically well-established: constant for stable and proportional to (z/L)1/3 for unstable stratification (Kaimal and Finnigan, 1994), therefore unlikely to cause model bias, either. The MOST parameterisation of ␧ =  ε z/u* 3 as a function of z/L is less certain. The expression used in WindTrax (Eq. A13 of Flesch et al., 2004) places the minimum of ␧ at z/L = 0, with a value of 1. This is in agreement with results of Wyngaard and Coté (1971), but not with the more recent studies of Kader and Yaglom (1990), Edson and Fairall (1998) and Smedman et al. (2007), which all show a minimum value of 0.7 (±0.15) near z/L≈–0.2(±0.1). It is thus possible that for moderately unstable stratification, WindTrax overestimates ε by up to 30%. A comparable overestimation of ␧ on the stable side is possible, too, where Pahlow et al. (2001) found ␧ = 0.61 + 5 z/L, rather than the conventional relationship ␧ = 1 + 5 z/L (used in WindTrax). The correct functional form of

J. Laubach / Agricultural and Forest Meteorology 150 (2010) 1428–1442

␧ is still unknown, and new high-quality experiments to determine it would be desirable. Such experiments should be designed to test the supposition of Laubach and McNaughton (2009) that in unstable stratification z/L is the wrong parameter for a functional dependence and should be replaced by z/zs , where zs is the depth of the surface friction layer. However, if WindTrax overestimated ε, this would lead to overestimated ˇ and that, in turn, to somewhat underestimated ˛k . The latter means that the modelled gas plume would disperse less rapidly than the true one, in all three dimensions. The effect on vertical dispersion would then be opposite to what was observed in this experiment. 6.5.3. The Kolmogorov constant Apart from ε, as discussed above, the model coefficient ˇ depends on the value of the universal constant C0 . For WindTrax, this was set to be consistent with the LS model of Wilson et al. (1981). In that model, a Lagrangian timescale TL = A z/ w was optimised for concentration profiles of the Project Prairie Grass Experiment (PPG), at neutral stability. This yielded A = 0.5, which implied Sc = 0.64. Equality of the gas diffusivities in the two models is obtained if: 2 ( 4 + u4∗ ) = Aw z C0 ε w

(14)

(Wilson et al., 2001). The value of A therefore sets the value of C0 ε, as well as that of Sc, both in inverse proportion. (Under the assumption that ␧ = 1 it then follows that C0 = 4.41 is equivalent to A = 0.5, yet as noted above, empirical evidence for this assumption is ambiguous.) Looking beyond PPG, the correct values of the turbulent Schmidt and Prandtl (Pr) numbers in the atmospheric surface layer are far from settled (Wilson et al., 2001; Flesch et al., 2002; Laubach and Kelliher, 2004). Reviews that include laboratory pipe and channel flows report that the most commonly found values of Sc and Pr are between 0.7 and 0.9 (Launder, 1978; Kays, 1994). The lower value obtained by Wilson et al. (1981) was a good fit for PPG, where “good” needs to be understood as demonstrably better than the two alternative values tested. However, revisiting Figs. 2, 5 to 7, and 11 of Wilson et al. (1981) suggests that such good fit, while correctly observed, is not very tightly constrained, and 10–20% decrease of A, and therefore increase of C0 ε, may still be compatible with both PPG and their own validation experiment. If C0 ε in WindTrax was altered accordingly, to simulate Sc > 0.7, then the vertical gas diffusivity would effectively be reduced and the vertical dispersion slowed down. Such a change would be of the right magnitude to explain the present results. On this basis it is suggested that too small C0 ε in WindTrax is the most probable explanation for the observed crossover points in the H profiles. It is cautioned that, in the light of the imperfections discussed in the previous sections, the present data do not constitute firm evidence. For a real test of the value of C0 , both better knowledge of ε and more data of the quality of PPG would be required. Since an imperfectly tuned C0 ε would have affected results of other experiments where WindTrax was used, some of these are reviewed next. 6.6. Evidence from other studies Laubach and Kelliher (2005a) analysed the Qcj /QcIHF ratios in three previous field experiments determining CH4 emissions from cattle, using 5-point profiles and similar instrumentation as in the present study. Some of their findings agree with the present data, others do not. In those three experiments the zj /x for the lowest four measurement heights were always less than 0.08, so from the present results one would expect the BLS method to overestimate Qc . Indeed Laubach and Kelliher (2005a) found Qcj /QcIHF ≈

1441

1.2(±0.2) for these heights, in excellent agreement with the expectation. By contrast, mean z5 /x > 0.06 in the majority of their setups, so one would expect Qc5 /QcIHF to be near or below 1, but it was found to be 1.3 (±0.2) instead. Also, for a given run one would generally expect Qcj to decrease with increasing height, but this was quite commonly not the case. In those three experiments the number of cattle was 9 times larger and the grazed paddock area typically 20 times larger than in the present study. Consequently, the source area contact distance, x, was typically between 100 and 300 m, 3–9 times larger than here. This means that for almost any mast location and measurement height used by Laubach and Kelliher (2005a), the nearer parts of the grazed area would have represented z/x above the crossover point, while the farther parts would have represented z/x below the crossover point, and the observed profiles would have been superpositions of quite wide z/x intervals, making it difficult to detect clear trends. McBain and Desjardins (2005) investigated the gas “recovery rate” of the BLS method (WindTrax) in a controlled-release experiment. They used five path-averaging concentration sensors at heights from 1.0 to 3.0 m (every 0.5 m), and they investigated if recovery rates depended on the average source–sensor distance (five classes of x, ranging from 13 to 38 m). As the path length (50 or 100 m) was much larger than the diameter of the release grid (2 m), the complete crosswind integral of the gas plume traversed the sensor paths, so inaccuracies of lateral dispersion modelling could not affect the accuracy of the results. They found at all x, except the longest, that the simulated emission rate decreased monotonically with increasing z. The decrease was steepest for shortest x (covering the largest absolute range of z/x, from 0.07 to 0.20). This trend agrees with the present study. For the class with longest x, McBain and Desjardins (2005) had only three valid runs, so the results of this class cannot be considered statistically reliable. Their data also fit excellently with the crossover point estimates from the present data, and with the notion that sensor placement at z/x above (below) the crossover point leads to under-(over-)estimates of Qc : for z =1.0 m, z/x was always below 0.07, and indeed they found QcBLS /Qc > 1 for the averages of all x classes; conversely, for z = 2.5 or 3.0 m, z/x was always above 0.08, and then QcBLS /Qc < 1 was true for the averages of all respective x classes. Flesch et al. (2004), also in a controlled-release experiment, used a single line concentration sensor at z = 0.9 (±0.1 m) and various x from 0 to 90 m. Ignoring x < 10 m because of large scatter, the bottom panel of their Fig. 6 shows an increasing (less than linear) trend of QcBLS /Qc with increasing x, from ca. 0.8 at 10 m to ca. 1.3 at 60 m, where the trend reverses towards ca. 0.7 at 90 m. The trend between 10 and 60 m suggests a crossover point (QcBLS /Qc = 1) at 16 (± 3) m, i.e. for zcr /x = 0.056 (±0.01), which is in good agreement with the estimates from the present vertical profile data. The figure of Flesch et al. (2004) is also in good qualitative agreement with the results from our line concentration measurements (Table 3), alleviating concerns about potential bias of either (CL − Cb ) or QcIHF BLS /Q . compromising the validity of QcL c A consequence of the crossover point being lower for stable than for unstable stratification is that a sensor at a given location is more likely to be above the crossover point for stable than for unstable stratification. Hence, stable runs should lead to underestimates of Qc more often than unstable runs. The present dataset is unsuitable to assess stability dependence since the range of stability variation is too small. Neither Flesch et al. (2004) nor McBain and Desjardins (2005) found a significant stability dependence. Gao et al. (2009b) found the opposite trend, an increase of mean QcBLS /Qc with increasing z/L. In their experiment, z = 1.2 m and x = 24 (±4) m, so z/x was close to where the above findings would suggest the crossover point. It is unclear how their result can be reconciled with the present data. It is possible that stability may affect the errors of other parameters (by either measurement error or incorrect

1442

J. Laubach / Agricultural and Forest Meteorology 150 (2010) 1428–1442

parameterisation), and that such indirect biases may overwhelm any direct stability dependence. Measurements of CH4 and NH3 emissions have more commonly shown daytime emissions to be larger than nighttime emissions (the author’s own unpublished observations as well as Flesch et al., 2007; McGinn et al., 2007; van Haarlem et al., 2008) – with the reverse reported only by Loh et al. (2008). However, in all these cases the observed temporal patterns fit well with the animals’ activity patterns, and are therefore interpreted as real, rather than as caused by a stability-dependent bias of the BLS method. 7. Conclusions In the present experiment, CH4 emission rates obtained with the BLS method differ typically by 10–20% to those from concurrent mass-budget measurements. This is in agreement with previous tests of the accuracy of WindTrax. It seems thus that the idealisation of a herd of animals as a homogeneous area source at ground level does not invalidate the inferred emission rates. The profile comparison suggests that WindTrax may overestimate the speed of vertical dispersion. As a consequence, for this experiment an ideal z/x ratio exists where the modelled emission rate is unbiased. Its value is about 0.080 in unstable and 0.067 in stable stratification. The generality of this observation should be tested with data from other experiments. Based on the above discussion, it is argued that the most likely cause of overestimated vertical dispersion is an underestimate of C0 ε. This is not really a “flaw” of the model, but reflects our current limited understanding of the TKE dissipation rate in the surface layer, as well as of the dissimilarity of momentum and gas transport due to the action of different kinds of eddies. If this interpretation is correct, then previously published WindTrax experiments should have been subject to similar trends as the observed one. Because this trend is small, it may in many experiments have been concealed in the statistical uncertainty of order 20%. The only other study comparing WindTrax results from several concentration sensor heights simultaneously is that of McBain and Desjardins (2005), which shows a trend compatible with the present one. Future verification of these conclusions is desirable, but requires careful layout with multiple concentration measurements at considerable accuracy, as well as a full and accurate characterisation of the surface layer flow. Acknowledgements This work was funded by New Zealand’s Foundation of Research, Science and Technology (FRST). The author thanks John Wilson and Tom Flesch for valuable comments on a draft manuscript and subsequent discussions. References Denmead, T., Chen, D., Turner, D., Li, Y., Edis.B., 2004. Micrometeorological measurements of ammonia emissions during phases of the grazing rotation of irrigated dairy pastures. In: SuperSoil 2004 – 3rd Australian and New Zealand Soils Conference, December 5–9, 2004, Sydney, University of Sydney, p. 8. Edson, J.B., Fairall, C.W., 1998. Similarity relationships in the marine atmospheric surface layer for terms in the TKE and scalar variance budgets. J. Atmos. Sci. 55, 2311–2327. Flesch, T.K., Wilson, J.D., Yee, E., 1995. Backward-time Lagrangian stochastic dispersion models, and their application to estimate gaseous emissions. J. Appl. Meteorol. 34, 1320–1332. Flesch, T.K., Prueger, J.H., Hatfield, J.L., 2002. Turbulent Schmidt number from a tracer experiment. Agric. Forest Meteorol. 111, 299–307. Flesch, T.K., Wilson, J.D., Harper, L.A., Crenna, B.P., Sharpe, R.R., 2004. Deducing ground-air emissions from observed trace gas concentrations: a field trial. J. Appl. Meteorol. 43, 487–502.

Flesch, T.K., Wilson, J.D., Harper, L.A., Todd, R.W., Cole, N.A., 2007. Determining ammonia emissions from a cattle feedlot with an inverse dispersion technique. Agric. Forest Meteorol. 144, 139–155. Gao, Z., Desjardins, R.L., Flesch, T.K., 2009a. Comparison of a simplified micrometeorological mass difference technique and an inverse dispersion technique for estimating methane emissions from small area sources. Agric. Forest Meteorol. 149, 891–898. Gao, Z., Mauder, M., Desjardins, R.L., Flesch, T.K., van Haarlem, R.P., 2009b. Assessment of the backward Lagrangian Stochastic dispersion technique for continuous measurements of CH4 emissions. Agric. Forest Meteorol. 149, 1516–1523. Harper, L.A., Flesch, T.K., Powell, J.M., Coblentz, W.K., Jokela, W.E., Martin, N.P., 2009. Ammonia emissions from dairy production in Wisconsin. J. Dairy Sci. 92, 2326–2337. Kader, B.A., Yaglom, A.M., 1990. Mean fields and fluctuation moments in unstably stratified turbulent boundary layers. J. Fluid Mech. 212, 637–662. Kaimal, J.C., Finnigan, J.J., 1994. Atmospheric Boundary Layer Flows. Oxford University Press, New York, 289 pp. Kays, W.M., 1994. Turbulent Prandtl number – where are we? ASME J. Heat Transfer 116, 284–295. Knight, T., Clark, H., Molano, G., Cavanagh, A., 2006. Comparison of the direct measurement of total methane production from a herd of cattle with measurements based on individual animals. Report prepared for Landcare Research. AgResearch, Palmerston North, New Zealand, 13 pp. Laubach, J., Kelliher, F.M., 2004. Measuring methane emission rates of a dairy cow herd by two micrometeorological techniques. Agric. Forest Meteorol. 125, 279–303. Laubach, J., Kelliher, F.M., 2005a. Measuring methane emission rates of a dairy cow herd (II): results from a backward-Lagrangian stochastic model. Agric. Forest Meteorol. 129, 137–150. Laubach, J., Kelliher, F.M., 2005b. Methane emissions from dairy cows: Comparing open-path laser measurements to profile-based techniques. Agric. Forest Meteorol. 135, 340–345. Laubach, J., Kelliher, F.M., Knight, T., Clark, H., Molano, G., Cavanagh, A., 2008. Methane emissions from beef cattle – a comparison of paddock- and animalscale measurements. Aust. J. Exp. Agric. 48, 132–137. Laubach, J., McNaughton, K.G., 2009. Scaling properties of temperature spectra and heat-flux cospectra in the surface friction layer beneath an unstable outer layer. Boundary-Layer Meteorol. 133, 219–252. Launder, B.E., 1978. Heat and mass transport. In: Bradshaw, P. (Ed.), Topics in Applied Physics, vol. 12, Turbulence. Springer, Berlin (Chapter 6). Loh, Z., Chen, D., Bai, M., Naylor, T., Griffith, D., Hill, J., Denmead, T., McGinn, S., Edis, R., 2008. Measurement of greenhouse gas emissions from Australian feedlot beef production using open-path spectroscopy and atmospheric dispersion modelling. Aust. J. Exp. Agric. 48, 244–247. McBain, M.C., Desjardins, R.L., 2005. The evaluation of a backward Lagrangian stochastic (bLS) model to estimate greenhouse gas emissions from agricultural sources using a synthetic tracer source. Agric. Forest Meteorol. 135, 61–72. McGinn, S.M., Flesch, T.K., Harper, L.A., Beauchemin, K.A., 2006. An approach for measuring methane emissions from whole farms. J. Environ. Qual. 35, 14–20. McGinn, S.M., Flesch, T.K., Crenna, B.P., Beauchemin, K.A., Coates, T., 2007. Quantifying ammonia emissions from a cattle feedlot using a dispersion model. J. Environ. Qual. 36, 1585–1590. Pahlow, M., Parlange, M.B., Porté-Agel, F., 2001. On Monin-Obukhov similarity in the stable atmospheric boundary layer. Boundary-Layer Meteorol. 99, 225–248. Schuepp, P.H., Leclerc, M.Y., MacPherson, J.I., Desjardins, R.L., 1990. Footprint prediction of scalar fluxes from analytical solutions of the diffusion equation. Boundary-Layer Meteorol. 50, 355–373. Smedman, A.-S., Högström, U., Hunt, J.C.R., Sahlée, E., 2007. Heat/mass transfer in the slightly unstable atmospheric surface layer. Quart. J. R. Meteorol. Soc. 133, 37–51. Sommer, S.G., McGinn, S.M., Flesch, T.K., 2005. Simple use of the backwards Lagrangian stochastic dispersion technique for measuring ammonia emission from small field-plots. Europ. J. Agronomy 23, 1–7. Stull, R.B., 1988. An Introduction to Boundary Layer Meteorology. Kluwer Acad. Publ., Dordrecht, 666 pp. Thomson, D.J., 1987. Criteria for the selection of stochastic models of particle trajectories in turbulent flows. J. Fluid Mech. 180, 529–556. van Haarlem, R.P., Desjardins, R.L., Gao, Z., Flesch, T.K., Li, X., 2008. Methane and ammonia emissions from a beef feedlot in western Canada for a twelve-day period in the fall. Can. J. Anim. Sci. 88, 641–649. Wilson, J.D., Sawford, B.L., 1996. Review of Lagrangian stochastic models for trajectories in the turbulent atmosphere. Boundary-Layer Meteorol. 78, 191–210. Wilson, J.D., Flesch, T.K., Harper, L.A., 2001. Micro-meteorological methods for estimating surface exchange with a disturbed windflow. Agric. Forest Meteorol. 107, 207–225. Wilson, J.D., Thurtell, G.W., Kidd, G.E., 1981. Numerical simulation of particle trajectories in inhomogeneous turbulence, III: Comparison of predictions with experimental data for the atmospheric surface layer. Boundary-Layer Meteorol. 21, 443–463. Wilson, J.D., Yee, E., Ek, N., d’Amours, R., 2009. Lagrangian simulation of wind transport in the urban environment. Q. J. R. Meteorol. Soc. 135, 1586–1602. Wyngaard, J.C., Coté, O.R., 1971. The budgets of turbulent kinetic energy and temperature variance in the atmospheric surface layer. J. Atmos. Sci. 28, 190–201.