FORest Canopy Atmosphere Transfer (FORCAsT) 2.0: model updates and evaluation with observations at a mixed forest site

The FORCAsT (FORest Canopy Atmosphere Transfer) model version 1.0 is updated to FORCAsT 2.0 by implementing five major changes, including (1) a change to the operator splitting, separating chemistry from emission and dry deposition, which reduces the run time of the gas-phase chemistry by 70 % and produces a more realistic in-canopy profile for isoprene; (2) a modification of the eddy diffusivity parameterization to produce greater and more realistic vertical mixing in the boundary layer, which ameliorates the unrealistic simulated end-of-day peaks in isoprene under well-mixed conditions and improves daytime air temperature; (3) updates to dry deposition velocities with available measurements; (4) implementation of the Reduced Caltech Isoprene Mechanism (RCIM) to reflect the current knowledge of isoprene oxidation; and (5) extension of the aerosol module to include isoprene-derived secondary organic aerosol (iSOA) formation. Along with the operator splitting, modified vertical mixing, and dry deposition, RCIM improves the estimation of first-generation isoprene oxidation products (methyl vinyl ketone and methacrolein) and some second-generation products (such as isoprene epoxydiols). Inclusion of isoprene in the aerosol module in FORCAsT 2.0 leads to a 7 % mass yield of iSOA. The most important iSOA precursors are IEPOX and tetrafunctionals, which together account for >86 % of total iSOA. The iSOA formed from organic nitrates is more important in the canopy, accounting for 11 % of the total iSOA. The tetrafunctionals compose up to 23 % of the total iSOA formation, highlighting the importance of the fate (i.e., dry deposition and gas-phase chemistry) of later-generation isoprene oxidation products in estimating iSOA formation.

1 Introduction

Forests cover 30 % of the land surface and play an important role in the Earth system through exchanges of energy, water, carbon dioxide, and reactive chemical species with the atmosphere (Bonan, 2008). Forest canopies emit large amounts of volatile organic compounds (VOCs) into the atmosphere (Guenther et al., 2006) that drive atmospheric chemistry (e.g., Chameides et al., 1992; Taraborrelli et al., 2012) and are precursors to climate-relevant species such as ozone (e.g., Wolfe et al., 2011) and particulate matter (e.g., Palm et al., 2016). In addition, forest canopies serve as a major sink of VOCs through dry deposition (e.g., Nguyen et al., 2015). These bidirectional exchanges and their influences on atmospheric chemistry are complicated by the three-dimensional structure of forest canopies, which creates turbulent flows significantly different from the overlying atmospheric boundary layer (e.g., Gao et al., 1993; Patton et al., 2001) and affects the vertical transport and the chemistry of trace gases (e.g., Kaser et al., 2015).

The complexity and interplay of these chemical and physical processes challenge our understanding of forest-driven climate impacts on local, regional, and global scales (e.g., Spracklen et al., 2008; Vanwalleghem and Meentemeyer, 2009). Improving our understanding of the chemical and physical processes governing forest–atmosphere interactions at a local scale is helpful to generalize the net impact of the terrestrial biosphere on chemistry and climate at broader scales. Canopy–chemistry models that explicitly represent the physical, chemical, and biological processes of an individual forest canopy are useful tools to investigate the chemically relevant interactions between forests and the atmosphere at local scales. These canopy–chemistry models calculate the environmental variables that drive emissions, dry deposition, turbulent mixing, and chemical reactions vertically throughout a canopy at very fine resolutions (e.g., on the order of meters), while atmospheric chemical transport models approximate canopy processes through parameterizations and operate on the 10–200 km scale. Several forest canopy–chemistry models (Stroud, 2005; Forkel et al., 2006; Boy et al., 2011; Wolfe and Thornton, 2011) have been developed to study chemically relevant forest–atmosphere exchanges with a focus on the gas-phase chemical processes. The FORCAsT model version 1.0 (FORest Canopy Atmosphere Transfer) (Ashworth et al., 2015) is one of the few canopy models currently capable of simulating the formation of secondary organic aerosol (SOA) from biogenic VOC oxidation.

Over the past decade, the understanding of isoprene chemistry under a wide range of NO_x conditions and their impact on atmospheric particles has greatly expanded. Specifically, understanding of isoprene oxidation under low-NO_x conditions has improved (Wennberg et al., 2018), and proper representation of isoprene oxidation and isoprene-derived SOA formation in canopy–chemistry models is now recognized to be important for a more accurate understanding of forest–atmosphere exchange. Isoprene-hydroxy-peroxy radicals (ISOPOO), produced by addition of a hydroxyl radical (OH) across one of the double bonds followed by the rapid addition of molecular oxygen (O₂), react with nitric oxide (NO), hydroperoxy radicals (HO₂), and other peroxy radicals or undergo unimolecular isomerization. Historically, the dominant fate of ISOPOO was thought to be reaction with NO, as mechanisms were developed for urban locations and the NO loss pathway dominates in polluted regions. Under low-NO_x conditions common in forested regions, unimolecular chemistry and reaction with HO₂ are also important. The first-generation product of the reaction of ISOPOO with HO₂, hydroxy hydroperoxide (ISOPOOH), is an important SOA precursor following their oxidation to epoxydiol products (Surratt et al., 2006; Paulot et al., 2009). The other novel low-NO_x pathway recently elucidated is ISOPOO isomerization, which can sustain elevated OH concentrations under low-NO conditions (Peeters and Müller, 2010; Crounse et al., 2011; Teng et al., 2017; Møller et al., 2019). These different branches of ISOPOO pathways produce a different ensemble of oxygenated compounds with low volatility and are thus crucial for accurate prediction of the environmental and climate impacts of isoprene chemistry. In addition to the gas-phase fate of isoprene, field studies found evidence of C5 compounds in ambient particles (Claeys, 2004; Kleindienst et al., 2007), and the modeling of isoprene-derived SOA has been significantly advanced in the past decade (e.g., Compernolle et al., 2009; Marais et al., 2016; Gaston et al., 2014). Wennberg et al. (2018) compile a comprehensive isoprene mechanism incorporating the current knowledge of isoprene chemistry and include the necessary isoprene SOA precursors. The Reduced Caltech Isoprene Mechanism (RCIM), a reduced version with the same product yields of known compounds and minimal simplifications beyond lumping of isomeric compounds and removal of minor (<2 % yield) pathways, is developed concurrently with the explicit mechanism and is more suitable for implementing in canopy–chemistry models.

Here, we develop FORCAsT 1.0 to FORCAsT 2.0 to incorporate important updates to gas-phase isoprene chemistry, isoprene-derived SOA formation, and two physical components of the model. Specifically, major updates include (i) separating the integration of chemistry from the emission and dry deposition (also called operator splitting) to provide more realistic representations of vertical gradients in the forest canopy and to make the chemical module more flexible with future chemical mechanism updates, (ii) improving the vertical mixing parameterization in the boundary layer, (iii) updating the dry deposition velocities for chemical species with available measurements (Nguyen et al., 2015), (iv) implementing the RCIM to reflect the current understanding of isoprene fate under low-NO_x conditions (Wennberg et al., 2018), and (v) extending the MPMPO (Model to Predict the Multiphase Partitioning of Organics) aerosol module (Griffin et al., 2005; Ashworth et al., 2015) to include isoprene-derived SOA formation. We evaluate FORCAsT 2.0's performance against FORCAsT 1.0 (Ashworth et al., 2015) and the observations from the AMOS (Atmospheric Measurements of Oxidants in Summer) field campaign conducted at the University of Michigan Biological Station (UMBS) during the summer of 2016.

2 Model description

The FORCAsT model, based on the CACHE canopy model by Forkel et al. (2006), is a one-dimensional model that couples atmospheric chemistry (gas-phase and gas–particle partitioning) and canopy processes. The vertical resolution of FORCAsT can be configured with a minimum of 20 and a maximum of 60 vertical layers. In this study, the number of model levels is set to 40, with 20 layers within the canopy, 8 layers representing the boundary layer (∼1 km), and the remaining layers extending to the lower troposphere (∼4 km). In addition to the aboveground layers, the model includes 15 soil layers for computing soil heat and moisture storage, as well as exchange with the atmosphere and root extraction (Forkel et al., 2006).

The FORCAsT model parameterizes the processes of radiative transfer, chemical species emission, advection, deposition, vertical exchange, and chemistry and then integrates the energy and mass balance equations at each vertical layer in the canopy. The parameterization of the radiative transfer, emission, advection, and turbulent mixing in the canopy layers has been described extensively in Bryan et al. (2012) and Ashworth et al. (2015) and remains unchanged in the updated version 2.0. Here we describe the major updates to the model, including the operator splitting, gas-phase chemistry, gas–particle partitioning, and some aspects of the dry deposition and turbulent mixing.

2.1 Operator splitting

Processes such as emission, turbulent mixing, dry deposition, and chemical reactions occur at the same time in the atmosphere. However, in numerical models, these processes (i.e., operators) are split and integrated over time and/or space in sequence, which is commonly referred to as operator splitting. It is generally faster to integrate the operators separately than to compute the solution when the operators are treated together (Lapointe et al., 2020). This computational efficiency comes at the cost of an error introduced by the splitting. In the context of atmospheric chemistry modeling, model accuracy is affected by the order in which the operators are applied (Santillana et al., 2016) and by the integration time steps of the operators (also called operator duration) (Philip et al., 2016). In prior studies of chemical transport models, the operator duration causes greater differences in concentrations of reactive emitted species such as nitrogen oxides (up to 5 times; Philip et al., 2016) than the order of operators (up to 10 %; Santillana et al., 2016).

In FORCAsT 1.0, the order of operation is from vertical transport to chemistry, with emissions and dry deposition integrated within the chemistry solver (Table 1). The chemistry solver typically dominates the computational cost of the simulations (Lapointe et al., 2020). To increase computational efficiency and to allow for flexible chemical mechanisms in the future, we separate the chemical solver from emission and dry deposition and integrate the latter two operators with the vertical transport (Table 1). The impacts of this operator splitting on FORCAsT 2.0's performance are discussed in Sect. 3.1.

Table 1 The order of operation in the FORCAsT model. $\frac{\partial C_i}{\partial t}$ denotes the time evolution of concentrations of chemical species i.

2.2 Turbulent mixing

In the surface layer (roughly 10 % of the boundary layer), the eddy diffusivity K is commonly defined as a simple function of height z of the form $K=\mathit{\kappa }{u}_{\ast }z$, where κ is the von Karman constant and u_∗ is the friction velocity (Stull, 1988). In the rest of the boundary layer, defining K is not as clear for numerical models. Several approaches have been used in the literature to define the K profile, such as linearly decreasing from the surface layer to the top of the boundary layer (Estoque, 1963), exponentially decreasing with height from the surface layer, and finding an interpolating polynomial passing through prescribed points with predefined slopes (O'Brien, 1970). The general approach for approximating K in the boundary layer has used a power-law dependence on $z/z_i$, where z_i is the top of the boundary layer and scale parameters are derived from similarity theory or empirically. A commonly used shape function is of the form (Troen and Mahrt, 1986)

$$\begin{array}{}\text{(1)}& K=\mathit{\kappa }{u}_{\ast }z{\left(\mathrm{1}-{\displaystyle \frac{z}{z_i}}\right)}^p.\end{array}$$

The exponent p=2 has been commonly used in the literature (Nissanka et al., 2018), while values between 2 and 3 agree with different observed profiles.

The eddy diffusivity (K) parameterization based on mixing-length theory in FORCAsT 1.0 is described as follows, with greater detail presented in Forkel et al. (2006) and Bryan et al. (2012):

$$\begin{array}{}\text{(2)}& {\displaystyle }K=l^{\mathrm{2}}{\displaystyle \frac{\partial u}{\partial z}},\text{(3)}& {\displaystyle }l={\displaystyle \frac{\mathit{\kappa }z}{\mathrm{1}+\frac{\mathit{\kappa }z}{\mathit{\lambda }}}},\text{(4)}& {\displaystyle \frac{\partial u}{\partial z}}={\displaystyle \frac{u_{\ast }}{\mathit{\kappa }z}},\end{array}$$

where l is the mixing length, u is the mean wind speed, κ is the von Karman constant (0.41), z is the height above the ground, and u_∗ is the friction velocity. The $\frac{\partial u}{\partial z}$ is derived from the common logarithmic expression for the boundary layer (Eq. 4). Combing Eqs. (2), (3), and (4) and rearranging the terms leads to

$$\begin{array}{}\text{(5)}& K=\mathit{\kappa }{u}_{\ast }{\displaystyle \frac{z}{(\mathrm{1}+\frac{\mathit{\kappa }z}{\mathit{\lambda }}{)}^{\mathrm{2}}}},\end{array}$$

where λ is a function of the height z as shown below.

$$\begin{array}{}\text{(6)}& \mathit{\lambda }=\left\{\begin{array}{ll}\mathrm{2.0},& z<h_{\mathrm{c}}\\ max(\mathrm{0.1}{z}_i,\mathrm{2.7}\times {\mathrm{10}}^{-\mathrm{4}}\frac{G}{f}),& h_{\mathrm{c}}\mathit{⩽}z\mathit{⩽}{z}_i\\ \mathrm{2.7}\times {\mathrm{10}}^{-\mathrm{4}}\frac{G}{f},& z\mathit{⩾}{z}_i\end{array}\right.\end{array}$$

Here, h_c is the canopy height, z_i is the height of the boundary layer, G is the geostrophic wind at the top of the boundary layer set to 17 m s⁻¹, and f is the Coriolis parameter. This parameterization yields a relatively small value of K (i.e., <10 m²s⁻¹) in the boundary layer and thus implicitly produces a low boundary layer height (∼250–300 m around 14:00 local time) in FORCAsT 1.0. In addition, this mixing parameterization in FORCAsT 1.0 produces an unrealistic end-of-day peak in isoprene under well-mixed conditions.

To produce more realistic K in the boundary layer and boundary layer height, we adopt the parameterization described in Eq. (1) to calculate the K in the boundary layer. In the present study, we use a cubic power of height z (i.e., p=2) for K in the boundary layer:

$$\begin{array}{}\text{(7)}& K_{\mathrm{new}}=\left\{\begin{array}{ll}\frac{\mathit{\kappa }{u}_{\ast }z}{\mathit{\varphi }(z/L)},& z<z_{\mathrm{sfc}}\\ \frac{a\mathit{\kappa }{u}_{\ast }z(\mathrm{1}-\frac{z}{zi}{)}^{\mathrm{2}}}{\mathit{\varphi }(z/L)},& z\mathit{⩾}{z}_{\mathrm{sfc}}\end{array}\right.,\end{array}$$

where z_sfc is the surface layer height (here assumed to be 10 % of the boundary layer height z_i), and the constant $a=\mathrm{1}/(\mathrm{1}-\frac{z_{\mathrm{sfc}}}{z_i}{)}^{\mathrm{2}}=\mathrm{1.23}$ is used to ensure a continuous transition of K_new from the surface layer to the boundary layer. The z_i is calculated as a function of the sensible heat flux at the top of the canopy (Stull, 1988):

$$\begin{array}{}\text{(8)}& z_i=\sqrt{{\displaystyle \frac{\mathrm{2}(\mathrm{2}c+\mathrm{1})}{\mathit{\gamma }}}\underset{\mathrm{0}}{\overset{t}{\int }}\stackrel{\mathrm{‾}}{(w^{\prime }{\mathit{\theta }}^{\prime }}{)}_{\mathrm{cpy}}},\end{array}$$

where c is a standard entrainment parameter (0.2), and γ is the lapse rate in the free atmosphere (0.0065 K m⁻¹). The integral ${\int }_{\mathrm{0}}^t{\stackrel{\mathrm{‾}}{\left(w^{\prime }{\mathit{\theta }}^{\prime }\right)}}_{\mathrm{cpy}}$ is approximated to ${\stackrel{\mathrm{‾}}{\left(w^{\prime }{\mathit{\theta }}^{\prime }\right)}}_{\mathrm{cpy}}\mathrm{\Delta }t$, where Δt is the elapsed time. The calculation for z_i starts when the heat flux ${\stackrel{\mathrm{‾}}{\left(w^{\prime }{\mathit{\theta }}^{\prime }\right)}}_{\mathrm{cpy}}$ first becomes positive in the morning and continues 1 h longer after ${\stackrel{\mathrm{‾}}{\left(w^{\prime }{\mathit{\theta }}^{\prime }\right)}}_{\mathrm{cpy}}$ becomes negative. The nighttime z_i is set to 200 m.

The stability function $\mathit{\varphi }(z/L)$ is given in Eq. (9) (Nissanka et al., 2018).

$$\begin{array}{}\text{(9)}& \mathit{\varphi }(z/L)=\left\{\begin{array}{ll}(\mathrm{1}-\mathrm{16}z/L{)}^{-\mathrm{1}/\mathrm{2}},& z/L<\mathrm{0}\\ \mathrm{1}+\mathrm{5}z/L,& z/L>\mathrm{0}\\ \mathrm{1},& z/L=\mathrm{0}\end{array}\right.\end{array}$$

Overall, the new parameterization produces a larger K, a more realistic boundary layer height, and thus a more realistic air temperature (see details in Sect. 3.2).

2.3 Dry deposition

The dry deposition velocity to the canopy foliage (v_d) is calculated using the resistance model (Meyers and Baldocchi, 1988; Wesely, 1989) in FORCAsT 1.0, with details provided in Bryan et al. (2012) and Ashworth et al. (2015). Recent work by Nguyen et al. (2015) reports dry deposition velocities based on measured fluxes and concentrations for 16 atmospheric compounds above a southeastern United States forest, suggesting that the parameterization in the FORCAsT 1.0 underestimates dry deposition velocities for these oxygenates (Table 2). We adopt this newer parameterization in FORCAsT 2.0, wherein the major revisions to the resistance model include (1) the addition of the aerodynamic resistance R_a, (2) the formulation of the molecular diffusion R_b, and (3) the addition of temperature dependence to the mesophyll resistance R_m and the cuticular resistance R_c. Details of the revised parameterization can be found in Nguyen et al. (2015). The revised resistance model is highly sensitive to Henry's law coefficient (Nguyen et al., 2015), which is unknown for many species. Here, we apply the revised resistance model to species whose Henry's law coefficients are available and measurement-based v_d exists to validate the estimates, including for H₂O₂ (hydrogen peroxide), HNO₃ (nitric acid), HCOOH (formic acid), HAC (hydroxyacetone), PROPNN (propanone nitrate), HMHP (hydroxymethyl hydroperoxide), INP (isoprene nitrooxy hydroperoxide), IHNs (isoprene hydroxy nitrates), HPALD (hydroperoxy aldehydes), IEPOX (epoxydiols), and ISOPOOH (isoprene hydroxy hydroperoxides) (Table 2). The revised deposition velocity parameterization is also applied to all the isomers for INP, IHN, ISOPOOH, IEPOX, and HPALD in RCIM. Note that the revised dry deposition scheme yields a much smaller v_d (0.35 cm s⁻¹) for methyl vinyl ketone (MVK) and methacrolein (MACR) than the observation-based estimate (up to 2.4 cm s⁻¹; Karl et al., 2010). Therefore, we prescribe the observation-based 2.4 cm s⁻¹ for MVK and MACR in FORCAsT. Because the resistance model assumes that each canopy layer is one “big leaf” perpendicular to the sunlight, the estimates of dry deposition velocities for all compounds are then scaled by the leaf area distribution in the canopy. Species other than those listed above use the old parameterization in FORCAsT 1.0. The comparison of dry deposition velocities between the old and new parameterizations is shown in Table 2, with the new parameterization generally increasing the deposition velocity with the exception of nitric and formic acids. As input data and measurements become available for other species, the revised parameterization can be evaluated against new observations and applied to other species.

Table 2 Estimates of the dry deposition velocities (cm s⁻¹) to the canopy foliage for relevant molecules at 12:00 local time. The measurement-based deposition velocities for a temperate forest are provided as a reference (Nguyen et al., 2015).

NA – not available.

2.4 Isoprene gas-phase mechanism: the Reduced Caltech Isoprene Mechanism (RCIM)

We replace the chemical reactions for isoprene oxidation in the Caltech Atmospheric Chemistry Mechanism (CACM) in FORCAsT 1.0 with RCIM. The original CACM was developed for application to urban conditions (CACM0.0) (Griffin et al., 2002, 2005), and the version incorporated into FORCAsT 1.0 was updated for low-NO_x conditions (Ashworth et al., 2015), hereinafter referred to as CACM1.0 in the following sections. RCIM is version 4.3 of the “reduced” mechanism in the Wennberg et al. (2018) mechanism repository (Bates and Wennberg, 2017), including the essential chemistry required to accurately simulate isoprene oxidation under remote conditions in the atmosphere. Compiled concurrently with the full explicit mechanism, RCIM groups isomers with similar reaction rates and products and lumps minor pathways (<2 % branching ratio) into the major channels to minimize the number of species and reactions while retaining an accurate description of the oxidative fate of those grouped species and lumped pathways (Wennberg et al., 2018). The reduced mechanism includes 119 species and 221 reactions, in contrast to the 385 species and 810 reactions in the explicit isoprene mechanism in Wennberg et al. (2018) and the 113 reactions for isoprene in CACM1.0.

RCIM treats the initial system of allylic and peroxy radicals formed following the addition of OH to isoprene dynamically. Older mechanisms, including CACM1.0 in FORCAsT 1.0, implicitly used fixed distributions of isoprene-hydroxy-peroxy radicals (ISOPOO) derived from experiments performed under high-NO conditions. Addition of O₂ to allylic radicals under ambient conditions is in fact a reversible process, resulting in a dynamic system with interconversion between the six major ISOPOO isomers (two subgroups of three defined by a common OH position) (Peeters et al., 2009; Teng et al., 2017). Wennberg et al. (2018) represent the reversibility of O₂ addition in the explicit mechanism by including 10 species (i.e., four allylic radicals and six major ISOPOOs; see Fig. 3 in Wennberg et al., 2018) and 69 reactions. The reduced mechanism, RCIM, retains this novel treatment of the ISOPOO system (Wennberg et al., 2018), although it simplifies the 10-species radical system to two major ISOPOO isomers, i.e., $(\mathrm{1}-\mathrm{OH},\mathrm{4}-\mathrm{OO})-\mathrm{ISOPOO}$ and $(\mathrm{4}-\mathrm{OH},\mathrm{1}-\mathrm{OO})-\mathrm{ISOPOO}$.

RCIM includes important updates to the formation and fates of isoprene hydroxyl nitrates (IHNs) through pressure- and temperature-dependent parameterizations of the branching ratios and a new structure–activity relationship for calculating the formation of nitrates from multifunctional peroxy radicals without measured yields (Wennberg et al., 2018). The dynamic representation of the ISOPOO isomers also contributes to higher production of IHN than previous mechanisms when included in global models (Bates and Jacob, 2019). In addition, RCIM includes 12 distinct C5 tetrafunctional compounds with unique combinations of functional groups. They are dihydroxy hydroperoxy nitrate, carbonyl hydroxy nitrooxy hydroperoxide, dihydroxy carbonyl nitrate, hydroxy hydroperoxy dialdehyde, hydroxy hydroperoxy dinitrate, dihydroxy dinitrate, carbonyl hydroperoxy-diol, carbonyl hydroxy dinitrate, dihydroxy hydroperxy epoxide, dihydroxy dihydroperoxide, hydroxy nitrooxy dihydroperoxide, and hydroxy nitrooxy hydroperoxy epoxide. Each C5 tetrafunctional compound represents a variety of isomers. Chamber experiments suggest that multifunctionals contribute to isoprene-derived SOA (iSOA) (Ng et al., 2008; Schwantes et al., 2019). Additional aspects of RCIM relative to CACM1.0 include decreased C5-hydroperoxy-aldehyde (HPALD) yields following the 1,6-H shifts of the $\mathrm{Z}-\mathit{\zeta }-\mathrm{OH}-\mathrm{peroxy}$ radicals in Teng et al. (2017) and additional intramolecular H shifts, including rapid peroxy–hydroperoxy shifts, resulting in higher OH recycling under low-NO conditions (Wennberg et al., 2018). These mechanism changes are manifested by changes in gas-phase isoprene oxidation products (Sect. 3.3).

2.5 Isoprene-derived secondary organic aerosol

FORCAsT simulates the partitioning of condensable species into the particle phase using the Model to Predict the Multiphase Partitioning of Organics (MPMPO; Griffin et al., 2005). In FORCAsT 1.0, 99 out of the 300 prognostic species in CACM1.0 are treated as condensable and lumped into 12 surrogate species according to their structures, sources (biogenic or anthropogenic), volatilities, and dissociative capabilities (Ashworth et al., 2015). Specifically, the 12 surrogate groups include 6 anthropogenic aromatic groups, 4 monoterpene-derived biogenic surrogates, and 1 group composed of non-volatile dimers of multifunctional acids from monoterpene oxidation (i.e., phthalic acid). Additionally, one surrogate group based on keto-propanoic acid and oxalic acid (formed from oxidation of methyl vinyl ketone and methacrolein) is considered condensable in CACM1.0-MPMPO. However, explicit formation of SOA from isoprene is missing in FORCAsT 1.0. We incorporate six new surrogate groups in MPMPO to account for the isoprene-derived SOA (iSOA) in FORCAsT 2.0 (Table 3).

Table 3 Properties of the six new iSOA surrogate groups in FORCAsT 2.0. The Henry's law constants below are from Sander (2015) and Safieddine et al. (2017) at 298 K. Temperature dependence of Henry's law constants is included in FORCAsT.

As the iSOA precursors are small organic molecules (number of carbon atoms ≤5) with numerous functional groups, they are expected to be highly hydrophilic. Under humid conditions representative of the summertime PROPHET (Program for Research on Oxidants: Photo-chemistry, Emissions, Transport) boundary layer, aqueous aerosol provides a medium for reactive uptake (Surratt et al., 2009), and thus iSOA is likely aqueous (Couvidat and Seigneur, 2011; Ervens et al., 2011). Marais et al. (2016) proposed a mechanism for irreversible reactive uptake of iSOA precursors by preexisting aqueous aerosols, dependent on the Henry's law coefficients, that is widely used in chemical transport models. Bates and Jacob (2019) estimate a global iSOA yield by mass of 25 % using this mechanism coupled with RCIM. On the other hand, dry chamber experiments (relative humidity <10 %) suggest up to 11 % of iSOA yield (Kroll et al., 2006), suggesting a partitioning between the gas phase and a non-aqueous organic phase. Thornton et al. (2020) estimate an upper bound of the non-aqueous iSOA yield of 3 % under atmospheric conditions using a volatility-driven gas–particle partitioning implemented in a box model.

The MPMPO in FORCAsT only considers organic aerosol species and assumes the partitioning of the gases into two particulate phases: a purely organic aerosol and an aqueous aerosol with associated molecular and ionic components. Equilibrium between the gas and aerosols is assumed for organic species. The equilibrium organic aerosol-phase mass concentration of an individual species i, O_i (µg m⁻³ air), is given by the following relationship (Griffin et al., 2005):

$$\begin{array}{}\text{(10)}& O_i=K_{\mathrm{om},i}{G}_i{M}_{\mathrm{o}}={\displaystyle \frac{G_i{M}_{\mathrm{o}}RT}{M_{\mathrm{om}}{\mathrm{10}}^{\mathrm{6}}{\mathit{\gamma }}_i{p}_{L,i}^{\mathrm{o}}}},\end{array}$$

where K_om,i (m³ air µg⁻¹) is the partitioning coefficient that describes the phase distribution of the condensing organic species (Pankow, 1994), G_i is its corresponding gas-phase concentration (µg m⁻³ air), and M_o is the total concentration (µg m⁻³ air) of organic aerosol mass available to act as the partitioning medium. R is the ideal gas constant ($\mathrm{8.2}\times {\mathrm{10}}^{-\mathrm{5}}$ m³ atm mol⁻¹ K⁻¹), T is temperature (K), M_om is the average molecular weight (g mol⁻¹) of the absorbing organics including both primary and secondary organic compounds, $p_{L,i}^{\mathrm{o}}$ is the pure component vapor pressure (atm) of species i, and γ_i is the activity coefficient of species i in the organic phase. The factor of 10⁶ converts grams (g) to micrograms (µg).

The aqueous-phase concentration of species i, A_i (µg m⁻³ air), is given by (Griffin et al., 2005):

$$\begin{array}{}\text{(11)}& A_i={\displaystyle \frac{G_i\cdot \text{LWC}\cdot H_i}{{\mathit{\gamma }}_{\mathrm{aq},i}}},\end{array}$$

where LWC is the aerosol liquid water content (µg H₂O m⁻³ air), H_i is the Henry’s law coefficient ($\left(\mathrm{\mu }\mathrm{g}\mathrm{\mu }{\mathrm{g}}^{-\mathrm{1}}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\right)/\left(\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}\text{air}\right)$), and γ_aq,i is the activity coefficient of organic species i in the aqueous phase normalized by that at infinite dilution. The aqueous-phase equilibrium is constrained by dissociation of the dissolved organic species, for which pH is needed to calculate the concentration of the charged ions (see Eqs. 13 and 14 in Ashworth et al., 2015).

The LWC in Eq. (11) is calculated offline as the sum of aerosol liquid water associated with inorganic species (LWC_inorg) and organics (LWC_org). The LWC_inorg and pH are calculated by the Extended AIM Aerosol Thermodynamics Model model II (E-AIM, http://www.aim.env.uea.ac.uk/aim/model2/model2d.php, last access: 14 February 2021) (Clegg et al., 1998) using measurements of NH₃(g) and PM_2.5 sulfate, nitrate, and ammonium (Table 4), following the method of Murphy et al. (2017). The LWC_org is calculated according to Petters and Kreidenweis (2007):

$$\begin{array}{}\text{(12)}& {\text{LWC}}_{\mathrm{org}}={\displaystyle \frac{m_{\mathrm{org}}{\mathit{\rho }}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{org}}}}{\displaystyle \frac{{\mathit{\kappa }}_{\mathrm{org}}}{\frac{\mathrm{1}}{\text{RH}}-\mathrm{1}}},\end{array}$$

where m_org is the organic mass concentration from high-resolution time-of-flight aerosol mass spectrometer (HR-ToF-AMS) measurements, ρ_org is the density of organics (1.4 g cm⁻³; Cerully et al., 2015), ρ_w is the density of water, RH is the relative humidity, and κ_org is the hygroscopicity growth for organics taken as 0.08 in this study. The κ_org is derived from the HR-ToF-AMS data using the method described in Cerully et al. (2015):

$$\begin{array}{}\text{(13)}& \mathit{\kappa }={\mathit{ϵ}}_{\mathrm{org}}{\mathit{\kappa }}_{\mathrm{org}}+{\mathit{ϵ}}_{\mathrm{inorg}}{\mathit{\kappa }}_{\mathrm{inorg}},\end{array}$$

where ϵ_org and ϵ_inorg are the volume fractions of organic and inorganic species calculated using AMS mass fraction data, and κ_inorg is the hygroscopicity growth determined from the speciated inorganic concentrations and κ_inorg for individual inorganic species from Padró et al. (2010). κ is the total hygroscopicity growth and has been shown to be insensitive to location and organic fraction (Padró et al., 2010). Here we use a value of 0.23 for κ based on measurements from previous studies (Padró et al., 2010; Cerully et al., 2015). We calculate LWC as the sum of LWC_inorg and LWC_org based on surface observations and vertically scale LWC by the relative humidity to extrapolate to other model heights.

The thermodynamic model UNIFAC (UNIversal Functional group Activity Coefficient) is used to calculate the activity coefficients γ in Eqs. (10) and (11) (Fredenslund et al., 1975). The standard UNIFAC parameters (e.g., alkane group) are found in Hansen et al. (1991), Balslev and Abildskov (2002), and Wittig et al. (2003). The UNIFAC parameters for the functional groups of nitrate and hydroperoxide are taken from Compernolle et al. (2009). The missing UNIFAC parameters (e.g., the unknown interaction parameters) for some functional groups (such as nitrate and hydroperoxide) are set to zero, which introduces uncertainties in estimating the activity coefficient for the new SOA surrogate groups.

3 Model evaluation with observations

The performance of FORCAsT 2.0 is evaluated against the observations obtained at the PROPHET (Program for Research on Oxidants: Photo-chemistry, Emissions, Transport) tower at UMBS during the 2016 AMOS field campaign. Full details of the PROPHET tower and the UMBS site can be found in Millet et al. (2018). Measurement details and references are listed in Table 4. The results presented in this study are based on a 2 d model simulation for the 2 sunny days of 22–23 July 2016; the first day (22 July) is a well-mixed day and the second day (23 July) is relatively stagnant based on micrometeorological analysis (Wei et al., 2020). The model spin-up time is 24 h. The 2 d period is relatively hot for the site, with a mean high of 30.0 ^∘C and a mean daily temperature of 23.9 ^∘C compared to the monthly averages of 25.4 and 21.0 ^∘C, respectively. The canopy structure in the model is identical to that used in Bryan et al. (2012) and Ashworth et al. (2015). The input data are based on the measurements during AMOS 2016, with the model driven by observed PAR (photosynthetically active radiation), standard deviation of the vertical velocity (σ_w), friction velocity (u_*), and the calculated aerosol liquid water content. The initial concentrations for the chemical species are taken from the measurements when available, including ozone (O₃), nitric oxide (NO), nitrogen dioxide (NO₂), formaldehyde (CH₂O), and methyl vinyl ketone and methacrolein (MVK+MACR). The isoprene emission factor is increased from prior FORCAsT 1.0 studies (Bryan et al., 2012; Ashworth et al., 2015), as the last measurement campaign (i.e., CABINEX 2009) occurred during a relatively cool period and utilized an emission factor roughly a factor of 2–3 lower than previously observed at the site (Pressley, 2005; Unger, 2013). The 2 d period evaluated in this study is warmer, and increasing the isoprene emission factor is consistent with the effects of prior day temperature (Guenther et al., 2006). The gas-phase chemical mechanism in FORCAsT includes explicit treatment of two monoterpene surrogate species: α-pinene (APIN) and d-limonene (DLMN). Light dependence of monoterpene emissions is included in the 2 d simulation by changing the synthesis emission factors from 0 to 0.4 nmol m⁻² s⁻¹ for the two monoterpene surrogates, as this has been observed by Ortega et al. (2007) at the site and this emission change improved the agreement with the measured monoterpene concentrations for a more realistic evaluation of the monoterpene-derived SOA. In the following sections, we describe the impact of model changes (Sect. 2.0) on the physical and chemical environment in the forest canopy during this 2 d period.

Millet et al. (2018)Geddes and Murphy (2014)Shi et al. (2021)Dusanter et al. (2009)Wood et al. (2016)Hottle et al. (2009)DiGangi et al. (2011)Cazorla et al. (2015)Vasquez et al. (2018b)Markovic et al. (2012)Bui et al. (2021)Bui et al. (2021)

Table 4 2016 UMBS PROPHET observations utilized in this study for model evaluation, including concentration measurements, measurement heights, instrumental techniques, and references. Chemical compound abbreviations are defined as follows: MVK+MACR (methyl vinyl ketone and methacrolein), NO_x (nitrogen oxides), O₃ (ozone), OH (hydroxyl radical), HO₂ (hydroperoxy radical), RO_x (peroxy radical), CH₂O (formaldehyde), ISOPOOH (isoprene hydroxy hydroperoxide), IEPOX (isoprene epoxydiol), IHN (isoprene hydroxy nitrate), MTN (monoterpene hydroxy nitrate). IEPOX-SOA and 91Fac represent the isoprene-derived epoxydiol organic aerosols and monoterpene-derived organic aerosols, respectively. The canopy height is 22.5 m.

3.1 Operator splitting

As described in Sect. 2.1, the chemistry is separated from the operators of emission and dry deposition in FORCAsT 2.0 (Table 1). The splitting of chemistry from the emission and dry deposition leads to weaker gradients in the vertical profiles for emitted species such as isoprene (Fig. 1). The old order of operators simulates higher concentrations between 0.7 and 1.0 $z/h$ where the emission occurs (Fig. 1a–c), likely due to the higher production rates of isoprene resulting from emission (on the order of 10⁻⁶–10⁻⁵ ppbv s⁻¹ depending on the leaf area distribution) compared to the reaction rates (on the order of 10⁻⁸ ppbv s⁻¹ for isoprene + OH) in the same solver, resulting in the vertical concentration profiles resembling the emission profiles. The new order of operators simulates more realistic in-canopy concentration gradients for the emitted species such as isoprene. Specifically, the vertical profile of isoprene before sunrise is better captured by the new order of operation with an RMSE (root mean square error) less than 0.1 (Fig. 1a, d). The midday in-canopy profiles are also improved by the new order, while simulated concentrations are higher above the canopy than observed (Fig. 1b, e). The overall RMSE for the midday case (Fig. 1e) suggests that the results for the new order agree slightly better with the observations. The differences between the two orders are up to 16 % around the height of 0.85 $z/h$ in the midday case and are in good agreement with Santillana et al. (2016), who report a 10 % difference. In the evening hours, the emission rates and chemistry quickly decrease due to the reduced radiation, leading to more fluctuations in the observed profile that are challenging for models to capture (Fig. 1c).

Figure 1 Impacts of operator splitting on the vertical distribution of isoprene. The modeled and observed isoprene vertical profiles (normalized by the concentration at the canopy) in the morning (a), midday (b), and evening (c) on 22 July 2016. Note that the observed isoprene at 34, 21, 17, 13, 9, and 5 m is interpolated into model levels for comparison. The corresponding root mean square error (RMSE) of the modeled isoprene profiles compared to the observations for the morning (d), midday (e), and evening (f) on 22 July 2016. The RMSEs are also calculated separately in the canopy (In-cpy) and above the canopy (Abv. cpy).

In addition, the time for running the model with the gas-phase chemistry alone (i.e., without running the aerosol module) is reduced from 10 to 3 min. In chemical models, the computational cost is mainly due to the chemical solver, which has a very small internal time step (about 1.2 s in our case). However, the transport solver has a time step of 1 min; therefore, moving the emission and dry deposition from the chemical solver to the transport solver reduces the run time for the gas-phase chemistry by 70 %, making the model more computationally efficient for adding increasingly complex chemical mechanisms in the future. In summary, the differences of the modeled gradients in isoprene concentrations between the two operator orders are relatively small (up to 16 %), with the new order having a more realistic in-canopy profile and a higher computational efficiency.

The impacts of the operator splitting on the gas concentrations are correlated with the chemical lifetime (Fig. 2). The OH and nitric oxide (NO) concentrations increase by 160 %–180 % and 130 % with the new order of operator, respectively, while CO (carbon monoxide) only differs by 13 % (Fig. 2). The increased NO, CH₂O (formaldehyde), and O₃ using the new order improve agreement with the observations (Fig. 4). In addition, the vertical gradients in concentrations between in- and above-canopy decrease with increasing lifetime, with 20 % difference for OH and almost zero for CO (Fig. 2). Overall, the results here draw attention to the influences of the operator splitting on reactive trace gases such as OH and NO, which are critically important for accurately predicting the gas-phase chemistry and aerosol formation. In summary, the advantages of the new operator order from the modeling perspective include (i) generally better agreement with the observations for the critical species NO and OH, (ii) improved in-canopy gradients for emitted species such as isoprene, and (iii) higher computational efficiency. Therefore, we implement the new order of operations in FORCAsT 2.0.

Figure 2 Impacts of operator splitting on species with chemical lifetimes ranging from seconds (OH; Di Carlo et al., 2004) to months (CO; Holloway et al., 2000). The relative differences between the two orders of operator (i.e., $\frac{C_{\mathrm{new}}-C_{\mathrm{old}}}{C_{\mathrm{old}}}$) for daytime (10:00–16:00 local time) average hydroxyl radicals (OH), nitric oxide (NO), formaldehyde (CH₂O), isoprene, ozone (O₃), and carbon monoxide (CO) at two heights (0.8 and 1.6 $z/h$).

3.2 Revised vertical mixing

The revised vertical mixing parameterization produces a larger eddy diffusivity (K) and a more realistic boundary layer height (Fig. 3a). The parameterization of K in the boundary layer remains challenging, as K is a derived parameter that is analogous to diffusion yet is not completely accurate for boundary layer turbulent mixing. Despite its limitations, it is still a useful approximation when a full solution for turbulence is computationally expensive to implement. Kumar and Sharan (2012) compiled estimated values for K in the boundary layer based on previous studies, suggesting that the magnitude of the K ranges from 60 to 200 m² s⁻¹ under weakly unstable conditions ($z/L=-\mathrm{2}$). In FORCAsT 1.0, K peaks at 16 m² s⁻¹ at 150 m above the ground during the daytime (Fig. 3a). This leads to weak mixing, resulting in an unrealistic end-of-day peak in isoprene concentrations around sunset for the well-mixed day of 22 July (Fig. 3c). The revised parameterization produces a larger K that falls within the lower range of previously reported values (Kumar and Sharan, 2012) and is sufficiently strong to produce a realistic isoprene diurnal cycle under well-mixed conditions (Fig. 3c). On the stagnant second day of the simulation (23 July), the revised parameterization also reproduces the end-of-day peak in isoprene but the modeled isoprene is much lower than the observed isoprene (Fig. 3c).

Figure 3 Impacts of the revised mixing parameterization on air temperature and isoprene. (a) Modeled vertical profile of the eddy diffusivity at 14:00 local time on 22 July 2016. (b) Modeled and observed air temperature (T_a) at 46 m. (c) Modeled and observed isoprene at 21 m.

FORCAsT predicts the air temperature based on mixing of surface heat instead of prescribing the measured temperature to nudge the model; therefore, the vertical mixing impacts the temperature profiles (Fig. 3b). During the daytime the canopy behaves as a heat source based on the leaf energy balance, and the stronger mixing distributes the heat more evenly throughout the model atmosphere. Therefore, a larger K reduces air temperature during the daytime in FORCAsT (Fig. 3b). At night the canopy is cooler than the air aloft due to the longwave radiation emission from the canopy, and the low nocturnal mixing (K<3.5 m² s⁻¹, not shown) fails to mix the warmer air down to the canopy, resulting in nighttime cooling of 8 ^∘C during 0–6 model hours and a smaller minimum temperature than the observation around 30 model hours (Fig. 3b). The unrealistic nighttime cooling during 0–6 model hours indicates that (i) heat capacity of leaves may be important as a heat source at night, and/or (ii) nighttime mixing is too low in FORCAsT. The idea of low nighttime mixing is also supported by the overestimated nighttime isoprene concentrations (Fig. 3c), suggesting the need for a better data-constrained nighttime mixing scheme. Note that air temperature (Thomas, 2011) and evening isoprene decay (Wei et al., 2020) also depend on horizontal advection, which is not considered in this study. Attention should be on these estimates when the homogeneity assumption is not met at the site, particularly under stable conditions.

3.3 Isoprene updates to the gas-phase chemistry mechanism

RCIM incorporates the current knowledge of isoprene chemistry under low-NO conditions, which is most notably manifested by the changes in the oxidation products of isoprene (Fig. 4). Adding the RCIM isoprene chemistry does not drive large changes in simulated isoprene concentrations (Fig. 4a), which respond to OH, the isoprene emission rate, and mixing, yet there is a substantial increase in MVK+MACR (Fig. 4c). Generally the RCIM simulates MVK+MACR that shows better agreement with the observations, although simulated concentrations are lower than observed on the second day, which is likely due to underestimated isoprene (Fig. 4a). The measured MVK+MACR peaks before isoprene on the first day, suggesting horizontal advection and/or vertical transport from the residual layer. Both RCIM and CACM1.0 underestimate formaldehyde (CH₂O), a high-yield product of isoprene oxidation, by over 50 % (Fig. 4j). Marvin et al. (2017) show that existing isoprene mechanisms generally underestimate CH₂O by 17 %–33 %. This suggests missing sources for CH₂O, possibly heterogeneous conversion of ISOPOOH on leaves, and/or unaccounted for VOC chemical oxidation in the chemical mechanisms (Canaval et al., 2020; DiGangi et al., 2011).

The modeled NO, NO₂, and O₃ are very similar between the two mechanisms (Fig. 4d–f). The timing of the early morning NO maxima, caused by photolysis of NO₂ transported downward during the morning breakup of the nocturnal boundary layer (Seok et al., 2013), is captured by both mechanisms. Both mechanisms overestimate NO₂ during 0–6 model hours, suggesting missing sinks for NO₂, possibly the aqueous-phase reaction of ${\mathrm{NO}}_{\mathrm{2}}+{\mathrm{NO}}_{\mathrm{3}}\to {\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{5}}$.

Figure 4 Comparison of FORCAsT 2.0 with the RCIM (red) and CACM1.0 (blue) mechanisms against observations in and above the canopy (heights noted). Gas-phase concentrations are evaluated versus observations (grey), including isoprene (34 m), monoterpenes (34 m), methyl vinyl ketone + methacrolein (MVK+MACR, 34 m), nitric oxide (NO, 29 m), nitrogen dioxide (NO₂, 29 m), ozone (O₃, 6 m), hydroxyl radicals (OH, 32 m), hydroperoxy radicals (HO₂, 32 m), total peroxy radicals (RO_x, 30 m), formaldehyde (CH₂O, 21 m), hydroxy hydroperoxide isomers (1,2-ISOPOOH + 4,3-ISOPOOH, 32 m), epoxydiol isomers (cis-IEPOX+rans-IEPOX, 32 m), isoprene hydroxy nitrates (1,2-IHN + 4,3-IHN, 32 m), IHN (19.5 m), and monoterpene hydroxy nitrates (MTN, 19.5 m) in sequence from (a) to (o). Grey shaded areas denote 1 standard deviation of the data when available. Instrumental information is in Table 4. The red dashed lines in (m) and (n) denote modeled IHN with a hydrolysis loss rate of 4×10⁵ M atm⁻¹ s⁻¹ for 1,2-IHN.

Adding RCIM also impacts HO_x, with little change in HO₂ (Fig. 4h) but a substantial increase in OH (by 50 % on the second day; Fig. 4g). Both mechanisms have HO_x regeneration from the H-shift isomerization of isoprene peroxy radicals (ISOPOO), which is important to sustain the OH concentrations under low-NO conditions (Peeters et al., 2009; Bates and Jacob, 2019). CACM1.0 recycles 1.0 HO₂ via the 1,6-H shift, while RCIM recycles 1.5 OH + 0.7 HO₂ through the 1,6-H shift pathway and 1.0 OH through the 1,5-H shift pathway (3.2 equivalents to HO_x). Note that the 1,6-H shift pathway dominates the 1,5-H shift pathway by a factor of 8 in RCIM. In addition, the ISOPOO concentrations in CACM1.0 are lower than in RCIM by 50 % (Fig. 5a). Therefore, the combination of higher ISOPOO and OH regeneration efficiency leads to larger OH concentrations in RCIM. Overall, the modeled OH is in good agreement with the campaign-average measurements (Fig. 4g).

Figure 5 Fate of isoprene hydroxy peroxy radicals (ISOPOO) in RCIM and in CACM1.0. (a) Diurnal profiles of ISOPOO concentrations in RCIM and CACM1.0. (b, c) Fraction of daytime-averaged (12:00–14:00 LT) production rates of ISOPOO through reactions with peroxy radicals (RO₂), hydroperoxy radicals (HO₂), nitric oxide (NO), and H-shift isomerization in RCIM and CACM1.0.

RCIM predicts higher RO_x (HO₂+RO₂) than CACM1.0 (Fig. 4i). Because both mechanisms predict similar HO₂, which is discussed in the previous paragraph, here we focus on RO₂ that is predominantly composed of ISOPOO and monoterpene-derived RO₂ (MT-RO₂). Daytime (10–16 and 34–40 model hours) RCIM-RO₂ is higher than CACM1.0-RO₂, mainly due to differences in ISOPOO (Fig. 5a), accounting for 60 % and 40 % of daytime RO₂ in RCIM and CACM1.0, respectively. Nighttime (0–6 and 24–30 model hours) RCIM-RO₂ comprises 16 % ISOPOO and 84 % MT-RO₂. The majority (87 %) of the nighttime MT-RO₂ is from NO₃-initiated reactions and the remaining (13 %) is O₃-initiated. The observed NO₃ is below the limit of detection of the instrument (1.4 pptv) on the 2 simulated days (not shown). While simulated nighttime NO₃ is <1.4 pptv (roughly 0.4 pptv) and monoterpenes agree with observations, the nighttime RO₂ is still overestimated, suggesting missing sinks for nighttime RO₂. One potential sink is reaction with NO. Using the observed NO to constrain the model reduces the nighttime RO₂ by roughly 30 % (not shown), but this is not sufficient to reproduce the observed RO_x. Another possible sink is the accretion reactions of monoterpene-derived RO₂. The most recent accretion reaction rates for OH-initiated and O₃-initiated monoterpene RO₂ are included in FORCAsT ($\mathrm{3.7}\times {\mathrm{10}}^{-\mathrm{11}}$ and $\mathrm{9.7}\times {\mathrm{10}}^{-\mathrm{12}}$ cm³ molecule⁻¹ s⁻¹, respectively; Berndt et al., 2018). However, laboratory updates on accretion rates for NO₃-initiated monoterpene RO₂ that dominate at night are not available. Finally, dry deposition for RO₂ is not included in the model due to lack of data, and this has the potential to be an additional nighttime sink. Overall, the results suggest that a better understanding of nighttime sinks of RO_x is needed, including chemical losses and dry deposition.

Because of the numerous mechanism changes for low-NO_x isoprene chemistry, large deviations between RCIM and CACM1.0 occur for the isoprene oxidation products, including ISOPOOH (isoprene hydroxy hydroperoxide), IEPOX (isoprene epoxydiol), and IHN (isoprene hydroxy nitrate) (Fig. 4k–n). ISOPOOH, formed by the reaction of HO₂ with ISOPOO, is increased by a factor of 2 in RCIM throughout most of the diurnal cycle (Fig. 4k) due to a combination of higher ISOPOO concentrations and a higher fraction of ISOPOO following the HO₂ pathway (Fig. 5). IEPOX, predominantly produced by the reaction of ISOPOOH+OH, is slightly higher in RCIM (Fig. 4l).

One important difference between the two mechanisms is that the ISOPOO−NO reaction pathway increases in RCIM (29 %) compared to CACM1.0 (11 %) (Fig. 5b, c). From a low-NO study in the Amazon (NO<0.1 ppbv), Liu et al. (2016) found that the ratio of the HO₂ pathway to NO pathway is about unity. For these UMBS simulations with low-NO conditions, RCIM shows a similar NO:HO₂ ratio of 1.3, in contrast to the ratio of 3 in CACM1.0 (Fig. 5b, c). The higher NO pathway percentage in RCIM consequently increases IHN concentrations (Fig. 4m, n). The in-canopy IHN (19.5 m) measured by an iodide-adduct chemical ionization mass spectrometer (CIMS) is well reproduced by RCIM, while the above-canopy IHN (32 m) by GC-HR-ToF-CIMS is overestimated by a factor of 2. Because we do not simulate strong vertical gradients in isoprene (supported by observations; e.g., Wei et al., 2020), the IHN differences in the observations are not likely due to vertical mixing or horizontal advection. Vasquez et al. (2020) show that the 1,2-IHN isomer undergoes a rapid hydrolysis loss not experienced by the 4,3-IHN isomer. The in-canopy (19.5 m) measurements by the iodide-adduct CIMS may be subject to uncertainty in the hydrolysis-related isomer sensitivity (Fig. 4m, n). Vasquez et al. (2020) estimate a condensed-phase hydrolysis coefficient of 4×10⁵ M atm⁻¹ s⁻¹ for 1,2-IHN to match their observed 1,2-IHN to 4,3-IHN isomer ratio. Using this hydrolysis loss coefficient, our model is able to reproduce the isomer ratio. However, the resulting total IHNs (1,2-IHN+4,3-IHN) are still significantly higher than the measurements (dashed line in Fig. 4m, n). A hydrolysis coefficient of 1×10⁸ M atm⁻¹ s⁻¹ is required for the model to match the measured 1,2-IHN concentrations at 32 m; however, the ratio is then underestimated because the model simulates more 4,3-IHN than observed (not shown), which may be due to the greater production of the 4,3-IHN driven by higher than observed HO₂. Overall, the differences in IHN observations and the measured–modeled discrepancies in the IHN isomers suggest that large uncertainties still exist in the simulated IHN concentrations.

The monoterpene hydroxy nitrates (MTNs) in RCIM are slightly lower than that in CACM1.0 at night (Fig. 4o), likely driven by lower nighttime monoterpene concentrations in RCIM (Fig. 4b). Overall, the MTN is underestimated by both RCIM and CACM1.0, and this could be affected by the rapid monoterpene OH RO₂ accretion reactions, which decrease the available monoterpene RO₂ for generating MTN during the daytime. Additionally, the dry deposition of MTN is prescribed to be 50 % of IHNs, which lack observational constraints and may result in an overconsumption of MTN. To further evaluate the representation of the low-NO chemistry of isoprene by the two mechanisms, we compare the ratio of ISOPOOH to MVK+MACR. For unpolluted regions, the reaction of ISOPOO with HO₂ is the dominant pathway, with ISOPOOH isomers as the major oxidation products. Reaction of ISOPOO with NO dominates in polluted regions, with major oxidation products MVK and MACR. MVK+MACR can also be produced through the isomerization and self-reaction of ISOPOO, with MVK dominating over MACR by approximately a factor of 2. Therefore, the ratio of ISOPOOH to MVK+MACR reflects the contribution of the HO₂ oxidation pathway relative to other oxidation pathways of ISOPOO. The consistency between the measured and modeled ISOPOOH fraction thus reflects the skill of the isoprene mechanisms in representing the isoprene chemistry under low-NO conditions. The most frequent (51 %) observed daytime NO levels when the ISOPOOH data are available range from 20 to 40 pptv. The measurement-based ratio of ISOPOOH to MVK+MACR decreases quickly and ranges from 0.15 to 0.05 with increasing NO when NO<40 pptv (Fig. 6), indicating that the fate of ISOPOO is highly sensitive to NO levels under low-NO conditions. The simulated dependence of the ratio on NO concentrations by RCIM is in good agreement with the observed, although the simulated ratio decays faster than the observations. In contrast, the ratios simulated by CACM1.0 are higher and decrease slower with increasing NO concentrations, partly due to lower OH (Fig. 6) that reacts faster with ISOPOOH than with MVK.

In summary, the fate of ISOPOO in RCIM is different than in CACM1.0 (Fig. 5), leading to differences in isoprene oxidation products. The HO₂ pathway ratio is higher in RCIM, leading to higher ISOPOOH concentrations, and the NO pathway of ISOPOO is also enhanced in RCIM, resulting in higher MVK+MACR and IHN concentrations. This comes at the expense of reducing the H-shift isomerization pathway in RCIM, yet OH concentrations are still enhanced due to higher HO_x recycling efficiency in RCIM. Overall, isoprene oxidation products in RCIM generally compare well with observations, and RCIM captures the decreasing trend of $\mathrm{ISOPOOH}/(\mathrm{MVK}+\mathrm{MACR})$ with increasing NO concentrations, indicating a better representation of low-NO chemistry. One aspect to note about the one-dimensional framework is that it prioritizes local processes (emissions, chemistry, and deposition) over that of advection. If we assume a horizontally homogeneous canopy as utilized in flux analysis, then advection should not be a factor. However, observations are clear that the site can be influenced by long-range transport, such as the advection of high-NO_x conditions from southern urban locations like Detroit, Milwaukee, and Chicago (e.g., Cooper et al., 2001; VanReken et al., 2015). One factor to consider is if horizontal advection also influences the diurnal cycle of isoprene and monoterpene oxidation products. For example, above-canopy short-lived oxidation products (such as ISOPOOH and IHN) are overestimated in the model compared to the observations, whereas longer-lived species such as IEPOX (Fig. 4k–m) show better agreement with the observations. Further analysis is needed that includes upwind sources of these oxidation products, as is potential analysis with wind direction and VOC sources.

Figure 6 Daytime (08:00–17:00 local time) measurement-based and modeled ratio of ISOPOOH (isoprene hydroxy hydroperoxides) to MVK+MACR (methyl vinyl ketone + methacrolein). Measurement-based ratios of ISOPOOH to MVK+MACR are from 23–27 July 2016 (Grey dots). The fitted exponential curve is for measurement-based ratios (black line), RCIM ratios (red line), CACM1.0 ratios (blue line), RCIM OH concentrations (red dash line), and CACM1.0 OH concentrations (blue dash line). The yellow patch denotes the most frequently observed nitric oxide (NO) levels (20–40 pptv) during the period of 23–27 July 2016.

3.4 Secondary organic aerosol formation

As discussed in Sect. 2.5, an isoprene-derived SOA (iSOA) parameterization based on Griffin et al. (2005) is implemented in FORCAsT 2.0, incorporating six new surrogate groups for the isoprene-derived precursors such as IEPOX (isoprene epoxydiol) and tetrafunctionals. As described in Ashworth et al. (2015), the partitioning of organic gases into the aerosol phase is a function of the liquid water content (Eq. 11).

Figure 7 Monoterpene-derived secondary organic aerosol (MNT-SOA) (a) and isoprene-derived secondary organic aerosol (iSOA) (b) as a function of time and height. (c, d) Comparison of modeled MNT-SOA and IEPOX-SOA with observations at 30 m. (e) Time series of calculated liquid water content (LWC) with different hygroscopicity growth κ. (f) Vertical profiles of modeled LWC during the daytime (14:00 local time) and in the early morning (6:00 local time).

Both monoterpene SOA (MNT-SOA) and iSOA mainly form in the boundary layer (∼1.5 km; Fig. 7a, b) in response to vertical mixing, vertical distributions of SOA precursors, and meteorological conditions such as relative humidity. The simulated boundary-layer-averaged MNT-SOA is roughly 0.15±0.1 µg m⁻³. The MNT-SOA above the canopy (30 m) is underestimated by a factor of 2–3 compared to the HR-ToF-AMS positive matrix factorization (PMF) concentrations (Fig. 3c), which is likely due in part to the underestimated precursors such as MTN (Fig. 4o). Another possible reason is that the partitioning efficiency for MNT-SOA is underestimated due to phase separation in MPMPO. MPMPO considers the two phases separately: the organic phase and the aqueous phase. In reality, water uptake to organics increases the partitioning medium, which in turn increases the partitioning efficiency of condensable organic species (Pye et al., 2017), and this could effectively reduce the simulated MNT-SOA in MPMPO.

The majority of simulated iSOA (>99 %) is formed through the aqueous phase is and therefore sensitive to the liquid water content (Fig. 7b, e, f). The boundary-layer-averaged iSOA is 1.0±1.0 µg m⁻³ (Fig. 7b). IEPOX-SOA dominates the modeled iSOA formation (Fig. 8a, b), and simulated concentrations agree well with the daytime HR-ToF-AMS PMF-derived IEPOX-SOA at 30 m (Fig. 7d). The modeled IEPOX-SOA shows a diurnal cycle that follows gas-phase IEPOX concentrations (Fig. 4l), and the rapid increase in IEPOX-SOA around 35 model hours results from the increase in liquid water content (Fig. 7e). No such increase in IEPOX-SOA is manifested around 6 model hours when liquid water content is also high due to the absence of iSOA precursors at the beginning of the simulation (Fig. 4l).

Figure 8 Composition of the average isoprene-derived secondary organic aerosol (iSOA) over the 2 d period at two different heights. The composition components are the surrogate groups of GLYX (glyoxal and methylglyoxal), IEPOX (isoprene-epoxydiol, isoprene-carbonyl-hydroxy-epoxide, and hydroxy methyl methyl-a-lactone), IHNs (isoprene hydroxy nitrates), INPs (nitrooxy hydroperoxides), C4 nitrate, and tetrafunctionals. Details about the surrogate groups can be found in Section 2.5. Composition of iSOA at 19 m (a) and 500 m (b). Isoprene-epoxydiol aerosols (IEPOX-SOA, c) and tetrafunctional compound aerosol (tetra-SOA, d) as a function of time and height.

FORCAsT 2.0 simulates an in-canopy iSOA mass yield of 7 %, which is slightly higher than values commonly used in global models (0.9 %–6.8 %), to represent the ambient atmosphere (Carlton et al., 2009). In the model, the two dominant iSOA surrogate groups are the tetrafunctionals (tetra-SOA) and IEPOX (IEPOX-SOA), which combined account for 86 % of the total iSOA in the canopy and 98 % in the boundary layer (Fig. 8a, b). RCIM has a low yield of glyoxal from isoprene, and the GLYX-SOA component is 4 % in the canopy and negligible in the boundary layer. The iSOA from organic nitrates (IHN-SOA and C4-SOA) is more important (11 %) in the canopy likely due to higher NO levels (Fig. 8a, b). Both IEPOX-SOA and tetra-SOA increase with height (Fig. 8c, d) following the gas-phase precursors (not shown). We note that dry deposition of tetrafunctionals lacks data constraints, and in this study, the dry deposition of these compounds is set to be the same as IHN (Table 2) due to the common nitrate functional group. This leads to the deposition of 60 % of the gas-phase tetrafunctionals in the canopy and is a major source of uncertainties for estimation of tetra-SOA. Using RCIM, Bates and Jacob (2019) report an even higher iSOA mass yield of 25 %, with IEPOX-SOA, organonitrates, and tetra-SOA each contributing about one-third to the total mass. The differences in the iSOA yield and composition between this study and Bates and Jacob (2019) may derive from (i) the NO_x levels, which are expected to significantly impact the iSOA precursors (i.e., gas-phase IEPOX, tetrafunctionals, and organic nitrates), and/or (ii) the different iSOA formation schemes.

There are a few points to note in our estimation for LWC. The LWC_org accounts for almost 40 % of total LWC on average (Fig. 7e), suggesting that water uptake to organics is important for iSOA formation at the site. However, LWC_org is highly sensitive to the hygroscopicity growth for organics (κ_org), ranging from 0.01 to 0.5 for slightly to very hygroscopic organic species (Petters and Kreidenweis, 2007). Our calculation gives a mean κ_org of 0.08 for the simulation period. Cerully et al. (2015) partitioned κ_org for different organic aerosol compositions with 0.08±0.02 for monoterpene SOA and 0.20±0.02 for isoprene SOA, respectively. To test the impact of κ_org on LWC estimation, we run a sensitivity simulation with a higher κ_org of 0.15, resulting in an increase of 0.4 µg m⁻³ in LWC_org on average and thus a relatively small increase in total LWC (Fig. 7e). Therefore, the higher κ_org of 0.15 leads to a negligible increase in iSOA formation. Slade et al. (2019) derived different κ_org values for daytime (0.02) and nighttime (0.15) due to diurnal changes in organic aerosol composition at the same study site, which would require a diurnal cycle of LWC_org with lower values at night. Because this result contrasts with current methods to derive LWC, we were unable to test whether this would impact our partitioning simulations, but we note that this large discrepancy could have implications for the simulation of biogenically derived organic aerosol over the diurnal cycle. Finally, we note that due to lack of measurements outside the canopy, the vertical distribution of LWC is scaled by relative humidity (RH) in the model. During the daytime (nighttime), LWC is relatively constant within the boundary (surface) layer and decreases rapidly with altitude above the boundary (surface) layer (Fig. 7f). In general, this trend is consistent with the meteorology soundings of RH from the nearest location (in Gaylord, Michigan), but the soundings suggest that the daytime boundary layer height and the nighttime surface layer height are underestimated in the model despite the improved vertical mixing formulation.

4 Conclusions

We update FORCAsT 1.0 to FORCAsT 2.0 by including (i) splitting the integration of chemistry from emission and dry deposition to provide more realistic representations of vertical gradients in the forest canopy and to make the chemical module more flexible for future chemical mechanism changes, (ii) updating eddy diffusivity in the boundary layer and the dry deposition velocity with available measurements, (iii) implementing the Reduced Caltech Isoprene Mechanism (RCIM) that reflects the current understanding of isoprene fate under low-NO_x conditions, and (iv) extending the aerosol module to include isoprene-derived SOA (iSOA) formation.

The new order of operations reduces the run time for the gas-phase chemistry by 70 %, making the model more computationally efficient for adding increasingly complex chemical mechanisms in the future. Meanwhile, the differences of the modeled gradients in isoprene concentrations between the two orders are relatively small (up to 16 %), with the new operator order having a more realistic in-canopy profile. The revised mixing parameterization yields greater eddy diffusivity in the boundary layer and improves the estimation of the diurnal cycles of isoprene. The end-of-day peak in isoprene is removed for the well-mixed conditions and retained for the stagnant conditions. The revised mixing also improves the simulated daytime air temperature. The revised dry deposition improves the deposition velocity for species when measurements are available.

The fate of isoprene peroxy radicals (ISOPOO) in RCIM is different than in CACM1.0. The key differences include a higher fraction of ISOPOO following the NO pathway (29 %), leading to higher concentrations of MVK+MACR (by a factor of 2) and in-canopy IHN (by a factor of 8) in RCIM. Despite a lower ratio of H-shift isomerization pathway (26 %), RCIM recycles more HO_x than CACM1.0. In addition, the ratio of ISOPOOH to MVK+MACR, which reflects the skill of the mechanisms in representing the isoprene chemistry under low-NO conditions, decreases from 0.15 to 0.02 for NO levels from 20 to 240 pptv. The decreasing trend of the ratio with increasing NO is generally reproduced by RCIM, indicating an improved representation of the chemistry under low-NO condition. FORCAsT 2.0 suggests a 7 % mass yield of isoprene-derived secondary organic aerosol (iSOA) in the canopy. The tetrafunctional-derived iSOA (tetra-SOA) and IEPOX-derived iSOA (IEPOX-SOA) combined account for >86 % of the total iSOA, with IEPOX-SOA dominant by a factor of 3.5. The iSOA from organic nitrates is more important in the canopy than above the canopy, accounting for 11 % of the total iSOA.

Generally, FORCAsT 2.0 shows better agreement with the observations. However, several limitations remain. Due to challenges in nighttime turbulence measurements, the poor simulation of nighttime mixing likely contributes to the cool bias in nighttime air temperature and overestimated isoprene concentrations. While there are many improvements in the evaluation of isoprene chemistry, there are still some notable discrepancies between the observations of isoprene oxidation products and both chemical mechanisms, including first-generation products such as ISOPOOH and CH₂O, as well as second-generation products such as IEPOX. Additionally, uncertainties in the biogenically derived nitrates (IHN and MTN) suggest that additional observations are necessary to understand the pathways of formation and removal. Heterogeneous reactions are not included in the current version, which could otherwise improve the simulation of nighttime NO₃, CH₂O, and IHN. For the formation of biogenically derived aerosols, some thermodynamic parameters (such as UNIFAC parameters and Henry's law constants) for the new surrogate groups (such as nitrate and hydroperoxide) carry uncertainties, and better constraints from laboratory data could improve the simulation of iSOA yields. The phase separation in MPMPO also adds uncertainty to SOA estimation. The diurnal cycle and vertical distribution of liquid water content lack data support and therefore represent additional sources of uncertainty for SOA estimation. Finally, dry deposition for the tetrafunctional compounds and RO_x is poorly constrained due to lack of information. More generally, this study assumes horizontal homogeneity at the study site, which might not be true due to the heterogeneity of the land cover caused by surrounding geography; this includes several large freshwater lakes. This problem is particularly critical under stable conditions. Further analysis of the role of horizontal advection at the site may be useful in understanding how long-range transport influences the local concentrations of long-lived biogenic VOC oxidation products. Taken together, these limitations carry certain levels of uncertainties which can be addressed in future development. Despite these limitations and uncertainties, the revised 1-D model provides a useful tool for understanding how local conditions can influence the oxidation of biogenic VOCs and subsequent aerosol processes.