Inferring Biogeochemical Sources and Sinks from Atmospheric Concentrations: General Considerations and Applications in Vegetation


M. R. Raupach 1. Introduction 41

CSIRO Land and Water. 2. Scalar and Isotopic Molar Balances 43

Australia Inverse Methods for Inferring Scalar Sources and Sinks in Canopies 47

4. Inverse Methods and Isotopes in Canopies 53

5. Summary and Conclusions 55

Appendix A: Single-Point Eulerian Molar Balance Equations 57

Appendix B: Lagrangian Molar Balance Equations and Green's Functions 57

References 58

This chapter is a review of the principles and application of atmospheric inverse methods in vegetation canopies. These methods enable the source-sink distributions of biogeochemically active entities (such as water, carbon dioxide, non-C02 greenhouse gases, and aerosols) to be inferred from measurements of their atmospheric concentrations. The chapter covers four topics. First, canopy-scale inverse methods are placed in the context of atmospheric inverse methods in general, at scales from small chambers to the globe. Next, because these methods depend on mass conservation, the balance equations for scalar mole fraction (c), iso-concentration (cS) and isotopic composition (<5) are analyzed in single-point Eulerian form and in Lagrangian form. This leads to the third topic, inverse Lagrangian methods for inferring source-sink profiles from concentration measurements in vegetation canopies. The theory of the approach is reviewed and results from several field experiments are summarised, showing that this approach is a practically useful tool for inferring the canopy source-sink distributions of scalars such as water vapor, heat, CO; and ammonia. However, there is a continuing need for improvement in the knowledge of the turbulence field in the canopy. The fourth topic is the extension of inverse Lagrangian analysis to describe the relationship between profiles of isotopic composition in canopy air and the profiles of isotopic sources and sinks. Lagrangian analysis is shown to provide a general basis for the Keeling plot and the Yakir-Wang expression for distinguishing assimilation and respiration, in the case where the isotopic composition of the exchanged scalar (<5P) is constant through the canopy. Using the Inverse Lagrangian approach, this analysis is extended to explain air isotopic composition profiles in circumstances where <5P is strongly nonuniform through the canopy.

1. Introduction

Exchange between the earth's surface and the atmosphere is a crucial part of the cycles of almost all biogeochemically active entities, including water, carbon dioxide, methane, oxides of nitrogen, volatile organic compounds, and others. As surface-atmosphere exchanges occur, the atmospheric concentrations of these entities are altered in both space and time, providing a transient imprint of both the magnitudes and the distributions of their sources and sinks at the surface. The information in this imprint can be used to infer the source-sink distribution of an entity from measurements of its atmospheric concentration field, by "inversion" of a "forward model" which specifies the concentration field in terms of the source-sink distribution. The forward model is based (in


Copyright V: 200! ■ Academic Preis. All rights ot'reproduction in ü;v form reserved.

different ways according to the application) on mass-balance principles augmented by knowledge of atmospheric advection and dispersion processes, which together determine the concentration field produced by a specified source-sink distribution. Such methods can be generically called atmospheric inverse methods. This chapter reviews some general principles underlying atmospheric inverse methods, and then applies these principles to the specific problem of inferring biogeochemical sources and sinks in vegetation canopies, from measurements of concentration and isotopic composition profiles in the air.

The plan of the chapter is as follows. This introductory section surveys atmospheric inverse methods from small-chamber to global scales, indicating the commonalities and differences among methods at various scales. Because of the fundamental reliance of atmospheric inverse methods on mass conservation, Section 2 discusses the generic balance equations for the mole fraction and isotopic composition of a scalar entity, considering both Eulerian (fixed) and Lagrangian (fluid-following) frameworks. The chapter then moves to its specific focus on vegetation canopies: Section 3 describes the methodology and application of methods for inferring scalar source-sink distributions from measured concentration profiles, and Section 4 extends this treatment to isotopic composition. Section 5 provides a summary and conclusions.

We begin with a general survey of atmospheric inverse methods. In all cases, the broad goal is to use concentration measurements in the air, together with information about atmospheric flow, to infer sources and sinks of entities at the earth's surface. Since the key concentration observations are remote from the surface sources and sinks, this entire class of methods relies explicitly or implicitly on an atmospheric mass or molar balance for the entity being measured, within a specified control volume. Such a balance can be either in an Eulerian framework, in which the control volume is fixed in space, or in a Lagrangian framework, in which the control volume moves with the flow. Considering the Eulerian framework first, the molar balance for a scalar entity can be written informally as dc Of c(%,) - c(x0

pz I

where c is the average mole fraction (specific concentration) of the entity in a control region extending from x0 to x, in the direction of the mean flow and from the surface to height z,; u is the mean flow velocity; p is air density; F is the vertical flux density of the scalar; and </> is the source (or sink) of the scalar in the control region, including contributions from fluxes at the surfacc. All quantities are suitably averaged in space (see next section and appendixes for details). The four terms respectively represent (I) storage change, (II) advection, (111) flux out of the top of the control region, and (IV) sources within the region. The aim is to infer the source lerm </> by choosing the control region so that most other terms in the balance arc small, leaving only one other major term (sometimes two) which can be measured to infer c/>. Such an

Eulerian approach is used at the scales of small chambers, small plots, vegetation canopies, the atmospheric surface layer, the atmospheric boundary layer (both convective and stable), large regions, and the globe. Reviews of micrometeorological methods are offered by Denmead and Raupach (1993), Denmead (1994, 1995), and Denmead et al. (1999), while aspects of methods at global scales are described by Enting et al. (1995), Ciais and Meijer

(1998), Bousquet et al (1999a, b), Enting (1999), and Rayner et al.

Commonalities and differences between these methods are highlighted by grouping the methods not by scale, but according to the term(s) in Eq. (1) used to infer </>. Table 1 shows such a grouping, in which four classes of Eulerian method are identified:

1. Methods based on spatial gradients in the direction of diffusion: These include the standard gradient and Bowen-ratio methods of surface-layer flux measurement in micrometeo-rology. The key assumption is that terms I and II are small, so that the source term <fi (in this case identifiable with the flux at the surface) equals the flux F at a measurement height z,.

2. Methods based on spatial gradients in the direction of advection: These include the small-plot techniques reviewed by Denmead (1994, 1995) and also chambers operated in a continuously ventilated, quasi-steady-state mode. The control region is designed so that all terms except II and IV are small.

3. Methods based on temporal gradients: This approach is used at all scales, including small chambers operated in fully enclosed mode for short periods, boundary-layer budget methods in both the daytime convective boundary layer (CBL) and the nocturnal stable boundary layer (SBL), and global atmospheric budget methods. The dominant terms are I and IV.

4. Flux-resolving, fluctuation-based methods: In these methods the vertical turbulent eddy flux F is measured as pw'c', where overbars and primes denote time means and fluctuations, respectively, and w is the vertical velocity component. They include eddy-covariance and eddy-accumulation methods for measuring fluxes from towers and aircraft. They are not normally regarded as atmospheric inverse methods, since that term usually refers to methods which relate a mean concentration field to a source density or surface flux distribution.

Lagrangian inverse methods also rely on a balance equation. In this case, the principle is to consider individual "marked fluid particles" which (with sufficient knowledge of the velocity field) can be followed as they move. The balance equation in a Lagrangian framework is simply

where D/Dtis, the material derivative, or time derivative following the motion of the fluid particle. This says that the concentration

TABLE 1 Atmospheric Inverse Methods, Grouped According to the Term(s) in the Scalar Conservation Equation Used as Surrogates for the Source-Sink Term at the Surface.


Quantity Measured

Key Assumptions and Constraints

Spatial Footprint (X) and Time Resolution ( T)

Eiilerinn methods based on spatial gradients in the direction of diffusion

Gradient, aerodynamic, Surface flux

Bowen ratio

Eulerian methods based on spatial gradients in the direction of advection

Ventilated chamber Surface flux

Small plot Surface flux

Eulerian methods based on temporal gradients

Enclosed chamber CBL or SBL budget

Surface flux

Surface flux (regional average)

Global budget Surface flux (global average)

Eulerian flux-resolving, fluctuation-based methods

Kddy covariance, eddy accumulation, relaxed eddy accumulation

Lagrangian methods

Inverse Lagrangian methods in canopies Global or regional synthesis inversion

Surface flux (sensor above canopy) or flux profile (sensor within canopy)

Source-sink profile <f>(z) Surface flux (regional distribution)

Steady, horizontally uniform

Chamber does not alter flux Control region is deep enough to neglect flux through top

Chamber does not alter flux Semi-Lagrangian and Eulerian averages are equivalent

Well-mixed global atmosphere

Steady, horizontally uniform

Steady, horizontally uniform

X = 1 m, T = 30 min T = 0.3 day X = 100 km (CBL) X = 10 km (SBL) X = globe, T > 1 year

X = 100 m, T= 30 min X = 3000 km, T = 1 month in the fluid particle changes in response to the sources or sinks it encounters along its path. Two problems can then be defined. The "forward" problem is to find the concentration field from the source-sink distribution and is solved by tracking the scalar added to (or removed from) each fluid particle by the sources or sinks and then calculating statistically the points to which the flow moves the scalar. To solve the opposite "inverse" problem, finding sources and sinks from concentrations, it is necessary to find the concentrations produced by a large number of point sources and then find the mix of point sources or sinks which best matches the observed concentration field.

Lagrangian inverse methods are listed as a fifth group in Table 1. There are several methods in this class, including the "synthesis inversion" methods for inferring the global distribution of atmospheric sources and sinks for entities such as C02 and its isotopes (Enting et al, 1995; Bousquet et al, 1999a, b; Enting, 1999a, b; Rayner et al., 1999) and the "inverse Lagrangian" methods for inferring scalar source-sink distributions in vegetation canopies which are the primary specific focus in this paper.

In this brief survey, atmospheric inverse methods are defined broadly to include both Lagrangian and Eulerian approaches because their common foundation is the inference of sources and sinks from atmospheric concentration measurements. The two broad streams offer different strengths: Eulerian approaches, by using tightly defined control volumes, can give quite precise estimates of average sources or sink strengths (ib) within the control volume, subject to the requirement that the neglected terms in Eq. (1) are indeed small. Their demands for information about the velocity field are usually modest, but they inherently produce results for (f> which are averaged through the control volume. By contrast, Lagrangian approaches offer much more resolution of the space-time distribution of </> at the expense of far greater requirements for information about the velocity field.

2. Scalar and Isotopic Molar Balances

2.1 General Principles

This section considers the scalar and isotopic balance equations in general, with attention to the source terms for scalars and isotopes which arise in vegetation canopies. After general principles are set out in Section 2.1, the Eulerian framework is described in Section

2.2 and the ensuing Lagrangian framework in Section 2.4. Section

2.3 discusses the source terms for isotopes.

We are concerned with the molar balance of a scalar entity with mole fraction c, with a minor isotopic constituent such as "C or lsO in CO,. The isotopic ratio R is the molar ratio of the minor (heavier) to the major (lighter) isotope, and the isotopic composition 8 is the normalized departure of R from its value R. in a standard reference material, so that 8 = R/R, — 1 and R = R,( 1 — 8). Henceforth c will denote the mole fraction of the major isotope unless otherwise stated, so the mole fraction of the minor isotope is cR. The molar balance applies in general to a region with volume V{t) enclosed by a bounding surface S(t) which may be moving, so the region may deform, expand, or contract. Usually, a part (say

So) of the surface 5 occurs in the open air, and the remainder (say SF) coincides with plant, soil, or water surfaces bounding V. Thus, S = So + ¿i>- It is useful to consider a general moving control volume in order to support Lagrangian or semi-Lagrangian applications in which V or some of its boundaries move with the flow, and also because major applications occur in the atmospheric convec-tive boundary layer (CBL) which grows through the day.

For the major isotope, the molar balance equation is d dt pc

where F is the flux density of the major isotope (mol m~2 s~')> p is the molar air density, u and v are respectively the vector velocities of the fluid and the moving surface S, and n is the inward unit normal vector on S (pointing into V). The term pc(u — v) is the advective scalar flux. The first integral covers the open-air part S0 of S and the second the surface part SP. The equation for the overall scalar entity is the same as Eq. (3), with appropriate redefinition of c and F.

Equation (3) takes a similar form for both time-averaged and instantaneous quantities. In the time-averaged case, F is a turbulent eddy flux pu'c', where overbars and primes denote time means and fluctuations, respectively. In the instantaneous case, F is a molecular diffusive flux which in practice can be neglected in the open air (on S0) relative to fluxes arising from fluid motion (the high Peclet number approximation). In this case the fluid-motion fluxes appear in the term pc(u — v), which is unaveraged and includes transport by turbulent fluctuations as well as by the mean flow. In contrast with the situation in the open air (on SQ), molecular fluxes can never be neglected at solid boundaries (SP), where they are responsible for all the scalar transport.

Turning to the minor isotope, the ratio of minor to major isotope fluxes in the open air (on SQ) is the same as the isotopic concentration ratio R, because there is negligible discrimination by fluid motion and mixing. Hence, on SQ, the advective flux for the minor isotope is pcR(u — v) and the molecular or turbulent flux is RF. However, transport across plant or other solid surfaces (Sp) does discriminate between minor and major isotopes. The flux of minor isotope across these surfaces is RPF, where is the isotopic ratio of the scalar exchanged (transported across Sp) by the flux F. Hence the molar balance for the minor isotope can be written in either of the forms d dt d dt pcR dx =

RpF- ndS

Vif J

where = R\JR, ~ 1 is the isotopic composition of the scalar exchanged across plant surfaces. The first form reverts to Eq. (3), when R = Rv = 1.

Balance equations are needed mainly for three quantities: c, cS, and 8. The quantity c8 = (cR — cR,)/R, is the "iso-concentration," a linear combination of the mole fractions for the minor and major isotopes which measures the departure of the minor isotope mole fraction cR from its value cR. when the isotopic composition is that of the reference material.

2.2 Single-Point Eulerian Equations

In the limit V —» 0 at a single fixed point x, Eqs. (3) and (4) reduce to conventional Eulerian differential conservation equations. Anticipating application in vegetation canopies, we consider that some of the region V contains plant leaves and stems. Provided that the canopy is sufficiently finely textured to be regarded as a porous continuum from the standpoint of the airflow, a scalar source density in Eq. (3) can be defined by


where ^(x, t) is the source density for absolute molar concentration and 4>{x, t) the source density for mole fraction. If there is no porous canopy at x, then <t> = = 0. The equivalent "iso-source" term for the iso-concentration, in Eq. (4), is

The "continuum porous canopy" assumption underlying Eqs. (5) and (6) is a simple approximation which gives the same result in the present case as the more formally rigorous procedure of taking finite volume averages over thin horizontal slabs with horizontal length scales large enough to average over the local heterogeneity of the canopy (Wilson and Shaw, 1977; Raupach and Shaw, 1982; Finnigan, 1985).

It is shown in Appendix A that in the limit V —* 0, Eqs. (2) to (4) yield the following conservation equations for scalar mole fraction (c), iso-concentration (câ), and isotopic composition (8) at a fixed point x:

Equations (7) and (8) (for c and c<5) include the usual time derivative and advection terms (the bracketed terms on the left-hand side), flux divergence terms, and source terms. Under steady, horizontally homogeneous conditions, both equations reduce to balances between vertical flux divergence and source density, so that (in the case of Eq. (7)) dF/dz = p</>, where z is height and F = (0, 0, F). Equation (9) (for 8) also includes familiar time derivative, advection, and source terms, but now the source term is nonzero only if ¿>P — 8 ^ 0, that is, if the scalar transferred across plant surfaces has an isotopic composition different from the air. The other significant term in Equation (9) is (F-VS)/(pc), the counterpart of the flux divergence term in Eqs. (7) and (8). It represents the contribution to local changes in 8 from the interaction of a scalar flux F with a gradient in <5 and is not a divergence. It balances the source term under steady, horizontally homogeneous conditions, so F(d5/dz) = (<5P — 8)p<f>.

2.3 Source Terms for C02

At this point it is useful to identify the isotopic source terms more precisely. If the scalar is C02, the source density for C02 mole fraction (</>) has an assimilation component ((/>A) and a respiration component (<f>R), so <f> = (f>A + </>R. The sign convention that is positive into V (that is, into the air) is retained, so t/>A < 0 and c/>K > 0. The source terms in Eqs. (8) and (9), for cS and 8, then become

4>(8P - 8) = 4>A(8P ~8) + <bR(8R - 8), where SA is the isotopic composition of the CO, assimilated into plants by current photosynthesis, and <5K is the isotopic composition of respired C02. The iso-source 8P(f> (the source density for the iso-concentration cS) is therefore the sum of contributions from assimilation and respiration, each the product of a C02 flux and the isotopic composition of the CO, exchanged (transported across the bounding surface SF) by the flux.

For the assimilation component, SA can be quantified in terms of the discrimination A or A*, defined by

^product i riv la 111 ^product

Krcactanl 1 + Srcaclanl where A is the definition of Farquhar and Richards (1984), convenient for calculating discrimination by a series of sequential processes, and A* is an alternative definition, convenient for processes acting in parallel (Farquhar et id., 1989a, Appendix Part III). The two arc related by A* = A/( 1 + A). For the present case, the product is the carbon fixed by current photosynthesis (composition 8a) and the rcactant is the CO, in the air at the point x (composition 8), so that

Neglecting second-order (A8, A2) and higher-order terms, the last of Eq. (12) simplifies to 8A~ 8 — A*. This approximate expression, or the exact Eq. (12), specifies SA in terms of A* and the air isotopic composition 8. The source terms in Eq. (10) now become

Sp4> = 4>a(8 ~ A* — SA*) + 4>r8r = cf>A(8 - a*) + </>rSr

0A(Sp - 8) = - 0A(A* + SA*) + 0R(SR - 8) = -</>aA* + <£r(Sr-S)

The discrimination A* (rather than A) is appropriate because surface elements or sources act in parallel rather than in series to influence ambient atmospheric isotopic composition.

It remains to specify SR and A*. For the respiration component, SR can be taken as the isotopic composition of the plant material, provided that discrimination during respiration can be ignored (Lloyd et al, 1996; Lin and Ehleringer, 1997), and discrimination on translocation of carbon within the plant (such as export of assimilated carbon from leaves) can also be ignored. For the assimilation component, net discrimination during photosynthesis depends on whether the photosynthetic pathway is C3 or C4 (Lloyd and Farquhar, 1994) and involves different mechanisms for discrimination against 13C in CO, (Farquhar et al, 1989a, b) and against ,sO in C02 (Farquhar et al, 1993). For C, photosynthetic discrimination against 13C, A is around about 20%o, decreasing with increasing water use efficiency (Farquhar and Richards, 1984; Farquhar et al, 1989b), and A* is about 0.4%o smaller. For C4 photosynthesis, discrimination against l3C is much lower with typical A values around 4-5%o.

2.4 Single-Point Lagrangian Equations

To identify the Lagrangian or fluid-following forms of the balance equations, we consider an arbitrary scalar entity with concentration a(x, f) and source density <p(x, t, a) which may depend on the scalar concentration a. Here a can stand for any of c, cR, or c8 (but not R or 8, for reasons given shortly). Table 2 gives the source density ip for each choice of a. Imagine an ensemble of realizations of the turbulent flow so that the instantaneous concentration a'", the instantaneous turbulent velocity uw, and the instantaneous source density ip"' = ip(x, t, a") are variables which differ randomly among realizations (distinguished by a superscript a>). The relationship between a'", uw, and <p'" is given by the instantaneous form of Eq (7) in which there is no Reynolds decomposition to produce eddy fluxes, so that the flux F'" accounts for scalar transfer by molecular diffusion only. In the body of the fluid, molecular fluxes are negligible relative to fluxes due to fluid motion and can safely be ignored (the high Peclet number approximation), while at source surfaces (SF), molecular fluxes are described by the source density. Under these conditions, the instantaneous equation (7) is

TABLE 2 Choices of the Arbitrary Concentration a and Source Density tp for the Scalars c, cR, cS, R, and 5

Single-Point Lagrangian

Kulerian Kquation Balance of a ¡p Methods Applicable?

TABLE 2 Choices of the Arbitrary Concentration a and Source Density tp for the Scalars c, cR, cS, R, and 5

Single-Point Lagrangian

Kulerian Kquation Balance of a ¡p Methods Applicable?


Scalar mole fraction




(8 i with S -

- R

Mole fraction of minor isotope









(9) with 5 -

-> R

Isotopic ratio





Isotopic composition



The last column indicates the applicability of the Lagrangian equation (16) and its canopy version, Eq. (17).

The last column indicates the applicability of the Lagrangian equation (16) and its canopy version, Eq. (17).

where D/Dt denotes the material derivative, or derivative following the motion.

The average scalar concentration field is found in the Lagrangian approach by tracking the motion of "marked fluid particles", connected parcels of fluid containing many molecules but small enough to be regarded as single points from the standpoint of the continuum turbulent flow. Hquation (14) can be integrated along the wandering path X"'(r) of a single fluid particlc to give a"'(f,) - a"'(t0

which says that the concentration in the fluid particle changes in response to the source densities it encounters along its path. A fluid passing through the space-time point (x,„f0) is therefore "marked" with a in proportion to the source density <p(x0, f0), and then contributes to the concentration a(x, t) at other points (x, t) according to the transition probability P(x, i|xo, t0), the probability that fluid motion transports the particle from (x0, t0) to (x, t). By considering a large number of fluid particles and taking an ensemble average, it can be shown (see Appendix B) that the mean concentration a(x, t) is given by a(x, t) =

The quantities a, <f>, and P are all ensemble-averaged, an operation which is ihe same as the more familiar lime average in a stationary (statistically steady) flow, but not in a nonstationary flow.

The transition probability P(x, f|xn, f0) carries all information about the velocity field that is needed to deduce a(x, t). In practice there are three main ways of determining P: first, one may obtain an exact or approximate analytic solution of the stochastic differential equations for the velocities of an ensemble of marked fluid particles. This approach was the basis of the Taylor (1921) solution for scalar dispersion in homogeneous turbulence and has been applied in vegetation canopies (Raupach, 1989a). A second approach is "random-flight" simulation, in which the stochastic differential equations for the marked particle velocities are solved by numerically constructing an ensemble of particle trajectories (Thomson, 1984, 1987; Sawford and Guest, 1987; Baldocchi, 1992; Katul et al, 1997). Third, one may obtain P from an Eulerian model for the velocity field, using higher-order closure or other methods. Examples of higher-order closure models in vegetation canopies include Wilson and Shaw (1977), Wilson (1988), Katul and Albertson (1998), Ayotte et al. (1999), and Massman and Weil (1999).

The arbitrary scalar entity a can be any conserved entity satisfying the superposition principle (that if source densities ip, and <p2 produce concentration fields a, and a,, then the source density <p, + <f>2 produces the concentration field a] + a2). This is true for the scalars c, cR, and c8, but not for R (isotopic ratio) or S (isotopic composition), because these scalars are ratios of the concentrations of the minor and major isotopes and so have a nonlinear dependence on sources. Hence, Eq. (16) cannot be applied to R or S (see Table 2).

The source density <p can depend on the concentration a at the point (x, t), provided that ip is statistically independent of the wind field (see Appendix B). Such a dependence, so that ip = ip(x, t, a), is indicated in Eq. (16). In this case Eq. (16) becomes an integral equation in a which must be solved by recursive or other means, rather than an explicit solution for a. This issue does not arise when the source density is specified independently, for instance, by the locations and strengths of point sources of air pollutants. However, it is crucial in most biogeochemical applications, because the source or sink densities for entities such as heat, water vapor, C02, and aerosols depend on the ambient concentrations of those entities. The inclusion of this dependence in Eq. (16) is considered in Appendix B, and the practical implications are further discussed in Section 3.5.

3. Inverse Methods for Inferring Scalar Sources and Sinks in Canopies

3.1 General Principles

Under steady conditions in a uniform, horizontally homogeneous vegetation canopy, the scalar source density <p(z) and concentration c(z) are functions only of height z. Equation (16) (with a = c and (p = <p, from Table 2) can then be written in the discrete form (Raupach, 1989b), ill ; i where ch; is the source density in canopy layer j, Azj is the thickness of layer j, m is the total number of source layers in the canopy, c-, is the concentration at height Z\, cr is the concentration at a reference height z, above the canopy, and Dti is the dispersion matrix, with n rows (i = 1 to n) corresponding to concentration measurement heights, and m columns (j = 1 to m) corresponding to source layers. The lowest source layer includes the ground, so that (/>,Az, is the sum of the scalar fluxes from the ground and the lowest canopy layer. Once I),, is known, inversion of Hq. (17) provides a solution to the inverse problem of calculating sources <f>\ from measured concentrations q. The linearity of Hq. (17) means that this solution is unique.

The dispersion matrix D-^ is a discrete form of the transition probability P in Eq. (16) and thus carries all the required information about the velocity field. Its elements have the dimension of aerodynamic resistance (s m '). The elements of column j of Dit are found by considering c/>|Azj to be a steady unit source, with sources in all other canopy layers set to zero. A theory of turbulent dispersion is used to calculate the concentration field c(z) resulting from this source distribution. The elements of column j of Dit are then given by

The heights of the sources (Zy) and the concentrations (z;) need not be the same in general.

Canopy-scale atmospheric inverse methods are all based explicitly or implicitly on Eq. (17). There are three main methods in this class, distinguished by the turbulent dispersion theory used to calculate D^ and aligning with the three means outlined in Section 2.4 for obtaining the transition probability P. Because all three rely on Eq. (17) and its Lagrangian foundation, the label "inverse Lagrangian" is appropriate in all three cases. The methods are (1) approximate analytic solution of Lagrangian equations for the velocities of an ensemble of marked fluid particles, leading to "localized near field" (LNF) theory (see below) and thence to an analytic specification for D,f, (2) random-flight numerical solution of stochastic differential equations for marked-

particle velocities, leading to numerical specification of DV] (Bal-docchi, 1992; Katul et at, 1997); and (3) definition of D,, by the solution of an Eulerian model for the velocity field and scalar transfer (Katul and Albertson, 1998; Massman and Weil, 1999; Katul et al., 2000). Specific calculations in the following use method (1) to find DYy

3.2 Localized Near Field Theory

LNF theory (Raupach, 1989a, b) is a semi-Lagrangian theory which provides an approximate means of calculating the concentration profile c(z) from a given source density profile 4>(z), recognizing the large-scale, coherent nature of turbulent eddies in vegetation canopies. Because of dominant role of these eddies, which have length scales on the same order as the strong shear layer at the top of the canopy (Raupach et al., 1996), the vertical turbulent transfer of scalars in the canopy cannot be described by gradientdiffusion theory. A clear indication of the failure of gradient-diffusion in canopies is provided by observations of counter-gradient or zero-gradient vertical fluxes of heat, water vapor, and CO, in forest canopies (Denmead and Bradley, 1987).

The theory centers on the evaluation of the transition probability P(x, i|x0, f0) in Eq. (16). The basic idea is that a cloud of marked fluid particles dispersing from an instantaneous point source at (x0, f0) can be regarded as undergoing random motions (in the sense that particle position is a Markov process) in the "far field" when the travel time (t — f0) is large compared with the Lagrangian time scale of the turbulence (TL). By contrast, in the "near field" when t — t0 is comparable with or smaller than TL, the dispersion of the cloud is governed by the persistence of the Lagrangian velocities of marked fluid particles close to the source (Taylor, 1921). Therefore, the dispersion of the cloud of marked particles (or of the scalar with which the particles are marked) is described by gradient-diffusion theory in the far field but not in the near field. Scalar sources in a vegetation canopy are spread throughout the canopy volume, so concentrations at any point x are made up of superposed contributions from sources at all travel times, including both far-field and near-field. Contributions from the latter are responsible for the observed failure of gradient-diffusion theory, including counter-gradient fluxes. A formal theory can be developed by splitting the transition probability P(x, f|xg, t0) into a diffusive "far-field" part PF and a nondiffusive "near-field" part Pn, so that P — Pp -F P^. The diffusive part Pp is described by a diffusion equation at all travel times. Two postulates are then made: first, that P PF at large travel times t — t0, and second, that PN = P — PF can be described by its value in locally homogeneous turbulence. The motivation for the second postulate is that for the small travel times for which PN is significant (t — t0 comparable with or less than TL) the particles are still quite close to the source and are in a turbulence field approximately the same as that at the source.

Corresponding to the partition P = PF + PN of the transition probability, the concentration c can be broken into two parts: c = cF + cN. Since the far-field part cF obeys a gradient-diffusion relationship between flux and concentration, it satisfies that

F(z) = - pKr(z)-^Kr(z) = oi(z)T,(z), (19) az where <xw(z) is the standard deviation, T] (z) is the Lagrangian time scale of the vertical velocity, Kf is the far-field eddy diffusiv-itv, F is the vertical eddy flux of the scalar, and z is height above the ground surface (z = 0). The expression for K¥ is the result of Lagrangian and Eulerian analyses in homogeneous turbulence (Taylor, 1921; Batchelor, 1949). The scalar flux F(z) is related to the source density by Eq. (7), which simplifies under steady, horizontally homogeneous conditions to dF_

where all elements are functions only of <rw(z) and Tl(z). These elements can be determined as follows: Let concentration measurement heights be Zj, so ct = c(Zj). For source layers, let z-} be the top of layer j and Z-} = (zj_, + Zj)/2 the center of layer j, with j = 0 at the ground so z0 = Z0 = 0. The far-field part of D^ is given by Eqs. (18) to (20):

where F0 is the flux at the ground surface. The flux at the top of the canopy is Fh = F(/z), and a concentration scale can be defined as c- = FJu..

The near-field part cN of the concentration accounts for the nondiffusive character of near-field dispersion. It is given by max(z,, Z,)

The near-field part needs to be calculated carefully because of the logarithmic singularity of /cN at £ = 0. If z-t is not between zM and Zj (that is, not within source layer j), then the singularity is not a problem. By replacing q — cx with cN(z|) — cN(zr) in Eq. (18), using Eq. (21) to find cN, and evaluating the integrals with rectangular approximations, D^'is found to be

0"wj TLi

where kN is the "near-field kernel," a weighted version of PN which can be calculated from the theory of dispersion in homogeneous turbulence. A good approximation (Raupach, 1989a) is

M£) = C]ln(l ~ e lfl) + c2e c, = - 1/^277 = -0.39894,

These equations provide a complete solution to the forward problem of calculating c(z) from a given <f>(z) in a uniform canopy. The required turbulence properties are canopy profiles of the standard deviation <7w(z) and the Lagrangian time scale TL(z) for vertical velocity.

3.3 The Dispersion Matrix

The specification of DVj is given here in more detail than in previous descriptions. From Eq. (18), D,j may be split into far-field term and near-field terms in the same way as c\

where crwj = crlv(Zj) and TLj = Tl(Zj). However, if z; is between Zj_, and Zj (within source layer j), then the first term in this expression must be replaced. We define /N(£) as the integral of fcN(£) from 0 to and evaluate it using a small-£ expansion (Raupach, 1989a, Eq. (A10)). This gives

[c, = c, = -0.39894, c4 = - 0.14127, c5 = 0.33333], (26) Then, Djj(,v) can be evaluated as Tu



Figure 1 illustrates the dispersion matrix by plotting the elements D-^, normalized as D^n, where u, is the friction velocity.

FIGURE 1 The normalised dispersion matrix D,j u,. Each profile is a column of the matrix, equal to the concentrations c-t - cr at heights z, (i = 1 to n) produced by a unit source in the layer centered at Z) (j = 1 to m). The number of concentration heights (») is 20, and the number of source layers (m) is 10. The reference height zr is 2h. Turbulence profiles aw(z)/u< and TL(z)uJh are from Eqs. (28) and (29), with csw = 1.6, «3(h) = 1.1, a3(l) = 1.25, cTL = 0.3, z,lh = 0.2, d/h = 0.75.

canopies, LJh is about 0.5. However, the ratio LJh shows more variation with canopy density than does LJ(h — d).

Several empirical forms have been proposed for <xw(z). A typical choice is

FIGURE 1 The normalised dispersion matrix D,j u,. Each profile is a column of the matrix, equal to the concentrations c-t - cr at heights z, (i = 1 to n) produced by a unit source in the layer centered at Z) (j = 1 to m). The number of concentration heights (») is 20, and the number of source layers (m) is 10. The reference height zr is 2h. Turbulence profiles aw(z)/u< and TL(z)uJh are from Eqs. (28) and (29), with csw = 1.6, «3(h) = 1.1, a3(l) = 1.25, cTL = 0.3, z,lh = 0.2, d/h = 0.75.

Each profile is the set of concentrations C; — cr (i = 1 to n) produced by a unit source in the layer centered at Zj (j = 1 to m). The near-field term D^ is responsible for the peak in each profile at the level of the unit source.

3.4 Turbulent Velocity Field

The turbulence properties o\v(z) and TL(z) can be specified with the aid of the observation that, for thermally neutral flow within and just above a vegetation canopy, the profiles of crw(z) and TL(z) are largely governed by a single velocity scale and a single length scale (Raupach, 1988; Kaimal and Finnigan, 1994; Raupach et al., 1996). Such a scaling is suggested by the mixing-layer analogy for flow within and just above the canopy (Raupach et al, 1996), which proposes that the turbulence structure in the strong shear layer near the top of the canopy is patterned on a plane mixing layer rather than on a boundary layer. This hypothesis implies that the velocity scale for the canopy turbulence is the friction velocity !/,, a measure of the turbulent momentum flux to the canopy, and that the length scale is the thickness of the mixing layer. This can be characterized by the length Ls = u(h)/u'(h) (where u is the mean flow velocity, h is the canopy height, and u'(z) = (du/dz), which is quite well determined by h — d (where d is the zero-plane displacement of the canopy). Data from a wide range of canopies show that LJ(h — d) is close to 2. Also, d/h is well constrained (Raupach, 1994),' being about 0.75 for typical canopies but depending slowly on the leaf area index and its profile with height in the canopy. Therefore Ls can also be related to h, and in typical

1 In Fig. 1 of Raupach (1994), the horizontal axis is incorrectly labeled as A

(Leaf Area Index). It should be A (Frontal Area Index). The relationship assumed in that paper between the two is A = 2A.

where csw is an empirical constant (typically around 1.5), ai|h| (typically about 1.1) is the ratio o\v/u, at the top of the canopy, n3(il (typically about 1.25) is the ratio crw/u, in the inertial sublayer above the roughness sublayer, and zwt-f is the height of the roughness sublayer (the layer just above the canopy within which the mixing-layer analogy applies). This profile assumes that crju, takes a slightly lower value at the top of the canopy than in the inertial sublayer well above the canopy, as implied by the mixing-layer analogy and observed in practice (Raupach et al., 1996). Within the canopy, erw/u, is assumed to have an exponential form.

For Tl(z), the principle of the parameterization is that u,TJLs is constant with height within the canopy (except very close to the ground) and in the roughness sublayer. This is a consequence of the mixing-layer analogy. Constancy of u.TJLs leads to two alternative parameterizations u,T, u.T{ h - d

depending on whether h or h — d is used as a surrogate for the shear length scale Ls. Appropriate values for the constants are cTL = 0.3 and c'JL = 1.2. As argued above, the mixing-layer analogy implies that the second form (involving (h — d)) is physically preferable, but it carries the penalty that d as well as h must be determined or estimated. Equation (29) is applicable in the range z, < z < zrutt, where z, is a low level in the canopy (typically about 0.25h) below which TL decreases with proximity to the ground. Above z,.lltt, Tl is given by standard inertial-sublayer expressions and increases linearly with z — d under thermally neutral conditions.

The height of the roughness sublayer, zrutt, can be estimated by noting that the far-field eddy diffusivity K¥ is given everywhere by ctw27'l (Eq. (19)). In particular, within the roughness sublayer (z< zrutt), Eq. (29) implies that K¥ = cnha^2/u.. Also, K¥ in the neutral inertial sublayer (z S' zrutt) is given by Kr = ku.(z — d), where k is the von Karman constant (assumed to be 0.4). Equating these two estimates for Ke at zrutf, it is seen that zrutt = d + ((¡:i(i)2/x)cTlh, where «3|I, = crw/u- in the inertial sublayer. Typically, zrutt is about 2h.

Recent works by Leuning et al. (2000) and Leuning (2000) have led to improvements in these parameterizations in two respects. First, smoothing of the slope discontinuities in crw(z) and TL(z) (implied by Eqs. (28) and (29)) produces more stable behavior in source distributions inferred from concentration profiles by inversion of Eq. (17). Second, Leuning (2000) introduced thermal stability into the parameterizations of trw(z) and TL(z) by assuming that the stability parameter throughout the roughness sublayer is h/LMQ (where LMO is the Monin-Obhukov length) and that the effect of this on TL below zrurt can be described by the factor [KH(h/LMO)aJ(Q)}/{KH (0)(V(WLmo)], where KH is a scalar eddy diffusivity. Standard, inertial-sublayer forms were used to describe the stability dependencies of KH and crw2.

3.5 Solutions for Forward, Inverse, and Implicit Problems

Having specified the turbulence properties <xw and TL and thence the dispersion matrix D;,, the apparatus is now in place to use Eq. (17) to solve three generic kinds of problem: the forward problem of determining the scalar concentration profile c(z) from a specified source density profile <f)(z), the inverse problem of determining 4>(z) from specified or measured information about c(z), and the implicit or coupled problem of determining both c(z) and 4>(z) together when cf> is a given function of c.

3.5.1 The Forward Problem

Of the three kinds of problem, this is the easiest theoretically but its main practical use is as a stepping stone in the solution of the other two problems. The solution in discrete form is given directly by Eq. (17). Sample results are shown in Figure 2, which shows the concentration profiles calculated from Eq. (17) (with JD;- specified using Eqs. (24) to (27)) for three assumed source profiles. A sharp, strongly peaked elevated source profile produces a peak in the concentration profile and hence a countergradient flux in the layer just below the peak where the flux is upward but dc/dz is positive, because of the strongly localized near-field contribution to the concentration field. By contrast, a more uniform source profile, in cluding significant contributions from low in the canopy, produces a concentration profile which decreases uniformly with height.

3.5.2 The Inverse Problem

This is the primary means of obtaining information about the canopy source distribution of a scalar from atmospheric concentration measurements. A formal discrete solution is found by matrix inversion of Eq. (17), choosing the number of source layers (m) to be equal to the number of concentration measurements (n) so that D|| is a square matrix. However, this solution provides no redundancy in concentration information, and therefore no possibility for smoothing measurement errors in the concentration profile, which can cause large errors in the inferred source profile. A simple means of overcoming this problem is to include redundant concentration information, and then find the sources 4>\ which produce the best fit to the measured concentrations q by maximum-likelihood estimation. By minimizing the squared error between measured values and concentrations predicted by Eq. (17), is found (Raupach, 1989b) to be the solution of m linear equations

This solution uses a nonsquare dispersion matrix in which n > m, to obtain the necessary redundancy in concentration information. The solution is valid whether or not <p is dependent on c, because it determines the values of cf> consistent with the current c field. The result is the source density profile <f)f in a (small) number of canopy layers, with the lowest layer including the ground source (F0). The total flux from the canopy is also obtained as the sum of </>,Azj over all layers.

Figure 3 shows a test of the sensitivity of this method to two key factors, the size of D¡j (determined by the number of concentration measurements, n, and source layers, m) and the presence or absence of the near-field term in in the dispersion matrix D^. Figure 3a shows the assumed concentration field, based on forward calculation of the concentration profile q (i = 1 to n) from a specified source profile using Eq. (17), with D^ as in Figure 1, and with n = 20, m = 10. This concentration field is then used to reconstruct the source profile <b{z) with Eq. (30), and the flux F(z) from Eq. (20), under three scenarios: (1) (ti,m) = (20,10), with D,/lV) included; (2) (n,m) = (10,5), with D,/N) included; and (3) (n,m) = (10,5), with -Dii(,V) omitted. Figures 3b and 3c respectively show the inferred profiles </>, and F-f at layers centered on heights

, , , , ... Z:. Scenario (1) represents optimum information, with (n, m)

(a) Assumed normalized source profiles <p(z)«/Fh and (b) pre- r .

. , ... t .. a i \i / u r i \ -.identical with that used in the forward calculation of c. Not sur-

dicted normalized concentration profiles (c - cr)/c. (where c. = Fhlu.), if-

lustrating a typical solution to the forward problem. Turbulence profM™111^'lhe inverse method exactly recovers the profile of </>, ini-(tw(z)/u> and TL(z)u,lh are from Eqs. (28) and (29), with csw = 1.6, am =tiallY assumed for the forward calculation, shown as a heavy line 1.1, am = 1.25, cTL = 0.3, z,/h = 0.2, d/h = 0.75. The dispersion matrix ftp Figures 3b and 3c. Scenarios (2) and (3) attempt to recover the


as in Figure 1.

initially assumed source profile with progressively more degraded

\ a



—•— (n,m) = (20,10)

—(n,m) = (10,5)

(n,m) = (10,5)

no near field


Source density

FIGURE 3 Sensitivity of the inverse Lagrangian (IL) method to the number of concentration measurements (n) and source points (m), and to the presence or absence of the near-field term in the dispersion matrix. The dispersion matrix is as in Figure 1. (a) Assumed normalized concentration field (c - c,.)/c,; (b) normalized source profile 0(z)/i/Fh; (c) normalized flux profile F(z)/Fh. Scenarios are: (1) (», m) = (20, 10), with D|j(M included; (2) (n, m) = (10,5), with Djm included; (3) (n, m) = (10, 5), with Dii(N) omitted.

I a








0 1 2 Source density 4>h/Fh

0 1 2 Source density 4>h/Fh

FIGURE 4 Sensitivity of the inverse Lagrangian (IL) method to errors in concentration measurements. To each c, is added a random error, uniformly distributed over the range [— ac c., ac c,]. Scenarios are: (1) (», m) = (20, 10) and ac = 0; (2) (», m) = (10, 5) and ac = 0.3; (3) (n, m) = (10, 5) and ac = 3. Other details are as for Figure 3.

information about the concentration field and the dispersion matrix. In scenario (2), alternate concentration points are omitted and the number of source layers is also reduced to achieve redundancy, so that (n, m) = (10, 5). Scenario (3) uses (n, m) = (10, 5) as in scenario (2), but also uses a degraded dispersion matrix in which the near-field component A,('v> is omitted, so that the inversion is done by solving Eq. (30) with D^ = D¡jiF). Figures 3b and 3c show the effects of these losses in information. Reduction of {n, m) from (20, 10) to (10, 5) in scenario (2) still produces a reasonable approximation to the initially assumed ("correct") source and flux profiles, but simultaneously omitting the near-field term in the dispersion matrix, in scenario (3), produces clearly unphysical source and flux profiles with large overshoots on both profiles.

Figure 4 shows a further test of the effect of imperfect information on the inferred source profiles. This time, random noise is added to the concentration profile before carrying out the inversion with Eq. (30), to test the ability of the method to cope with random errors in concentration measurement. Again, three scenarios are used: the first is identical with scenario (1) above. The other two use (n, m) = (10, 5) and the full D^ as in scenario (2), but also add to each C\ a random error, uniformly distributed over the interval [— crcc., <rcc.], where crc is either 0.3 or 3 (and c- = FJu.). The smaller error produces little effect on the inferred profiles and Fj, but the larger one leads to significant degradation in both profiles.

These illustrative calculations have been carried out with a very simple inversion procedure based on Eq. (30). Rapid recent developments in the application inverse theory offer several possibilities for improvement. First, the use of singular-value decomposition as a formal inversion framework (e. g., Press et al, 1992)

provides means for assessing the effect of measurement errors and guarding against degeneracies in the solution caused by the inability of the measurements to distinguish between alternative different solutions to the inverse problem. Second, Bayesian synthesis approaches (Tarantola, 1987) offer means for making the inversion process more robust through the application of prior constraints. Third, there are possibilities for further improving the robustness of the method by fitting for parameters in an empirical or process-based model for the source profile, rather than for the layer-discretized source-sink profile directly.

3.5.3 The Implicit or Coupled Problem

This problem is important in modeling the canopy concentration (c) and source-sink (<f>) profiles of entities for which <fi depends on c, including most biogeochemically significant entities exchanged between vegetation and the atmosphere: heat, water vapor, aerosols, C02, other trace gases which react with stomatal or other canopy surfaces, and isotopes (to be discussed in more detail in the next section). In these cases both <f)(z) and c(z) are unknowns, and the prescribed constraints in the problem are the concentration cr at a reference height above the canopy and forcing or surface parameters determining the relationship 0(c). The biogeochemical importance of such situations is the reason for careful consideration of concentration-dependent sources in the Lagrangian analysis (Section 2.4 and Appendix B). An Eulerian version of an analogous problem for diffusive flow in the atmospheric surface layer is defined by the scalar conservation equation (7) and mixed or radiation boundary conditions (Philip, 1959). A simple solution for the implicit problem in a canopy can be obtained when the relationship (f>(c) is linear (or can be linearized) so that <p = q — rc, where q is a specified concentration-independent scalar source density and r is a rate constant describing a sink of scalar proportional to c (that is, governed by first-order reaction kinetics). This describes a concentration-independent source when r = 0, and a pure first-order sink (</>= — rc) when q = 0. Putting this form for 4> into Eq. (17), we obtain

It is necessary to make D^ square (n = m) so that c is available at all heights for which sources must be evaluated. Then Eq. (31) is a set of m coupled linear equations in q — cr.

This approach has been used in several multilayer canopy models for coupled heat and water vapor transfer (Dolman and Wallace, 1991; van den Hurk and McNaughton, 1995; McNaughton and van den Hurk, 1995). The solution in this case can be obtained with a linear

0 0

Post a comment