Investigating a Hierarchy of Eulerian Closure 
Models for Scalar Transfer Inside Forested 
Canopies
Authors: Jehn-Yih Juang, Gabriel G. Katul, Mario B. 
Siqueira, Paul C. Stoy, and Heather R. McCarthy
The final publication is available at Springer via https://dx.doi.org/10.1007/s10546-008-9273-2.
Juang, Jehn-Yih, Gabriel G. Katul, Mario B. Siqueira, Paul C. Stoy, and Heather R. McCarthy. 
“Investigating a Hierarchy of Eulerian Closure Models for Scalar Transfer Inside Forested 
Canopies.” Boundary-Layer Meteorology 128, no. 1 (April 25, 2008): 1–32. doi:10.1007/
s10546-008-9273-2.
Made available through Montana State University’s ScholarWorks 
scholarworks.montana.edu 
Investigating a Hierarchy of Eulerian Closure Models
for Scalar Transfer Inside Forested Canopies
Jehn-Yih Juang · Gabriel G. Katul ·
Mario B. Siqueira · Paul C. Stoy ·
Heather R. McCarthy
Received: 21 February 2007 / Accepted: 27 March 2008 / Published online: 25 April 2008
© Springer Science+Business Media B.V. 2008
Abstract Modelling the transfer of heat, water vapour, and CO2 between the biosphere
and the atmosphere is made difficult by the complex two-way interaction between leaves and
their immediate microclimate. When simulating scalar sources and sinks inside canopies on
seasonal, inter-annual, or forest development time scales, the so-called well-mixed assump-
tion (WMA) of mean concentration (i.e. vertically constant inside the canopy but dynamically
evolving in time) is often employed. The WMA eliminates the need to model how vegeta-
tion alters its immediate microclimate, which necessitates formulations that utilize turbulent
transport theories. Here, two inter-related questions pertinent to the WMA for modelling
scalar sources, sinks, and fluxes at seasonal to inter-annual time scales are explored: (1) if
the WMA is to be replaced so as to resolve this two-way interaction, how detailed must
the turbulent transport model be? And (2) what are the added predictive skills gained by
resolving the two-way interaction vis-à-vis other uncertainties such as seasonal variations in
physiological parameters. These two questions are addressed by simulating multi-year mean
J.-Y. Juang · G. G. Katul · M. B. Siqueira · P. C. Stoy · H. R. McCarthy
Nicholas School of the Environment and Earth Sciences, Duke University, Durham, NC, USA
J.-Y. Juang (B)
Department of Geography, National Taiwan University, No. 1, Sec. 4, Roosevelt Road, Taipei 10617,
Taiwan
e-mail: jjuang@ntu.edu.tw
G. G. Katul
Department of Civil and Environmental Engineering, Duke University, Durham, NC, USA
Present Address:
M. B. Siqueira
Departamento de Engenharia Mecânica, Universidade de Brasília, Brasilia, Brazil
Present Address:
P. C. Stoy
School of GeoSciences, Department of Atmospheric and Environmental Science, University of Edinburgh,
Edinburgh, UK
Present Address:
H. R. McCarthy
Department of Earth System Science, University of California at Irvine, Irvine, CA, USA
scalar concentration and eddy-covariance scalar flux measurements collected in a Loblolly
pine (P. taeda L.) plantation near Durham, North Carolina, U.S.A. using turbulent transport
models ranging from K-theory (or first-order closure) to third-order closure schemes. The
multi-layer model calculations with these closure schemes were contrasted with model calcu-
lations employing the WMA. These comparisons suggested that (i) among the three scalars,
sensible heat flux predictions are most biased with respect to eddy-covariance measurements
when using the WMA, (ii) first-order closure schemes are sufficient to reproduce the seasonal
to inter-annual variations in scalar fluxes provided the canonical length scale of turbulence is
properly specified, (iii) second-order closure models best agree with measured mean scalar
concentration (and temperature) profiles inside the canopy as well as scalar fluxes above
the canopy, (iv) there are no clear gains in predictive skills when using third-order closure
schemes over their second-order closure counterparts. At inter-annual time scales, biases in
modelled scalar fluxes incurred by using the WMA exceed those incurred when correcting
for the seasonal amplitude in the maximum carboxylation capacity (Vcmax,25) provided its
mean value is unbiased. The role of local thermal stratification inside the canopy and possible
computational simplifications in decoupling scalar transfer from the generation of the flow
statistics are also discussed.
Keywords Canopy turbulence · Closure models · Net ecosystem CO2 exchange ·
Scalar transfer · Stratified canopy flows
“The tree, tilting its leaves to capture bullets of light; inhaling,
exhaling; its many thousand stomata breathing, creating the air”.
Ruth Stone, 2002, In the Next Galaxy
1 Introduction
Quantifying the exchange of scalars (e.g., carbon dioxide CO2, water vapour H2O, tem-
perature, or other chemical species) between leaves and their local environment (hereaf-
ter referred to as microenvironment) is made difficult by a two-way interaction in which 
the microenvironment exerts controls over scalar exchange at the leaf surface, and leaves 
have some capacity to regulate their own microenvironment through stomatal opening and 
closure. This two-way interaction is further complicated by the vertical distribution of 
foliage within the canopy, leading to significant vertical gradients in the radiation envi-
ronment and airflow regimes. The intrinsic non-linearity in leaf physiological responses 
(e.g. leaf-level photosynthesis and transpiration) to radiation further exasperates this 
difficulty.
To date, most eco-physiological approaches to modelling annual and inter-annual 
ecosystem-scale carbon and water fluxes have focused on radiative transfer and the non-
linearity in the physiological response to incident radiation at the leaf surface (Aber et al. 
1996; De Pury and Farquhar 1997; Kirschbaum et al. 1998; Leuning et al. 1995; Luo et 
al. 2001; Naumburg et al. 2001; Nikolov et al. 1995; Wang and Leuning 1998; Williams 
et al. 1996, 1998). These models did not specifically resolve the vertical distribution of the 
scalar concentration (CO2 and H2O) and temperature, and assume that within the canopy 
volume, mean scalar concentration is constant and identical to its state above the canopy 
(hereafter, referred to as the well-mixed assumption, WMA). However, the time evolution 
of the mean scalar concentrations (and temperature) is generally included given the avail-
ability of such measurements above the canopy. The basic premise behind the WMA is that
temporal variation of the mean scalar concentration (e.g. diurnally or seasonally) far exceeds
their vertical variation inside the canopy. Besides, resolving such two-way interaction adds
significant computational burden and model complexity because of the need to compute the
flow field with perhaps modest gains in predictive skills, though the degree of improvement in
resolving sources, sinks, and fluxes remains largely unexplored. Furthermore, uncertainties in
describing the non-stationarity and vertical inhomogeneity in physiological parameters (e.g.
in photosynthesis calculations) may overshadow any improvements gained by resolving this
two-way interaction.
While the WMA may be defensible for some canopy types, it is too simplistic for tall-
forested ecosystems, especially when such an assumption is confronted with experimen-
tal evidence that vertical variations in excess of 50 ppm for CO2 concentration and 3 K or
more for air temperature occur inside the canopy volume even during daytime conditions
(Lai et al. 2002b; Siqueira and Katul 2002). Because the vertical variations in mean scalar
concentration profiles are not random within the canopy, the WMA may inject systematic
biases in modelling scalar sources, sinks, and fluxes. Hence, future developments in ecosys-
tem carbon–water source–sink, and flux modelling are likely to benefit from answering two
inter-related questions:
(1) If the WMA in mean scalar concentration is to be revised, then how detailed must the
turbulent transport model be to resolve the two-way interaction between the canopy
and its microclimate?
(2) When incorporating a detailed turbulent transport model, are the improvements in the
overall prediction of scalar sources, sinks, and fluxes sufficiently large when com-
pared to improvements gained by resolving the seasonality in physiological parameters
such as the ones most pertinent to leaf photosynthesis? The variability in physiological
parameters is rarely included in biogeochemical cycle models.
These two questions frame the study objectives here and are explored using the multi-year
record available at the AmeriFlux Duke Forest loblolly pine site as a case study.
Several multi-layer one-dimensional models have been developed to resolve the two-way
interaction between leaf and microclimate using turbulent transport theories in conjunction
with detailed physiological and radiative transfer principles (Baldocchi 1992; Baldocchi et
al. 1997; Baldocchi and Meyers 1998; Leuning et al. 1995; Meyers and Paw 1986, 1987;
Raupach 1988, 1989; Simon et al. 2005a,b). The term ‘CANVEG’ (for ‘canopy vegetation’)
was coined for such multi-layer models that employ time and planar-averaged equations
(Baldocchi 1992; Baldocchi et al. 1997; Baldocchi and Meyers 1998). CANVEG models
do not resolve explicitly rapid transient physiological responses occurring on turbulent time
scales (Gu et al. 1999; Su et al. 1996) or sun-fleck events (Naumburg et al. 2001). The net
effect of these phenomena on half-hourly time scales may be parameterized (e.g. by lags
or stochastic jumps) but are not considered here. In the original CANVEG (e.g. Baldocchi
1992), linkages between scalar sources and mean concentration relied on the principles of
Lagrangian fluid mechanics for characterizing turbulent dispersion (Raupach 1989), which
was a major theoretical improvement over classical Eulerian first-order flux-gradient closure
models (or K-theory). However, the Lagrangian framework in CANVEG suffered from two
fundamental limitations: (i) it cannot explicitly treat thermal stratification inside the canopy
volume, known to be significant in tall-forested ecosystems (Malhi et al. 1998; Siqueira and
Katul 2002), and (ii) it requires detailed higher-order velocity statistics profiles (e.g., the
vertical velocity standard deviation profile), which either have to be measured or simulated.
To circumvent these limitations, Lai et al. (2000a,b) and Siqueira et al. (2002) modified the
CANVEG approach in Baldocchi (1992) and Baldocchi and Meyers (1998) by employing a
Eulerian-Lagrangian hybrid framework for turbulent dispersion. The basic premise is to uti-
lize higher order Eulerian closure approaches (Meyers and Paw 1986; Wilson 1989; Wilson 
and Shaw 1977) to compute the velocity statistics within the canopy, and then use Lagrangian 
transport theory to couple scalar sources to their mean concentration, while retaining similar 
physiological and radiative transfer schemes as in CANVEG. The approach in Lai et al.
(2002a,b) eliminates one of the two limitations earlier mentioned, but thermal stratification 
inside the canopy was treated as a modifier to the upper boundary conditions of the modelled 
velocity statistics and on the Lagrangian integral time scale, an approach also employed by 
Leuning (2000).
Here, a multilayer biosphere–atmosphere model is developed that retains the detailed 
eco-physiological parameterization and radiative transfer principles in CANVEG, but the 
complexity in the turbulent transport scheme needed to capture this two-way interaction 
between the canopy and its microclimate is varied in a hierarchical manner. Different closure 
approximations ranging from first-order to third-order schemes are employed to parameterize 
higher order turbulent moments in the governing conservation equations for both momentum 
and scalar transfer. As a reference, these model calculations are contrasted with scalar sources 
(and subsequently fluxes) calculated by the WMA (i.e. a zeroth order closure in the model 
hierarchy). In addition, these model calculations are repeated with and without resolving the 
effects of atmospheric stability inside the canopy. Resolving atmospheric stability locally 
inside the canopy necessitates a computationally demanding iterative scheme for solving 
simultaneously the Reynolds stress budget and heat flux budget equations, which may not 
be feasible for inter-annual water and carbon flux calculations.
Because the interest here is in ‘ecological’ (i.e. seasonal to inter-annual) time scales, clo-
sure models with increasing closure complexity are compared against a multi-year record 
of eddy-covariance sensible and latent heat fluxes and CO2 flux collected in a uniform Lob-
lolly pine (Pinus taeda L.) stand in the Blackwood Division of Duke Forest, described next. 
Because the CANVEG model does not explicitly account for soil moisture stress (i.e. soil–
plant hydrodynamics are not considered), the data-model intercomparison is restricted to 
years with well-watered soil moisture states.
2 Experiment
2.1 Study Site
All observations were made at a pine stand in the Blackwood Division of the Duke Forest 
near Durham, North Carolina, U.S.A. (site location: 35◦58′N, 79◦05′W, 163 m above sea 
level) as part of a long-term H2O/CO2 flux monitoring initiative (Baldocchi et al. 2001). 
The study site is a uniformly planted Loblolly pine forest (planted in 1983 at 2 m× 2.4 m 
spacing) extending some 300–600 m in the east-west direction and 1,000 m in the south-
north direction. The long-term mean annual precipitation is 1,185±177 mm, and the annual 
mean air temperature is 14.9 ± 0.9◦C. The local topographic variations around the 
micrometeo-rological tower are small enough (slope < 5%) to ignore the effect of complex 
terrain on the flow statistics (Siqueira et al. 2002). In addition to the dominant pine 
overstorey, the sub-canopy (roughly ranging in height from 0 to 0.4 h, where  h is the 
mean height of the main canopy) contains some 40 woody species, of which Liquidambar 
styraciflua L., Acer rubrum L., Ulmus alata Michx., and  Cornus florida L. are the most 
prevalent (Palmroth et al. 2005).
Fig. 1 Top: Leaf area index (LAI, m2 m−2) variation from 2001 to 2005 for the pine canopy and the hardwood
understorey. The light grey area and the darker grey area respectively, represent the severe drought event during
the growing season and the ice storm event in December 2002. Bottom: the leaf area density (LAD, m2 m−3)
profiles at different times of the year based on different LAI values in 2005
While the forest is evergreen, leaf area index (LAI) variations can be significant, as
evidenced by the measurements in Fig. 1 (taken from a model-data analysis described by
McCarthy et al. (2007) in which multiple data sources were used). During the growing season
in 2002, the ecosystem experienced a severe drought lasting several months, and on December
4 and 5 2002, an ice storm occurred causing severe damage to the canopy structure (McCarthy
et al. 2006b). The lower panel of Fig. 1 shows the vertical distribution of the foliage at differ-
ent times of the year. These complex spatial and temporal variations in foliage distribution
will have appreciable effects on both the radiation/energy environment and the attenuation
of turbulent flow statistics inside the canopy volume, thereby significantly influencing the
two-way interaction.
2.2 Measurements
The velocity statistics (wind velocity and Reynolds stresses), sensible heat flux (HS), latent
heat flux (Lv E , where Lv is the latent heat of vaporization (J kg−1K−1)), and CO2 fluxes
were measured at 22.5 m using an eddy-covariance (hereafter referred to as EC) system com-
posed of a LI-7500 open-path CO2/H2O gas analyzer (Li-Cor Inc., Lincoln, NE, USA) and a
CSAT3 tri-axial sonic anemometer (Campbell Scientific Inc., Logan, UT, USA). The EC data
were sampled at 10 Hz and averaged every 30 min. The effects of air density fluctuations on
CO2 and H2O flux measurements were corrected after Webb et al. (1980), while corrections
to scalar statistics beyond the vertical turbulent fluxes are discussed elsewhere (Detto and
Katul 2007).
The net radiation (Rn) was measured using a Fritschen-type net radiometer at the tower
top (at 24 m) from 1997 until 2003, and then this radiometer was replaced with a Kipp &
Zonen CNR1 net radiometer (Kipp & Zonen USA Inc., Bohemia, NY, USA) in 2004. The
photosynthetically active radiation (PAR) was measured with a Li-Cor LI-190SA quantum
sensor at the canopy top, and mean air temperature and relative humidity were measured using
a Campbell Scientific HMP35C temperature/relative humidity probe at (2/3)h. To sample the
vertical distribution of scalar concentrations, a multi-port system was installed at 10 different
levels (0.1, 0.75, 1.5, 3.5, 5.5, 7.5, 9.5, 11.5, 13.5, and 15.5 m) and a Li-Cor 6262 CO2/H2O
infrared gas analyzer was used to measure the mean concentrations. Another profiling system
was employed to measure the mean air temperature at eight levels (1.5, 3.5, 5.5, 7.5, 9.5,
11.5, 13.5, and 15.5 m) using shielded copper-constantan thermocouple sensors.
The forest floor CO2 efflux was measured with an automated carbon efflux system (ACES,
US Patent 6692970) developed by the USDA Forest Service, Southern Research Station
Laboratory in Research Triangle Park, NC (Butnor et al. 2003; Butnor and Johnsen 2004;
Palmroth et al. 2005), whose configuration details, quality checks, and spatial sampling area
are described in Palmroth et al. (2005). These measurements were used to drive the lower
boundary conditions for the mean CO2 continuity equation (see Appendix A).
The plant area index (PAI, m2 m−2) was measured several times per year using a pair
of Li-Cor LAI 2000 optical sensors. The plant area density (hereafter PAD, m2 m−3) mea-
surements were conducted at 1-m intervals from the forest floor up to the canopy top (Lai
et al. 2000a; Pataki et al. 1998; Schafer et al. 2003). The dynamics of leaf area were recon-
structed for broadleaf species, using data on leaf litterfall mass and timing, specific leaf area
and allometry, and for Pinus taeda, needle litterfall (lagged by 2 years to account for foliage
longevity) and timing, combined with needle elongation rates, and fascicle and shoot counts
(McCarthy et al. 2006a). Further details about the site and data processing are described
elsewhere (Juang et al. 2006; Oren et al. 1998; Palmroth et al. 2005; Pataki and Oren 2003;
Stoy et al. 2005).
3 Theory
3.1 Governing Equation
For completeness, the models employed in the turbulent closure hierarchy are briefly reviewed.
For a steady-state, incompressible, high Reynolds number and Peclet number flow, the time-
and spatially-averaged budget equation for the mean longitudinal momentum, in the absence
of Coriolis effects and for the case of a planar-homogeneous canopy, can be expressed as
(Katul and Albertson 1998; Meyers and Paw 1987; Wilson and Shaw 1977)
∂ 〈u〉
∂t
= 0 = −Cd A(z) 〈u¯〉2 −
∂
〈
u′w′
〉
∂z
, (1)
where u and w are the instantaneous longitudinal and vertical velocity components, Cd is the
drag coefficient (Wilson and Shaw 1977), and A(z) is the bulk PAD at level z. Overbars and
brackets represent temporal and spatial averaging of a flow variable, and primed quantities
denote departures from this average. Using multiple sonic anemometer measurements, Katul
and Albertson (1998) determined Cd = 0.2 for this stand, a value used throughout our study.
Using similar simplifications and averaging procedures, the steady state budget equation
for an arbitrary scalar ξ is given by
∂
〈
ξ
〉
∂t
= 0 = Sξ −
∂
〈
w′ξ ′
〉
∂z
, (2)
where Sξ are local sources and sinks related to the vegetation elements.
Assuming Fickian diffusion (or its resistor analogy) for mass transfer from the leaf to the
surrounding atmosphere, the difference in scalar quantities between the leaves and the sur-
rounding air becomes the driver for Sξ . Upon scaling with the local leaf area density (LAD),
a(z), leads to a spatially-averaged source given by
Sζ = a(z)Gζ
(〈
ξL
〉 − 〈ξ 〉) , (3)
where Gζ (m s−1) is the bulk conductance for scalar ξ between the leaves and the surrounding
air, and ξL is the mean quantity of a specific scalar (ξ) inside the leaf, which may be one
of the three considered here: the leaf surface temperature (at the skin) 〈TL
〉
, the intercellu-
lar saturated water vapour concentration 〈qL〉, and the intercellular CO2 concentration
〈
CL
〉
.
For water vapour and CO2 concentrations, Gζ includes both stomatal conductance (gs,ξ ,
mol m−2 s−1) and boundary-layer conductance (gb,ξ , mol m−2 s−1), while only boundary-
layer conductance contributes to heat transfer. Note that PAD is responsible for momentum
extraction and light interception, but only LAD is responsible for the simultaneous water
vapour release andCO2 uptake by the leaves, though the two quantities are considered to be
approximately the same for this evergreen stand.
Using a flat-plate boundary-layer analogy at the leaf surface, the boundary-layer conduc-
tance is given by (Campbell and Norman 1998):
gb,ξ = υξ
( 〈u〉
0.7d
)0.5
, (4)
where υξ is a molecular transfer coefficient for scalar ξ , and d is the characteristic length
of the leaf. Here, d = 0.0015 m for the pine foliage and d = 0.15 m for the understorey
broadleaf foliage.
3.2 Closure Approximation
In Eqs. 1 and 2, the turbulent momentum
〈
u′w′
〉
and scalar
〈
w′ξ ′
〉
fluxes are additional
unknowns that require parameterization so that sources, sinks, fluxes, and concentration can
be solved. In the hierarchy of models, either parameterizations or full budget equations for
these higher-order terms are derived.
3.2.1 First-order Closure Model
The first-order closure approximation (or K-theory) simply parameterizes the turbulent fluxes
in Eqs. 1 and 2 into
〈
u′w′
〉
= −Km ∂ 〈u〉
∂z
, (5a)
〈
w′ξ ′
〉
= −Kξ ∂
〈
ξ
〉
∂z
, (5b)
where K is the eddy diffusivity, parameterized based on the Prandtl–von Karman mixing
length hypothesis as
K = L(z)2
∣∣∣∣
∂ 〈u〉
∂z
∣∣∣∣ , (6)
where L(z) is a canonical mixing length at height z. Katul et al. (2004) recently summarized
previous work on specifying time or length scales inside canopies (Liu et al. 1996; Massman
and Weil 1999; Poggi et al. 2004) and concluded that L(z) may be assumed approximately
constant inside the canopy volume for a dense canopy and can be estimated from:
L(z) =
{
αml h; z < h
k(z − d0) z ≥ h, (7)
where k = 0.4 is the von Karman constant, and d0 is the zero-plane displacement height
determined from the centroid of the drag force (Jackson 1981; Thom 1971). The parameter
αml is calculated by assuming that L(z) is continuous at the canopy top (Katul et al. 2004).
3.2.2 Second-order Closure Model
Instead of using K-theory to characterize turbulent fluxes in Eqs. 1 and 2, budget equations
can be derived for the turbulent flux though triple correlation terms must be subsequently
parameterized. For the momentum and Reynolds stress equations, standard closure approxi-
mations result in (Wilson and Shaw 1977; Katul and Albertson 1998; Katul and Chang 1999)
∂
〈
u′w′
〉
∂t
= 0 = −
〈
w′w′
〉 ∂ 〈u〉
∂z
−
Q
〈
u′w′
〉
3λ2
+ Cw Q2 ∂ 〈u〉
∂z
− ∂
∂z
〈
w′u′w′
〉
, (8a)
∂
〈
u′u′
〉
∂t
= 0 = −2
⎛
⎝
〈
u′w′
〉 ∂ 〈u〉
∂z
+ 〈u〉
∂
〈
u′w′
〉
∂z
⎞
⎠ − Q
3λ2
(〈
u′u′
〉
− Q
2
3
)
−2
3
Q3
λ3
− ∂
∂z
〈
w′u′u′
〉
, (8b)
∂
〈
v′v′
〉
∂t
= 0 = − Q
3λ2
(〈
v′v′
〉
− Q
2
3
)
− 2
3
Q3
λ3
− ∂
∂z
〈
w′v′v′
〉
, (8c)
∂
〈
w′w′
〉
∂t
= 0 = − Q
3λ2
(〈
w′w′
〉
− Q
2
3
)
− 2
3
Q3
λ3
+ 2g〈
T
〉
〈
w′T ′
〉
− ∂
∂z
〈
w′w′w′
〉
, (8d)
where g is the gravitational acceleration, Q is the characteristic turbulent velocity given by√〈
u′i u′i
〉
, and λ1, λ2, and λ3 are the characteristic length scales for the triple velocity correla-
tions, the pressure–velocity gradient correlations, and the viscous dissipation, respectively.
These three length scales are determined by λi = ai × L , where ai and Cw are similarity
constants.
Similarly, the closure approximation of the turbulent scalar flux for ξ (
〈
w′ξ ′
〉
) can be
derived using:
∂
〈
w′ξ ′
〉
∂t
= 0 = −
〈
w′w′
〉 ∂ 〈ξ 〉
∂z
− Q
3λ2
〈
w′ξ ′
〉
+ g〈
T
〉
〈
T ′ξ ′
〉
− 2 Q
λ3
〈
w′ξ ′
〉
+ ∂
∂z
〈
w′w′ξ ′
〉
.
(9)
Note that the third terms on the right-hand side of equations (8d) and (9) are the buoyant
production/consumption that necessitate coupling between momentum and scalar exchange.
The steady-state budget equation for the correlation term
〈
T ′ξ ′
〉
is given as
∂
〈
T ′ξ ′
〉
∂t
= 0 = −
〈
w′ξ ′
〉 ∂ 〈ξ 〉
∂z
−
〈
w′T ′
〉 ∂ 〈ξ 〉
∂z
− 2 Q
λ3
〈
T ′ξ ′
〉
+ ∂
∂z
〈
w′T ′ξ ′
〉
. (10)
In Eqs. 8–10, the triple correlation terms are extra unknowns that require closure assump-
tions. The following second-order closure approximation proposed by Mellor (1973) and
modified by Wilson and Shaw (1977) is employed to parameterize the turbulent fluxes of the
Reynolds stresses in Eq. 8,
〈
w′u′w′
〉
= −2Qλ1
∂
〈
u′w′
〉
∂z
, (11a)
〈
w′u′u′
〉
= −Qλ1
∂
〈
u′u′
〉
∂z
, (11b)
〈
w′v′v′
〉
= −Qλ1
∂
〈
v′v′
〉
∂z
, (11c)
〈
w′w′w′
〉
= −3Qλ1
∂
〈
w′w′
〉
∂z
. (11d)
Similarly, the turbulent fluxes of the scalars and buoyant terms can be expressed as
〈
w′w′ξ ′
〉
= −2Qλ1
∂
〈
w′ξ ′
〉
∂z
, (11e)
〈
w′T ′ξ ′
〉
= −Qλ1
∂
〈
T ′ξ ′
〉
∂z
. (11f)
Details of the closure approximations for Eqs. 8–11 are discussed in Appendix A.
3.2.3 Third-order Closure Model
Unlike the second-order closure model, third-order closure schemes do not parameterize the
triple correlation terms but employ full budget equations. Here, the equations derived and
discussed in previous studies (Katul and Albertson 1998; Meyers and Paw 1986, 1987) are
used. To contrast with second-order closure modelling (Eq. 11b), their outcome for the triple
moments are summarized below (Katul and Albertson 1998; Meyers and Paw 1986, 1987):
〈
w′u′u′
〉
= − τ
C1
⎛
⎝2
〈
w′u′w′
〉 ∂ 〈u〉
∂z
+
〈
w′w′
〉 ∂
〈
u′u′
〉
∂z
+ 4
〈
u′w′
〉 ∂
〈
u′w′
〉
∂z
⎞
⎠, (12a)
〈
w′u′w′
〉
= − τ
C1
⎛
⎝
〈
w′w′w′
〉 ∂ 〈u〉
∂z
+
〈
u′w′
〉 ∂
〈
w′w′
〉
∂z
+ 3
〈
u′w′
〉 ∂
〈
u′w′
〉
∂z
⎞
⎠, (12b)
〈
w′w′w′
〉
= − τ
C1
⎛
⎝3
〈
w′w′
〉 ∂
〈
w′w′
〉
∂z
⎞
⎠, (12c)
〈
w′v′v′
〉
= − τ
C1
⎛
⎝
〈
w′w′
〉 ∂
〈
v′v′
〉
∂z
⎞
⎠, (12d)
〈
w′w′ξ ′
〉
= − τ
C1
⎛
⎝
〈
w′w′w′
〉 ∂ 〈ξ 〉
∂z
+
〈
w′ξ ′
〉 ∂
〈
w′w′
〉
∂z
+ 2
〈
w′w′
〉 ∂
〈
w′ξ ′
〉
∂z
⎞
⎠, (12e)
〈
w′T ′ξ ′
〉
= − τ
C1
(〈
w′w′ξ ′
〉 ∂ 〈T 〉
∂z
+
〈
w′w′T ′
〉 ∂ 〈ξ 〉
∂z
+
〈
w′w′
〉 ∂
〈
T ′ξ ′
〉
∂z
+
〈
w′T ′
〉 ∂
〈
w′ξ ′
〉
∂z
+
〈
w′ξ ′
〉 ∂
〈
w′T ′
〉
∂z
⎞
⎠. (12f)
In Eq. 12, the coefficient C1 is a similarity constant and τ is a relaxation time scale defined
as Q2/ε, where ε = 2
3
Q3
λ3
is the mean turbulent kinetic energy dissipation rate. More details
are provided in Appendix A.
3.3 Radiation Budget within the Canopy Volume
The radiative distribution within the canopy volume uses a multi-layer light attenuation model
described in Leuning et al. (1995) and Schafer et al. (2003). Basic concepts and applications
can be found elsewhere (Campbell and Norman 1998; Erbs et al. 1982; Goudriaan and van
Laar 1994; Pataki et al. 1998; Spitters et al. 1986; Spitters 1986; Stenberg 1998), but for
completeness, salient features of the model are summarised below.
The total incoming solar radiation contains two major components with different spectra:
the visible (VIS) and near-infrared radiation (NIR). Because the attenuation properties of VIS
and NIR within the canopy are different, the distributions of these two portions are treated
separately (Weiss and Norman 1985).
The isothermal form of net radiation (Rn) absorbed by the leaves at level z is given by
(Leuning et al. 1995)
Rn(z) = SA(z) − RL(z), (13)
where SA is the shortwave radiation absorbed by the leaves, and RL is the isothermal net
longwave radiation, discussed in Appendix B.
Within the canopy volume, sunlit leaves absorb all the radiative components, the direct
beam SA,b(z), diffuse SA,d(z), and scattered radiation SA,s(z), given by SA,sl(z) = SA,b(z)+
SA,d(z) + SA,s(z), while shaded leaves only absorb the diffuse and scattered components,
i.e. SA,sh(z) = SA,d(z)+ SA,s(z) (Campbell and Norman 1998; Leuning et al. 1995; Schafer
et al. 2003), where SA,sl(z) and SA,sh(z) are the absorbed shortwave radiation fluxes on sunlit
and shaded leaves, respectively (see Appendix B for a detailed discussion).
The vertical distribution of absorbed shortwave radiation for the sunlit (SA,sl) and shaded
leaves (SA,sh) can be calculated as
SA,sl(z) = fsl(z)SA,sl,VIS(z) + fsl(z)SA,sl,NIR(z), (14a)
SA,sh(z) = fsh(z)SA,sh,VIS(z) + fsh(z)SA,sh,NIR(z), (14b)
where fsl(z) and fsh(z) are the fractions of leaf area for sunlit and shaded leaves respectively,
with fsl(z)+ fsh(z) = 1. These two fractions are determined by the transmission coefficient
of the beam component, τb (see Appendix B):
fsl(z) = τb, (15a)
fsh(z) = 1 − τb. (15b)
3.4 Leaf-Level Energy Balance
The energy balance at the leaf scale is used to calculate leaf surface temperature TL (Campbell
and Norman 1998), and using the linearization technique of Penman (1948), TL is approxi-
mated as (Campbell and Norman 1998; Lai et al. 2000b):
TL = T +
SA − εLσ T 4 − Lvρ gbgs,vgb + gs,v D
Cpρgb + Lvρ gbgs,vgb + gs,v + Cpρgr
, (16)
where D is the water vapour deficit (kg kg−1), Cp is the specific heat of air (J mol−1 K−1), gr
is the radiative conductance (m s−1), σ is the Stefan–Boltzman constant, and gs,v and gb,v are
the stomatal and boundary-layer conductances (mol m−2 s−1) of water vapour, respectively,
and Lv is the latent heat of vaporization.
3.5 Physiological Model and Photosynthesis
The stomatal conductance for CO2, gs,c, can be estimated from the physiological model of
Collatz et al. (1991), originally proposed by Ball et al. (1987), given by
gs,c = m An RHsCs + b (17)
where m and b are parameters determined from gas-exchange measurements, RHs and Cs
are the relative humidity and CO2 concentration on the leaf surface, and An is the leaf net
assimilation rate, which can be calculated from the An −Ci curve after Farquhar et al. (1980),
and given by
An = κ1 (Ci − 
∗)
Ci + κ2 − Rd. (18)
Here κ1 = αpem Ip and κ2 = 2 ∗ for light-limited photosynthesis, κ1 = Vcmax and κ2 =
KC (1 + Oi/KO) for rubisco limited photosynthesis, Vcmax is the maximum catalytic capacity
of rubisco per unit leaf area, Rd is the dark respiration rate, αp is the leaf absorptivity for
PAR, em is the maximum quantum efficiency, Ip is the incident PAR flux density on the leaf
surface, ∗ is the CO2 concentration in the absence of Rd, KC and KO are Michaelis–Menten
constants for CO2 fixation and for oxygen inhibition with respect to CO2, and Oi is the oxy-
gen concentration (Campbell and Norman 1998; Katul et al. 2003; Lai et al. 2000a). Rd can
be estimated from Vcmax after Farquhar et al. (1980):
Rd = χVcmax (19)
where χ = 0.015 for this stand (Campbell and Norman 1998; Ellsworth 2000; Juang et al.
2006). The temperature dependency of Vcmax(t, z) can be expressed as:
Vcmax = Vcmax,25 exp[a1(TL − 25)]1 + exp[a2(TL − 41)] (20)
where a1 and a2 are the species-specific adjustment coefficients, which are obtained via
porometry measurement and Vcmax,25 is the value of Vcmax at 25◦C (Campbell and Norman
1998; Collatz et al. 1991; Farquhar et al. 1980). From previous studies conducted at the site
(Ellsworth 1999; Naumburg and Ellsworth 2000; Naumburg et al. 2001), the average Vcmax,25,
a1 and a2 are 59 µmol m−2 s−1, 0.051 and 0.205, respectively, for the upper canopy pine
foliage, and are 30µmol m−2 s−1, 0.088 and 0.290, respectively, for the sub-canopy broad-
leaved plants (Lai et al. 2002b). The seasonal acclimation of Vcmax,25, derived from leaf-gas
exchange measurements collected over a 3-year period, was reported by Ellsworth (1999) to
have an amplitude of about 18%. The importance of this seasonal acclimation vis-à-vis resolv-
ing the two-way interaction between the leaves and their microclimate is discussed below.
3.6 Combining the Sub-models: The CANVEG Framework
Summing up, the CANVEG equations for an arbitrary layer within the canopy volume are
as follows:
(1) The three mean scalar continuity equations for CO2, H2O, and temperature provide
three equations with nine unknowns (two mean scalar concentrations and one mean
temperature, three turbulent scalar fluxes, and three scalar sources or sinks).
(2) The turbulent transport models (at any closure level) establish three additional equations
linking turbulent fluxes to mean concentration (albeit the complexity of this expression
varies with the closure scheme and can be affected by thermal stratification). Hence,
the physical system alone combining mean scalar continuity and turbulent transport
theories provides six equations with nine unknowns.
(3) The assumption that mass and heat transfer from the leaf surface (or stomatal pores) to
the atmosphere is Fickian provides three additional equations that mathematically close
the original system (nine equations with nine unknowns) at the expense of introducing
four additional unknowns: the leaf stomatal conductance, and three internal or leaf-
level mean concentration states (surface temperature, surface humidity, and internal
CO2 concentration).
(4) The fact that the three scalars considered here are heat, water vapour, and CO2, permits
the formulation of three independent equations based on leaf energy balance (for surface
temperature), the Farquhar photosynthesis model (for inter-cellular CO2 concentration),
the hydration state of the leaves (for linking internal water vapour concentration with sur-
face temperature). Finally, the system is mathematically ‘closed’ via a semi-empirical
stomatal conductance model for well-hydrated leaves (i.e. Eq. 17).
The WMA eliminates steps (1) and (2) in the above procedure because the scalar con-
centration is assumed constant inside the canopy (as is the case in virtually all forest growth
and biogeochemical cycling models). Again, because leaf stomatal pores are assumed to be
saturated (in step 4), drought conditions known to induce stomatal closure are excluded from
the model-data comparisons here.
3.7 Boundary Conditions
Boundary conditions are needed to solve the governing equations. The long-term EC mea-
surements at the canopy top are used to derive similarity expressions for the normalized
standard deviation of the velocity components u, v, w, and the scalar quantities, θ , q, c as
a function of the atmospheric stability parameter ς = z−d0L , where Lis the Obukhov length.
The similarity functions for different quantities and for different stability classes are shown
in Fig. 2; the regression parameters obtained from a least squares analysis on these data
are summarized in Table 1. Despite the fact that the measurements were collected in the
roughness sublayer just above the canopy, they appear to be consistent with a number of field
experiments in the atmospheric surface layer summarized by Stull (1988) and Kaimal and
Finnigan (1994). The boundary conditions at the forest floor are given in Appendix C.
4 Results and Discussion
To address the two main questions here, the results and discussion are structured as follows:
the seasonal variation in leaf area density and its effect on the main driver for scalar pro-
duction and uptake inside the canopy (i.e, light levels) is first presented. Next, the modelled
space–time variations of the scalar sources, sinks, and fluxes for the three closure models
are presented to illustrate their dissimilarity inside the canopy at various times of the day
and different days of the year. These dissimilarities are needed as a bridge to explain the
differences in annual uptake or release of scalars by the canopy. With this background, the
mean scalar concentration profile calculations using the three closure schemes are compared
to the measurements (and to the calculations that employ the WMA). The computed scalar
fluxes using the WMA and the three closure models are also compared with the EC-based flux
measurements above the canopy. Statistical analyses of these model-model and model-data
intercomparisons provide quantitative answers to the first question. The second question is
addressed by contrasting two sets of model calculations conducted using the higher order
closure version: one set accounts for the seasonal variations in Vc max,25, while the other does
not (though the annual mean Vc max,25 is the same in both sets). The scalar flux differences
at the canopy top between these two sets of calculations are then compared to their counter-
parts when employing model calculations that use the WMA. Discussions about the need to
resolve thermal stratification inside the canopy are presented last.
4.1 Radiation Environment Within the Canopy Volume
The radiative components (PAR and net) directly influence leaf photosynthesis (described
later) and leaf energy exchange (latent heat and sensible heat). Incident shortwave radiation
Fig. 2 Normalized standard deviation of the velocity components u, u, w, the scalar quantities, θ , q, and c for 
unstable (left column) and stable (right column) atmospheric stability conditions. All quantities are analyzed 
when the friction velocity at the canopy top exceeds 0.1 m s−1. The solid lines represent the regression curves 
whose parameters are summarized in Table 1
at the top of the canopy was not directly measured before 2004, so to estimate the incident 
shortwave radiation for earlier periods, a relationship between PAR (in µmol m−2 s−1) and 
shortwave radiation (in W m−2) measured by the Kipp & Zonen CNR1 net radiometer 
installed at the top of the canopy (in the summer of 2003) was derived. This relationship 
is affected by zenith angle, location, and sky properties because the wavebands for PAR (in 
the spectral range of VIS, 400–700 nm) and shortwave radiation (VIS and NIR) are different. 
To convert units of PAR from µmol m−2 s−1 to W m−2, a mean waveband = 550 nm and a 
photon energy = 2.30×105 J mol−1 were used. The coefficient of determination (R2 = 0.99) 
between PAR and shortwave radiation measurements, shown in Fig. 3, suggests an excellent 
linear relationship although minor scatter remains due to cloudy conditions.
To illustrate how the canopy attenuates PAR throughout the day and across seasons, model 
results for PAR penetration ratios (i.e. the ratio of PAR at a given level zinside the canopy to 
the incident PAR at the top of the canopy) at midday (1200 LT) and in the morning (0800 LT) 
are shown in Fig. 4. The results are presented as 7-day moving averages for the entire year 
of 2005 for illustration. The normalized leaf area density profiles are shown for reference
Table 1 The similarity function expressions for normalized standard deviation of the wind components u, v,
w, and the scalar quantities, θ , q, c, for different atmospheric stability conditions (ς = (z − do)/L)
Stability Variable Form Constants
Unstable σu/u∗ a (b + cζ )d a = 1.95; b = 1.0; c = −1.5; d = 1/3
σv/u∗ a = 1.94; b = 1.0; c = −1.5; d = 1/3
σw/u∗ a = 1.15; b = 1.0; c = −1.5; d = 1/3
σθ /θ∗ a = 1.18; b = 0.05; c = −1.0; d = −1/3
σq/q∗ a = 1.28; b = 0.05; c = −1.0; d = −1/3
σc/c∗ a = 1.84; b = 0.05; c = −1.0; d = −1/3
Stable σu/u∗ a + bζ d a = 1.95; b = 0.39; d = 0.5
σv/u∗ a = 1.94; b = 0.37; d = 0.5
σw/u∗ a = 1.15; b = 0.23; d = 0.6
σθ /θ∗ C c = 3.20
σq/q∗ c = 3.47
σc/c∗ c = 5.00
Neutral σu/u∗ C c = 1.95
σv/u∗ c = 1.94
σw/u∗ c = 1.15
σθ /θ∗ c = 3.20
σq/q∗ c = 3.47
σc/c∗ c = 5.00
Fig. 3 The relationship between half-hourly measured photosynthetically active radiation (PAR) and mea-
sured incident shortwave radiation at the top of the canopy during the summer of 2003 (symbols). The
regression equation (solid line), derived from 1235 half-hourly data points, is given by y = 0.51x − 1.62, and
the coefficient of determination (R2) is 0.99. The theoretical slope for clear-sky conditions is 0.5
to illustrate their seasonal variations (and the contribution of the under-story dynamics). As
expected, more PAR (≈ 40% of the canopy top value) reaches the forest floor at midday
in winter months (when LAI is minimum). However, less than 15% of the canopy top PAR
Fig. 4 The vertical variation of modelled PAR penetration ratio throughout the year at different times during 
the day (1200 LT and 0800 LT) during 2005. The upper panel shows the normalized leaf area density profiles 
(LAD) for reference
penetrates the canopy volume during the growing season, when LAI is at its maximum value. 
For integration of scalar fluxes to annual time scales, it is important to note that the under-
storey dynamics and PAR attenuation play an important role. Because of the low LAI in the 
pine stand during the spring, significant PAR arrives at the understorey levels for at least 
2 months per year (e.g., see Julian days 130–180). Whether this understorey PAR contributes 
significantly to scalar fluxes is explored next.
4.2 Scalar Source Strength and Flux Distribution
Figure 5 shows the source (or sink) strength for different scalars at different times of day 
throughout the entire year of 2005 using the second-order closure model with atmospheric 
stability corrections as an example. For the CO2 source (or sink) strength, it is clear that while 
forest floor respiration dominates during nighttime conditions, the upper canopy becomes a 
strong sink during the day, especially during the growing season. The co-existence of both 
net carbon sinks and sources within the various canopy layers during the day is also evi-
dent during the winter season, and when the zenith angle is large (i.e., early morning and 
early evening). The model calculations suggest that the ecosystem is a large source of both 
latent and sensible heat during the growing season (as expected). Furthermore, it appears 
that canopy heating primarily occurs in the top one-third of the canopy where much of the 
radiation is intercepted by the main canopy foliage. The understorey sources and sinks of 
heat and water vapour appear much smaller than their pine counterpart. At noon, the under-
storey contributes to CO2 sequestration during the spring and early summer period, though 
its contribution during the remainder of the summer period appears minor.
Fig. 5 The time-depth distributions of modelled CO2 sink/source strength (mgC m−3 s−1, left column), latent
heat source strength (W m−3, middle column), and sensible heat sink/source strength (W m−3, right column)
at different times of day and different days of the year for 2005. All the model results shown here use the
second-order closure scheme with atmospheric stability corrections
Figure 6 presents the simulated time-depth variations of the turbulent fluxes for the same
time period shown in Fig. 5. These fluxes are computed by depth integrating the scalar
sources (or sinks) and adding to this sum the scalar flux emitted (or absorbed) by the forest
floor. Again, the CO2 flux distribution demonstrates that, while forest floor respiration is
large, this ecosystem remains a net carbon sink throughout the year as verified by eddy-
covariance measurements. Also, Fig. 6 shows how the understorey plays a significant role
in rapidly increasing CO2 uptake from the atmosphere during the early spring period, but its
relative impact on water vapour and sensible heat fluxes during this same period appears less
significant.
4.3 Comparison Between Different Closure Approximations
The ensemble-averaged scalar profiles (mean air temperature, mean water vapour concen-
tration, and mean CO2 concentration) within the canopy volume for all the three closure
schemes (i.e. first-order to third-order) and over different time periods in year 2005 are
shown in Fig. 7. The WMA and the scalar profile measurements are also presented for refer-
ence. The comparisons in Fig. 7 demonstrate that second- and third-order closure schemes
best agree with the mean concentration (or temperature) profile measurements for all three
scalars. Furthermore, the analysis here suggests that there is no gain in predictive skills by
increasing the closure from second-order to third-order as these two model results are almost
indistinguishable (Fig. 7). The main scalar concentration differences emerge when com-
paring the results from these two higher-order closure schemes with the first-order closure
model calculations. Table 2 summarizes the relative standard deviation (RSD, in percentage)
when comparing different closure approximations with the scalar profile measurements.
Fig. 6 The time-depth distributions of CO2 fluxes (mgC m−2 s−1, left column), latent heat fluxes (W m−2, 
middle column), and sensible heat fluxes (W m−2, right column) at different times of day and different days 
of the year for 2005. All the model results shown here use the second-order closure scheme with atmospheric 
stability corrections
These data-model comparisons, shown in Fig. 7, suggest that higher-order closure approxi-
mations (order 2 or more) are needed if accurate mean scalar profiles inside the canopy are 
desired. The follow-up question is whether the predictive skills in modelling scalar sources 
(or sinks) and fluxes significantly improve by increasing the order of the closure scheme (as 
is the case for the scalar concentration and temperature). Recall that the two-way interaction 
between the canopy and its micro-environment is better resolved with higher-order closure 
schemes.
Figure 8 presents the ensemble-averaged hourly comparisons for the scalar fluxes (sen-
sible and latent heat fluxes, and CO2 fluxes) at the top of canopy for the all closure models 
(in year 2005). Differences in nocturnal fluxes are very small among all the models, but 
during daytime conditions, differences between the models (and data) become significant. 
The reason the modelled CO2 fluxes above the canopy all agree with each other at night is 
due to the fact that these fluxes are dominated by the forest floor respiration (≈70% of the 
flux at the canopy top). This forest floor respiration is a lower boundary condition identically 
imposed on all the models. However, as noted earlier, linking this flux boundary condition at 
the ground to CO2 concentration becomes sensitive to the turbulence representation, which 
is significantly different among the models.
All closure models overestimate the EC measured sensible heat and CO2 fluxes at night. 
For the sensible heat, it is difficult to separate all the possible factors. These factors may 
include parameterization of the longwave radiation inside the canopy, the small mean wind 
speeds inside the canopy that lead to small boundary-layer conductance, and hence, small 
errors in this conductance leads to large errors in leaf temperature, and the fact that the 
sensible heat affects both the fluid and scalar transport equations so that errors in the flow 
parameterization can also affect the sensible heat flux. It is difficult to attribute the data-
model
Fig. 7 Ensemble-averaged profiles of mean air temperature (K, upper panel), mean water vapour concentration
(g kg−1, middle panel), and mean CO2 concentration (ppm, lower panel) for different closure approximations
over different times of day in 2005. Model results include the atmospheric stability corrections. Note that
the model results for the second- and third-order closure schemes are almost indistinguishable. The WMA is
shown as vertical dotted lines
difference to a single cause. For CO2, the overestimation by the model is not too different
from the underestimation expected by the EC measurements at night (even when the storage
flux is small). The fact that EC measurements tend to underestimate nighttime CO2 fluxes,
as discussed extensively for this site by Juang et al. (2006), suggests that the overestimation
by the model is reasonable.
The relative standard deviations in Table 2 demonstrate that higher-order closure schemes
perform significantly better than zeroth- and first-order closure models, especially for sen-
sible and latent heat fluxes. Note that the value of the modelled sensible heat flux changes
sign earlier in the afternoon when using the well-mixed assumption. This can have significant
impact on calculations of convective boundary-layer height and triggers of convective rainfall
(Juang et al. 2007a,b). Whether such differences in model performances inject flux biases at
ecologically relevant time scales (i.e. annual and longer) is explored next.
Table 2 The average relative standard deviation (in %) for all model results when compared with the half-
hourly eddy-covariance and canopy scalar profile measurements. ‘With stability’ refers to model results that
resolve the effects of atmospheric stability within the canopy volume
First-order
closure
Second-order
closure
Third-order
closure
With
stability
Without
stability
With
stability
Without
stability
Air temperature profiles 16.3 11.8 14.3 11.3 14.0
Water vapour concentration profiles 10.1 5.6 7.2 6.0 7.7
CO2 concentration profiles 21.5 14.0 16.7 13.6 16.0
Sensible heat flux at top of canopy 12.3 7.4 10.0 7.8 10.3
Latent heat flux at top of canopy 16.7 11.2 15.8 11.8 16.1
Carbon flux at top of canopy 23.0 17.7 20.2 17.4 19.8
Fig. 8 Hourly ensemble-averaged measured (open circles) and modelled (lines) sensible heat flux (HS), latent 
heat flux (Lv E) and carbon dioxide flux (FC ) at the top of the canopy in 2005. The model results were derived 
using different closure approximations (first-order to third-order) as well as the WMA for scalars. Ensemble 
averaging was conduced across each hour of day in 2005. Vertical bars represent one standard deviation around 
the ensemble-averaged eddy covariance measurements
4.4 Estimations of Annual NEE and ET
The main objective is to explore how the two-way interaction between leaves and their
microenvironment affects ecosystem annual net ecosystem CO2 exchange (NEE) and evapo-
transpiration (ET), and whether resolving this two-way interaction is as important as resolving
the seasonality in key physiological parameters. From previous discussions, the second-order
closure scheme was shown to optimally reproduce the mean scalar concentration (and tem-
perature) profiles and fluxes; its predictive skills were also comparable to the computationally
more-expensive third-order closure scheme. Hereafter, the second-order closure model with
atmospheric stability corrections is used as a ‘reference model’ for resolving the two-way
interaction between leaves and their microenvironment.
As earlier stated, Vcmax,25 of the pine foliage is known to vary seasonally for a number of
reasons including temperature acclimation, age, and leaf nitrogen changes. Detailed leaf-level
gas-exchange measurements reported by Ellsworth (1999) and inverse-model calculations
reported in Juang et al. (2006) suggest that seasonal variations in Vcmax,25 of pine foliage can
be approximated by a sinusoidal function given as
Vcmax,25[
Vcmax,25
]
annual
= 0.18 sin
[
2π
(
x − 120
365
)]
. (21)
where x is the day of the year. The annual NEE (from 2001 to 2005) derived from model
calculations with a fixed Vcmax,25 (= to the mean value) and the sinusoidally varying Vcmax,25
are compared in Fig. 9. Although differences in annual NEE values between these two model
results are small (< 5%), some differences in the patterns at different times of year are
discernable. For our purposes, the effect of non-stationarity (i.e. sinusoidal variations) in
Vcmax,25 on annual NEE fluxes appears relatively minor (because of partial cancellations at
the seasonal scales), at least when compared to biases in annual NEE introduced by using
the WMA.
Figure 9 also compares model results derived from the WMA with the EC measurements
for the years 2001–2005. While the comparison for all 5 years is shown, the discussion is
limited to the wet years (2001, 2004, and 2005). Recall that year 2002 experienced a severe
drought during the entire growing season (resulting in leaf area reduction in the pine stand),
and the forest floor respiration in 2003 was dramatically different because of excess branches
and litter produced by the damage from the ice storm in December of 2002. It suffices to note
that annual NEE was overestimated by the reference model in years 2002 and 2003 for these
two reasons. However, because the forest floor evaporation is much less sensitive to increased
respiring biomass in 2003, good agreement between ET predictions and EC measurements
were observed, in contrast to the poor agreement between EC-measured and modelled NEE.
Figure 9 also suggests that (i) the WMA model always predicts higher annual NEE
(≈6.6%) and ET (≈8.3%) when compared to the second-order closure model, and this over-
estimation appears larger than the annual flux adjustments obtained by resolving the 18%
variation in Vcmax,25.
4.5 The Effect of Atmospheric Stability
Having demonstrated that the two-way coupling between the vegetation and its microenvi-
ronment is well reproduced by the second-order closure model, we explore next whether
non-adiabatic simplifications may still yield good agreements with the EC measurements.
Resolving the adiabatic conditions inside the canopy at every level requires an iterative
scheme for every 30-min simulation run because the velocity and scalar transfer equations
Fig. 9 Measured and modelled annual net ecosystem exchange (NEE, upper panel) and annual evapotranspira-
tion (ET, lower panel). Model calculations are presented for the second-order closure scheme with atmospheric
stability but (i) with a constant Vc max,25, (ii) a sinusoidal-varied Vc max,25. Model calculations assuming well-
mixed scalar profiles inside the canopy are also shown. The eddy-covariance measurements from 2001 to 2005 
are presented as well
must be simultaneously solved. A non-adiabatic simplification permits us to ‘de-couple’ 
the solution of the flow field from the scalar transfer calculations, as was done earlier by 
Baldocchi (1992), Baldocchi and Meyers (1998), and Lai et al.  (2000a,b) in their Lagrangian 
CANVEG simulations. A previous field experiment in a deciduous forest at Camp Borden 
in Canada described by Shaw et al. (1988) reported that the effect of thermal stratification 
inside the canopy on the flow statistics can be as important as the effect of leaf area index 
variations. Whether this result holds for the pine forest here remains to be explored.
Table 2 summarizes the model-data comparisons when considering the effects of atmo-
spheric stability and when neglecting them. Accounting for atmospheric stability improves 
the modelled relative standard deviations by 10–25% for mean scalar concentration (and 
temperature) profiles and turbulent fluxes at the top of the canopy. The effect of atmospheric 
stability on annual NEE and ET is also explored in Fig. 10. The annual NEE comparisons 
suggest that assuming a neutral atmosphere degrades the modelled NEE by roughly 12–15%. 
For the annual ET, the model degradation is between 9 and 13%. Hence, atmospheric stability 
corrections can improve model performance by about 10% or more, on average.
Fig. 10 Annual net ecosystem exchange (NEE, upper panel) and annual evapotranspiration (ET, lower panel)
computed from the second-order closure model with and without atmospheric stability corrections. The EC
measurements collected from 2001 to 2005 are repeated here for reference
5 Summary and Conclusions
The study objectives sought answers to two interrelated questions via model calculations and
long-term measurements of mean scalar concentration and temperature profiles and eddy-
covariance scalar fluxes above the canopy. These questions were, (i) how detailed must the
turbulent transport model be to resolve the two-way interaction between leaves and their
microclimate? and (ii) what were the predictive skills gained by resolving this two-way
interaction vis-à-vis correcting for seasonality in key physiological parameters (chosen here
as Vcmax,25)? To address the first question, a number of turbulent models were used ranging
from zeroth-order (well-mixed assumption) to third-order closure. To address the second
question, model calculations using a non-stationary Vcmax,25 were compared to model cal-
culations with a constant Vcmax,25 case (though both have the same annual mean Vcmax,25).
We found that:
(1) First-order closure modelling captures much of the scalar source–sink variations within
the canopy as well as scalar fluxes above the canopy (to within 15%); however, second-
order closure models were necessary to reproduce the mean scalar concentration (and
temperature) profiles within the canopy (to within 8–15%).
(2) Third-order closure models perform no better than their second-order counterparts in
terms of overall predictive skills for scalar concentration profiles and fluxes, though
they were computationally much more demanding.
(3) The vertical flux most sensitive to the well-mixed assumption is sensible heat with a
clear bias in the afternoon periods. This bias can have important consequences when
modelling boundary-layer development and estimating the triggers of convective rain-
fall (Juang et al. 2007a,b).
(4) The comparison between measured and modelled annual NEE and ET using static and
dynamic Vcmax,25 do not show significant differences (≈5%). However, for shorter inte-
gration periods (e.g., 1 month or one season), the CO2 flux can be sensitive to variations
in Vcmax,25. This conclusion may not be general for all forest species, but indicative that
seasonal perturbations of ca. 20% around the annual mean in physiological parameters
do not bias mean annual fluxes. However, assuming a well-mixed concentration (and
temperature) can bias annual fluxes in this forested ecosystem by more than 5%.
(5) Changes in near-surface atmospheric stability have significant effects (10–20%) on the
scalar concentration profile distributions inside the canopy and scalar flux estimates
at the top of the canopy. The calculation of annual NEE and ET suggests that assum-
ing a non-neutral atmosphere can be more important than resolving the 18% amplitude
variation in Vcmax,25 for this evergreen forest, as long as the mean Vcmax,25 is not biased.
Acknowledgements The authors thank Ross McMurtrie for the helpful discussions on the use of WMA
in forest growth models. This study was supported, in part, by US Department of Energy (DOE) through
the Office of Biological and Environmental Research (BER), Terrestrial Carbon Processes (TCP) program
(Grants # 10509-0152, DE-FG02-00ER53015, and DE-FG02-95ER62083), the National Science Foundation
(NSF-EAR 0628342 and 0635787), and the Bi-national Agricultural Research and Development (BARD)
Fund (IS-3861-96).
Appendix A: Closure Approximation
A.1 Second-order Closure Parameterization for Triple Correlations
The triple correlation terms of the turbulent velocity components for the second-order closure
model are represented as (Katul and Albertson 1998; Mellor 1973):
〈
u′iu′ju′k
〉
= −Qλ1
⎛
⎝∂
〈
u′iu′j
〉
∂xk
+
∂
〈
u′ju′k
〉
∂xi
+
∂
〈
u′ku′i
〉
∂xj
⎞
⎠ . (22)
Similarly, the second-order closure approximation for the turbulent transport of scalar
fluxes and buoyant terms can be expressed as
〈
u′iu′jζ ′
〉
= −Qλ1
⎛
⎝∂
〈
u′jζ ′
〉
∂xi
+
∂
〈
u′iζ ′
〉
∂xj
⎞
⎠, (23a)
〈
u′iT ′ζ ′
〉
= −Qλ1
⎛
⎝∂
〈
T ′ζ ′
〉
∂xi
⎞
⎠. (23b)
A.2 Third-order Closure Parameterization
The closure expression for the triple correlation terms in Eq. 12 is derived from a simplified
steady-state budget described in Meyers and Paw (1986):
∂
〈
u′iu′jw′
〉
∂t
= −
〈
u′iw′w′
〉 ∂ 〈uj
〉
∂z
−
〈
u′jw′w′
〉 ∂ 〈ui〉
∂z
−
〈
w′w′
〉 ∂
〈
u′iu′j
〉
∂z
−
〈
u′iw′
〉 ∂
〈
w′u′j
〉
∂z
−
〈
u′jw′
〉 ∂
〈
w′u′i
〉
∂z
−
〈
u′iu′j
∂p′
∂z
〉
−
〈
u′iw′
∂p′
∂z
〉
−
〈
u′jw′
∂p′
∂z
〉
−
〈
u′iw′
〉 〈∂p′
∂xj
〉
−
〈
u′jw′
〉 〈∂p′
∂xi
〉
− Mij3 = 0, (24)
where Mij3 is the dissipation rate and is expressed as:
Mij3 = 2ν
〈
w′
∂u′i
∂xk
∂u′j
∂xk
+ u′i
∂w′
∂xk
∂u′j
∂xk
+ u′j
∂u′i
∂xk
∂w′
∂xk
〉
. (25)
A simplified return-to-isotropic model introduced by Andre et al. (1981) is applied to
parameterize the pressure and molecular terms in Eq. 24
〈
u′iu′j
∂p′
∂z
〉
+
〈
u′iw′
∂p′
∂z
〉
+
〈
u′jw′
∂p′
∂z
〉
+ Mij3 =
C8
〈
u′iu′jw′
〉
τ
. (26)
After rearranging equations 24 and 26, the general form of the triple correlation terms can
be represented as:
〈
u′iu′jw′
〉
= − τ
C8
⎛
⎝
〈
u′iw′w′
〉 ∂ 〈uj
〉
∂z
+
〈
u′jw′w′
〉 ∂ 〈ui〉
∂z
+
〈
w′w′
〉 ∂
〈
u′i u′j
〉
∂z
+
〈
u′iw′
〉 ∂
〈
w′u′j
〉
∂z
+
〈
u′jw′
〉 ∂
〈
w′u′i
〉
∂z
⎞
⎠ . (27)
Likewise, the steady-state expression for the triple correlation terms in the scalar term and
temperature fluctuation can be derived following the same procedure described in
Meyers and Paw (1987) to yield:
〈
u′iu′jζ ′
〉
= − τ
C8
⎛
⎝
〈
u′iu′jw′
〉 ∂ 〈ζ 〉
∂z
+
〈
w′u′jζ ′
〉 ∂ 〈ui〉
∂z
+
〈
w′u′iζ ′
〉 ∂ 〈uj
〉
∂z
+
〈
w′u′i
〉 ∂
〈
u′jζ ′
〉
∂z
+
〈
w′u′j
〉 ∂
〈
u′iζ ′
〉
∂z
+
〈
w′ζ ′
〉 ∂
〈
u′iu′j
〉
∂z
⎞
⎠ , (28)
and
〈
u′iT ′ζ ′
〉
= − τ
C8
⎛
⎝
〈
w′u′iζ ′
〉 ∂ 〈T 〉
∂z
+
〈
w′u′iT ′
〉 ∂ 〈ζ 〉
∂z
+
〈
w′T ′ζ ′
〉 ∂ 〈ui〉
∂z
+
〈
w′u′i
〉 ∂
〈
T ′ζ ′
〉
∂z
+
〈
w′T ′
〉 ∂
〈
u′iζ ′
〉
∂z
+
〈
w′ζ ′
〉 ∂
〈
u′iT ′
〉
∂z
⎞
⎠ . (29)
For lower boundary conditions, vertical gradients in the triple moments were assumed to be
zero at the forest floor.
Appendix B: Light Attenuation Model
B.1 Partitioning of Incoming Shortwave Radiation
The incoming shortwave radiation received by the surface is composed of two portions: beam
Sb0 and diffuse Sd0 components. To quantify the partitions of Sb0 and Sd0, the atmospheric
transmissivity (or clearness index) Ka needs to be determined, and estimated as
Ka = S0Ssc sin β , (30)
where S0 is the measured incident radiation on the surface, Ssc is the solar radiative flux
density normal to the solar beam outside the atmosphere, and sin β is the solar elevation
(Erbs et al. 1982; Goudriaan and van Laar 1994; Leuning et al. 1995).
Once Ka is determined from S0, the fraction of diffuse component fd = Sd0/S0 can be
quantified based on the experimental relationship by Erbs et al. (1982):
fd = 1.0 − 0.09Ka if Ka ≤ 0.22
fd = 0.95 − 0.16Ka + 4.38K 2a − 16.63K 3a + 12.33K 4a if 0.22 < Ka ≤ 0.8
fd = 0.16 if Ka > 0.8
(31)
Then Sb0 and Sd0 at the top of the canopy are calculated from fd using:
Sb0 = S0 (1 − fd), (32a)
Sd0 = S0 fd. (32b)
B.2 Radiation Attenuation Calculation
B.2.1 Beam Component
The extinction coefficient of the beam component, Kb, for an ellipsoidal leaf distribution can
be determined as (Campbell 1986; Campbell and Norman 1998):
Kb (ψ) = (x
2 + tan2 ψ)0.5
x + 1.774 (x + 1.182)−0.733 , (33)
where ψ is the solar zenith angle, x is the leaf angle distribution index. Note that Eq. 33 is the
general expression for Kb, which includes the case of ellipsoidal leaf distribution (x 	= 1)
and spherical leaf angle distribution (x = 1). In this study, x = 1 was assumed, which
is reasonable for coniferous foliage (Schafer et al. 2003). The fraction τb (or transmission
coefficient) for the incoming beam from a given ψ at different levels is given as:
τb(ψ) = exp [−Kb(ψ) L t(z) ], (34)
where L t(z) is the cumulative PAD above level z, and  is the clumping factor due to the
shading effect of other leaves (Stenberg 1998). The direct beam absorbed at level z can be
expressed as:
Sabs,b(z) = Sbo Kbα, (35)
where α is the absorptivity of the leaves for a given radiative waveband, which is ≈ 0.2 for
VIS, and ≈ 0.8 for NIR (Campbell and Norman 1998).
B.2.2 Diffuse Component
To determine the extinction coefficient of diffuse radiation (Kd) we approximated the fitted
curve as a function of L t(z) for x = 1 after Campbell and Norman (1998) to give:
Kd(L t) = −0.035 log210 L t(z) − 0.16 log10 L t(z) + 0.82 if L t > 0.01
Kd(L t) = 1 if L t ≤ 0.01 . (36)
The absorbed diffuse radiation Sabs,d(z) can be expressed as:
Sabs,d(z) = Sdo α0.5 Kd(1 − ρcd)τd, (37)
whereρcd is the canopy reflection coefficient for diffuse component, and τd is the transmission
coefficient for the diffuse component and is given as:
τd = exp
[
−α0.5 Kd(L t)L t(z)
]
. (38)
B.2.3 Scattered Component Inside the Canopy
The total beam transmitted and scattered through the canopy to a given depth L t (z) can be
derived using the beam extinction coefficient (Eq. 31):
τs (ψ) = exp
[
−α0.5 Kbe(ψ)L t(z)
]
. (39)
Therefore, the absorbed transmitted/scattered radiation at a given depth Lt (z) can be
expressed as:
Sabs,s(z) = Sb0
[
α0.5 Kb(1 − ρcb)τs − Kbατb
]
. (40)
B.3 Longwave Radiation Inside the Canopy Volume
The isothermal net longwave radiation (Leuning et al. 1995) at the top of the canopy can be
represented as:
RL0 = (εc − εatm) σ T 40 , (41)
where εc (≈ 0.97) is the bulk canopy emissivity, σ is the Stephan Boltzmann constant, T0
is the air temperature at the canopy top, and εatm is the sky emissivity empirically given by
(Brutsaert 1975):
εatm = 0.642 (e0/T0)1/7 , (42)
where e0 is water vapour pressure at the top of the canopy.
Once RL0 is obtained, the isothermal net longwave radiation through the canopy volume
is estimated as (Leuning et al. 1995):
RL (z) = RL0 Kd exp [−Kd L t(z)]. (43)
Appendix C: Forest Floor Boundary Conditions
For the mean velocity, the no-slip boundary condition was imposed at the forest floor. For
all higher-order flow statistics, a free-slip condition was assumed (i.e. ∂(·)/∂z = 0). As for
CO2, the forest floor efflux (Fff ), derived by fitting chamber respiration measurements to soil
temperature (Tsoil, in ◦C) and soil moisture (θ , in m3 m−3), was used (Palmroth et al. 2005).
The canonical form of the regression equation used to fit the chamber measurements is:
Fff (Tsoil, θ) = Rbea1Tsoil
[
1 − e(−bθ+c)
]
(44)
where Rb is the base respiration, a1 is the temperature sensitivity (Q10 = exp(10a1)) when
soil moisture content is not limiting, and the constants b and c are soil moisture reduction
parameters. The respiration parameters for 2001–2003 were given in Palmroth et al. (2005)
and summarized by Juang et al. (2006). For 2004–2005, these parameters were Rb = 0.817,
a1 = 0.103, b = −31.121, c = 4.628 for 2004, and Rb = 0.548, a1 = 0.118, b = −29.511,
c = 3.301 for 2005.
As for the latent heat flux from the forest floor (Lv Eff ), it is expressed as a function of
the Bowen ratio βo, and given by
Lv Eff = 11 + βo
(
Rn,0 − G
) (45)
where Rn,0 is the forest floor net radiation estimated from the radiative attenuation model in
Sect. 3.3 and Appendix B, and G is the soil heat flux simply calculated from the measured
soil temperature Tsoil at 0.1 m depth, air temperature just above the ground, and the mean
thermal conductivity for soil (≈0.5 W m−2 K−1). Because the soil heat flux is small, heat
storage in the upper soil layer is neglected. The value of βo was computed from the surface
gradients of air temperature ((∂T /∂z)
s
and humidity (∂q/∂z)s) using
βo =
(
C p
) (
∂T /∂z
)
s
+ d
, (46)Lv (∂q/∂z)s
where d is the dry adiabatic lapse rate. The surface gradients above were extrapolated 
from the measured mean water vapour concentration and temperature data. After Lv Eff is
determined, the forest floor sensible heat flux (HSff ) can be estimated as a residual using
HSff = Rn,0 − G − Lv Eff . Note that unlike the CO2 fluxes at the forest floor (specified
using the chamber data), only extrapolated mean gradients from the profiles measurements
were used to determine βo(i.e. the partitioning between sensible and latent heat fluxes). The
radiative model was used to compute Rn,0 in Eq. 45.
References
Aber JD, Reich PB, Goulden ML (1996) Extrapolating leaf CO2 exchange to the canopy: a generalized model
of forest photosynthesis compared with measurements by eddy correlation. Oecologia 106:257–265
Andre JC, Lacarrere P, Traore K (1981) Pressure effects on triple correlations in turbulent convective flows.
Turbulent shear flows III. Springer-Verlag, Berlin, pp 243–252
Baldocchi D (1992) A Lagrangian random-walk model for simulating water vapour, CO2 and sensible heat flux
densities and scalar profiles over and within a soybean canopy. Boundary-Layer Meteorol 61:113–144
Baldocchi D, Meyers T (1998) On using eco-physiological, micrometeorological and biogeochemical theory
to evaluate carbon dioxide, water vapour and trace gas fluxes over vegetation: a perspetive. Agric Forest
Meteorol 90:1–25
Baldocchi DJ, Vogel C, Hall B (1997) Seasonal variation of carbon dioxide exchange rates above and below
a boreal jack pine forest. Agric Forest Meteorol 83:147–170
Baldocchi DD, Falge E, Gu L, Olson R, Hollinger D, Running S, Anthoni P, Bernhofer C, Davis K, Fuentes J,
Goldstein A, Katul G, Law B, Lee X, Malhi Y, Meyers T, Munger JW, Oechel W, Pilegaard K, Schmid
HP, Valentini R, Verma S, Vesala T, Wilson K, Wofsy S (2001) FLUXNET: a new tool to study the tem-
poral and spatial variability of ecosystem-scale carbon dioxide, water vapour and energy flux densities.
Bull Amer Meteorol Soc 82:2415–2435
Ball JT, Woodrow IE, Berry JA (1987) A model predicting stomatal conductance and its contribution to the
control of photosynthesis under different environmental conditions. In: Biggins J (ed) Progress in pho-
tosynthesis research. Martinus Nijhoff, Amsterdam Netherlands pp 221–224
Brutsaert W (1975) On a derivable formula for long-wave radiation from clear skies. Water Resour Res 11:742–
744
Butnor JR, Johnsen KH (2004) Calibrating soil respiration measures with a dynamic flux apparatus using
artificial soil media of varying porosity. Eur J Soil Sci 55:639–647
Butnor JR, Johnsen KH, Oren R (2003) Reduction of forest floor respiration by fertilization on both carbon
dioxide-enriched and reference 17-year-old loblolly pine stands. Global Change Biol 9:849–861
Campbell GS (1986) Extinction coefficients for radiation in plant canopies calculated using an ellipsoidal
inclination angle distribution. Agric Forest Meteorol 36:317–321
Campbell GS, Norman JM (1998) An introduction to environmental biophysics. Springer-Verlag, New York
307 pp
Collatz GJ, Ball JT, Crivet C, Berry JA (1991) Physiological and environmental regulation of stomatal con-
ductance, photosynthesis and transpiration: a model that Includes a laminar boundary layer. Agric Forest
Meteorol 54:107–136
De Pury DGG, Farquhar GD (1997) Simple scaling of photosynthesis from leaves to canopies without the
errors of big-leaf models. Plant Cell Environ 20:537–557
Detto M, Katul GG (2007) Simplified expressions for adjusting higher-order turbulent statistics obtained from
open path gas analyzers. Boundary-Layer Meteorol 122:205–216
Ellsworth DS (1999) CO2 enrichment in a maturing pine forest: are CO2 exchange and water status in the
canopy affected. Plant Cell Environ 22:461–472
Ellsworth DS (2000) Seasonal CO2 assimilation and stomatal limitations in a Pinus taeda canopy. Tree Physiol
20:435–445
Erbs DG, Klein SA, Duffie JA (1982) Estimation of the diffuse radiation fraction for hourly, daily and monthly
average global radiation. Sol Energy 28:293
Farquhar GD, Von Caemmerer S, Berry JA (1980) A biochemical model of photosynthetic CO2 assimilation
in leaves of C3 species. Planta 149:79–90
Goudriaan J, van Laar HH (1994) Modelling potential crop growth processes. Kluwer Academic Publishers,
Dordrecht 256 pp
Gu L, Shugart HH, Fuentes JD, Black TA, Shewchuk SR (1999) Micrometeorology, biophysical exchanges
and NEE decomposition in a two-story boreal forest—development and test of an intergrated model.
Agric Forest Meteorol 94:123–148
Jackson PS (1981) On the displacement height in the logarithmic velocity profile. J Fluid Mech 111:15–25
Juang J-Y, Katul GG, Siqueira MBS, Stoy PC, Palmroth S, McCarthy HR, Kim H, Oren R (2006) Modeling
nighttime ecosystem respiration from measured CO2 concentration and air temperature profiles using
inverse methods. J Geophys Res-Atmos 111:D08S05
Juang J-Y, Katul GG, Porporato A, Stoy PC, Siqueira M, Detto M, Kim H-S, Oren R (2007a) Eco-hydrological
controls on summertime convective rainfall triggers. Global Change Biol 13(4):887–896
Juang J-Y, Porporato A, Stoy PC, Siqueira MBS, Oishi CA, Detto M, Kim H-S, Katul GG (2007b) Hydrologic
and atmospheric controls on initiation of convective precipitation events. Water Resour Res 43:W03421
Kaimal JC, Finnigan JJ (1994) Atmospheric boundary layer flows: their structure and measurements. Oxford
University Press, 289 pp
Katul GG, Albertson JD (1998) An investigation of higher-order closure model for a forested canopy. Bound-
ary-Layer Meteorol 89:47–74
Katul GG, Chang WH (1999) Principal length scales in second-order closure models for canopy turbulence.
J Appl Meteorol 38:1631–1643
Katul GG, Leuning R, Oren R (2003) Relationship between plant hydraulic and biochemical properties derived
from a steady-state coupled water and carbon transport model: Plant Cell Environ 26:339–350
Katul GG, Mahrt L, Poggi D, Sanz C (2004) One and two equation models for canopy turbulence. Boundary-
Layer Meteorol 113:91–109
Kirschbaum MUF, Küppers M, Schneider HGC, Noe S (1998) Modelling photosynthesis in fluctuating light
with inclusion of stomatal conductance, biochemical activation and pools of key photosynthetic inter-
mediates. Planta 204:16–26
Lai C-T, Katul GG, Oren R, Ellsworth D, Schafer K (2000a) Modeling CO2 and water vapour turbulent flux
distributions within a forest canopy. J Geophys Res Atmos 105(D21):26333–26351
Lai C-T, Katul GG, Ellsworth D, Oren R (2000b) Modelling vegetation-atmosphere CO2 exchange by a
coupled Eulerian–Langrangian approach. Boundary-Layer Meteorol 95:91–122
Lai C-T, Katul GG, Butnor JR, Siqueira M, Ellsworth DS, Maier C, Johnsen KH, McKeand S, Oren
R (2002a) Modeling the limits on the response of net carbon exchange to fertilization in a southeastern
pine forest. Plant Cell Environ 25:1095–1119
Lai C-T, Katul GG, Butnor J, Ellsworth D, Oren R (2002b) Modelling night-time ecosystem respiration by a
constrained source optimization method. Global Change Biol 8:124–141
Leuning R, Kelliher FM, De Pury DGG, Schulze E-D (1995) Leaf nitrogen, photosynthesis, conductance and
transpiration: scaling from leaves to canopies. Plant Cell Environ 18:1183–1200
Leuning R (2000) Estimation of scalar sources/sink distributions in plant canopies using Lagrangian disper-
sion analysis: corrections for atmospheric stability and comparison with a multilayer canopy model.
Boundary-Layer Meteorol 96:293–314
Liu J, Black TA, Novak MD (1996) E-epsilon modeling of turbulent air flow downwind of a model forest
edge. Boundary-Layer Meteorol 77:21–44
Luo Y, Medlyn B, Hui D, Ellsworth DS, Reynolds JF, Katul G (2001) Gross primary productivity in the Duke
forest: modeling synthesis of the free-air CO2 enrichment experiment and eddy-covariance measure-
ments. Ecol Appl 11:239–252
Malhi Y, Nobre AD, Grace J, Kruijt B, Pereira MGP, Culf A, Scott S (1998) Carbon dioxide transfer over a
Central Amazonian rain forest. J Geophys Res Atmos 103:31593–31612
Massman WJ, Weil JC (1999) An analytical one-dimensional second-order closure model of turbulence statis-
tics and the Lagrangian time scale within and above plant canopies of arbitrary structure. Boundary-Layer
Meteorol 91:81–107
McCarthy HR, Oren R, Finzi AC, Johnsen KH (2006a) Canopy leaf area constrains [CO2]-induced enhance-
ment of productivity and partitioning among aboveground carbon pools. Proc Natl Acad Sci USA
103:19356–19361
McCarthy HR, Oren R, Kim H-S, Johnsen KH, Maier C, Pritchard SG, Davis MA (2006b) Interaction of ice
storms and management practices on current carbon sequestration in forests with potential mitigation
under future CO2 atmosphere. J Geophys Res Atmos 111:D15103
McCarthy HR, Oren R, Finzi AC, Ellsworth DS, Kim H-S, Johnsen KH, Millar B (2007) Temporal dynamics
and spatial variability in the enhancement of canopy leaf area under elevated atmospheric CO2. Global
Change Biol 13:2479–2497
Mellor G (1973) Analytic prediction of the properties of stratified planetary boundary layer. J Atmos Sci
30:1061–1069
Meyers T, Paw UKT (1986) Testing of a higher-order closure model for modeling airflow within and above
plant canopies. Boundary-Layer Meteorol 37:297–311
Meyers T, Paw UKT (1987) Modeling the plant canopy micrometeorology with higher-order closure princi-
ples. Agric Forest Meteorol 41:143–163
Naumburg E, Ellsworth DS (2000) Photosynthetic sunfleck utilization potential of understory saplings grow-
ing under elevated CO2 in FACE. Oecologia 122:163–174
Naumburg E, Ellsworth DS, Katul GG (2001) Modeling daily understory photosynthesis of species with dif-
fering photosynthetic light dynamics in ambient and elevated CO2. Oecologia 126:487–499
Nikolov NT, Massman WJ, Schoettle AW (1995) Coupling biochemical and biophysical processes at the leaf
level: an equilibrium photosynthesis model for leaves of C3 plants. Ecol Model 80:205–235
Oren R, Ewers B, Todd P, Phillips N, Katul GG (1998) Water balance delineates the soil layer in which
moisture affects canopy conductance. Ecol Appl 8:990–1002
Palmroth S, Maier CA, McCarthy HR, Oishi AC, Kim H-S, Johnsen KH, Katul GG, Oren R (2005) Contrast-
ing responses to drought of forest floor CO2 efflux in a loblolly pine plantation and a nearby oak-hickory
forest. Global Change Biol 11:1–14
Pataki DE, Oren R (2003) Species differences in stomatal control of water loss at the canopy scale in a mature
bottomland deciduous forest. Adv Water Resour 26(12):1267–1278
Pataki DE, Oren R, Tissue D (1998) Elevated carbon dioxide does not affect average canopy stomatal con-
ductance of Pinus taeda L. Oecologia 127:549–559
Penman HL (1948) Natural evaporation from open water, bare soil and grass. Proc Roy Soc London Ser A
193:120–146
Poggi D, Katul GG, Albertson JD (2004) Momentum transfer and turbulent kinetic energy budgets within a
dense model canopy. Boundary-Layer Meteorol 111:589–614
Raupach MR (1988) Canopy transport processes. In: Steffen WL, Denmead OT (eds) Flow and transport in
the natural environment: advances and applications. Springer-Verlag, New York, pp 95–127
Raupach MR (1989) A practical Lagrangian method for relating scalar concentrations to source distributions
in vegetation canopies. Quart J Roy Meteorol Soc 115:609–632
Schafer KVR, Oren R, Ellsworth DS, Lai C-T, Herrick JD, Finzi AC, Richter DD, Katul GG (2003) Exposure
to an enriched CO2 atmosphere alters carbon assimilation and allocation in a pine forest ecosystem.
Global Change Biol 9:1378–1400
Shaw RH, Hartog GD, Neumann HH (1988) Influence of foliar density and thermal stability on profiles of
Reynolds stress and turbulence intensity in a deciduous forest. Boundary-Layer Meteorol 45(4):391–409
Simon E, Meixner FX, Ganzeveld L, Kesselmeier J (2005a) Coupled carbon-water exchange of the Amazon
rain forest, I. Model description, parameterization and sensitivity analysis. Biogeosciences 2:231–253
Simon E, Meixner FX, Rummel U, Ganzeveld L, Ammann C, Kesselmeier J (2005b) Coupled carbon-water
exchange of the Amazon rain forest, II. Comparison of predicted and observed seasonal exchange of
energy, CO2, isoprene and ozone at a remote site in Rondônia. Biogeosciences 2:255–275
Siqueira M, Katul G (2002) Estimating heat sources and fluxes in thermally stratified canopy flows using
higher-order closure models. Boundary-Layer Meteorol 103:125–142
Siqueira M, Kutal GG, Lai C-T (2002) Quantifying net ecosystem exchange by multilevel ecophysiological
and turbulent transport models. Adv Water Resour 25:1357–1366
Spitters CJT, Toussaint HAJM, Goudriaan J (1986) Separating the diffuse and direct component of global
radiation and its implications for modelling canopy photosynthesis. Part I. Components of incoming
radiation. Agric Forest Meteorol 38:225–237
Spitters CJT (1986) Separating the diffuse and direct component of global radiation and its implications for
modelling canopy photosynthesis. Part II. Calculation of canopy photosynthesis. Agric Forest Meteorol
38:238–241
Stenberg P (1998) Implications of shoot structure on the rate of photosynthesis at different levels in a conif-
erous canopy using a model incorporating grouping and penumbra. Funct Ecol. 12: 82–91
Stoy P, Katul GG, Siqueira MBS, Juang J-Y, McCarthy HR, Kim H-S, Oishi C, Oren R (2005) Variability in
net ecosystem exchange from hourly to inter-annual time scales at adjacent pine and hardwood forests:
a wavelet analysis. Tree Physiol 25:887–902
Stull RB (1988) An introduction to boundary layer meteorology. Kluwer Academic Publishers Group, Dordr-
echt, 666 pp
Su H-B, Paw UKT, Shaw RH (1996) Development of a coupled leaf and canopy model for the simulation of
plant-atmosphere interaction. J Appl Meteorol 35:733–748
Thom AS (1971) Momentum absorption by vegetation. Quart J Roy Meteorol Soc 97: 414–428
Wang YP, Leuning R (1998) A two-leaf model for canopy conductance, photosynthesis and partitioning of
available energy I: Model description and comparison with a multi-layered model. Agric Forest Meteorol
91:89–111
Webb EK, Pearman GI, Leuning R (1980) Correction of flux measurements for density effects due to heat and
water vapour transfer. Quart J Roy Meteorol Soc 106: 85–100
Weiss A, Norman JM (1985) Partitioning solar radiation into direct and diffuse, visible and near-infrared
components. Agric Forest Meteorol 34:205–213
Williams M, Rastetter EB, Fernades DN, Goulden MLC, Wofsy S, Shaver GR, Meillo JM, Munger JW, Fan
SM, Nadelhoffer KJ (1996) Modeling the soil–plant–atmosphere continuum in a Quercus-Acerstand at
Harvard forest: the regulation of stomatal conductance by light, nitrogen and soil/plant hydraulic prop-
erties. Plant Cell Environ 19:911–927
Williams M, Malhi Y, Nobre AD, Rastetter EB, Grace J, Pereira MGP (1998) Seasonal variation in net car-
bon exchange and evapotranspiration in a Brazilian rain forest: a modelling analysis. Plant Cell Environ
21:953–968
Wilson JD (1989) Turbulent transport within the plant canopy. In: Blackman TA et al (eds) Estimation of areal
evapotranspiration. IAHS Pub., Gentbrugge, Belgium, pp 43–80
Wilson NR, Shaw RH (1977) A higher order closure model for canopy flow. J Appl Meteorol 16:1198–1205