The Chemistry of Atmosphere-Forest Exchange ( CAFE ) Model – Part 1 : Model description and characterization

We present the Chemistry of Atmosphere-Forest Exchange (CAFE) model, a vertically-resolved 1-D chemical transport model designed to probe the details of nearsurface reactive gas exchange. CAFE integrates all key processes, including turbulent diffusion, emission, deposition and chemistry, throughout the forest canopy and mixed layer. CAFE utilizes the Master Chemical Mechanism (MCM) and is the first model of its kind to incorporate a suite of reactions for the oxidation of monoterpenes and sesquiterpenes, providing a more comprehensive description of the oxidative chemistry occurring within and above the forest. We use CAFE to simulate a young Ponderosa pine forest in the Sierra Nevada, CA. Utilizing meteorological constraints from the BEARPEX-2007 field campaign, we assess the sensitivity of modeled fluxes to parameterizations of diffusion, laminar sublayer resistance and radiation extinction. To characterize the general chemical environment of this forest, we also present modeled mixing ratio profiles of biogenic hydrocarbons, hydrogen oxides and reactive nitrogen. The vertical profiles of these species demonstrate a range of structures and gradients that reflect the interplay of physical and chemical processes within the forest canopy, which can influence net exchange.


Introduction
At the forest-atmosphere interface, biogenic emissions, surface deposition and anthropogenic pollutants interact with significant impacts on atmospheric composition and ecosystem function.Globally, biogenic activity accounts for more than 80% of volatile organic compound (VOC) emissions (Guenther et al., 1995).Forest emissions include iso-Correspondence to: J. A. Thornton (thornton@atmos.washington.edu)prene (C 5 H 8 ), 2-methyl-3-buten-2-ol (MBO), monoterpenes, sesquiterpenes and a host of small-chain organic compounds.Oxidation of these VOC influences ozone (O 3 ) production and secondary organic aerosol (SOA) formation (Hallquist et al., 2009), with broad implications for air quality and climate (Goldstein et al., 2009;Isaksen et al., 2009).Forests also provide ample surface area that facilitates dry deposition of reactive nitrogen, O 3 , sulfur dioxide (SO 2 ), aerosols and other atmospheric constituents (Fowler et al., 2009).Atmospheric deposition, in turn, influences forest productivity.For example, wet and dry deposition of atmospheric ammonia, nitric acid and other oxidized nitrogen compounds represents an important source of bio-available nitrogen that can modulate carbon sequestration by forests (Magnani et al., 2007;Thomas et al., 2010).In contrast, ecosystem stress from O 3 deposition can reduce carbon uptake (Sitch et al., 2007) and even alter emissions (Schade and Goldstein, 2002).The complexity of this system lends itself to an array of biosphere-atmosphere feedbacks that are sensitive to temperature, radiation, atmospheric composition and other parameters (Carslaw et al., 2010;Fuentes et al., 2001).
Detailed chemical-transport models that resolve the vertical structure of processes throughout the forest canopy are powerful tools for assessing and documenting chemical contributions to the fluxes of reactive species.Several such models have been developed to explore forest-atmosphere exchange of SO 2 , O 3 , VOC and nitrogen oxides (Forkel et al., 2006;Ganzeveld et al., 2002;Gao et al., 1993Gao et al., , 1991;;Makar et al., 1999;Meyers, 1987;Stroud et al., 2005;Walton et al., 1997;Baldocchi, 1988;Boy et al., 2010).We present a new model optimized for investigating the effects of canopy-atmosphere interactions on atmospheric composition and reactive gas exchange, the Chemistry of Atmosphere-Forest Exchange (CAFE) model.CAFE represents a unique extension of many previous approaches in that it incorporates the extensive Master Chemical Mechanism (MCM), permitting a more comprehensive description of the chemical Published by Copernicus Publications on behalf of the European Geosciences Union.
G. M. Wolfe and J. A. Thornton: The (CAFE) Model -Part 1 contributions to fluxes over forested regions.Model design has focused on simulating results from the Biosphere Effects on Aerosols and Photochemistry Experiment (BEARPEX) 2007 field campaign, which we review briefly in Sect. 2. In Sect. 3 we detail the various aspects of CAFE, including canopy structure, meteorology, vertical transport, emissions, deposition, advection and chemistry.Section 4 examines the sensitivity of modeled concentrations and fluxes to several key parameterizations that influence or control physical and chemical processes in the forest.We conclude in Sect. 5 with a characterization of the chemical environment within and immediately above the forest.This approach to model characterization is similar to several analogous model descriptions (Baldocchi, 1988;Gao et al., 1993).Detailed comparisons with BEARPEX-2007 data are presented in a companion paper (Wolfe et al., 2010).

BEARPEX-2007
Though adaptable to any forest, the CAFE model was originally designed to simulate a young Ponderosa pine plantation located in the western foothills of the Sierra Nevada Mountains, CA (38 • 58 42.9 N, 120 • 57 57.9 W, elevation 1315 m).This site is managed by Sierra Pacific Industries and is adjacent to the University of California's Blodgett Forest Research Station (BFRS) as described in detail elsewhere (Goldstein et al., 2000).Meteorological and chemical observations have been ongoing at BFRS since 1997 and have included two major collaborative field intensives, the Biosphere Effects on Aerosols and Photochemistry Experiment (BEARPEX) campaigns in August-October 2007 and June-July 2009.Below, we describe the configuration of CAFE for BEARPEX-2007.For evaluating key physical parameterizations in the model design, we restrict our simulations to mean mid-day (11:30-12:30 PST) conditions from 17 September 2007 (day of year 260).

Model description
The CAFE model domain consists of 86 layers in the vertical ranging from 0.1 m to 800 m.Layer spacing is constant at 0.1 m up to 1 m (the first 10 layers), after which layer heights are given by the exponential formulation of Gao et al. (1993), where z is the layer height, i is the layer index (11 to 86), a = −48.82and b = 11.41.This formulation yields a fineresolution grid of 36 layers within our (10 m tall) forest canopy and 50 layers within the atmospheric boundary layer (ABL).The ABL height of 800 m is typical of summer midday conditions at BFRS (Choi et al., 2010).Figure 1 summarizes the key processes included in the model.Within each layer, the 1-D time-dependent continuity equation is solved to determine the rate of change for all chemical species: Terms on the right respectively represent rates of chemical production, chemical loss, emission, deposition, advection (horizontal mixing) and vertical turbulent flux divergence.Below we outline considerations for modeling each of these processes in the context of BEARPEX-2007.

Canopy structure
The BFRS overstory is primarily Ponderosa pine, with a few interspersed White fir, Douglas fir, Incense cedar, Black oak and Sugar pine.A tree survey conducted in early October 2007 gave an average tree height of 7.9 ± 2.8 m, with 75% of the trees below 10 m.For the model, we choose an overstory height (h) of 10 m, a one-sided leaf area index (LAI os ) of 3.2 m 2 m −2 and a leaf area dry mass (d os ) of 219 g m −2 (Table 1).LAI os and d os were estimated from tree survey data via the allometric equations developed by Xu et al. (2001).Vertically-resolved canopy structure is described by the leaf area density function (LADF os ), which is parameterized following a modified Weibull distribution (Teske and Thistle, 2004).
The constants b = 0.6 and c = 3.6 are chosen to fit the parameterized LADF to results from destructive harvesting measurements (L.Misson, personal communication, 2008).The understory at BFRS is primarily Manzanita and Ceanothus shrubs and comprises a significant portion (∼40%) of  . Misson, personal communication, 2008).LADF us is parameterized with a parabolic shape to mimic a similar Manzanita stand described by Law et al. (2001): where z 1 is the height of the first layer (0.01 m). Figure 2a displays the modeled vertical LADF profiles, as well as the Ponderosa pine LADF as determined by destructive harvesting in 2003 and scaled to the 2007 canopy height and LAI.By definition, integration under these curves yields the total one-sided LAI.

Meteorology
Meteorological constraints, including pressure, temperature, water vapor concentration, radiation and friction velocity, are derived from a combination of measurements and parameterizations.Table 2 lists the mean noontime meteorological observations for day 260 during BEARPEX-2007.Atmospheric pressure is measured at a single height (12.5 m) and assumed to decay exponentially with increasing altitude.Measured air temperature profiles are fit using a spline interpolation below 12.5 m, above which an adiabatic profile is assumed (Fig. 2c).Water vapor concentration is set constant at the 12.5 m measurement.Above the canopy, total incoming solar radiation and photosynthetically-active radiation (PAR) are set constant to the 12.5 m measured values.
In-canopy radiation extinction is parameterized according to a modified Beer-Lambert law (Makar et al., 1999), where k rad is the radiation extinction coefficient, LAI cum (z) is the top-down cumulative LAI and SZA is the solar zenith angle as calculated by the Troposphere Ultraviolet and Visible (TUV) Radiation Model (http://cprm.acd.ucar.edu/Models/TUV/).The radiation extinction profile calculated from this method is similar to that derived from more complex parameterizations (Guenther et al., 1995;Gao et al., 1993) but is more flexible due to the tunable parameter k rad .
For k rad = 0.4, radiation at the lowest model level (0.01 m) is attenuated by a factor of 10 compared to above the canopy (Fig. 2b).We explore model sensitivity to k rad in Sect.4.3.
Friction velocity (u * ) is intimately related to the nearsurface diffusion profile.Since within-canopy micrometeorological observations are not available for BFRS, abovecanopy u* measurements are attenuated within the canopy by applying an LAI-dependent exponential decay (Yi, 2008): Figure 2d illustrates the in-canopy u* gradient, which closely tracks the LADF and has decreased by more than a factor of 10 at the lowest level.A typical measured horizontal wind speed profile is also shown, demonstrating a qualitatively similar decay.

Vertical transport
Turbulent mixing between layers is represented as a purely diffusive process using the flux-gradient or "K-theory" approach: K-theory is known to have limitations within dense forest canopies, where "near-field" effects of individual canopy elements and large-scale coherent processes like intermittent sweep-ejections influence turbulence structure (Arya, 1988;Harman and Finnigan, 2007).We neglect such intricacies in CAFE, as reproducing the full mixing structure in a forest requires a more refined numerical approach, such as largeeddy simulation (Patton et al., 2001).Flux-gradient relationships may also fail when reaction timescales are similar to or shorter than mixing timescales, leading to segregation of reactants that can reduce the effective rate of second-order chemical reactions (Dlugi et al., 2010;Pugh et al., 2010).Previous efforts to account for such effects include modifying eddy diffusivities to account for chemistry or implementation of higher-order closure schemes (Hamba, 1993;Gao et al., 1991).Such corrections can be computationally expensive in a model with thousands of reactions, like CAFE.Despite these drawbacks, K-theory continues to be a standard approach for this type of model (Gao et al., 1993;Makar et al., 1999;Stroud et al., 2005) and is likely sufficient for representing vertical diffusion in an average sense.
Figure 3 shows the vertical profile of the eddy diffusion coefficient, K(z).Above 12.5 m, K(z) is based on the values used by Gao et al. (1993), scaled to an ABL height of 800 m.Below 12.5 m, K(z) is given by Here, the vertical wind speed standard variance (σ 2 w ) and Lagrangian timescale (T L ) are defined as in Raupach (1989): The constant r represents a "near-field" correction, as derived in Makar et al. (1999) from the work of Raupach (1989), that accounts for the influence of canopy elements on eddy diffusivity: Here, τ (z) represents the "time since emission" for a theoretical diffusing cloud.In practice, this equation is tuned by choosing an appropriate value for τ /T L , typically ranging from 1 to 5 and constant throughout the canopy.For the current study, we choose τ /T L = 4.The effects of this choice on in-canopy diffusion are explored in Sect.4.1.

BVOC emissions
Biogenic VOC (BVOC) emissions are modeled in each canopy layer as a function of leaf density, light, temperature and vegetation type (overstory and understory).For each emitted compound and in each layer, the emission rate is calculated in units of molecules cm −3 s −1 as E b is the basal emission rate in molecules per gram of leaf per second and C L (z) and C T (z) are dimensionless correction factors for light and temperature (Guenther et al., 1995).The rightmost terms collectively represent the verticallydistributed leaf dry mass (grams of leaf per cm 3 air), determined by the leaf area dry mass (d, grams of leaf per cm 2  ground), the LADF (cm 2 ground per cm 3 air) and the 1sided LAI (cm 2 leaf per cm 2 ground).Modeled BVOC emissions include MBO, isoprene, methyl chavicol (MCHAV, also known as estragole), and a suite of speciated MT and SQT.
Table 3 lists all parameters used for emission calculations.Basal emission rates, defined at T = 30 • C and PAR = 1000 µmol m 2 s −1 , are largely based on values reported previously for this forest (Bouvier-Brown et al., 2009b, c;Harley et al., 1998;Schade et al., 2000) but are slightly adjusted in some cases to improve model agreement with BEARPEX-2007 observations (Wolfe et al., 2010).Basal emission rates for MBO include a 27% reduction for weaker emissions in 2 and 3-year-old needles (Schade et al., 2000).MBO emissions from Ponderosa pine are independent of soil moisture content (Gray et al., 2003), thus we do not include corrections for drought stress.
The temperature dependence for MCHAV, MT and SQT takes the exponential form with temperature coefficients (β) as reported by Bouvier-Brown et al. (2009c).Temperature corrections for MBO and isoprene follow functionally similar forms to one another (Guenther et al., 1995;Harley et al., 1998): where R is the ideal gas constant (8.314J mol −1 K −1 ) and all other coefficients are empirical constants (Table 3).
The light correction factor (C L ) is set to unity for MT and SQT, since their emissions are not strongly radiationdependent.Bouvier-Brown et al. (2009b) note that MCHAV mixing ratios at BFRS track well with those of MBO, and thus MCHAV emissions may be light-sensitive.As a strict characterization of the MCHAV light-dependence has yet to be performed, we set C L = 1 for MCHAV in the current study.MBO and isoprene share the same light correction factor (Guenther et al., 1995(Guenther et al., , 1999)): Figure 4a illustrates C L and C T calculated for MBO emissions.The temperature factor dominates the environmental correction, though light extinction becomes the limiting factor below ∼3 meters.At higher temperatures (not shown), C T for MBO approaches unity and C L becomes the controlling factor.(b) BVOC emission rates for the full canopy (overstory and understory) during the hot period, including MBO, isoprene, monoterpenes (MT), methyl chavicol (MCHAV) and sesquiterpenes (SQT).Note that bulk MT and SQT emissions are further partitioned into individual species (Table 4).
Figure 4b shows the calculated emission profile at BFRS for our meteorological conditions (Table 2 and Fig. 2).MBO is the dominant emission of Ponderosa pine, while understory emissions are primarily monoterpenes.By design, the canopy-integrated emission rates of isoprene are consistent with values used in larger-scale models for this region (Steiner et al., 2007).Local isoprene emissions at BFRS are likely smaller than our estimates here, as isoprene is not emitted from Ponderosa pine.Isoprene can, however, originate from other trees within and upwind of the plantation, particularly Black Oak.While likely compensating for an underestimate of isoprene advection from upwind sources (Dreyfus et al., 2002), modeling the isoprene source as an emission provides the best agreement with BFRS observations (Wolfe et al., 2010).Representing the isoprene source as strictly due to advection (Sect.3.7) results in a large overestimate of the first-generation isoprene oxidation products methyl vinyl ketone and methacrolein.In the emission parameterization, MT and SQT are given a bulk emission rate and then partitioned into individual species (Table 4) based on leaflevel cuvette measurements .(Bouvier-Brown et al., 2009c).Speciated MT include α-pinene, β-pinene, limonene, 3carene, myrcene, camphene, terpinolene, α-terpinene and γ -terpinene.SQT include α-bergamotene (ABERG), βcaryophyllene (BCARY), α-farnesene (AFARN) and unspeciated SQT (USQT).USQT are a proxy for the non-speciated SQT observed by Bouvier-Brown et al. (2009a, c).

Soil NO Emission
Emission of nitric oxide (NO) from soils can be an important source of oxidized nitrogen to the atmosphere, especially in rural and remote regions.NO emission rates are sensitive to soil and vegetation type, temperature, soil moisture and nitrogen deposition (Herman et al., 2003;Williams et al., 1992;Yienger and Levy, 1995).We model soil NO emissions as a function of temperature only, following the parameterization of Yienger and Levy (1995) for dry soil: NO is the basal emission flux (defined at 0 • C), z 1 is the height of the first layer and the function for calculating soil temperature (T soil ) from the air temperature in the lowest layer (T 1 ) is taken from the Williams et al. 1992) estimate for forests.A recent inventory of soil NO emission observations suggest a range for the basal NO emission flux of 9.88 +46.66  −8.17 ngN m −2 s −1 for dry-soil coniferous forests (Steinkamp and Lawrence, 2010).We choose a basal NO emission flux of 3 ngN m −2 s −1 , which gives a temperaturecorrected NO emission flux of 2.4 ngN m −2 s −1 .This value is at the low end of estimates derived from above-canopy NO 2 fluxes at BFRS (Farmer and Cohen, 2008) but is within the range of observations from California oak woodlands at lower elevations (Herman et al., 2003) and from recent soil chamber measurements at BFRS (E.Browne, personal communication, 2010).e Terpene and isoprene oxidation products (see Sect. 3.7), assumed to have the same deposition velocity as HNO 3 .f Adjusted to give an above-canopy deposition velocity at the aerodynamic limit (Paulot et al., 2009a).g Adjusted to give an above-canopy deposition velocity of ∼ 1.6 cm s −1 as in Hall and Claiborn (1997).h Adjusted to give an above-canopy deposition velocity of ∼ 2.7 cm s −1 as in Farmer and Cohen (2008).i Modified from Zhang et al. (2003) recommended values (see text).

Deposition
Leaf-level deposition is calculated for 35 species (Table 5) using a standard resistance parameterization (Wesely, 1989;Wesely and Hicks, 2000).In this approach, individual processes controlling the deposition of a molecule are assigned a characteristic resistance inversely related to the rate of that process.Summation of these resistances, in a manner analogous to the treatment of electrical circuits, yields the total deposition rate.All constants for deposition parameterizations are listed in Table 5.
Transfer of material between the canopy airspace and a leaf starts with molecular diffusion across a thin boundary layer, known as the laminar sublayer, which develops adjacent to the leaf surface.Our parameterization of the laminar sublayer resistance (R b ) is based on the derivation of Jensen andHummelshøj (1995, 1997), modified slightly to represent a resolved canopy rather than a "bulk canopy": Here, ν = 0.146 cm 2 s −1 is the kinematic viscosity of air, D is the species-dependent molecular diffusion coefficient, l w is the "aerodynamic leaf width," and c is a canopy-dependent tunable constant, set to 1 for the current study.Little information is available regarding the appropriate values of l w for different canopies, especially conifers.For the understory, we choose l w = 2 cm, which is the typical width of a shrub leaf at BFRS.For the overstory, we choose l w = 1 cm, which is the square root of the projected leaf area of a Ponderosa pine needle.The sensitivity of modeled fluxes to the choice of l w is explored in Sect.4.2.
G. M. Wolfe and J. A. Thornton: The (CAFE) Model -Part 1 Once past the laminar sublayer, a molecule may deposit through either stomatal or non-stomatal pathways.Resistance to passage through the leaf stomata (R s ) is primarily a function of in-canopy PAR, with additional corrections for temperature (f (T )) and water vapor pressure deficit (f(VPD)) as described in Zhang et al. (2003): In the above, R s,min is the minimum stomatal resistance, β PAR is a light correction coefficient, D H 2 O and D are the molecular diffusion coefficients for water and the molecule of interest, and the remaining variables are empirical coefficients that vary by canopy type (Table 5).Stomatal resistance is thus calculated for water vapor and subsequently scaled via molecular diffusion coefficients for other gases.Originally, all parameters for the above equations were adopted from the Zhang et al. (2003) recommendations for the land-use category "evergreen needleleaf trees."These values, however, do not adequately predict the stomatal resistance at BFRS when compared to independent canopy-scale calculations using the Penman-Monteith (Monteith and Unsworth, 1990) or leaf temperature (Thom, 1975) methods, both of which rely on observed above-canopy latent heat fluxes.In particular, the VPD dependence was found to be too strong with the default R s,min and b VPD values.Thus, we adopt the Zhang et al. (2003) values for β PAR , T min , T max and T opt , but we set R s,min = 7.4 s cm −1 and b VPD = 0.10 kPa −1 (Table 5) to optimize agreement with the observationally-constrained stomatal resistance calculations from the full BEARPEX-2007 dataset.The correction for leaf water potential stress from Zhang et al. (2003) is not included, as it was found to be negligible for the conditions of our study.Since f (T ) and f (VPD) were originally developed for a "big-leaf" model, we compute these using meteorological measurements at a single height (12.5 m) and apply them as constant corrections throughout the canopy.Thus, vertical variations in R s (z) are driven solely by PAR, while seasonal changes depend on both PAR and meteorology.
After passing through the stomata, molecules are assimilated by the mesophyll.Resistance to mesophyll uptake (R m ) is generally treated as a function of the effective Henry's law constant, H * , and a "reactivity" factor, f 0 , assigned for each depositing molecule (Table 5): In this formulation, H * and f 0 are meant to account for the effects of solubility and reactivity on uptake (Wesely, 1989).The sum, R s and R m defines the total resistance to stomatal uptake.
In addition to stomatal uptake, molecules can be lost via adhesion to -or reaction on -non-stomatal surfaces such as the waxy leaf cuticle.The cuticular resistance (R cut ) for individual molecules is calculated by assuming a characteristic cuticular resistance for ozone that varies by canopy type, R cut (O 3 ), and scaling this value by H * and f 0 (Wesely, 1989;Zhang et al., 2003): Figure 5a illustrates example sub-laminar, stomatal, mesophyll and cuticular resistances for O 3 deposition to the overstory.R b and R s increase in the canopy due to decreasing friction velocity and PAR, respectively.In general, the relative importance of stomatal and non-stomatal uptake will depend on the magnitude of R cut , which can vary by orders of magnitude between molecules.We caution that the choice of f 0 is somewhat arbitrary, as few flux or laboratory uptake measurements are available to validate modeled deposition rates for many species (Zhang et al., 2003(Zhang et al., , 2002b)).Constraints on R cut (O 3 ) are also tenuous, as they are typically derived from or validated against above-canopy flux observations that do not account for potential intra-canopy gasphase chemistry contributions to the net forest -atmosphere exchange (Zhang et al., 2002a, b).Such uncertainties may limit assessments of the relative contributions of stomatal uptake, non-stomatal deposition and chemistry to above-canopy fluxes.
The total resistance to deposition in each model layer, R dep (z), is the sum of the laminar sublayer and total surface resistances, where the latter consists of contributions from stomatal and nonstomatal uptake added in parallel: R dep (z) is calculated separately for the overstory and understory and scaled by LADF to give a first-order loss rate constant: Multiplication of k dep (z) by a concentration yields the firstorder loss rate due to deposition in each layer.Figure 5b illustrates deposition rates for several key species.HNO 3 and H 2 O 2 deposition are controlled primarily by R b , as their Henry's Law constants are set high to force R cut ∼ = 0.The effective Henry's law coefficient for H 2 O 2 is set much higher than typical values (∼10 5 M atm −1 ) to force H 2 O 2 to deposit at the "aerodynamic limit," meaning that its deposition is limited only by turbulent and molecular diffusion, in accordance with recent observations at BFRS (Paulot et al., 2009a).Formic acid (HC(O)OH) and acetic acid (CH 3 C(O)OH, not shown) also deposit quickly due to their relatively high solubility, as do the alkyl nitrates (ANs).AN deposition velocities are tuned to match the above-canopy value of 2.7 cm s −1 suggested in Farmer and Cohe (2008) by increasing H * to 1 × 10 8 M atm −1 , effectively lowering the cuticular resistance.A similar modification is made for the most abundant organic peroxides to yield an above-canopy deposition velocity of 1.6 cm s −1 (Hall and Claiborn, 1997).We note that setting f 0 to 1 for these molecules -following a recent study showing fast uptake of oxidized VOC (Karl et al., 2010) would not be sufficient to generate the required deposition velocity in this resistance framework.Ozone and acyl peroxy nitrates (APNs) deposit more slowly through both stomatal and non-stomatal pathways.The model includes deposition of eight APNs, six AN species and six ROOH species (Table 5), with a single k dep for each class.Most other species exhibit deposition rates smaller than O 3 , the exception being a set of terpene oxidation tracers (MTOX and SQTOX) and the isoprene epoxide IEPOX, which are assigned the same deposition velocity as HNO 3 (see Sect. 3.7).We do not include deposition of other oxidized VOC as suggested by Karl et al. (2010), partly due to a lack of observational constraints and partly because the pathway of metabolic consumption is likely limited by the stomatal conductance at this forest during the late summer.
Ground deposition includes resistances to aerodynamic transfer (R a0 ) between the first model layer and the ground and to actual loss at the surface (R g ): R a0 (20 s cm −1 ) and ground resistances for O 3 and SO 2 (both 2 s cm −1 ) are taken from Wesely (1989).Ground resistances for other species are scaled to O 3 and SO 2 according to H * and f 0 , akin to cuticular and mesophyll resistances.Ground deposition is treated separately from canopy deposition in the numerical integration scheme (see Sect. 3.9).

Chemistry
Chemistry in CAFE is based on the Master Chemical Mechanism (MCM) version 3.1 (http://mcm.leeds.ac.uk/MCM/).The MCM is a nearly-explicit reaction set that aims to track all oxidation processes and products throughout the photochemical degradation of primary VOC (Jenkin et al., 1997;Saunders et al., 2003).For the current study, a subset of the MCM is used that includes all reactions stemming from oxidation of MBO, isoprene, α-  grated) and that all species with lifetimes shorter than 0.01 s instantaneously proceed to products.The latter includes decomposition of alkoxy radicals (RO) and Criegee biradicals and is a standard simplification (Taraborrelli et al., 2009;Whitehouse et al., 2004).
(ii) Product branching for the reaction of methyl glyoxal (MGLYOX) with OH is updated following Baeza-Romero et al. (2007).Product branching for reaction of   the methacrolein-derived peroxy radical MACO 3 with NO is updated following Orlando et al. (1999).
(iii) An OH-production channel for reaction of HO 2 with peroxyacetyl radical (CH 3 CO 3 ) is added (Hasson et al., 2004).Branching ratios for the two pathways already included in the MCM are updated accordingly.Analogous changes are incorporated for reaction of HO 2 with peroxypropionyl (C 2 H 5 CO 3 ) and peroxymethacryloyl (MACO 3 ) radicals.
(iv) The oxidation of glycoaldehyde (HOCH 2 CHO) is updated following (Butkovskaya et al., 2006).First, product yields for the abstraction of the central hydrogen by OH are adjusted to include a 19% yield of OH and HCOOH.Next, a fast unimolecular decomposition is implemented for HOCH 2 CO 3 , the acyl peroxy radical formed after abstraction of the aldehydic hydrogen, with a rate constant of 10 s −1 (Butkovskaya et al. (2006) do not give a rate constant for this reaction but they note that it should be fast).A similar reaction is added for the IPRHOCO 3 radical, an analogous β-hydroxy acyl peroxy radical produced during reaction of IBUTALOH with OH.These changes are subject to verification with observations and, to our knowledge, these processes have not been independently verified.
(v) Alkyl nitrate chemistry for MVK and MACR is updated following Paulot et al. (2009b).This includes increasing the yield of HMVKANO 3 from 0.017 to 0.11 and adding the formation of a new compound, MACRNO 3 .We also include a reaction for the oxidation of MACRNO 3 by OH.
(vi) A set of 5 reactions are included for the formation and oxidation of isoprene dihydroxyepoxides (IEPOX) during reactions of OH with the four first-generation isoprene hydroxyhydroperoxides (ISOPAOOH, ISOPBOOH, ISOPCOOH and ISOP-DOOH) following Paulot et al. (2009c).Our IEPOX mechanism deviates from the original formulation in that reaction of IEPOX with OH produces two isoprenederived peroxy radicals already in the MCM, C57O 2 and C58O 2 .Though this is only strictly correct for IEPOX derived from ISOPBOOH and ISOPDOOH (Archibald et al., 2010), these two are the most dominant of the four MCM isomers.Moreover, the fate of IEPOX from ISOPAOOH and ISOPCOOH is not known (Paulot et al., 2009c).IEPOX is forced to deposit at the aerodynamic limit by assigning it a deposition rate equal to that of HNO 3 .
(vii) 36 reactions are added to model oxidation of MCHAV and all MT and SQT (listed in Sect.3.4) by OH, O 3 and NO 3 .This excludes α and β-pinene, which are already included in the MCM.Rate constants and reaction products are taken from the literature (Atkinson and Arey, 2003b;Bouvier-Brown et al., 2009c;Lee et al., 2006a, b;Hites and Turner, 2009).Unspeciated SQT are given the same reaction rate constants as βcaryophyllene.Though these reactions do not follow the explicit oxidation scheme of the MCM, we attempt to track oxidation products by including formation of the generic peroxy radicals MTO 2 and SQTO 2 .Oxidation of MCHAV also makes MTO 2 (though MCHAV is not strictly a monoterpene).These radicals react with NO, HO 2 and RO 2 to form the species MTOX and SQ-TOX, which represent closed-shell oxidized products.Reaction rate constants and product yields are estimated from the analogous reactions of β-pinene-derived RO 2 (Table 6).As further reaction or decomposition of MT and SQT oxidation products would also lead to the formation of smaller VOC (e.g.HCHO, CH3COCH3, etc.), we have included non-zero yields for these later generation products in the initial oxidation step using experimentally-determined values when available (Lee et al., 2006a, b).The compounds represented by MTOX and SQTOX are likely semi-volatile and thus should deposit to canopy surfaces or particles.For simplification, MTOX and SQTOX do not currently react further within CAFE but are allowed to deposit with a deposition rate constant equal to that of nitric acid.The behavior of these oxidation products is discussed further in Sect.5.1.
(viii) Rate constants are updated for reactions of OH with MPAN, PAN, MBO, IBUTALOH and the four isomers of ISOPOOH following IUPAC recommendations or literature-reported measurements.These are listed in Table 6.

Advection
Atmospheric composition at BFRS is influenced by upwind emissions, chemistry and mixing (Dillon et al., 2002).In particular, horizontal transport and mixing of localized plumes with regional background air can sustain or otherwise affect concentrations of species with large upwind emissions, such as NO x (= NO + NO 2 ) and isoprene (Dillon et al., 2002;Murphy et al., 2007).For simplicity, we refer to such processes collectively as advection.Advection is treated as a simple mixing process within each model layer (Dillon et al., 2002): The mixing rate constant is taken as k mix = 0.3 h −1 (Perez et al., 2009).Advection concentrations (C a ) are set constant throughout the model domain (no vertical profile) and are chosen such that final model concentrations represent typical conditions at BFRS for this period (Table 7); thus, in a sense, CAFE is constrained by measured chemical mixing ratios.Strictly speaking, it would be more appropriate to include a horizontal wind speed dependence in our advection terms; however, we do not have the necessary measurements to constrain such a parameterization.Advection is necessary to maintain reasonable concentrations for species that would otherwise build up to unreasonable levels or decay below measured values during integration.As parameterized here, advection has a negligible effect on modeled vertical exchange in CAFE, since the timescale for advective mixing is long (>3 h) compared to the timescales for turbulent transport and chemistry.

Numerical considerations
During model initialization, meteorological inputs are used to calculate diffusion parameters, emission rates, deposition velocities and chemical rate constants.These are held constant throughout a model run, as are advection concentrations.The latter are also used to initialize all chemical concentrations (Table 7).Integration is accomplished via operator splitting (Jacobson, 2005): Changes in concentration due to diffusion and chemistry are thus calculated separately and incrementally over each integration interval ( t = 2 s).For a single integration interval, each operator in Eq. ( 33) includes 20 operations, each with a time step of 0.05 s, giving a total of 40 chemistry and 40 diffusion operations for each 2 s interval.Intervals shorter than 2 s or time steps shorter than 0.05 s do not noticeably improve model performance.Diffusion (Eq.7) is solved numerically using a Crank-Nicolson scheme (Jacobson, 2005).At the upper boundary (the ABL) we assume a net-zero flux divergence (∂F ∂z = 0), corresponding to a constant concentration gradient through the last layer.This boundary condition implicitly assumes some level of entrainment across the ABL and was deemed better than simply fixing the concentrations at the top of the model, as it allows concentration and flux profiles to remain smooth.For most species, we assume no downward flux at the lower boundary (i.e.C 0 = C 1 at the ground).For species that are emitted or deposited at the ground, this boundary condition is replaced by a ground flux, calculated as the ground emission or deposition rate multiplied by the height of the lowest layer.The chemistry operator is executed with a simple forward Euler scheme.Canopy emission and deposition are represented in the chemistry operator as 0th-order production and 1st-order loss processes, respectively.Advection is also included in the chemistry operator.The entire model is written and executed in MAT-LAB.
The model is run for two hours, after which fluxes and exchange velocities are calculated from concentration profiles via Figure 6 depicts the time evolution of O 3 , PAN and MBO fluxes during a typical model run.A total integration time of two hours is sufficient for relaxation of flux and exchange velocity profiles.

Sensitivity to parameterizations
Several of the parameterizations implemented in CAFE are not fully constrained by observations or are highly simplified representations of complex processes.To better characterize sources of uncertainty, we have performed a series of sensitivity simulations that demonstrate the contributions of incanopy diffusion, laminar sublayer resistance and radiation extinction to modeled concentration and flux profiles.Model runs for these studies consisted of taking the base setup (described above) and modifying a single parameter: the diffusion timescale ratio (τ /T L ), the aerodynamic leaf width (l w ) or the radiation extinction coefficient (k rad ).The default values for these parameters are listed in Tables 1 and 5.In the following discussion, we focus primarily on concentrations and fluxes of HNO 3 , O 3 , MBO and the nitrogen oxides NO and NO 2 .Nitric acid and ozone are lost primarily via deposition to the canopy, though their controlling mechanisms are different.Turbulent and molecular diffusion limits HNO 3 deposition, while stomatal and cuticular resistances control O 3 uptake.Thus, O 3 and HNO 3 should respond differently to changes in diffusion, deposition and radiation parameterizations.MBO is the dominant emission of Ponderosa pine and its flux and profile will be sensitive to diffusion and radiation parameterizations.The NO/NO 2 ratio will be sensitive to the radiation profile, with implications for hydrogen oxide radical partitioning and the formation and loss of higher oxides of nitrogen.Our goals for these sensitivity tests are to demonstrate the effects of our choices on model output and to provide a more thorough, but still general, picture of the mechanics that underlie forest-atmosphere exchange.

In-canopy turbulent diffusion
As noted in Sect.3.3, turbulent transport within forest canopies is difficult to accurately simulate.In our K-theory approach, we compensate for "near field" effects of canopy elements by inclusion of a correction factor, r, that is linked to the canopy structure through our choice of τ /T L .For τ /T L ≤ 1, r is undefined.Increasing τ /T L leads to a nonlinear increase in r, corresponding to faster diffusion and shorter residence times within the canopy.For high enough values of τ /T L (> 5), r is approximately unity and K is defined by the "far field" limit (Makar et al., 1999).
Figure 7 shows results for model runs with τ /T L values of 1.1, 1.5, 2 and 4.These correspond to r = 0.074,0.45,0.71and 0.97.At the lowest value, diffusion is slow enough that substantial gradients develop in depositing and emitting species.HNO 3 mixing ratios decrease within the canopy to ∼20% of the canopy-top value.O 3 does not deposit as readily, thus the O 3 gradient is less steep, reaching 75% of the canopy-top value.In contrast, MBO pools near the base of the canopy, increasing by more than a factor of 5 relative to the base case.Decreasing in-canopy K-values also reduce the magnitude of above-canopy fluxes, though this effect is somewhat offset by steeper concentration gradients (see Eq. 34).As noted above, the effects of τ /T L are nonlinear, and there is very little difference in mixing ratio profiles between τ /T L = 2 and 4.
For BFRS, we find τ /T L = 4 (giving r = 0.97) to be the best choice, based on several metrics.First, this value provides the best agreement between modeled and measured vertical concentration profiles of strongly emitting and depositing species such as BVOC and ozone (Wolfe et al., 2010).Second, modeled above-canopy K-values are consistent with those calculated independently from fluxes and gradients of carbon dioxide (CO 2 ), water vapor, sensible heat and momentum obtained during BEARPEX-2007.Table 8 shows our derivation of K-values from these observations, which is accomplished by rearranging Eq. ( 34) and substituting in the observed fluxes (at 12.5 m) and vertical gradients (from observations at 12.5 m and 8.75 m).We broaden the data selection window to include all noontime observations from the BEARPEX-2007 dataset for this analysis, as calculated K-values can be quite variable for any single point in time.Observationally-derived average K-values range from 2.4-3.9 m 2 s −1 , which bracket the canopy-top model value of 2.8 m 2 s −1 .Within the canopy -where the uncertainty in K is largest -data is not available for such a calculation.Instead, we estimate an average canopy residence time by integrating over the "turbulent aerodynamic resistance" within each layer, in analogy to the resistance formulation for deposition parameterizations (e.g.Baldocchi, 1988): For our conditions, this calculation yields τ can = 134 s, within the range of observation-based estimates for BFRS (Farmer and Cohen, 2008;Holzinger et al., 2005).For comparison, values of τ /T L = 1.1,1.5 and 2 give canopy residence times of 1700, 290 and 182 s.

Laminar sublayer resistance
The laminar sublayer is a thin region of laminar flow that develops above individual canopy surfaces in response to surface wind drag.This layer is a barrier to deposition, as molecules in the canopy airspace must diffuse across it before uptake can occur.The thickness of the laminar sublayer is related to the size and shape of canopy elements (e.g.leaves and needles) and is represented in our resistance parameterization by the aerodynamic leaf width, l w .According to Eq. ( 21), R b should scale with the square root of l w .The impact of changing l w on fluxes will be strongest for species whose deposition is not limited by stomatal or cuticular resistances, such as HNO 3 and H 2 O 2 .Sievering et al. (2001) note that HNO 3 deposition can be much faster over coniferous canopies than over broad-leaf canopies due to the smaller R b in the former.Figure 8 displays modeled HNO 3 and O 3 exchange velocities for a range of l w values.In these model runs, the base l w for both the overstory (1 cm) and understory (2 cm) were multiplied by factors of 0.1, 0.5 and 2. At the lower limit, HNO 3 exchange velocities are quite fast (−7 to −10 cm s −1 ) and are essentially limited by turbulent diffusion and the available leaf area (Eq.29).Exchange velocities follow the expected trend with increasing l w .The effects are dampened in the case of O 3 fluxes as their uptake is limited by stomatal and cuticular resistances.Default values of l w give an abovecanopy HNO 3 exchange velocity of −3.5 cm s −1 .This value is consistent with the results of Farmer and Cohen (2008), who scaled wintertime HNO 3 flux measurements with friction velocity to derive a noontime HNO 3 deposition-only exchange velocity of −3.4 cm s −1 for BFRS in August 2005.

Radiation extinction
Though solar radiation is heterogeneously distributed in a forest canopy, we represent the average vertical profile by a Beer's Law parameterization (Eq.5).Without in-canopy radiation measurements, the radiation extinction coefficient (k rad ) must be estimated from measurements at other forests; typical values of k rad range from 0.4-0.65 for conifers and understory shrubs (Law and Waring, 1994;Runyon et al., 1994).As solar radiation can influence stomatal uptake, BVOC emissions and photochemistry, it is prudent to examine model sensitivity to k rad .Figure 9 illustrates modeled concentration and flux profiles of O 3 and MBO for k rad ranging from 0 to 0.6.For both species, the effects of increasing k rad on both mixing ratios and fluxes are fairly linear.As light extinction increases, O 3 concentrations increase slightly (by ∼2%) due to higher stomatal resistance in a darker canopy.This effect is more apparent in O 3 fluxes, which decrease by ∼30% on going from low to high k rad .Ozone mixing ratio trends may also include a photochemical effect from shifts in the NO-NO 2 -O 3 equilibrium (see below), but such behavior does not significantly influence O 3 fluxes.MBO displays a much stronger dependence on radiation, stemming from its emission parameterization (Eq.16).MBO mixing ratios are ∼50% lower for our base case (k rad = 0.4) compared to the case with no extinction (k rad = 0), and the gradient is slightly steeper near the ground.Within CAFE, changes in radiation will influence MBO  will not affect those of MT, SQT and MCHAV, though there is some evidence that the emissions of the latter three are also sensitive to photosynthetically-active radiation in certain ecosystems (Bouvier-Brown et al., 2009b;Fowler et al., 2009;Geron and Arnts, 2010).
In addition to affecting surface processes like deposition and emission, radiation extinction also influences in-canopy photochemistry, as demonstrated by the partitioning of nitrogen oxides (NO x = NO + NO 2 ).Chemistry generally controls the relative balance of NO and NO 2 , which interconvert rapidly via reaction of NO and O 3 to form NO 2 , followed by photolysis of NO 2 in the presence of molecular oxygen to reform NO and O 3 .
Figure 10 shows the sensitivity of modeled NO/NO 2 profiles to changes in k rad .In these model runs, total NO x changes by <2%.The equilibrium defined by reactions (Eqs.37-39) is light-dependent; thus, NO/NO 2 shifts towards lower values as radiation extinction increases and less NO 2 is photolyzed.Near the ground, soil NO emissions enhance NO/NO 2 , but the effect is limited to z/h < 0.2 for even moderate values of k rad .The influence of in-canopy radiation extinction persists above the canopy, with NO/NO 2 ratios for the various runs converging only above 30 m (z/h = 3).Variations in NO/NO 2 will feed back into other chemical processes, including those controlling the cycling of hydrogen oxide radicals and the fate of acyl peroxy nitrates; thus, it is important to have an accurate estimate of radiation in the canopy airspace.Our choice of k rad = 0.4 gives near-surface NO/NO 2 ratios of 0.3, in good agreement with recent measurements at BFRS (E.Browne, personal communication, 2010).Furthermore, modeled gradients and fluxes of MBO with k rad = 0.4 are consistent with previous observations at this site (Baker et al., 1999;Holzinger et al., 2005;Schade et al., 2000).

Mixing ratio profiles
The most unique aspect of the CAFE model is the incorporation of the extensive MCM reaction scheme.CAFE is also the first 1-D canopy model to be implemented for BFRS despite many years of chemical measurements at this site.To characterize the chemical environment in this forest, we present a brief overview of modeled near-surface mixing ratio profiles of biogenic VOC (BVOC), hydrogen oxides (RO x ) and reactive nitrogen (NO y ).All results presented below are taken from the base run, and all profiles have been normalized to their canopy top (z/h = 1) values.
In a companion paper, we explore concentrations and fluxes in more detail through an extensive comparison with results from BEARPEX-2007.

BVOC
Oxidation of biogenic VOC is integral to the formation of tropospheric ozone, organic nitrates and secondary organic aerosol (SOA) and can potentially influence forestatmosphere exchange of O 3 (Kurpius and Goldstein, 2003), acyl peroxy nitrates (Wolfe et al., 2009)  Figure 11a shows modeled profiles of MBO and isoprene, the two dominant BVOC emissions in CAFE.MBO and isoprene exhibit the same gradient, with ∼20% higher concentrations at the ground than at the top of the canopy.Averaged over the full 800 meter model domain, MBO and isoprene mixing ratios are only ∼25% of those at canopy top, illustrating that surface layer concentrations may not be entirely representative of the overlying regional mixed layer.The major 1st-generation oxidation products of MBO are 2-methyl-2-hydroxypropanal (IBUTALOH) and glycoaldehyde (HOCH 2 CHO), while those of isoprene are methyl vinyl ketone (MVK) and methacrolein (MACR); Appendix A lists the structures of these compounds.Glycoaldehyde is the product of a number of VOC reactions, but MBO is by far the dominant source for our model conditions.Incanopy gradients in these oxidation products are <5% of the canopy top concentrations (Fig. 11b), much smaller than for their parent VOC.This reflects the relative timescales of chemistry and turbulent diffusion; that is, vertical mixing of these species is faster than their chemical production and/or loss.IBUTALOH demonstrates a slightly steeper gradient than the other three because its production rate is slightly faster, thus it is more sensitive to the gradient of MBO.Little is known about the deposition of oxidized VOC, though recent work suggests that it may be faster than currently assumed in regional and global models (Karl et al., 2010).The shape of these profiles might be different than our model suggests if these molecules are taken up by the canopy.Particularly, we might expect that the hydroxyl functionality on IBUTALOH and glycoaldehyde would enhance the deposition rate of these molecules relative to that of MVK and MACR, leading to variations in their respective gradients.Furthermore, in-canopy gradients of isoprene, MVK and MACR at BFRS could be less pronounced than predicted by CAFE given the primary isoprene source is believed to be advection and not emission.
Figure 12 displays a similar set of predictions for terpenes.In this plot, the composite MT' includes MCHAV and all monoterpenes other than α and β-pinene (Sect.3.4).All monoterpenes (α-pinene, β-pinene and MT') exhibit similar profiles, with enhancements of 40-70% at the ground relative to the canopy top, consistent with early measurements of monoterpene profiles at BFRS (Holzinger et al., 2005).We expect this enhancement to be stronger for terpenes than for MBO and isoprene because the understory is a major terpene source (Fig. 4).The dominant species in the MT' group are 3-carene (24% of total mixing ratio), MCHAV (20%) and camphene (19%).Boundary layer-averaged monoterpene concentrations are generally 20% of their canopy top values, consistent with their high reactivity.Sesquiterpene (SQT) gradients are even stronger, with a of 2 enhancement at the ground and domain-averaged concentrations that are <10% of the canopy-top mixing ratio.Dominant resolved SQT include α-bergamotene (44% of total mixing ratio) and α-farnescene (Table 6), relative terpene mixing ratios will not necessarily reflect their relative emission rates (Table 4).Each of the four terpene classes shown in Fig. 12a has a distinct oxidation product in our chemical mechanism.α and β-pinene are the precursors to pinonaldehyde (PINAL) and nopinone, respectively, while oxidation of other monoterpenes, MCHAV and sesquiterpenes yields MTOX or SQ-TOX (see Sect. 3.7 and Table 6).PINAL and nopinone display in-canopy gradients of ∼5% (Fig. 12b), similar to the 1st-generation products of MBO and isoprene.These compounds are semivolatile (Kavouras et al., 1999), and partitioning to particles or deposition to surfaces -neither of which are represented in the current version of CAFE -may change their profiles near the surface.The MTOX profile represents an upper limit for such effects, as we choose the deposition velocity for MTOX to be near the aerodynamic limit.MTOX decreases by ∼10% between the top of the canopy and the ground.SQTOX does not decay as much in the canopy, despite having the same deposition rate; indeed, there even appears to be a slight bulge in the SQTOX profile at z/h = 1.5.The primary contributors to SQTOX formation (through SQTO 2 ) are reactions of O 3 with β-caryophyllene and unspeciated sesquiterpenes, the rates of which are more than double those of other sesquiterpene oxidation reactions and more than 10 times faster than monoterpene oxidation.Combined with the enhanced SQT profile near the ground, this leads to a suppressed depositional decay in the canopy and a slight bulge over the canopy.Above z/h = 1.5, turbulent mixing and a lack of precursors causes the SQTOX profile to again decrease.
Previous observations at BFRS have suggested the presence of BVOC oxidation products that maximize above the canopy (Holzinger et al., 2005).In that study, the profile of unidentified oxidation products showed a pronounced increase of a factor of 2 between the ground and directly above the canopy.Holzinger et al. (2005) note that the observed oxidation product mixing ratios require emissions of highly reactive BVOC at a rate 6-30 times higher than the total monoterpene emission rate.SQT were not considered by Holzinger et al. (2005) but are a significant contributor to total terpene reactivity.In CAFE, the total SQT emission rate is 33% of the total MT emission rate on a per-molecule basis, while the reactivity with ozone, averaged over the full model domain, is 10 times higher for SQT than MT.Reactivities with OH are similar for both classes.It is intriguing that our modeled SQTOX profile is qualitatively consistent with the results of Holzinger et al. (2005); however, it appears that we would require significant additional emissions of highly reactive BVOC to fully reproduce their results.

Hydrogen oxides
As evident from Eqs. (40-43), the hydrogen oxide and organic peroxy radicals OH, HO 2 and RO 2 are central to photochemical cycling and hydrocarbon degradation in the troposphere.The lifetimes of these radicals are fairly short (<1 s for OH and <100 s for peroxy radicals), and the primary sources in CAFE are through photolysis of O 3 and HCHO: Thus, we might expect their concentrations to show significant in-canopy gradients.Figure 13 displays the normalized mixing ratios for all three radicals (RO 2 represents the sum of 154 speciated organic peroxy radicals).OH decreases by 38% between canopy top and the ground due to decaying sources (O 3 photolysis and reaction of HO 2 with NO) and faster losses via reactions with VOC.OH continues to increase above the canopy as VOC mixing ratios continue to decrease (Fig. 11) and HO 2 and NO increase.Our OH gradient is consistent with other model estimates (Gao et al., 1993;Stroud et al., 2005), though these results should be viewed with caution.Growing evidence suggests that OH chemistry in high-BVOC environments is not well understood (Hofzumahaus et al., 2009;Lelieveld et al., 2008;Thornton et al., 2002) and thus is likely not accurately represented by the current MCM.Indeed, elucidation of this chemistry was a key scientific target of the BEARPEX-2007 campaign.Modeled HO 2 and RO 2 profiles are more vertical, with increasing HO 2 and decreasing RO 2 as a function of height.Gradients of these two species mirror one another above z/h = 0.2.This behavior is linked to the profile of NO (Fig. 14), which drives the conversion of RO 2 to HO 2 via Eqs.(41-42) and is lower in the canopy than above.The cycling of OH, HO 2 and RO 2 via Eqs.(40-43) is relatively fast compared to the source and sinks of these molecules, thus it is often convenient to define the sum of these as the chemical family RO x .RO x (also shown in Fig. 13) is mostly conserved throughout the surface layer, demonstrating that gradients in RO x components are driven by their interchange through radical-propagating reactions.This result might change for a denser forest than BFRS, where substantial light attenuation could significantly limit radical production and NO-NO 2 cycling.

Reactive nitrogen
The reactive nitrogen family (NO y ) encompasses a wide spectrum of atmospheric oxidized nitrogen compounds, including NO x (= NO + NO 2 ), acyl peroxy nitrates (APNs), alkyl nitrates (ANs) and nitric acid (HNO 3 ).NO x enters the atmosphere via anthropogenic (e.g.combustion) or biogenic (e.g.soil) emissions, while the higher oxides of nitrogen are formed via reactions of NO x and RO x : The quantity and partitioning of NO y in the troposphere influences O 3 production, pollution transport and deposition of bio-available nitrogen to ecosystems.Near the surface, NO y gradients are determined by both deposition and photochemistry, with potentially important consequences for net forestatmosphere exchange (Farmer and Cohen, 2008;Wolfe et al., 2009;Horii et al., 2006).Figure 14 illustrates the modeled vertical profiles for all major forms of NO y .NO and NO 2 follow the behavior discussed in Sect.4.3, with profiles controlled by NO emissions and the radiation-dependent NO x -O 3 equilibrium.Total APNs, composed primarily (77%) of peroxyacetyl nitrate (PAN), do not exhibit a strong concentration gradient under our model conditions.APN deposition is not particularly fast (Fig. 5), thus chemical production and loss can also influence the APN gradient.APN chemistry is a function of both temperature (through the thermal equilibrium, Eq. 47) and the NO/NO 2 ratio, both of which can exhibit marked gradients in a forest canopy.In a detailed analysis of APN flux observations collected during BEARPEX-2007, Wolfe et al. (2009) demonstrated that the considerable increase in air temperature near the ground at BFRS (Fig. 2c) can enhance APN losses, which adds a negative (downward) contribution to the net above-canopy flux (or, equivalently, a positive contribution to the in-canopy concentration gradient).Observations at BFRS have also suggested that in-canopy chemical production can alter APN gradients and fluxes.From observations of positive sum peroxy nitrate ( PN) fluxes in summer 2005, Farmer and Cohen (2008) deduced a PN concentration gradient of −2% between 7.3 m and 14.3 m (z/h = 1 and 2).This inferred slope is opposite in sign but similar in magnitude to our modeled APN gradient, which we find to be driven by a combination of deposition and enhanced thermal decomposition near the ground.Comparison to measured APN gradients suggests intra-canopy losses are under-estimated (Wolfe et al., 2010).
Total ANs and HNO 3 demonstrate profiles characteristic of strong deposition, with in-canopy gradients of 6% and 12%, respectively.Little information is available on AN deposition, thus we set their Henry's Law coefficient ( to give an above-canopy deposition velocity of ∼2.5 cm s −1 , inline with AN flux measurements at BFRS (Farmer and Cohen, 2008).Our modeled AN gradient agrees with that predicted by Farmer and Cohen (2008), while our modeled HNO 3 gradient is somewhat steeper than their "depositiononly" gradient despite having similar deposition velocities.
Our HNO 3 profile is also within the range of above-canopy HNO 3 gradient observations from a similar forest (Sievering et al., 2001).Total NO y is mostly conserved throughout the vertical.The NO y gradient is most similar to that of APNs, which comprise 42% of total NO y averaged over 0-30 m, followed by NO 2 (24%), HNO 3 (19%), NO (9%) and AN (6%).Other components (e.g.HONO, NO 3 , N 2 O 5 , ClNO 2 , etc.) are negligibly small for the current study, either because photochemistry prevents significant concentrations or, in the case of HONO, because we neglect sources for which quantitative theoretical descriptions are lacking (Kleffmann et al., 2005;Zhang et al., 2009;Ren et al., 2010).The NO y gradient is not affected by interchange among individual nitrogencontaining species, thus its shape is determined by deposition and soil NO emissions (note the negative slope near the ground).Our modeled NO y profile looks similar -both in shape and magnitude -to the model results of Gao et al. (1993), despite the fact that the latter study was focused on a different forest and used a fairly condensed chemical mechanism.NO, NO 2 and NO x profiles also compare well between our results and those of Gao et al. (1993).Our modeled NO x /NO y ratio of 0.33 (averaged over 0-30 m, profile not shown) is ∼50% of that reported by Gao et al. (1993), which may reflect greater production of organic nitrates (APNs and ANs) in our detailed chemical mechanism.Gao et al. (1993) do not report the details of their modeled NO y speciation.

Conclusions
We have introduced and assessed the initial performance of the CAFE model, a vertically-resolved 1-D chemical transport model designed to examine forest-atmosphere exchange of reactive gases, by simulating a young Ponderosa pine forest with meteorological constraints from the BEARPEX-2007 field campaign.We characterized the sensitivity of modeled profiles to key parameterizations: diffusion, laminar sublayer resistance and radiation extinction.Slowing incanopy diffusion leads to steep gradients in depositing and emitting species, leading to a net decreases in the magnitude of deposition-driven modeled fluxes.Changes to the laminar sublayer resistance produce expected variations in exchange velocities of depositing species, with the strongest effect on diffusion-limited HNO 3 deposition.In-canopy radiation decay affects both stomatal uptake and light-dependent BVOC emissions, though the effects are relatively minor over a realistic range of radiation extinction coefficients.In addition, canopy radiation extinction has a substantial effect on the NO/NO 2 ratio, the influence of which persists up to ∼3 times the canopy height.Generally, these sensitivity tests confirm that our chosen parameterizations are reasonable in that they are within ranges inferred from observations made previously at the site or represent a midpoint between possible extreme values where constraints are lacking.They also provide some insight into how fluxes of reactive species might vary between forests having different canopy architectures and/or meteorology.
To illustrate the chemical environment in this forest, we examined modeled mixing ratio profiles of biogenic VOC, RO x and NO y within and immediately above the canopy.Directly emitted VOC demonstrate decaying profiles characteristic of turbulent mixing, though chemical loss also affects the gradients of more reactive VOC such as sesquiterpenes.First-generation oxidation products of these VOC are more evenly distributed in the vertical but can show interesting structure, such as an above-canopy local maximum in the sesquiterpene oxidation tracer SQTOX.Photochemistry, VOC abundance and radical-propagating cycles control the gradients of individual RO x components, while total RO x is conserved throughout the surface layer.NO x partitioning reflects the radiation-dependent photostationary state between NO, NO 2 and O 3 , which is only influenced by soil NO emissions in the lower canopy.The APN gradient is fairly weak but reflects net in-canopy sinks, consistent with BEARPEX-2007 observations.Gradients in AN and HNO 3 are primarily deposition-driven, while the vertical profile of total NO y is determined by deposition and, near the ground, soil NO emissions.Quantitatively, these results are likely specific to BFRS, though similar relative gradients should be expected for analogous environments, e.g.Ponderosa pine forests in Mediterranean climates.
Detailed canopy modeling is a useful tool for investigating the physics and chemistry underlying forest-atmosphere exchange.Though CAFE is not currently viable as a fully prognostic model, it can be used to examine the present state of knowledge on reactive gas exchange and related hypotheses, thus allowing us to identify potential avenues for future research.In a companion paper (Wolfe et al., 2010), we utilize CAFE in conjunction with the comprehensive chemical dataset obtained during BEARPEX-2007 to evaluate model processes and chemical contributions to forest-atmosphere exchange of reactive gases at this site.Though CAFE is currently optimized to simulate Blodgett Forest, we note that it could be adapted to another forest with sufficient information about canopy structure, meteorology and chemical concentrations.
Fig. 2. (a) Modeled leaf area distribution functions (LADF) for overstory (solid line) and understory (dashed line).Red circles represent the LADF for Ponderosa pine as measured by destructive harvesting and scaled to the BEARPEX 2007 canopy height.(b) Modeled in-canopy radiation extinction ratio, shown as a fraction of above-canopy incoming radiation.(c) Measured (red circles) and modeled (blue lines) near-surface temperature profile.(d) Measured above-canopy friction velocity (red circle) and modeled profile (blue line).The measured wind speed profile for the same period is also shown (black asterisks).
Fig. 4. (a) Environmental adjustment factors for MBO emissions: radiation (C L ), temperature (C T ) and the product of these two.(b)BVOC emission rates for the full canopy (overstory and understory) during the hot period, including MBO, isoprene, monoterpenes (MT), methyl chavicol (MCHAV) and sesquiterpenes (SQT).Note that bulk MT and SQT emissions are further partitioned into individual species (Table4).

Fig. 6 .
Fig. 6.Evolution of modeled fluxes for (a) ozone, (b) PAN and (c) MBO over the course of a model run.Profiles have been normalized to their respective canopy-top (z/h = 1) values at the end of the model run (7200 s).Arrows denote the general progression of profiles as a function of model runtime.

Fig. 7 .
Fig. 7. Vertical profiles of mixing ratios and fluxes of nitric acid (ab), ozone (c-d) and MBO (e-f) for the diffusion sensitivity runs.In each plot, curves are normalized by the canopy top values from the base run (τ /T L = 4).Profiles correspond to runs with τ /T L values of 1.1, 1.5, 2 and 4.

Fig. 8 .
Fig. 8. Exchange velocity profiles of (a) nitric acid and (b) ozone for the laminar sublayer sensitivity runs.Profiles correspond to runs with the base l w values multiplied by 0.1, 0.5, 1 and 2.

Fig. 9 .
Fig. 9. Vertical profiles of mixing ratios and fluxes of ozone (ab) and MBO (c-d) for the radiation sensitivity runs.In each plot, curves are normalized by the canopy top values from the base run (k rad = 0.4).Profiles correspond to runs with k rad values of 0, 0.2, 0.4 and 0.6.

Fig. 10 .
Fig.10.Vertical profiles the NO/NO 2 ratio for the radiation sensitivity runs.Profiles correspond to runs with k rad values of 0, 0.2, 0.4 and 0.6.

Fig. 11 .
Fig. 11.Mixing ratio vertical profiles for (a) MBO and isoprene, and (b) IBUTALOH, HOCH 2 CHO, MVK and MACR.Each profile is normalized to its value at the top of the canopy (z/h = 1).
Figure12displays a similar set of predictions for terpenes.In this plot, the composite MT' includes MCHAV and all monoterpenes other than α and β-pinene (Sect.3.4).All monoterpenes (α-pinene, β-pinene and MT') exhibit similar profiles, with enhancements of 40-70% at the ground relative to the canopy top, consistent with early measurements of monoterpene profiles at BFRS(Holzinger et al., 2005).We expect this enhancement to be stronger for terpenes than for MBO and isoprene because the understory is a major terpene source (Fig.4).The dominant species in the MT' group are 3-carene (24% of total mixing ratio), MCHAV (20%) and camphene (19%).Boundary layer-averaged monoterpene concentrations are generally 20% of their canopy top values, consistent with their high reactivity.Sesquiterpene (SQT) gradients are even stronger, with a of 2 enhancement at the ground and domain-averaged concentrations that are <10% of the canopy-top mixing ratio.Dominant resolved SQT include α-bergamotene (44% of total mixing ratio) and α-farnescene (39%).Due to different oxidation lifetimes

Fig. 13 .Fig. 14 .
Fig. 13.Mixing ratio vertical profiles for OH, HO 2 , RO 2 and total RO x .Each profile is normalized to its value at the top of the canopy (z/h = 1).
the total leaf area as of 2007.Understory height (h us ) is estimated at 2 m, while LAI us (1.9 m 2 m −2 , one-sided) and d us (377 g m −2 ) are approximated by scaling 2003 destructive harvesting measurements to the 2007 understory height (L

Table 2 .
Mean noontime meteorological parameters for day 260.
b Measured at 12.5 m.c From TUV model.d Divide by 2.92 to convert to W m −2 .

Table 6 .
Additional and modified reactions.

Table 7 .
Initial/advection concentrations.Species not listed have initial/advection concentrations set to 0.
d Gradient in air temperature.e Gradient in horizontal wind speed.f Model value at 12 m.