A single-point modeling approach for the intercomparison and evaluation of ozone dry deposition across chemical transport models (Activity 2 of AQMEII4)

A primary sink of air pollutants and their precursors is dry deposition. Dry deposition estimates differ across chemical transport models, yet an understanding of the model spread is incomplete. Here, we introduce Activity 2 of the Air Quality Model Evaluation International Initiative Phase 4 (AQMEII4). We examine 18 dry deposition schemes from regional and global chemical transport models as well as standalone models used for impact assessments or process understanding. We configure the schemes as single-point models at eight Northern Hemisphere locations with observed ozone fluxes. Single-point models are driven by a common set of site-specific meteorological and environmental conditions. Five of eight sites have at least 3 years and up to 12 years of ozone fluxes. The interquartile range across models in multiyear mean ozone deposition velocities ranges from a factor of 1.2 to 1.9 annually across sites and tends to be highest during winter compared with summer. No model is within 50 % of observed multiyear averages across all sites and seasons, but some models perform well for some sites and seasons. For the first time, we demonstrate how contributions from depositional pathways vary across models. Models can disagree with respect to relative contributions from the pathways, even when they predict similar deposition velocities, or agree with respect to the relative contributions but predict different deposition velocities. Both stomatal and nonstomatal uptake contribute to the large model spread across sites. Our findings are the beginning of results from AQMEII4 Activity 2, which brings scientists who model air quality and dry deposition together with scientists who measure ozone fluxes to evaluate and improve dry deposition schemes in the chemical transport models used for research, planning, and regulatory purposes.

1 Introduction

Dry deposition is a sink of many air pollutants and their precursors, removing compounds from the atmosphere after turbulence transports them to the surface and the compounds stick to or react with surfaces. Dry deposition may be a key influence on air pollution levels, including during high-pollution episodes (Vautard et al., 2005; Solberg et al., 2008; Emberson et al., 2013; Huang et al., 2016; Anav et al., 2018; Baublitz et al., 2020; Clifton et al., 2020b; Lin et al., 2020; Gong et al., 2021). Dry deposition can also harm plants when gases diffuse through stomata (Krupa, 2003; Ainsworth et al., 2012; Lombardozzi et al., 2013; Grulke and Heath, 2019; Emberson, 2020). In particular, stomatal uptake of ozone adversely impacts crop yields (Mauzerall and Wang, 2001; McGrath et al., 2015; Guarin et al., 2019; Hong et al., 2020; U.S. EPA, 2020a, b; Tai et al., 2021) and alters terrestrial carbon and water cycles (Ren et al., 2007; Sitch et al., 2007; Lombardozzi et al., 2015; Oliver et al., 2018).

Chemical transport models are key tools for research, planning, and regulatory purposes, including quantifying the influence of meteorology and emissions on air pollution. Accurate estimates of sinks like dry deposition are needed for source attribution, and simulated tropospheric and near-surface abundances of air pollutants are highly sensitive to dry deposition (Wild, 2007; Tang et al., 2011; Walker, 2014; Bela et al., 2015; Beddows et al., 2017; Hogrefe et al., 2018; Baublitz et al., 2020; Sharma et al., 2020; Ryan and Wild, 2021; Liu et al., 2022). However, chemical transport models do not always reproduce the observed variability in dry deposition nor in the near-surface abundances of air pollutants expected to be influenced strongly by dry deposition (Hardacre et al., 2015; Clifton et al., 2017; Kavassalis and Murphy, 2017; Silva and Heald, 2018; Travis and Jacob, 2019; Visser et al., 2021; Wong et al., 2022; Ye et al., 2022; Lam et al., 2023).

Previous work has shown that dry deposition rates differ across chemical transport models (Dentener et al., 2006; Flechard et al., 2011; Hardacre et al., 2015; Li et al., 2016; Vivanco et al., 2018). Differences can stem from the dry deposition scheme (Le Morvan-Quéméner et al., 2018; Wu et al., 2018; Wong et al., 2019; Otu-Larbi et al., 2021; Sun et al., 2022) as well as from the near-surface concentrations of the air pollutant and model-specific forcing related to meteorology and land use/land cover (LULC) (Hardacre et al., 2015; Tan et al., 2018; Zhao et al., 2018; Huang et al., 2022). Even with the same forcing, deposition velocities, or the strength of the dry deposition independent of near-surface concentrations, can vary by 2- to 3-fold across models (Flechard et al., 2011; Schwede et al., 2011; Wu et al., 2018; Wong et al., 2019; Cao et al., 2022; Sun et al., 2022), highlighting the roles of process representation and parameter choice. Minimizing uncertainties in dry deposition schemes is not only important for the chemical transport models used for forecasting and regulatory applications but also for improved understanding of long-term trends and variability in air pollution and impacts on humans, ecosystems, and resources, and building the related predictive ability in global Earth system and chemistry–climate models (Archibald et al., 2020; Clifton et al., 2020a).

In addition to occurring after diffusion through stomata, dry deposition occurs via nonstomatal pathways, including soil and leaf cuticles, as well as snow and water (Wesely and Hicks, 2000; Helmig et al., 2007; Fowler et al., 2009; Hardacre et al., 2015; Clifton et al., 2020a). For ozone, a recent review estimates that nonstomatal uptake is 45 % on average of dry deposition over physiologically active vegetation (Clifton et al., 2020a). For highly soluble gases, nonstomatal uptake may dominate dry deposition (e.g., Karl et al., 2010; Nguyen et al., 2015; Clifton et al., 2022). Observations show strong unexpected spatiotemporal variations in nonstomatal uptake (Lenschow et al., 1981; Godowitch, 1990; Fuentes et al., 1992; Rondón et al., 1993; Coe et al., 1995; Mahrt et al., 1995; Fowler et al., 2001; Coyle et al., 2009; Helmig et al., 2009; Stella et al., 2011; Rannik et al., 2012; Potier et al., 2015; Wolfe et al., 2015; Fumagalli et al., 2016; Clifton et al., 2017, 2019; Stella et al., 2019). In general, a dearth of common process-oriented diagnostics has prevented a clear picture of the stomatal versus nonstomatal deposition pathways driving differences in past model intercomparisons.

Measured turbulent fluxes are the best existing observational constraints on dry deposition, but they are limited with respect to providing information on the relative roles of individual deposition pathways (Fares et al., 2018; Clifton et al., 2020a; He et al., 2021). While we can build a mechanistic understanding of individual processes with laboratory and field chamber measurements (Fuentes and Gillespie, 1992; Cape et al., 2009; Fares et al., 2014; Fumagalli et al., 2016; Sun et al., 2016a, b; Potier et al., 2017; Finco et al., 2018), the dry deposition models that are used to scale processes to the ecosystem level, often the same models used in dry deposition schemes in chemical transport models, are highly empirical and poorly constrained. For example, a recent synthesis found that, while we have a basic knowledge of processes controlling ozone dry deposition, the relative importance of various processes remains uncertain, and we lack the ability to predict spatiotemporal changes well (Clifton et al., 2020a).

Launched in 2009, the Air Quality Model Evaluation International Initiative (AQMEII) has organized several activities (Rao et al., 2011). The fourth phase of AQMEII emphasizes process-oriented investigation of deposition in a common framework (Galmarini et al., 2021). AQMEII4 has two main activities. Activity 1 evaluates both wet and dry deposition across regional air quality models (Galmarini et al., 2021). Here, we introduce Activity 2, which examines dry deposition schemes as standalone single-point models at eight sites with ozone flux observations. Importantly, single-point models are forced with the same site-specific observational datasets of meteorology and ecosystem characteristics; thus, the intercomparison and evaluation can focus on deposition processes and parameters, as recommended by a recent review (Clifton et al., 2020a).

The four aims of Activity 2 are as follows:

1. to quantify the performance of a variety of dry deposition schemes under identical conditions,

2. to understand how different deposition pathways contribute to the inter-model spread,

3. to probe the sensitivity of schemes to environmental factors and to examine variability in the sensitivities across schemes, and

4. to understand differences in dry deposition simulated in regional models in Activity 1.

Our effort builds on recent work using observation-driven single-point modeling of dry deposition schemes at Borden Forest (Wu et al., 2018), Ispra and Hyytiälä (Visser et al., 2021), and two sites in China (Cao et al., 2022), but it is designed to test more sites and schemes as well as to gain a better understanding of inter-model differences. For example, the sites examined represent a range of ecosystems in North America, Europe, and Israel, and single-point models are required to archive process-level diagnostics to facilitate an understanding of simulated variations. Although our fourth aim is to contextualize differences among regional air quality models in Activity 1, we also include additional schemes in Activity 2 (e.g., from global chemical transport models and schemes that are always used as standalone models) to allow for a more comprehensive range of inter-model variation.

Below, we describe the single-point modeling approach (Sect. 2) and fully document the individual single-point models using consistent language, units, and variable names (when appropriate) (Sect. 3). We also describe the Northern Hemisphere locations and site-specific meteorological and environmental datasets used to drive and evaluate the single-point models and the post-processing of observed and simulated values (Sect. 4). Our focus on ozone dry deposition reflects the availability of long-term ozone flux measurements. In the results (Sect. 5), we present how models differ with respect to capturing the observed seasonality in ozone deposition velocities, including the contribution of different deposition pathways and how some environmental factors drive changes. We focus on multiyear averages and, thus, climatological evaluation but examine some aspects of interannual variability for sites with ozone flux records with 3 or more years. We then present a summary of our findings (Sect. 6). To our knowledge, this is the first model intercomparison demonstrating how the contribution of different pathways varies across dry deposition schemes and contributes to the model spread in ozone deposition velocities.

2 Single-point modeling approach

The single-point models used here are standalone dry deposition schemes driven by a consistent set of meteorological and environmental inputs from observations at sites with ozone fluxes. The single-point models were extracted from regional models used in AQMEII4 Activity 1 as well as other chemical transport models or have always been configured as single-point models. In general, dry deposition schemes vary in structure and level of detail in terms of the processes represented. Because there is limited documentation in the peer-reviewed literature of dry deposition schemes (especially as the schemes are configured in chemical transport models) and as complete and consistent model descriptions aid our effort here, we fully describe the participating single-point models using consistent language, units, and variable names (when appropriate). Due to our focus on ozone, we limit our description to dry deposition of ozone. For brevity, we also limit our description to the implementation of the schemes in the single-point models at the eight sites examined, as opposed to how the schemes work as embedded within the chemical transport models (hereinafter, “host models”).

We note that surface- and soil-dependent variable choices (e.g., volumetric soil water content at wilting point) in the host model implementation of the schemes have likely been optimized for generalized LULC and soil classification schemes as well as environmental conditions and meteorology generated or used by the host model. Thus, our prescription of common site-specific variables across the single-point models in this study may create potential inconsistencies with the performance of the schemes inside host models. However, this separation and unification of variables that describe the surface and soil states is key for realistic estimates of the model spread due to structural uncertainty with respect to the processes and parameters directly related to dry deposition.

Table 1 gives the measured and inferred variables used to force single-point models as well as other common variables used in the models. The meaning and units of variables listed in Table 1 are consistent throughout the paper. If a variable is not listed in Table 1, that variable's meaning and units cannot be assumed to be consistent across models nor within the paper. The first time that we mention the variables included in Table 1, we refer to Table 1.

Table 1 Variables related to forcing datasets for single-point models.

The forcing variables provide inputs to drive models with detailed dependencies on biophysics, such as coupled photosynthesis–stomatal conductance models, as well as models that depend mainly on atmospheric conditions. Not every model uses every forcing variable. In general, the input variables used by each single-point model should reflect the operation of the dry deposition scheme. For example, if the scheme in the host model ingests precipitation to calculate canopy wetness, rather than ingesting canopy wetness, the single-point model should ingest precipitation to calculate canopy wetness.

We note that dry deposition schemes in many chemical transport models use methods derived from classic schemes like Wesely (1989). Implementations of classic schemes may deviate from original parameterization description papers in ways that can affect simulated rates (e.g., Hardacre et al., 2015) but may not be well documented. For example, there may be changes to LULC-specific parameters or the use of different LULC categories. In addition, implementations may tie processes to variables like leaf area index to capture seasonal changes rather than relying on season-specific parameters. To foster understanding of how adaptations from original schemes influence simulated dry deposition rates, we encouraged participation in Activity 2 from models using schemes based on classic parameterizations, in addition to models with different approaches.

Like many model intercomparisons, our effort is an “ensemble of opportunity” (e.g., Galmarini et al., 2004; Tebaldi and Knutti, 2007; Potempski and Galmarini, 2009; Solazzo and Galmarini, 2015; Young et al., 2018) and may underestimate structural uncertainty due to process and parameter differences across models. Nonetheless, the design of our effort, with emphasis on processes, parameters, and sensitivities, is designed to explore uncertainty more systematically than past attempts.

The first set of Activity 2 simulations is driven by inputs from observations, and those simulations are examined here. Future work will examine sensitivity tests in which dry deposition is calculated with perturbed values of input variables (e.g., air temperature and leaf area index). We will also design tests that isolate the influence of input parameters (e.g., initial resistance to stomatal uptake and the field capacity of soil).

Diagnostic outputs required from single-point models follow the requirements of Activity 1 (see Table 4 in Galmarini et al., 2021). Among the required outputs are effective conductances (Paulot et al., 2018; Clifton et al., 2020b) for dry deposition to plant stomata, leaf cuticles, the lower canopy, and soil. (Note that not all single-point models simulate deposition to the lower canopy.) As explained and defined in Galmarini et al. (2021), an effective conductance [m s⁻¹] represents the portion of v_d that occurs via a single pathway. An effective conductance is distinct from an absolute conductance, which represents an individual process. (Note that a conductance is the inverse of a resistance.) The sum of the effective conductances across all pathways represented is v_d. In contrast, calculating v_d with absolute conductances requires considering the resistance framework. Archiving effective conductances facilitates comparison of the contribution of each pathway across dry deposition schemes with varying resistance frameworks and differing resistances to transport. Previous model comparisons examine absolute conductances and suggest that differences in pathways or processes lead to differences in v_d (Wu et al., 2018; Huang et al., 2022). Our approach with effective conductances offers a more apples-to-apples comparison across models, allowing us to definitively say whether a given pathway leads to inter-model differences in v_d.

3 Documentation of single-point models

The classic big-leaf resistance network for ozone deposition velocity (v_d) [m s⁻¹] (Table 1) is based on three resistances, which are added in series, as follows:

$$\begin{array}{}\text{(1)}& v_{\mathrm{d}}={\left(r_{\mathrm{a}}+r_{\mathrm{b}}+r_c\right)}^{-\mathrm{1}}.\end{array}$$

Here, the variable r_a is aerodynamic resistance, r_b is quasi-laminar boundary layer resistance around the bulk surface, and r_c is surface resistance. Throughout the paper, all resistances (denoted by r) are in units of seconds per meter [s m⁻¹]. The single-point models examined here employ Eq. (1), with two exceptions. The exceptions are the Multi-layer Canopy and Chemistry Exchange Model (MLC-CHEM), which is a multilayer canopy model that simulates the ozone concentration gradient within the canopy, and the Community Multiscale Air Quality Model (CMAQ) Surface Tiled Aerosol and Gaseous Exchange (STAGE), which uses surface-specific quasi-laminar resistances. In this section, we describe methods for r_a and r_b across models (Tables S1, S2, and S3 in the Supplement) and the ozone-specific dry deposition parameters (Table S4). Equations for r_c and the v_d equation for CMAQ STAGE, which deviates from Eq. (1), are in the individual model subsections below. In the model subsection for MLC-CHEM, we describe how the model diagnoses v_d from the canopy-top ozone flux and the resistances associated with dry deposition.

With one exception (CMAQ STAGE), the single-point models use r_a equations based on Monin–Obukhov similarity theory (Table S1). However, the exact forms of the Monin–Obukhov similarity theory equations vary across the models.

Obukhov length (L) [m] (Table 1) is often used in r_a equations but is not observed. Most model L equations are similar, apart from whether models use virtual or ambient temperature and whether they include bounds on L (and what the bounds are) (Table S2).

Models are configured to accept inputs and return predicted values at the specified ozone flux measurement height at the given site (i.e., reference height z_r [m]; Table 1). Roughness length (z₀) [m] (Table 1) and displacement height (d) [m] (Table 1) are also often used in r_a equations; however, they are not observed and are especially important for estimating fluxes at z_r rather than the lowest atmospheric level of the host model. We supply estimates of z₀ and d for the models that employ them. Estimates follow Meyers et al. (1998):

$$\begin{array}{}\text{(2)}& {\displaystyle }{z}_{\mathrm{0}}=h\left(\mathrm{0.23}-{\displaystyle \frac{{\mathrm{LAI}}^{\mathrm{0.25}}}{\mathrm{10}}}-{\displaystyle \frac{a-\mathrm{1}}{\mathrm{10}}}\right),\text{(3)}& {\displaystyle }d=h\left(\mathrm{0.05}+{\displaystyle \frac{{\mathrm{LAI}}^{\mathrm{0.2}}}{\mathrm{2}}}+{\displaystyle \frac{a-\mathrm{1}}{\mathrm{20}}}\right).\end{array}$$

Here, the variable h [m] is canopy height (Table 1), LAI [m² m⁻²] is leaf area index (Table 1), and a [unitless] is a parameter based on LULC. Meyers et al. (1998) suggest a correction for z₀ if LAI is less than 1, but we do not employ this correction given that it creates discontinuities in the time series.

Table S3 provides the quasi-laminar boundary layer resistance equations. Most models treat this resistance for the bulk surface (i.e., r_b in Eq. 1) and most use r_b from Wesely and Hicks (1977). A key part of r_b parameterizations is the ratio scaling the quasi-laminar boundary layer resistance for heat to ozone (R_diff,b) (Table S4). $R_{\mathrm{diff},\mathrm{b}}=Sc/Pr$, where Sc [unitless] is the Schmidt number (Table 1) and Pr [unitless] is the Prandtl number (Table 1). All but one employ $R_{\mathrm{diff},\mathrm{b}}=\mathit{Sc}/\mathit{Pr}=\mathit{\kappa }/D_{{\mathrm{O}}_{\mathrm{3}}}$ where κ [m² s⁻¹] is thermal diffusivity of air (Table 1) and $D_{{\mathrm{O}}_{\mathrm{3}}}$ [m² s⁻¹] is ozone diffusivity in air (Table 1); however, values of κ and $D_{{\mathrm{O}}_{\mathrm{3}}}$ vary across models (Table S4).

Table S4 presents model prescriptions for ozone-specific dry deposition parameters: the ratio that scales stomatal resistance from water vapor to ozone (R_diff,st), the reactivity factor for ozone (f₀) [unitless] (Table 1), and the Henry's law constant for ozone (H) [M atm⁻¹] (Table 1). Where used, values of f₀ and H are very similar across models. Some models employ temperature dependencies on H. Notably, values of R_diff,st vary from 1.2 to 1.7 across models. (The current estimate of this ratio is 1.51 according to Massman, 1998.) The Global Environmental Multiscale – Modelling Air quality and Chemistry model (GEM-MACH) Zhang and models based on GEOS-Chem are the models that prescribe lower R_diff,st values.

3.1 WRF-Chem Wesely

The Weather Research and Forecasting (WRF) model coupled with Chemistry (WRF-Chem) uses a scheme based on Wesely (1989). Parameters in Table S5 are site- and season-specific. WRF-Chem has two seasons: midsummer with lush vegetation [day of year between 90 and 270] and autumn with unharvested croplands [day of year less than 90 or greater than 270].

3.1.1 Surface resistance

Surface resistance (r_c) is calculated as follows:

$$\begin{array}{}\text{(4)}& r_c={\left({\displaystyle \frac{\mathrm{1}}{r_{\mathrm{st}}+r_{\mathrm{m}}}}+{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{cut}}}}+{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{dc}}+r_{cl}+r_T}}+{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{ac}}+r_g+r_T}}\right)}^{-\mathrm{1}}.\end{array}$$

To consider the effects of air temperature (T_a) [^∘C] (Table 1), resistance r_T (Walmsley and Wesely, 1996) is calculated as follows:

$$\begin{array}{}\text{(5)}& r_T=\mathrm{1000}{e}^{-T_{\mathrm{a}}-\mathrm{4}}.\end{array}$$

In addition to the use of r_T in Eq. (4), r_T is used in the equation for cuticular resistance below.

3.1.2 Stomatal and mesophyll resistances

Stomatal resistance (r_st) is expressed as follows:

$$\begin{array}{}\text{(6)}& r_{\mathrm{st}}=R_{\mathrm{diff},\mathrm{st}}{\displaystyle \frac{r_{\mathrm{i}}}{f\left(T_{\mathrm{a}}\right)f\left(G\right)}}.\end{array}$$

The parameter r_i is initial resistance for stomatal uptake (Table S5).

The effects of T_a are calculated as follows:

$$\begin{array}{}\text{(7)}& f\left(T_{\mathrm{a}}\right)=T_{\mathrm{a}}{\displaystyle \frac{\left(\mathrm{40}-T_{\mathrm{a}}\right)}{\mathrm{400}}}.\end{array}$$

The effects of incoming shortwave radiation (G) [W m⁻²] (Table 1) are expressed as follows:

$$\begin{array}{}\text{(8)}& f\left(G\right)={\left(\mathrm{1}+{\left({\displaystyle \frac{\mathrm{200}}{G+\mathrm{0.1}}}\right)}^{\mathrm{2}}\right)}^{-\mathrm{1}}.\end{array}$$

Mesophyll resistance (r_m) is calculated as follows:

$$\begin{array}{}\text{(9)}& r_{\mathrm{m}}={\left({\displaystyle \frac{H}{\mathrm{3000}}}+\mathrm{100}{f}_{\mathrm{0}}\right)}^{-\mathrm{1}}.\end{array}$$

3.1.3 Cuticular resistance

Cuticular resistance (r_cut) is expressed as follows:

$$\begin{array}{}\text{(10)}& r_{\mathrm{cut}}=\left\{\begin{array}{l}\frac{r_{\mathrm{lu}}+r_T}{\frac{H}{{\mathrm{10}}^{\mathrm{5}}}+f_{\mathrm{0}}},\mathrm{RH}\le \mathrm{0.95}\mathrm{and}P=\mathrm{0}\\ {\left(\frac{\mathrm{1}}{W}+\frac{\mathrm{3}}{r_{\mathrm{lu}}+r_T}\right)}^{-\mathrm{1}},\mathrm{RH}>\mathrm{0.95}\mathrm{or}P>\mathrm{0}.\end{array}\right.\end{array}$$

Here, the parameter r_lu is initial resistance for cuticular uptake (Table S5), RH is relative humidity [fractional] (Table 1), and P is precipitation rate [mm h⁻¹] (Table 1). The parameter W is used to account for leaf wetness:

$$\begin{array}{}\text{(11)}& W=\left\{\begin{array}{l}\mathrm{3000},P=\mathrm{0}\\ \mathrm{1000},P>\mathrm{0}.\end{array}\right.\end{array}$$

3.1.4 Resistances to the lower canopy and ground (and associated resistances to transport)

The resistance associated with within-canopy convection (r_dc) is calculated as follows:

$$\begin{array}{}\text{(12)}& r_{\mathrm{dc}}=\mathrm{100}\left(\mathrm{1}+{\displaystyle \frac{\mathrm{1000}}{G}}\right).\end{array}$$

Resistances to the lower canopy (r_cl), in-canopy turbulence (r_ac), and the ground (r_g) are prescribed (Table S5).

3.2 GEOS-Chem Wesely

GEOS-Chem is based on Wesely (1989). Wang et al. (1998) describe the initial implementation. We examine the scheme from GEOS-Chem v13.3. Parameters in Table S6 are site-specific. If there is snow, surface resistance (r_c) is calculated with the snow parameters in Table S6.

3.2.1 Surface resistance

Surface resistance (r_c) is expressed as follows:

$$\begin{array}{}\text{(13)}& r_c={\left({\displaystyle \frac{\mathrm{1}}{r_{\mathrm{st}}+r_m}}+{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{cut}}}}+{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{dc}}+r_{\mathrm{cl}}}}+{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{ac}}+r_{\mathrm{g}}}}\right)}^{-\mathrm{1}}.\end{array}$$

To consider effects of T_a, resistance r_T is calculated as follows:

$$\begin{array}{}\text{(14)}& r_T=\mathrm{1000}{e}^{-T_{\mathrm{a}}-\mathrm{4}}.\end{array}$$

The variable r_T is used in the below equations for the resistances to cuticles, the lower canopy, and the ground.

3.2.2 Stomatal and mesophyll resistances

Stomatal resistance (r_st) is expressed as follows:

$$\begin{array}{}\text{(15)}& r_{\mathrm{st}}=R_{\mathrm{diff},\mathrm{st}}{\displaystyle \frac{r_{\mathrm{i}}}{{\mathrm{LAI}}_{\mathrm{eff}}f\left(T_{\mathrm{a}}\right)}}.\end{array}$$

Here, the parameter r_i is the initial resistance to stomatal uptake (Table S6), and LAI_eff [m² m⁻²] is the effective LAI, which is the surface area of actively transpiring leaves per ground surface area. The variable LAI_eff is calculated as a function of LAI and solar zenith angle (θ) [^∘] (Table 1), and the cloud fraction using a parameterization developed by Wang et al. (1998). In GEOS-Chem, if G is 0, LAI_eff equals 0.01. For the single-point model, we set G to 0 when θ is greater than 95^∘ so that nighttime r_st values in the single-point model are more similar to GEOS-Chem. GEOS-Chem almost never has nonzero G values at night, but measured values are frequently small and nonzero. Here, the cloud fraction is assumed to be zero.

The effects of T_a are expressed as follows:

$$\begin{array}{}\text{(16)}& f\left(T_{\mathrm{a}}\right)=\left\{T_{\mathrm{a}}\begin{array}{l}\mathrm{0.01},T_{\mathrm{a}}\le \mathrm{0}\\ \frac{\left(\mathrm{40}-T_{\mathrm{a}}\right)}{\mathrm{400}},\mathrm{0}<T_{\mathrm{a}}<\mathrm{40}\\ \mathrm{0.01},\mathrm{40}\le T_{\mathrm{a}}.\end{array}\right.\end{array}$$

Mesophyll resistance (r_m) is calculated as follows:

$$\begin{array}{}\text{(17)}& r_{\mathrm{m}}={\left({\displaystyle \frac{H}{\mathrm{3000}}}+\mathrm{100}{f}_{\mathrm{0}}\right)}^{-\mathrm{1}}.\end{array}$$

3.2.3 Cuticular resistance

Cuticular resistance (r_cut) is expressed as follows:

$$\begin{array}{}\text{(18)}& r_{\mathrm{cut}}=\left\{\begin{array}{l}\frac{r_{\mathrm{lu}}+min\mathit{\{}{r}_T,r_{\mathrm{lu}}\mathit{\}}}{\mathrm{LAI}}{\left(\frac{H}{{\mathrm{10}}^{\mathrm{5}}}+f_{\mathrm{0}}\right)}^{-\mathrm{1}},\frac{r_{\mathrm{lu}}+min\mathit{\{}{r}_T,r_{\mathrm{lu}}\mathit{\}}}{\mathrm{LAI}}<\mathrm{9999}\\ {\mathrm{10}}^{\mathrm{12}},\frac{r_{\mathrm{lu}}+min\mathit{\{}{r}_T,r_{\mathrm{lu}}\mathit{\}}}{\mathrm{LAI}}\ge \mathrm{9999}.\end{array}\right.\end{array}$$

The parameter r_lu is initial resistance for cuticular uptake (Table S6).

3.2.4 Resistances to the lower canopy and ground (and associated resistances to transport)

The resistance associated with in-canopy convection (r_dc) is expressed as follows:

$$\begin{array}{}\text{(19)}& r_{\mathrm{dc}}=\mathrm{100}\left(\mathrm{1}+{\displaystyle \frac{\mathrm{1000}}{G+\mathrm{10}}}\right).\end{array}$$

The resistance to surfaces in the lower canopy (r_cl) is calculated as follows:

$$\begin{array}{}\text{(20)}& r_{\mathrm{cl}}={\left({\displaystyle \frac{H}{{\mathrm{10}}^{\mathrm{5}}\left(r_{\mathrm{cl},\mathrm{S}}+\min \left\{r_T,r_{\mathrm{cl},\mathrm{S}}\right\}\right)}}+{\displaystyle \frac{f_{\mathrm{0}}}{r_{\mathrm{cl},\mathrm{O}}+min\left\{r_T,r_{\mathrm{cl},\mathrm{O}}\right\}}}\right)}^{-\mathrm{1}}.\end{array}$$

Here, the parameters r_cl,S and r_cl,O are initial resistances to the lower canopy (Table S6).

The resistance to turbulent transport to the ground (r_ac) is constant (Table S6).

Resistance to the ground (r_g) is expressed as follows:

$$\begin{array}{}\text{(21)}& r_{\mathrm{g}}={\left({\displaystyle \frac{H}{{\mathrm{10}}^{\mathrm{5}}\left(r_{\mathrm{g},\mathrm{S}}+min\left\{r_T,r_{\mathrm{g},\mathrm{S}}\right\}\right)}}+{\displaystyle \frac{f_{\mathrm{0}}}{r_{\mathrm{g},\mathrm{O}}+min\left\{r_T,r_{\mathrm{g},\mathrm{O}}\right\}}}\right)}^{-\mathrm{1}}.\end{array}$$

Here, the parameters r_g,S and r_g,O are initial resistances to uptake on the ground (Table S6).

3.3 IFS

The ECMWF Integrated Forecasting System (IFS) uses two schemes based on Wesely (1989): Météo-France's SUMO (Michou et al., 2004) (“IFS SUMO Wesely”) and GEOS-Chem 12.7.2 (“IFS GEOS-Chem Wesely”). Unless stated otherwise, the components are the same between schemes. The IFS SUMO Wesely parameters in Table S7 are site- and season-specific. Seasons are defined as follows: “transitional spring” (March, April, and May), “midsummer” (June, July, and August), “autumn” (September, October, and November), and “late autumn” (December, January, and February). Otherwise, if there is snow, the model employs the “winter, snow” parameter values. The IFS GEOS-Chem Wesely parameters in Table S8 are site-specific. If there is snow, the model employs the snow type. For snow type, only the resistance to surfaces in the lower canopy (r_cl) is defined [1000 s m⁻¹].

3.3.1 Surface resistance

Surface resistance (r_c) is expressed as follows:

$$\begin{array}{}\text{(22)}& r_c={\left({\displaystyle \frac{\mathrm{1}}{r_{\mathrm{st}}+r_m}}+{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{cut}}}}+{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{dc}}+r_{\mathrm{cl}}}}+{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{ac}}+r_{\mathrm{g}}+r_T}}\right)}^{-\mathrm{1}}.\end{array}$$

To consider the effects of T_a, resistance r_T is calculated as follows:

$$\begin{array}{}\text{(23)}& r_T=\mathrm{1000}{e}^{-T_{\mathrm{a}}-\mathrm{4}}.\end{array}$$

In addition to the use of r_T in Eq. (22), r_T is included in cuticular resistance equations below.

3.3.2 Stomatal and mesophyll resistances

For IFS SUMO Wesely, stomatal resistance (r_st) is expressed as follows:

$$\begin{array}{}\text{(24)}& r_{\mathrm{st}}=R_{\mathrm{diff},\mathrm{st}}{\displaystyle \frac{r_{\mathrm{i}}}{\mathrm{LAI}f\left(G\right)f\left(\mathrm{VPD}\right)f\left(w_{\mathrm{2}}\right)}}.\end{array}$$

The parameter r_i is initial resistance to stomatal uptake (Table S7).

The effects of G are calculated as follows:

$$\begin{array}{}\text{(25)}& f\left(G\right)=min\left\{{\displaystyle \frac{\mathrm{0.004}G+\mathrm{0.5}}{\mathrm{0.81}\left(\mathrm{0.004}G+\mathrm{1}\right)}},\mathrm{1}\right\}.\end{array}$$

The effects of vapor pressure deficit (VPD) [kPa] (Table 1) are expressed as follows:

$$\begin{array}{}\text{(26)}& f\left(\mathrm{VPD}\right)=\left\{\begin{array}{l}{e}^{\mathrm{0.3}\mathrm{VPD}},\mathrm{forests},\\ \mathrm{1},\mathrm{otherwise}.\end{array}\right.\end{array}$$

The effects of root-zone soil water content (w₂) [m³ m⁻³] (Table 1) are calculated as follows:

$$\begin{array}{}\text{(27)}& f\left(w_{\mathrm{2}}\right)=\left\{\begin{array}{l}\mathrm{0},w_{\mathrm{2}}<w_{\mathrm{wlt}}\\ \frac{w_{\mathrm{2}}-w_{\mathrm{wlt}}}{w_{\mathrm{fc}}-w_{\mathrm{wlt}}},w_{\mathrm{wlt}}<w_{\mathrm{2}}<w_{\mathrm{fc}}\\ \mathrm{1},w_{\mathrm{2}}>w_{\mathrm{fc}}.\end{array}\right.\end{array}$$

Here, the parameter w_wlt is the soil water content at wilting point [m³ m⁻³] (Table 1) and w_fc is the soil water content at field capacity [m³ m⁻³] (Table 1).

For IFS GEOS-Chem Wesely, stomatal resistance (r_st) is expressed as follows:

$$\begin{array}{}\text{(28)}& r_{\mathrm{st}}=R_{\mathrm{diff},\mathrm{st}}{\displaystyle \frac{r_{\mathrm{i}}}{{\mathrm{LAI}}_{\mathrm{eff}}f\left(T_{\mathrm{a}}\right)}}.\end{array}$$

Here, the parameter r_i is initial resistance to stomatal uptake (Table S8), and LAI_eff [m² m⁻²] is the effective LAI, which is the surface area of actively transpiring leaves per ground surface area of actively transpiring leaves. The variable LAI_eff is calculated as a function of the LAI, θ, and the cloud fraction using a parameterization developed by Wang et al. (1998). In GEOS-Chem, if G is 0, LAI_eff is equal to 0.01. For the single-point model, we set G to 0 when θ is greater than 95^∘. GEOS-Chem almost never has nonzero G at night, but measured values are frequently small and nonzero. This change makes nighttime r_st values in the single-point model more similar GEOS-Chem. Here, the cloud fraction is assumed to be zero.

The effects of T_a are calculated as follows:

$$\begin{array}{}\text{(29)}& f\left(T_{\mathrm{a}}\right)=T_{\mathrm{a}}{\displaystyle \frac{\mathrm{40}-T_{\mathrm{a}}}{\mathrm{400}}}.\end{array}$$

For both configurations, mesophyll resistance (r_m) as expressed as follows:

$$\begin{array}{}\text{(30)}& r_{\mathrm{m}}={\left({\displaystyle \frac{H}{\mathrm{3000}}}+\mathrm{100}{f}_{\mathrm{0}}\right)}^{-\mathrm{1}}.\end{array}$$

3.3.3 Cuticular resistance

For IFS SUMO Wesely,

$$\begin{array}{}\text{(31)}& r_{\mathrm{cut}}=(r_{\mathrm{lu}}+r_T){\left({\displaystyle \frac{H}{{\mathrm{10}}^{\mathrm{5}}}}+f_{\mathrm{0}}\right)}^{-\mathrm{1}}.\end{array}$$

The parameter r_lu is initial resistance for cuticular uptake (Table S7).

For IFS GEOS-Chem Wesely,

$$\begin{array}{}\text{(32)}& r_{\mathrm{cut}}={\displaystyle \frac{(r_{\mathrm{lu}}+r_T)}{\mathrm{LAI}}}{\left({\displaystyle \frac{H}{{\mathrm{10}}^{\mathrm{5}}}}+f_{\mathrm{0}}\right)}^{-\mathrm{1}}.\end{array}$$

The parameter r_lu is initial resistance to cuticular uptake (Table S8).

3.3.4 Resistances to the lower canopy and ground (and associated resistances to transport)

The resistance associated with in-canopy convection (r_dc) is expressed as follows:

$$\begin{array}{}\text{(33)}& r_{\mathrm{dc}}=\mathrm{100}\left(\mathrm{1}+{\displaystyle \frac{\mathrm{1000}}{G}}\right).\end{array}$$

Resistances to surfaces in the lower canopy (r_cl), in-canopy turbulence (r_ac), and ground (r_g) are prescribed (Tables S7 and S8).

3.4 GEM-MACH Wesely

Operationally, GEM-MACH uses a dry deposition scheme based on Wesely (1989) (Makar et al., 2018). Parameters defined in Table S9 are site- and sometimes season-specific. Table S10 describes how seasons are distributed as a function of month and latitude.

3.4.1 Surface resistance

Surface resistance (r_c) is calculated as follows:

$$\begin{array}{}\text{(34)}& r_c={\left({\displaystyle \frac{\mathrm{1}-W}{r_{\mathrm{st}}+r_{\mathrm{m}}}}+{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{cut}}}}+{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{dc}}+r_{\mathrm{cl}}}}+{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{ac}}+r_{\mathrm{g}}}}\right)}^{-\mathrm{1}}.\end{array}$$

The parameter W [fractional] is used to account for leaf wetness:

$$\begin{array}{}\text{(35)}& W=\left\{\begin{array}{l}\mathrm{0.5},P>\mathrm{1}\mathrm{mm}{\mathrm{h}}^{-\mathrm{1}}\mathrm{or}\mathrm{RH}>\mathrm{0.95}\\ \mathrm{0},\mathrm{otherwise}.\end{array}\right.\end{array}$$

3.4.2 Stomatal resistance and mesophyll resistance

Stomatal resistance (r_st) is based on Jarvis (1976), Zhang et al. (2002a, 2003), and Baldocchi et al. (1987):

$$\begin{array}{}\text{(36)}& r_{\mathrm{st}}=R_{\mathrm{diff},\mathrm{st}}{\displaystyle \frac{r_{\mathrm{i}}}{\mathrm{LAI}max\left\{f\left(G\right)f\left(\mathrm{VPD}\right)f\left(T_{\mathrm{a}}\right)f\left(c_a\right),\mathrm{0.0001}\right\}}}.\end{array}$$

Here, the parameter r_i is initial resistance to stomatal uptake (Table S9).

Curve-fitting of data from Jarvis (1976) and Ellsworth and Reich (1993) was used to infer the following:

$$\begin{array}{}\text{(37)}& f\left(G\right)=\max \mathit{\{}\mathrm{0.206}\ln \left(G\right)-\mathrm{0.605},\mathrm{0}\mathit{\}}.\end{array}$$

The effects of VPD are expressed as follows:

$$\begin{array}{}\text{(38)}& \begin{array}{rl}f& \left(\mathrm{VPD}\right)=\max \left\{\mathrm{0.0},max\left\{\mathrm{1.0},\right.\right.\\ & \left.\left.\left(\mathrm{1.0}-\mathrm{0.03}\left(\mathrm{1}-\mathrm{RH}\right){\mathrm{10}}^{\frac{\mathrm{0.7859}+\mathrm{0.03477}{T}_{\mathrm{a}}}{\mathrm{1}+\mathrm{0.00412}{T}_{\mathrm{a}}}}\right)\right\}\right\}.\end{array}\end{array}$$

The effects of T_a are calculated as follows:

$$\begin{array}{}\text{(39)}& f\left(T_{\mathrm{a}}\right)={\left({\displaystyle \frac{\left(T_{\mathrm{a}}-T_{\min }\right)\left(T_{\max }-T_{\mathrm{a}}\right)}{\left(T_{\mathrm{opt}}-T_{\min }\right)\left(T_{\max }-T_{\mathrm{opt}}\right)}}\right)}^{\mathrm{0.62}}.\end{array}$$

Here, the parameters T_min, T_max, and T_opt [^∘C] are the minimum, maximum, and optimum temperatures, respectively (Table S9).

The effects of the ambient carbon dioxide mixing ratio ([CO₂]) [ppmv] (Table 1) are expressed as follows:

$$\begin{array}{}\text{(40)}& f\left(c_a\right)=\left\{\begin{array}{l}\mathrm{1},\left[{\mathrm{CO}}_{\mathrm{2}}\right]\le \mathrm{100}\\ \mathrm{1}-(\mathrm{7.35}\times {\mathrm{10}}^{-\mathrm{4}}\ln (\ln \left(G\right))-\mathrm{8.75}\\ \times {\mathrm{10}}^{-\mathrm{4}}\left)\right[{\mathrm{CO}}_{\mathrm{2}}],\mathrm{100}<[{\mathrm{CO}}_{\mathrm{2}}]<\mathrm{1000}\\ \mathrm{0},\left[{\mathrm{CO}}_{\mathrm{2}}\right]\ge \mathrm{1000}.\end{array}\right.\end{array}$$

Mesophyll resistance (r_m) is calculated as follows:

$$\begin{array}{}\text{(41)}& r_{\mathrm{m}}={\left(\mathrm{LAI}\left({\displaystyle \frac{H}{\mathrm{3000}}}+\mathrm{100}{f}_{\mathrm{0}}\right)\right)}^{-\mathrm{1}}.\end{array}$$

3.4.3 Cuticular resistance

Cuticular resistance (r_cut) is expressed as follows:

$$\begin{array}{}\text{(42)}& r_{\mathrm{cut}}={\displaystyle \frac{r_{\mathrm{lu}}}{\mathrm{LAI}}}{\left({\displaystyle \frac{H}{{\mathrm{10}}^{\mathrm{5}}}}+f_{\mathrm{0}}\right)}^{-\mathrm{1}}.\end{array}$$

The parameter r_lu is initial resistance to cuticular uptake (Table S9).

3.4.4 Resistances to the lower canopy and ground (and associated resistances to transport)

The resistance associated with in-canopy convection (r_dc) is calculated as follows:

$$\begin{array}{}\text{(43)}& r_{\mathrm{dc}}=\mathrm{100}+\left(\mathrm{1}+{\displaystyle \frac{\mathrm{1000}}{G+\mathrm{10}}}\right).\end{array}$$

The resistance posed by uptake to the lower canopy (r_cl) is expressed as follows:

$$\begin{array}{}\text{(44)}& r_{\mathrm{cl}}={\left({\displaystyle \frac{H}{{\mathrm{10}}^{\mathrm{5}}{r}_{\mathrm{cl},\mathrm{S}}}}+{\displaystyle \frac{f_{\mathrm{0}}}{r_{\mathrm{cl},\mathrm{O}}}}\right)}^{-\mathrm{1}}.\end{array}$$

Here, the parameters r_cl,S and r_cl,O are initial resistances to uptake by surfaces in the lower canopy (Table S9).

The parameter r_ac is resistance to in-canopy turbulence and r_g is resistance to the ground; both are prescribed (Table S9).

3.5 GEM-MACH Zhang

GEM-MACH also has an implementation of Zhang et al. (2002b). Parameters in Table S11 are site-specific.

3.5.1 Surface resistance

Surface resistance (r_c) is expressed as follows:

$$\begin{array}{}\text{(45)}& r_c=min\left\{\mathrm{10},{\left({\displaystyle \frac{\mathrm{1}-W}{r_{\mathrm{st}}}}+{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{cut}}}}+{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{ac}}+r_{\mathrm{g}}}}\right)}^{-\mathrm{1}}\right\}.\end{array}$$

The variable W [fractional] is used to account for leaf wetness:

$$\begin{array}{}\text{(46)}& W=\left\{\begin{array}{l}\min \left\{\mathrm{0.5},\frac{G-\mathrm{200}}{\mathrm{800}}\right\},\mathrm{precipitation}\mathrm{or}\mathrm{dew},\\ T_{\mathrm{a}}>\mathrm{1},G>\mathrm{200}\\ \mathrm{0},\mathrm{otherwise}.\end{array}\right.\end{array}$$

Precipitation is assumed to occur if P is greater than 0.20 mm h⁻¹. Dew is assumed to occur if P is less than 0.20 mm h⁻¹ and

$$\begin{array}{}\text{(47)}& u^{*}<c_{\mathrm{dew}}{\displaystyle \frac{\mathrm{1.5}}{\max \left\{\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}},\frac{\mathrm{0.622}{e}_{\mathrm{sat}}\left(\mathrm{1}-\mathrm{RH}\right)}{p_{\mathrm{a}}}\right\}}}.\end{array}$$

Here, the variable e_sat [Pa] is saturation vapor pressure (Table 1), p_a [Pa] is air pressure (Table 1), and c_dew is the dew coefficient [0.3].

3.5.2 Stomatal resistance

Stomatal resistance (r_st) is expressed as follows:

$$\begin{array}{}\text{(48)}& r_{\mathrm{st}}=R_{\mathrm{diff},\mathrm{st}}{\displaystyle \frac{r_{\mathrm{i}}(\mathrm{LAI},\mathrm{PAR})}{f\left(T_{\mathrm{a}}\right)f\left(\mathrm{VPD}\right)f\left({\mathit{\psi }}_{\mathrm{leaf}}\right)}}.\end{array}$$

Here, the variable r_i (LAI, PAR) is initial resistance to stomatal uptake that varies with LAI and PAR, based on Norman (1982) and Zhang et al. (2001):

$$\begin{array}{}\text{(49)}& r_{\mathrm{i}}(\mathrm{LAI},\mathrm{PAR})={\left({\displaystyle \frac{{\mathrm{LAI}}_{\mathrm{sun}}}{r_{\mathrm{i}}\left(\mathrm{1}+\frac{b_{\mathrm{rs}}}{{\mathrm{PAR}}_{\mathrm{sun}}}\right)}}+{\displaystyle \frac{{\mathrm{LAI}}_{\mathrm{shd}}}{r_{\mathrm{i}}\left(\mathrm{1}+\frac{b_{\mathrm{rs}}}{{\mathrm{PAR}}_{\mathrm{shd}}}\right)}}\right)}^{-\mathrm{1}}.\end{array}$$

The parameter r_i is initial resistance to stomatal uptake (Table S11), b_rs [W m⁻²] is empirical (Table S11), and LAI_sun and LAI_shd [m² m⁻²] are sunlit and shaded LAI, respectively. The latter two parameters are calculated as follows:

$$\begin{array}{}\text{(50)}& {\displaystyle }{\mathrm{LAI}}_{\mathrm{sun}}={\displaystyle \frac{\mathrm{1}-e^{-K_b\mathrm{LAI}}}{K_b}},\text{(51)}& {\displaystyle }{\mathrm{LAI}}_{\mathrm{shd}}=\mathrm{LAI}-{\mathrm{LAI}}_{\mathrm{sun}}.\end{array}$$

The variable K_b is the canopy light extinction coefficient [unitless]:

$$\begin{array}{}\text{(52)}& K_b={\displaystyle \frac{\mathrm{0.5}}{\cos \left(\frac{\mathit{\pi }}{\mathrm{180}}\mathit{\theta }\right)}}.\end{array}$$

The variables PAR_sun and PAR_shd [W m⁻²] are photosynthetically active radiation reaching sunlit and shaded leaves, respectively:

$$\begin{array}{}\text{(53)}& {\displaystyle }\begin{array}{rl}{\mathrm{PAR}}_{\mathrm{shd}}& ={\mathrm{PAR}}_{\mathrm{diff}}{e}^{-\mathrm{0.5}{\mathrm{LAI}}^a}\\ & +\mathrm{0.07}{\mathrm{PAR}}_{\mathrm{dir}}\left(\mathrm{1}-\mathrm{0.1}\mathrm{LAI}\right)e^{-\cos \left(\frac{\mathit{\pi }}{\mathrm{180}}\mathit{\theta }\right),}\end{array}\text{(54)}& {\displaystyle }{\mathrm{PAR}}_{\mathrm{sun}}={\mathrm{PAR}}_{\mathrm{shd}}+{\displaystyle \frac{\mathrm{0.5}{\mathrm{PAR}}_{\mathrm{dir}}^b}{\cos \left(\frac{\mathit{\pi }}{\mathrm{180}}\mathit{\theta }\right).}}\end{array}$$

If LAI is greater than 2.5 m² m⁻² and G is less than 200 W m⁻², the empirical parameter a equals 0.8 and b equals 0.8. Otherwise, a equals 0.07 and b equals 1. The calculation of direct and diffuse components of PAR (PAR_dir and PAR_diff, respectively) has been updated from Zhang et al. (2001) to follow Iqbal (1983):

$$\begin{array}{}\text{(55)}& {\displaystyle }{\mathrm{PAR}}_{\mathrm{dir}}=G{\mathrm{FRAD}}_V{\mathrm{FD}}_V,\text{(56)}& {\displaystyle }{\mathrm{PAR}}_{\mathrm{diff}}=G{\mathrm{FRAD}}_V\left(\mathrm{1}-{\mathrm{FD}}_V\right).\end{array}$$

The variable FRAD_V is calculated as follows:

$$\begin{array}{}\text{(57)}& {\mathrm{FRAD}}_V={\displaystyle \frac{R_V}{R_V+R_N}}.\end{array}$$

The variables R_V and R_N are expressed as follows:

$$\begin{array}{}\text{(58)}& {\displaystyle }{R}_N={\mathrm{RD}}_M+{\mathrm{RD}}_N,\text{(59)}& {\displaystyle }{R}_V={\mathrm{RD}}_U+{\mathrm{RD}}_V.\end{array}$$

The variable RD_U is calculated as follows:

$$\begin{array}{}\text{(60)}& {\mathrm{RD}}_U=\mathrm{600}\cos \left({\displaystyle \frac{\mathit{\pi }}{\mathrm{180}}}\mathit{\theta }\right)e^{\frac{-\mathrm{0.185}{p}_{\mathrm{a}}}{p_{\mathrm{std}}\cos \left(\frac{\mathit{\pi }}{\mathrm{180}}\mathit{\theta }\right)}}.\end{array}$$

The variable p_std is standard air pressure [1.0132×10⁵ Pa].

The variable RD_V is expressed as follows:

$$\begin{array}{}\text{(61)}& {\mathrm{RD}}_V=\mathrm{0.42}\left(\mathrm{600}-{\mathrm{RD}}_U\right)\cos \left({\displaystyle \frac{\mathit{\pi }}{\mathrm{180}}}\mathit{\theta }\right).\end{array}$$

The variable RD_M is calculated as follows:

$$\begin{array}{}\text{(62)}& \begin{array}{rl}{\mathrm{RD}}_M& =\cos \left({\displaystyle \frac{\mathit{\pi }}{\mathrm{180}}}\mathit{\theta }\right)\left(\mathrm{720}{e}^{\left(-\frac{\mathrm{0.06}{p}_{\mathrm{a}}}{p_{\mathrm{std}}\cos \left(\frac{\mathit{\pi }}{\mathrm{180}}\mathit{\theta }\right)}\right)}\right.\\ & \left.-\left(\mathrm{1320}\cdot \mathrm{0.077}{\left({\displaystyle \frac{\mathrm{2}{p}_{\mathrm{a}}}{p_{\mathrm{std}}\cos \left(\frac{\mathit{\pi }}{\mathrm{180}}\mathit{\theta }\right)}}\right)}^{\mathrm{0.3}}\right)\right).\end{array}\end{array}$$

The variable RD_N is expressed as follows:

$$\begin{array}{}\text{(63)}& \begin{array}{rl}{\mathrm{RD}}_N& =\mathrm{0.65}\cos \left({\displaystyle \frac{\mathit{\pi }}{\mathrm{180}}}\mathit{\theta }\right)\left(\mathrm{720}\right.\\ & \left.-{\mathrm{RD}}_M-\left(\mathrm{1320}\cdot \mathrm{0.077}{\left({\displaystyle \frac{\mathrm{2}{p}_{\mathrm{a}}}{p_{\mathrm{std}}\cos \left(\frac{\mathit{\pi }}{\mathrm{180}}\mathit{\theta }\right)}}\right)}^{\mathrm{0.3}}\right)\right).\end{array}\end{array}$$

The variable FD_V is calculated as follows:

$$\begin{array}{}\text{(64)}& {\mathrm{FD}}_V=\left\{\begin{array}{ll}\mathrm{0.941124}{\mathrm{RD}}_U/R_V,& \frac{G}{R_V+R_N}\ge \mathrm{0.89}\\ \left(\mathrm{1}-{\left(\frac{\left(\mathrm{0.9}-\frac{G}{R_V+R_N}\right)}{\mathrm{0.7}}\right)}^{\frac{\mathrm{2}}{\mathrm{3}}}\right){\mathrm{RD}}_U/R_V,& \mathrm{0.21}\ge \frac{G}{R_V+R_N}<\mathrm{0.89}\\ \mathrm{0.00955}{\mathrm{RD}}_U/R_V,& \frac{G}{R_V+R_N}<\mathrm{0.21}.\end{array}\right.\end{array}$$

The effects of T_a are as follows:

$$\begin{array}{}\text{(65)}& f\left(T_{\mathrm{a}}\right)=\left({\displaystyle \frac{T_{\mathrm{a}}-T_{\min }}{T_{\mathrm{opt}}-T_{\min }}}\right){\left({\displaystyle \frac{T_{\max }-T_{\mathrm{a}}}{T_{\max }-T_{\mathrm{opt}}}}\right)}^{\frac{T_{\max }-T_{\mathrm{opt}}}{T_{\max }-T_{\min }}}.\end{array}$$

Here, the parameters T_min, T_max, and T_opt [^∘C] are the minimum, maximum, and optimum temperatures, respectively (Table S11).

The effects of VPD are expressed as follows:

$$\begin{array}{}\text{(66)}& f\left(\mathrm{VPD}\right)=\min \left\{\max \left\{\mathrm{1}-b_{\mathrm{vpd}}\mathrm{VPD},\mathrm{0}\right\},\mathrm{1}\right\}.\end{array}$$

The parameter b_vpd [kPa⁻¹] is empirical (Table S11).

The effects of the leaf water potential (ψ_leaf) [MPa] (Table 1) are calculated as follows:

$$\begin{array}{}\text{(67)}& f\left({\mathit{\psi }}_{\mathrm{leaf}}\right)=\min \left\{\max \left\{{\displaystyle \frac{{\mathit{\psi }}_{\mathrm{leaf}}-{\mathit{\psi }}_{\mathrm{leaf},\mathrm{2}}}{{\mathit{\psi }}_{\mathrm{leaf},\mathrm{1}}-{\mathit{\psi }}_{\mathrm{leaf},\mathrm{2}}}},\mathrm{0}\right\},\mathrm{1}\right\}.\end{array}$$

The variable ψ_leaf is approximated as follows:

$$\begin{array}{}\text{(68)}& {\mathit{\psi }}_{\mathrm{leaf}}=-\mathrm{0.72}-\mathrm{0.0013}G.\end{array}$$

The parameters ψ_leaf,1 and ψ_leaf,2 [MPa] are empirical (Table S11).

3.5.3 Cuticular resistance

Cuticular resistance (r_cut) is expressed as follows:

$$\begin{array}{}\text{(69)}& r_{\mathrm{cut}}=\left\{\begin{array}{l}\max \left\{\mathrm{100},\frac{c_{\mathrm{cut},\mathrm{dry}}}{u^{*}{\mathrm{LAI}}^{\mathrm{0.25}}{e}^{\mathrm{3}\mathrm{RH}}}\right\},T_{\mathrm{a}}\ge -\mathrm{1},\mathrm{neither}\mathrm{precipitation}\mathrm{nor}\mathrm{dew}\\ \frac{c_{\mathrm{cut},\mathrm{wet}}}{u^{*}\sqrt{\mathrm{LAI}}},T_{\mathrm{a}}\ge -\mathrm{1},\mathrm{precipitation}\mathrm{or}\mathrm{dew}\mathrm{occurring}\\ \max \left\{\mathrm{100},\frac{C_{\mathrm{cut},\mathrm{dry}}}{u^{*}{\mathrm{LAI}}^{\mathrm{0.25}}{e}^{\mathrm{3}\mathrm{RH}}}\min \left\{\mathrm{2},e^{\mathrm{0.2}\left(-\mathrm{1}-T_{\mathrm{a}}\right)}\right\}\right\},T_{\mathrm{a}}<-\mathrm{1}.\end{array}\right.\end{array}$$

The variable u^* [m s⁻¹] is friction velocity (Table 1) and c_cut,dry [unitless] is a coefficient related to dry cuticular uptake (Table S11).

If the fraction of snow coverage (f_snow) is greater than 10⁻⁴, a correction is applied:

$$\begin{array}{}\text{(70)}& r_{\mathrm{cut}}={\left({\displaystyle \frac{\mathrm{1}-f_{\mathrm{snow}}}{r_{\mathrm{cut}}}}+{\displaystyle \frac{f_{\mathrm{snow}}}{\mathrm{2000}}}\right)}^{-\mathrm{1}}.\end{array}$$

If LAI is less than $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{6}}$ m² m⁻², r_cut is very large.

The fraction of snow coverage (f_snow) is calculated as follows:

$$\begin{array}{}\text{(71)}& f_{\mathrm{snow}}=\min \left\{\mathrm{1},{\displaystyle \frac{\mathrm{SD}}{{\mathrm{SD}}_{\max }}}\right\}.\end{array}$$

The variable SD [cm] is snow depth (Table 1) and SD_max [cm] is maximum snow depth (Table S11).

3.5.4 Resistance to the ground (and associated resistance to transport)

The resistance to in-canopy turbulence (r_ac) is expressed as follows:

$$\begin{array}{}\text{(72)}& r_{\mathrm{ac}}=r_{\mathrm{ac}\mathrm{0}}{\displaystyle \frac{{\mathrm{LAI}}^{\mathrm{0.25}}}{{\left(u^{*}\right)}^{\mathrm{2}}}}.\end{array}$$

The variable r_ac0 is calculated as follows:

$$\begin{array}{}\text{(73)}& r_{\mathrm{ac}\mathrm{0}}=r_{\mathrm{ac}\mathrm{0},\min }+{\displaystyle \frac{\mathrm{LAI}-{\mathrm{LAI}}_{\min }}{{\mathrm{LAI}}_{\max }-{\mathrm{LAI}}_{\min }}}\left(r_{\mathrm{ac}\mathrm{0},\max }-r_{\mathrm{ac}\mathrm{0},\min }\right).\end{array}$$

Here, the parameters LAI_min and LAI_max [m² m⁻²] are the minimum and maximum LAI values across the site's observational record, respectively, and r_ac0,min and r_ac0,max are initial resistances (Table S11).

Ground resistance (r_g) is prescribed but modified under certain conditions. If T_s is less than −1 ^∘C,

$$\begin{array}{}\text{(74)}& r_{\mathrm{g}}=r_{\mathrm{g}}\min \left\{\mathrm{2},e^{-\mathrm{0.2}\left(T_{\mathrm{s}}+\mathrm{1}\right)}\right\}.\end{array}$$

The near-surface air temperature (T_s) is approximated from a linear interpolation between T_a and T_g to a height of 1.5 m.

If f_snow (see Eq. 71) is greater than or equal to 10⁻⁴,

$$\begin{array}{}\text{(75)}& r_{\mathrm{g}}={\left({\displaystyle \frac{\mathrm{1}-min\left\{\mathrm{1},\mathrm{2}{f}_{\mathrm{snow}}\right\}}{r_{\mathrm{g}}}}+{\displaystyle \frac{min\left\{\mathrm{1},\mathrm{2}{f}_{\mathrm{snow}}\right\}}{\mathrm{2000}}}\right)}^{-\mathrm{1}}.\end{array}$$

3.6 CMAQ M3Dry

M3Dry (Pleim and Ran, 2011) is designed to couple with the Pleim–Xiu land surface model (PX LSM; Pleim and Xiu, 1995) in the Weather Research and Forecasting (WRF) model and is used operationally in CMAQ. There is also M3Dry-psn, which follows M3Dry but uses a coupled photosynthesis–stomatal conductance model. M3Dry-psn was developed and evaluated with the intention to supplement PX LSM and M3Dry in CMAQ (Ran et al., 2017). To date, however, M3Dry-psn has not been implemented in CMAQ. The parameters in Table S12 are site-specific.

3.6.1 Surface resistance

Surface resistance (r_c) is expressed as follows:

$$\begin{array}{}\text{(76)}& r_c={\left(\begin{array}{l}{f}_{\mathrm{veg}}\left(\frac{\mathrm{1}}{r_{\mathrm{st}}+r_{\mathrm{m}}}+\frac{\left(\mathrm{1}-f_{\mathrm{wet}}\right)\mathrm{LAI}}{r_{\mathrm{cut},\mathrm{dry}}}+\frac{f_{\mathrm{wet}}LAI}{r_{\mathrm{cut},\mathrm{wet}}}+\frac{\mathrm{1}}{r_{\mathrm{ac}}+r_{\mathrm{g}}}\right)\\ +\frac{\mathrm{1}-f_{\mathrm{veg}}}{r_{\mathrm{g}}}\end{array}\right)}^{-\mathrm{1}}.\end{array}$$

Here, the parameter f_veg is the fraction of the site covered by the vegetation canopy (Table S12) and f_wet is the fraction of canopy that is wet (Table 1).

3.6.2 Stomatal and mesophyll resistances

For M3Dry, stomatal resistance (r_st) follows Xiu and Pleim (2001):

$$\begin{array}{}\text{(77)}& r_{\mathrm{st}}=R_{\mathrm{diff},\mathrm{st}}{\displaystyle \frac{r_{\mathrm{i}}}{\mathrm{LAI}f\left(\mathrm{PAR}\right)f\left(w_{\mathrm{2}}\right)f\left({\mathrm{RH}}_l\right)f\left(T_{\mathrm{a}}\right)}}.\end{array}$$

The parameter r_i is initial resistance to stomatal uptake (Table S12).

The effects of photosynthetically active radiation (PAR) [µmol m⁻² s⁻¹] (Table 1) follow Echer and Rosolem (2015):

$$\begin{array}{}\text{(78)}& f\left(\mathrm{PAR}\right)=\left(\mathrm{1}-a\mathrm{LAI}\right)\left(\mathrm{1}-e^{-\mathrm{0.0017}\mathrm{PAR}}\right).\end{array}$$

The parameter a [unitless] is empirical (Table S12).

The effects of w₂ follow Xiu and Pleim (2001):

$$\begin{array}{}\text{(79)}& f\left(w_{\mathrm{2}}\right)={\left(\mathrm{1}+e^{-\mathrm{5}\left(\frac{w_{\mathrm{2}}-w_{\mathrm{wlt}}}{w_{\mathrm{fc}}-w_{\mathrm{wlt}}}-\left(\frac{w_{\mathrm{fc}}-w_{\mathrm{wlt}}}{\mathrm{3}}+w_{\mathrm{wlt}}\right)\right)}\right)}^{-\mathrm{1}}.\end{array}$$

The effects of leaf-level RH (RH_l) [fractional] are expressed as follows:

$$\begin{array}{}\text{(80)}& f\left({\mathrm{RH}}_l\right)={\mathrm{RH}}_l={\displaystyle \frac{q_{\mathrm{a}}{\left(r_{\mathrm{a}}+r_{\mathrm{b},\mathrm{v}}\right)}^{-\mathrm{1}}+q_{\mathrm{s}}{r}_{\mathrm{st},\mathrm{v}}^{-\mathrm{1}}}{\left(r_{\mathrm{st},\mathrm{v}}^{-\mathrm{1}}+{\left(r_a+r_{\mathrm{b},\mathrm{v}}\right)}^{-\mathrm{1}}\right)q_{\mathrm{s}}}}.\end{array}$$

Here, the variable q_a is the ambient air humidity mixing ratio, q_s is the saturation mixing ratio at leaf temperature (T_leaf), r_b,v is the quasi-laminar boundary layer resistance for water vapor, and r_st,v is the stomatal resistance for water vapor. M3Dry assumes that, when the sensible heat flux (SH) [W m⁻²] (Table 1) is greater than zero, the T_leaf equals $T_{\mathrm{a}}-\frac{\mathrm{SH}}{(r_a+r_{\mathrm{b},\mathrm{h}})\mathit{\rho }{c}_p}$, where r_b,h is quasi-laminar boundary layer resistance for heat. Otherwise, T_leaf equals T_a. Equation (80) is computed using an implicit quadratic solution as described by Xiu and Pleim (2001).

The effects of T_a are expressed as follows:

$$\begin{array}{}\text{(81)}& f\left(T_{\mathrm{a}}\right)=\left\{\begin{array}{l}{\left(\mathrm{1}+e^{-\mathrm{0.41}\left(T_{\mathrm{a}}-\mathrm{8.9}\right)}\right)}^{-\mathrm{1}},T_{\mathrm{a}}\le \mathrm{29}\\ {\left(\mathrm{1}+e^{\mathrm{0.5}\left(T_{\mathrm{a}}-\mathrm{40.85}\right)}\right)}^{-\mathrm{1}},T_{\mathrm{a}}>\mathrm{29}.\end{array}\right.\end{array}$$

For M3Dry-psn, r_st is simulated at the leaf level using the Ball–Woodrow–Berry approach (Ball et al., 1987), as described by Collatz et al. (1991, 1992) and Bonan et al. (2011):

$$\begin{array}{}\text{(82)}& r_{\mathrm{st}}={\left(g_{\mathrm{0}}+g_{\mathrm{1}}{\displaystyle \frac{A_{\mathrm{n}}}{\frac{p_{{\mathrm{CO}}_{\mathrm{2}},l}}{p_{\mathrm{a}}}}}{\mathrm{RH}}_{\mathrm{l}}\right)}^{-\mathrm{1}}{\displaystyle \frac{D_{{\mathrm{CO}}_{\mathrm{2}}}}{D_{{\mathrm{O}}_{\mathrm{3}}}}}{\displaystyle \frac{\mathrm{1000.0}\mathit{\rho }}{M_{\mathrm{air}}}}.\end{array}$$

Here, the parameter g₀ equals 0.01 mol CO₂ m⁻² s⁻¹ for C₃ plants; g₁ equals 9 [unitless]; A_n is the leaf-level net photosynthesis [mol CO₂ m⁻² s⁻¹]; $p_{{\mathrm{CO}}_{\mathrm{2}},l}$ is the carbon dioxide partial pressure at the leaf surface [Pa]; RH_l is the leaf-level RH [fractional], which follows Eq. (80) as described for M3Dry; $D_{{\mathrm{CO}}_{\mathrm{2}}}$ [m² s⁻¹] is the carbon dioxide diffusivity in air (Table 1); ρ [kg m⁻³] is the air density (Table 1); and M_air [g mol⁻¹] is the molar mass of air (Table 1). Leaf-level A_n is estimated based on Farquhar et al. (1980) as described by Ran et al. (2017), based on co-limitation among three potential assimilation rates, limited by Rubisco, light, and transport of photosynthetic products. The maximum rate of the carboxylation of Rubisco (V_cmax) [µmol m² s⁻¹] is key for A_n; thus, we include values at 25 ^∘C in Table S12.

Leaf-level A_n and r_st are calculated separately for sunlit versus shaded leaves in M3Dry-psn. Sunlit and shaded portions of the LAI (LAI_sun and LAI_shd, respectively) follow Campbell and Norman (1998) and Song et al. (2009). Canopy-scale r_st is expressed as follows:

$$\begin{array}{}\text{(83)}& r_{\mathrm{st}}={\left(\left({\displaystyle \frac{{\mathrm{LAI}}_{\mathrm{sun}}}{r_{\mathrm{st},\mathrm{sun}}}}+{\displaystyle \frac{{\mathrm{LAI}}_{\mathrm{shd}}}{r_{\mathrm{st},\mathrm{shd}}}}\right)f\left(w_{\mathrm{2}}\right)\right)}^{-\mathrm{1}}.\end{array}$$

The variables r_st,sun and r_st,shd are leaf-level stomatal resistances for sunlit and shaded leaves, respectively, calculated via Eq. (82). The function f(w₂) follows Eq. (79).

For both M3Dry and M3Dry-psn, mesophyll resistance (r_m) is expressed as follows:

$$\begin{array}{}\text{(84)}& r_{\mathrm{m}}={\displaystyle \frac{\mathrm{0.01}}{\mathrm{LAI}}}.\end{array}$$

3.6.3 Cuticular resistances

The variable r_cut,wet is the resistance to wet cuticles:

$$\begin{array}{}\text{(85)}& r_{\mathrm{cut},\mathrm{wet}}=\left\{\begin{array}{l}\mathrm{1250},T_g>\mathrm{0}\\ \mathrm{6667},T_g<\mathrm{0}.\end{array}\right.\end{array}$$

The variable T_g [^∘C] is ground temperature near the surface (Table 1).

The variable r_cut,dry is resistance to dry cuticles:

$$\begin{array}{}\text{(86)}& r_{\mathrm{cut},\mathrm{dry}}=r_{\mathrm{cut},\mathrm{dry},\mathrm{0}}\left(\mathrm{1}-f\left(\mathrm{RH}\right)\right)+r_{\mathrm{cut},\mathrm{wet}}f\left(\mathrm{RH}\right).\end{array}$$

The parameter $r_{\mathrm{cut},\mathrm{dry},\mathrm{0}}$ equals 2000 s m⁻¹.

The effects of RH are expressed as follows:

$$\begin{array}{}\text{(87)}& f\left(\mathrm{RH}\right)=max\left\{\mathrm{100}{\displaystyle \frac{\mathrm{RH}-\mathrm{0.7}}{\mathrm{0.3}}},\mathrm{0}\right\}.\end{array}$$

3.6.4 Resistance to the ground (and associated resistance to transport)

The resistance to in-canopy turbulence (r_ac) follows Erisman et al. (1994):

$$\begin{array}{}\text{(88)}& r_{\mathrm{ac}}=\mathrm{14}{\displaystyle \frac{h\mathrm{LAI}}{u_{*}}}.\end{array}$$

Ground resistance (r_g) is calculated as follows:

$$\begin{array}{}\text{(89)}& r_{\mathrm{g}}=\left\{\begin{array}{l}{\left(\frac{\mathrm{1}-f_{\mathrm{wet}}}{r_{\mathrm{g},\mathrm{dry}}}+\frac{f_{\mathrm{wet}}}{r_{\mathrm{g},\mathrm{wet}}}\right)}^{-\mathrm{1}},\mathrm{no}\mathrm{snow},\\ {\left(\frac{\mathrm{1}-X_{\mathrm{m}}}{r_{\mathrm{snow}}}+\frac{X_{\mathrm{m}}}{r_{\mathrm{sndiff}}+r_{\mathrm{g},\mathrm{wet}}}\right)}^{-\mathrm{1}},\mathrm{snow}.\end{array}\right.\end{array}$$

The variable r_g,wet is expressed as

$$\begin{array}{}\text{(90)}& r_{\mathrm{g},\mathrm{wet}}=\left\{\begin{array}{l}\mathrm{500},T_g>\mathrm{0}\\ \mathrm{6667},T_g<\mathrm{0}.\end{array}\right.\end{array}$$

The variable r_g,dry is expressed as follows (Massman, 2004; Mészáros et al., 2009):

$$\begin{array}{}\text{(91)}& r_{\mathrm{g},\mathrm{dry}}=\mathrm{200}+\left(r_{\mathrm{g},\mathrm{wet}}-\mathrm{200}\right){\displaystyle \frac{w_g}{w_{\mathrm{fc}}}}.\end{array}$$

If the near-surface soil water content (w_g) [m³ m⁻³] (Table 1) is greater than w_fc, the soil is wet (i.e., r_g,dry equals r_g,wet). The parameter r_snow is resistance to snow or ice [6667 s m⁻¹] and r_sndiff is resistance to diffusion through the snowpack [10 s m⁻¹]. Parallel pathways to frozen snow/ice and diffusion through the snowpack to liquid water follow Bales et al. (1987). Snow liquid water mass (X_m) is calculated as follows:

$$\begin{array}{}\text{(92)}& X_{\mathrm{m}}=\left\{\begin{array}{l}max\left\{\mathrm{0.02}{\left(T_{\mathrm{a}}+\mathrm{1}\right)}^{\mathrm{2}},\mathrm{0.5}\right\},T_{\mathrm{a}}>-\mathrm{1}\\ \mathrm{0},T_{\mathrm{a}}<-\mathrm{1}.\end{array}\right.\end{array}$$

3.7 CMAQ STAGE

The Surface Tiled Aerosol and Gaseous Exchange (STAGE) parameterization is an option in CMAQ. Parameters in Table S13 are site-specific.

3.7.1 Deposition velocity

The ozone deposition velocity (v_d) follows:

$$\begin{array}{}\text{(93)}& \begin{array}{rl}{v}_{\mathrm{d}}& =f_{\mathrm{veg}}{\left(r_{\mathrm{a}}+{\displaystyle \frac{\mathrm{1}}{\frac{\mathrm{1}}{r_{\mathrm{b},\mathrm{v}}+\frac{\mathrm{1}}{\frac{\mathrm{1}}{r_{\mathrm{st}}+r_m}+\frac{\mathrm{1}}{r_{\mathrm{cut}}}}}+\frac{\mathrm{1}}{r_{\mathrm{ac}}+r_{\mathrm{b},\mathrm{g}}+r_{\mathrm{g}}}}}\right)}^{-\mathrm{1}}\\ & +\left(\mathrm{1}-f_{\mathrm{veg}}\right){\left(r_{\mathrm{a}}+r_{\mathrm{b},\mathrm{g}}+r_{\mathrm{g}}\right)}^{-\mathrm{1}}\end{array}\end{array}$$

CMAQ STAGE considers separate quasi-laminar boundary layer resistances around vegetation versus the ground (r_b,v and r_b,g, respectively) (Table S3). The parameter f_veg is the vegetated fraction of the site; the M3Dry value is used (Table S12).

3.7.2 Stomatal and mesophyll resistances

Stomatal resistance (r_st) follows Pleim and Ran (2011):

$$\begin{array}{}\text{(94)}& r_{\mathrm{st}}=R_{\mathrm{diff},\mathrm{st}}{\displaystyle \frac{r_{\mathrm{i}}}{\mathrm{LAI}f\left(\mathrm{PAR}\right)f\left(w_{\mathrm{2}}\right)f\left({\mathrm{RH}}_l\right)f\left(T_{\mathrm{a}}\right)}}.\end{array}$$

The parameter r_i is initial resistance to stomatal uptake (Table S13). The functions follow M3Dry (Eqs. 78–81).

Mesophyll resistance (r_m) follows Wesely (1989):

$$\begin{array}{}\text{(95)}& r_{\mathrm{m}}={\left({\displaystyle \frac{H}{\mathrm{3000}}}+\mathrm{100}{f}_{\mathrm{0}}\right)}^{-\mathrm{1}}.\end{array}$$

3.7.3 Cuticular resistance

Cuticular resistance (r_cut) is expressed as follows:

$$\begin{array}{}\text{(96)}& r_{\mathrm{cut}}={\left(\mathrm{LAI}\left({\displaystyle \frac{f_{\mathrm{wet}}}{\mathrm{1250}}}+{\displaystyle \frac{\mathrm{1}-f_{\mathrm{wet}}}{\mathrm{2000}}}\right)\right)}^{-\mathrm{1}}.\end{array}$$

3.7.4 Resistance to the ground (and associated resistance to transport)

The resistance to in-canopy turbulence (r_ac) is similar to Shuttleworth and Wallace (1985):

$$\begin{array}{}\text{(97)}& r_{\mathrm{ac}}={\int }_{\mathrm{0}}^h{\displaystyle \frac{\mathrm{d}z}{K_t}}.\end{array}$$

The variable K_t is in-canopy eddy diffusivity [m² s⁻¹]. By applying the drag coefficient ($C_{\mathrm{d}}=\frac{u_{*}^{\mathrm{2}}}{u^{\mathrm{2}}}$), assuming a uniform vertical distribution of leaves, and using an in-canopy attenuation coefficient of momentum following Yi (2008) $\left[\frac{\mathrm{LAI}}{\mathrm{2}}\right]$, we obtain the following:

$$\begin{array}{}\text{(98)}& r_{\mathrm{ac}}=\mathit{Pr}{\displaystyle \frac{u}{u_{*}^{\mathrm{2}}}}\left(e^{\frac{\mathrm{LAI}}{\mathrm{2}}}-\mathrm{1}\right)=r_{\mathrm{a}}\left(e^{\frac{\mathrm{LAI}}{\mathrm{2}}}-\mathrm{1}\right).\end{array}$$

The variable u [m s⁻¹] is wind speed (Table 1).

The resistance to the ground (r_g) changes whether the ground is snow-covered, dry, or wet (wet is w_g greater than or equal to w_sat, where w_sat [m³ m⁻³] is soil water content at saturation; Table 1). For dry ground, r_g follows Fares et al. (2014) and Fumagalli et al. (2016). An asymptotic function bounds the resistance, following observations reported in Fumagalli et al. (2016):

$$\begin{array}{}\text{(99)}& r_{\mathrm{g}}=\left\{\begin{array}{l}\mathrm{250}+\mathrm{2000}\mathrm{atan}\left(\frac{{\left(\frac{w_g-w_{\mathrm{wlt}}}{w_{\mathrm{fc}}}\right)}^B}{\mathit{\pi }}\right),w<w_{\mathrm{sat}}\\ \frac{\mathrm{62}\mathrm{500}}{HR(T_g+\mathrm{273.15})},w\ge w_{\mathrm{sat}}\\ \frac{\mathrm{1}-X_{\mathrm{m}}}{r_{\mathrm{snow}}}+\frac{X_{\mathrm{m}}}{r_{\mathrm{sndiff}}+\frac{\mathrm{62}\mathrm{500}}{HR(T_g+\mathrm{273.15})}},\mathrm{snow}.\end{array}\right.\end{array}$$

Here, the parameter R [L atm K⁻¹ mol⁻¹] is the universal gas constant, B [unitless] is an empirical parameter related to soil moisture (Table 1), r_snow is resistance to snow or ice [6667 s m⁻¹], and r_sndiff is resistance to diffusion through the snowpack [10 s m⁻¹]. The liquid fraction of the quasi-liquid layer in snow (X_m) is modeled as a system dominated by van der Waals forces using the temperature parameterization following Huthwelker et al. (2006) and assuming a maximum of 20 % to match gas–liquid partitioning findings in Conklin et al. (1993):

$$\begin{array}{}\text{(100)}& X_{\mathrm{m}}=\left\{\begin{array}{cc}\frac{\mathrm{0.025}}{{\left(\mathrm{273.15}-T_g\right)}^{\mathrm{1}/\mathrm{3}}},& \mathrm{0.002}<\mathrm{273.15}-T_g<\mathrm{10}\\ \mathrm{0.2},& \mathrm{273.15}-T_g<\mathrm{0.002}.\end{array}\right.\end{array}$$

3.8 TEMIR

The Terrestrial Ecosystem Model in R (TEMIR) (Tai et al., 2023) provides two dry deposition schemes (Sun et al., 2022): Wesely and Zhang. Wesely in TEMIR largely follows GEOS-Chem version 12.0.0, while Zhang follows Zhang et al. (2003). In both schemes, the default stomatal resistance is highly empirical. TEMIR can also use two photosynthesis-based stomatal conductance models (hereinafter, “psn”): the Farquhar–Ball–Berry model (hereinafter, “BB”; Farquhar et al., 1980; Ball et al., 1987) and the Medlyn et al. (2011) model (hereinafter, “Medlyn”). Thus, three stomatal conductance models are used for TEMIR Wesely and TEMIR Zhang, respectively. TEMIR Zhang parameters in Table S14 and TEMIR psn parameters in Table S15 are site-specific.

3.8.1 Surface resistance

For Wesely, surface resistance (r_c) is calculated as follows:

$$\begin{array}{}\text{(101)}& r_c={\left({\displaystyle \frac{\mathrm{1}}{r_{\mathrm{st}}}}+{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{cut}}}}+{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{dc}}+r_{\mathrm{cl}}}}+{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{ac}}+r_{\mathrm{g}}}}\right)}^{-\mathrm{1}}.\end{array}$$

For Zhang, surface resistance (r_c) is expressed as follows:

$$\begin{array}{}\text{(102)}& r_c={\left({\displaystyle \frac{\mathrm{1}-W}{r_{\mathrm{st}}}}+{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{cut}}}}+{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{ac}}+r_{\mathrm{g}}}}\right)}^{-\mathrm{1}}.\end{array}$$

The parameter W [fractional] is used to account for leaf wetness. If P is greater than 0.2 mm h⁻¹,

$$\begin{array}{}\text{(103)}& W=\left\{\begin{array}{l}\mathrm{0},G\le \mathrm{200}\\ \frac{G-\mathrm{200}}{\mathrm{800}},\mathrm{200}\le G\le \mathrm{600}\\ \mathrm{0.5},G>\mathrm{600}.\end{array}\right.\end{array}$$

3.8.2 Stomatal resistance

For Wesely, stomatal resistance (r_st) is expressed as follows:

$$\begin{array}{}\text{(104)}& r_{\mathrm{st}}=R_{\mathrm{diff},\mathrm{st}}{\displaystyle \frac{r_{\mathrm{i}}}{{\mathrm{LAI}}_{\mathrm{eff}}f\left(T_{\mathrm{a}}\right)}}.\end{array}$$

Here, the parameter r_i is initial resistance to stomatal uptake (same for GEOS-Chem Wesely; Table S6), and LAI_eff [m² m⁻²] is the effective LAI, which is the surface area of actively transpiring leaves per ground surface area. The variable LAI_eff is calculated as a function of the LAI, θ, and the cloud fraction using a parameterization developed by Wang et al. (1998). In GEOS-Chem, if G is 0, LAI_eff equals 0.01. For the single-point model, we set G to 0 when θ is greater than 95^∘ so that nighttime r_st values in the single-point model are more similar GEOS-Chem. GEOS-Chem almost never has nonzero G at night, but measured values are frequently small and nonzero. Here, the cloud fraction is assumed to be zero.

The effects of T_a are expressed as follows:

$$\begin{array}{}\text{(105)}& f\left(T_{\mathrm{a}}\right)=\left\{T_{\mathrm{a}}\begin{array}{l}\mathrm{0.01},T_{\mathrm{a}}\le \mathrm{0}\\ \frac{\left(\mathrm{40}-T_{\mathrm{a}}\right)}{\mathrm{400}},\mathrm{0}<T_{\mathrm{a}}<\mathrm{40}\\ \mathrm{0.01},\mathrm{40}\le T_{\mathrm{a}}.\end{array}\right.\end{array}$$

For Zhang, stomatal resistance (r_st) is calculated as follows:

$$\begin{array}{}\text{(106)}& r_{\mathrm{st}}=R_{\mathrm{diff},\mathrm{st}}{\displaystyle \frac{r_{\mathrm{i}}(\mathrm{LAI},\mathrm{PAR})}{f\left(T_{\mathrm{a}}\right)f\left(\mathrm{VPD}\right)f\left({\mathit{\psi }}_{\mathrm{leaf}}\right).}}\end{array}$$

Dependencies on T_a, VPD, and ψ_leaf are as described in Brook et al. (1999).

The variable r_i (LAI, PAR) is expressed as follows:

$$\begin{array}{}\text{(107)}& r_{\mathrm{i}}(\mathrm{LAI},\mathrm{PAR})={\left({\displaystyle \frac{{\mathrm{LAI}}_{\mathrm{sun}}}{r_{\mathrm{i}}\left(\mathrm{1}+\frac{b_{\mathrm{rs}}}{{\mathrm{PAR}}_{\mathrm{sun}}}\right)}}+{\displaystyle \frac{{\mathrm{LAI}}_{\mathrm{shd}}}{r_{\mathrm{i}}\left(\mathrm{1}+\frac{b_{\mathrm{rs}}}{{\mathrm{PAR}}_{\mathrm{shd}}}\right)}}\right)}^{-\mathrm{1}}.\end{array}$$

Here, the parameter r_i is the initial resistance to stomatal uptake (Table S14), b_rs [W m⁻²] is empirical (Table S14), and LAI_sun and LAI_shd [m² m⁻²] are sunlit and shaded LAI, respectively. The latter two parameters are expressed as follows:

$$\begin{array}{}\text{(108)}& {\displaystyle }{\mathrm{LAI}}_{\mathrm{sun}}={\displaystyle \frac{\mathrm{1}-e^{-K_b\mathrm{LAI}}}{K_b}},\text{(109)}& {\displaystyle }{\mathrm{LAI}}_{\mathrm{shd}}=\mathrm{LAI}-{\mathrm{LAI}}_{\mathrm{sun}}.\end{array}$$

The variable K_b is the canopy light extinction coefficient [unitless]:

$$\begin{array}{}\text{(110)}& K_b={\displaystyle \frac{\mathrm{0.5}}{\cos \left(\frac{\mathit{\pi }}{\mathrm{180}}\mathit{\theta }\right)}}.\end{array}$$

The variables PAR_sun and PAR_shd [W m⁻²] are PAR reaching sunlit and shaded leaves, respectively:

$$\begin{array}{}\text{(111)}& {\displaystyle }\begin{array}{rl}{\mathrm{PAR}}_{\mathrm{shd}}& =R_{\mathrm{diff}}{e}^{-\mathrm{0.5}{\mathrm{LAI}}^a}\\ & +\mathrm{0.07}{R}_{\mathrm{dir}}\left(\mathrm{1.1}-\mathrm{0.1}\mathrm{LAI}\right)e^{-\cos \left(\frac{\mathit{\pi }}{\mathrm{180}}\mathit{\theta }\right),}\end{array}\text{(112)}& {\displaystyle }{\mathrm{PAR}}_{\mathrm{sun}}={\mathrm{PAR}}_{\mathrm{shd}}+{\displaystyle \frac{R_{\mathrm{dir}}^b\cos \left(\frac{\mathit{\pi }}{\mathrm{180}}\mathit{\alpha }\right)}{\cos \left(\frac{\mathit{\pi }}{\mathrm{180}}\mathit{\theta }\right).}}\end{array}$$

Here, the parameter α is the angle between the leaf and the sun [60^∘], and R_diff and R_dir are downward visible radiation fluxes from diffuse and direct-beam radiation above the canopy, respectively. We use the diffuse fraction from the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) reanalysis product (GMAO, 2015) to separate R_diff and R_dir from observed PAR. If the LAI is less than 2.5 m² m⁻² or G is less than 200 W m⁻², a equals 0.7 and b equals 1. Otherwise, a equals 0.8 and b equals 0.8.

The effects of T_a are expressed as follows:

$$\begin{array}{}\text{(113)}& f\left(T_{\mathrm{a}}\right)=\left({\displaystyle \frac{T_{\mathrm{a}}-T_{\min }}{T_{\mathrm{opt}}-T_{\min }}}\right){\left({\displaystyle \frac{T_{\max }-T_{\mathrm{a}}}{T_{\max }-T_{\mathrm{opt}}}}\right)}^{\frac{T_{\max }-T_{\mathrm{opt}}}{T_{\mathrm{opt}}-T_{\min }}}.\end{array}$$

The parameters T_min, T_max, and T_opt [^∘C] are the minimum, maximum, and optimum temperatures, respectively (Table S14).

The effects of VPD are expressed as follows:

$$\begin{array}{}\text{(114)}& f\left(\mathrm{VPD}\right)=\mathrm{1}-b_{\mathrm{VPD}}\mathrm{VPD}.\end{array}$$

The parameter b_VPD [kPa⁻¹] is empirical (Table S14).

The effects of ψ_leaf are calculated as follows:

$$\begin{array}{}\text{(115)}& f\left({\mathit{\psi }}_{\mathrm{leaf}}\right)={\displaystyle \frac{{\mathit{\psi }}_{\mathrm{leaf}}-{\mathit{\psi }}_{\mathrm{leaf},\mathrm{2}}}{{\mathit{\psi }}_{\mathrm{leaf},\mathrm{1}}-{\mathit{\psi }}_{\mathrm{leaf},\mathrm{2}}}}.\end{array}$$

The parameters ψ_leaf,1 and ψ_leaf,2 [MPa] are empirical (Table S14), whereas ψ_leaf is parameterized as follows:

$$\begin{array}{}\text{(116)}& {\mathit{\psi }}_{\mathrm{leaf}}=-\mathrm{0.72}-\mathrm{0.0013}G.\end{array}$$

We now describe psn options for TEMIR Wesely and TEMIR Zhang. For BB (Ball et al., 1987; Farquhar et al., 1980; von Caemmerer and Farquhar, 1981; Collatz et al., 1991, 1992),

$$\begin{array}{}\text{(117)}& r_{\mathrm{st}}={\left({\mathit{\beta }}_t{g}_{\mathrm{0}}+g_{\mathrm{1}}{\displaystyle \frac{A_{\mathrm{n}}\mathrm{RH}}{\frac{p_{{\mathrm{CO}}_{\mathrm{2}},l}}{p_{\mathrm{a}}}}}\right)}^{-\mathrm{1}}{\displaystyle \frac{p_{\mathrm{a}}}{R{\mathit{\theta }}_a}}.\end{array}$$

Here, the parameter g₀ equals 0.01 mol m⁻² s⁻¹, g₁ equals 9, A_n is net photosynthesis [mol m⁻² s⁻¹], β_t is a soil water stress factor [unitless], $p_{{\mathrm{CO}}_{\mathrm{2}},l}$ is the carbon dioxide partial pressure at the leaf surface [Pa], R is the universal gas constant [J mol⁻¹ K⁻¹], and θ_a is potential air temperature [K].

For Medlyn (Medlyn et al., 2011),

$$\begin{array}{}\text{(118)}& r_{\mathrm{st}}={\left({\mathit{\beta }}_t{g}_{\mathrm{0}}+{\displaystyle \frac{D_{\mathrm{w}}}{D_{{\mathrm{CO}}_{\mathrm{2}}}}}\left(\mathrm{1}+{\displaystyle \frac{g_{\mathrm{1}M}}{\sqrt{\mathrm{VPD}}}}\right){\displaystyle \frac{A_{\mathrm{n}}}{\frac{p_{{\mathrm{CO}}_{\mathrm{2},l}}}{p_{\mathrm{a}}}}}\right)}^{-\mathrm{1}}{\displaystyle \frac{p_{\mathrm{a}}}{R{\mathit{\theta }}_a}}.\end{array}$$

Here, the parameter g_1M [kPa^0.5] is empirical (Table S15), g₀ equals 0.0001 mol m⁻² s⁻¹, D_w [m² s⁻¹] is the diffusivity of water vapor in air (Table 1), and the ratio of diffusivities is 1.6.

A single-layer bulk soil formulation considering the root zone (0–100 cm) is used to calculate β_t:

$$\begin{array}{}\text{(119)}& {\mathit{\beta }}_t=\left\{\begin{array}{l}\mathrm{1},{\mathit{\psi }}_{\mathrm{soil}}>{\mathit{\psi }}_{\mathrm{soil},\mathrm{fc}}\\ \frac{{\mathit{\psi }}_{\mathrm{soil},\mathrm{wlt}}-{\mathit{\psi }}_{\mathrm{soil}}}{{\mathit{\psi }}_{\mathrm{soil},\mathrm{wlt}}-{\mathit{\psi }}_{\mathrm{soil},\mathrm{fc}}},{\mathit{\psi }}_{\mathrm{soil},\mathrm{wlt}}\le {\mathit{\psi }}_{\mathrm{soil}}\le {\mathit{\psi }}_{\mathrm{soil},\mathrm{fc}}\\ \mathrm{0},{\mathit{\psi }}_{\mathrm{soil}}<{\mathit{\psi }}_{\mathrm{soil},\mathrm{fc}}.\end{array}\right.\end{array}$$

The variable ψ_soil [kPa] is soil matric potential (Table 1):

$$\begin{array}{}\text{(120)}& {\mathit{\psi }}_{\mathrm{soil}}={\mathit{\psi }}_{\mathrm{soil},\mathrm{sat}}{w}_{\mathrm{2}}^{-B}.\end{array}$$

For both Medlyn and BB, leaf-level r_st is calculated individually for sunlit and shaded leaves and is then scaled up:

$$\begin{array}{}\text{(121)}& r_{\mathrm{st}}=R_{\mathrm{diff},\mathrm{st}}{\left({\displaystyle \frac{{\mathrm{LAI}}_{\mathrm{sun}}}{r_{\mathrm{b},\mathrm{leaf}}+r_{\mathrm{st},\mathrm{sun}}}}+{\displaystyle \frac{{\mathrm{LAI}}_{\mathrm{shd}}}{r_{\mathrm{b},\mathrm{leaf}}+r_{\mathrm{st},\mathrm{shd}}}}\right)}^{-\mathrm{1}}.\end{array}$$

The variables r_st,sun and r_st,shd are leaf-level stomatal resistances for sunlit and shaded leaves, respectively; LAI_sun and LAI_shd are sunlit and shaded LAI, respectively; and r_b,leaf is leaf boundary layer resistance:

$$\begin{array}{}\text{(122)}& r_{\mathrm{b},\mathrm{leaf}}={\displaystyle \frac{\mathrm{1}}{c_v}}\sqrt{{\displaystyle \frac{u_{*}}{l}}}.\end{array}$$

The parameter c_v [0.01 m s^−0.5] is the turbulent transfer coefficient and l [0.04 m] is the characteristic dimension of leaves.

The variables LAI_sun and LAI_shd are expressed as follows:

$$\begin{array}{}\text{(123)}& {\displaystyle }{\mathrm{LAI}}_{\mathrm{sun}}={\mathrm{PAI}}_{\mathrm{sun}}{\displaystyle \frac{\mathrm{LAI}}{\mathrm{LAI}+\mathrm{SAI}}},\text{(124)}& {\displaystyle }{\mathrm{LAI}}_{\mathrm{shd}}={\mathrm{PAI}}_{\mathrm{shd}}{\displaystyle \frac{\mathrm{LAI}}{\mathrm{LAI}+\mathrm{SAI}}}.\end{array}$$

Here, the variable SAI [m² m⁻²] is the stem area index, and PAI_sun and PAI_shd [m² m⁻²] are the sunlit and shaded plant area index, respectively. The latter two variables are expressed as follows:

$$\begin{array}{}\text{(125)}& {\displaystyle }{\mathrm{PAI}}_{\mathrm{sun}}={\displaystyle \frac{\mathrm{1}-e^{-K_b\left(\mathrm{LAI}+\mathrm{SAI}\right)}}{K_b}},\text{(126)}& {\displaystyle }{\mathrm{PAI}}_{\mathrm{shd}}=\mathrm{LAI}+\mathrm{SAI}-{\mathrm{PAI}}_{\mathrm{sun}}.\end{array}$$

The variable SAI follows Zeng et al. (2002):

$$\begin{array}{}\text{(127)}& {\mathrm{SAI}}_n=\max \mathit{\{}\mathrm{0.5}{\mathrm{SAI}}_{n-\mathrm{1}}+max\left\{{\mathrm{LAI}}_{n-\mathrm{1}}-{\mathrm{LAI}}_n,\mathrm{0}\right\},\mathrm{1}\mathit{\}}.\end{array}$$

The parameter n is nth month of the year.

Leaf-level photosynthesis of C₃ plants is represented by the formulation that relates to Michaelis–Menten enzyme kinetics and photosynthetic biochemical pathways, as in the Community Land Model 4.5 (CLM4.5) (Oleson et al., 2013) and following Collatz et al. (1992):

$$\begin{array}{}\text{(128)}& A_{\mathrm{n}}=min\left\{A_c,A_j,A_{\mathrm{p}}\right\}-R_{\mathrm{d}}.\end{array}$$

The Rubisco-limited photosynthetic rate (A_c) [mol m⁻² s⁻¹] is expressed as follows:

$$\begin{array}{}\text{(129)}& A_c=V_{\mathrm{cmax}}{\displaystyle \frac{c_{\mathrm{i}}-{\mathrm{\Gamma }}_{*}}{c_{\mathrm{i}}+K_{\mathrm{c}}\left(\mathrm{1}+\frac{o_{\mathrm{i}}}{K_{\mathrm{o}}}\right)}}.\end{array}$$

Here, the variable c_i is the intercellular carbon dioxide partial pressure [Pa]; K_c and K_o are Michaelis–Menten constants for carboxylation and oxygenation [Pa], respectively; o_i is the intercellular oxygen partial pressure [0.029 p_a Pa]; Γ_* is the carbon dioxide compensation point [Pa]; and V_cmax is the maximum rate of carboxylation [mol m⁻² s⁻¹] adjusted for leaf temperature. The latter variable is calculated as follows:

$$\begin{array}{}\text{(130)}& V_{\mathrm{cmax}}=V_{\mathrm{cmax},\mathrm{25}}f\left(T_{\mathrm{l}}\right)f_H\left(T_{\mathrm{l}}\right){\mathit{\beta }}_t.\end{array}$$

The parameter V_cmax,25 is the value of V_cmax at 25 ^∘C (Table S15).

The function of leaf temperature (T_l) [K] is expressed as follows:

$$\begin{array}{}\text{(131)}& f\left(T_{\mathrm{l}}\right)=e^{\frac{\mathrm{\Delta }{H}_a}{\mathrm{298.15}\cdot \mathrm{0.001}R}\left(\mathrm{1}-\frac{\mathrm{298.15}}{T_{\mathrm{l}}}\right)}.\end{array}$$

The parameter R is the universal gas constant [J kg⁻¹ K⁻¹]. The high-temperature function of T_l is calculated as follows:

$$\begin{array}{}\text{(132)}& f_H\left(T_{\mathrm{l}}\right)={\displaystyle \frac{\mathrm{1}+e^{\frac{\mathrm{298.15}\mathrm{\Delta }S-\mathrm{\Delta }{H}_{\mathrm{d}}}{\mathrm{298.15}\cdot \mathrm{0.001}R}}}{\mathrm{1}+e^{\frac{\mathrm{\Delta }{\mathrm{ST}}_v-\mathrm{\Delta }{H}_{\mathrm{d}}}{\mathrm{0.001}R{T}_{\mathrm{l}}}}}}.\end{array}$$

The variables ΔH_a [J mol⁻¹], ΔS [J mol⁻¹ K⁻¹], and ΔH_d [J mol⁻¹] are temperature-dependent and follow the definitions in CLM4.5 (see Table S15 for the CLM4.5 plant functional types used for each site).

The ribulose-1,5-bisphosphate (RuBP)-limited photosynthetic rate (A_j) [mol m⁻² s⁻¹] is expressed as follows:

$$\begin{array}{}\text{(133)}& A_j={\displaystyle \frac{J}{\mathrm{4}}}{\displaystyle \frac{c_{\mathrm{i}}-{\mathrm{\Gamma }}_{*}}{c_{\mathrm{i}}+\mathrm{2}{\mathrm{\Gamma }}_{*}}}.\end{array}$$

The parameter J is the electron transport rate [mol m⁻² s⁻¹], taken as the smaller of the two roots of the equation below:

$$\begin{array}{}\text{(134)}& {\displaystyle }{\mathit{\theta }}_{\mathrm{PSII}}{J}^{\mathrm{2}}-\left(I_{\mathrm{PSII}}+J_{\max }\right)J+I_{\mathrm{PSII}}{J}_{\max }=\mathrm{0},\text{(135)}& {\displaystyle }{J}_{\max }=\mathrm{1.97}{V}_{\mathrm{cmax},\mathrm{25}}f\left(T_{\mathrm{l}}\right)f_H\left(T_{\mathrm{l}}\right),\text{(136)}& {\displaystyle }{I}_{\mathrm{PSII}}=\mathrm{0.5}{\mathrm{\Phi }}_{\mathrm{PSII}}\mathrm{4.6}\times {\mathrm{10}}^{-\mathrm{6}}\mathit{\varphi }.\end{array}$$

Here, the parameter θ_PSII [unitless] represents curvature; I_PSII [mol m⁻² s⁻¹] is the light utilization in electron transport by photosystem II; J_max [mol m⁻² s⁻¹] is the potential maximum electron transport rate; Φ_PSII [unitless] is the quantum yield of photosystem II; and ϕ [W m⁻²] is the photosynthetically active radiation absorbed by leaves, converted to photosynthetic photon flux density with $\mathrm{4.6}\times {\mathrm{10}}^{-\mathrm{6}}$ mol J⁻¹.

The product-limited photosynthetic rate (A_p) [mol m⁻² s⁻¹] is calculated as follows:

$$\begin{array}{}\text{(137)}& A_{\mathrm{p}}=\mathrm{3}{T}_{\mathrm{p}}.\end{array}$$

The parameter T_p is the triose phosphate utilization rate [mol m⁻² s⁻¹]:

$$\begin{array}{}\text{(138)}& T_{\mathrm{p}}=\mathrm{0.167}{V}_{\mathrm{cmax},\mathrm{25}}f\left(T_{\mathrm{l}}\right)f_H\left(T_{\mathrm{l}}\right).\end{array}$$

Dark respiration (R_d) [mol m⁻² s⁻¹] is expressed as follows:

$$\begin{array}{}\text{(139)}& R_{\mathrm{d}}=\mathrm{0.015}{V}_{\mathrm{cmax},\mathrm{25}}f\left(T_{\mathrm{l}}\right)f_H\left(T_{\mathrm{l}}\right){\mathit{\beta }}_t.\end{array}$$

The calculation for A_n and r_st involves a coupled set of equations that are solved iteratively at each time step until c_i converges (see Sect. 8.5 of Oleson et al., 2013):

$$\begin{array}{}\text{(140)}& \begin{array}{rl}{A}_{\mathrm{n}}& ={\displaystyle \frac{p_{{\mathrm{CO}}_{\mathrm{2}},a}-p_{{\mathrm{CO}}_{\mathrm{2}},i}}{\left(\mathrm{1.4}{r}_{\mathrm{b},\mathrm{leaf}}+\frac{D_{\mathrm{w}}}{D_{{\mathrm{CO}}_{\mathrm{2}}}}{r}_{\mathrm{st}}\right)p_{\mathrm{a}}}}={\displaystyle \frac{p_{{\mathrm{CO}}_{\mathrm{2}},\mathrm{a}}-p_{{\mathrm{CO}}_{\mathrm{2}},\mathrm{l}}}{\mathrm{1.4}{r}_{\mathrm{b},\mathrm{leaf}}{p}_{\mathrm{a}}}}\\ & ={\displaystyle \frac{p_{{\mathrm{CO}}_{\mathrm{2}},\mathrm{l}}-p_{{\mathrm{CO}}_{\mathrm{2}},\mathrm{i}}}{\frac{D_{\mathrm{w}}}{D_{{\mathrm{CO}}_{\mathrm{2}}}}{r}_{\mathrm{st}}{p}_{\mathrm{a}}}}.\end{array}\end{array}$$

Here, the variables $p_{{\mathrm{CO}}_{\mathrm{2}},\mathrm{a}}$, $p_{{\mathrm{CO}}_{\mathrm{2}},\mathrm{l}}$, and $p_{{\mathrm{CO}}_{\mathrm{2}},\mathrm{i}}$ are the carbon dioxide partial pressure [Pa] in air, at the leaf level, and in the intercellular space, respectively.

3.8.3 Cuticular resistance

For Wesely, cuticular resistance (r_cut) is expressed as follows:

$$\begin{array}{}\text{(141)}& r_{\mathrm{cut}}=\left\{\begin{array}{l}{r}_{\mathrm{lu}}min\mathit{\{}\mathrm{2},e^{\mathrm{0.2}(-\mathrm{1}-T_{\mathrm{a}})}\mathit{\}}{\left(\frac{H}{{\mathrm{10}}^{\mathrm{5}}}+f_{\mathrm{0}}\right)}^{-\mathrm{1}},T_{\mathrm{a}}<-\mathrm{1}\\ \left(\frac{r_{\mathrm{lu}}}{\mathrm{LAI}}+\mathrm{1000}{e}^{-T_{\mathrm{a}}-\mathrm{4}}\right){\left(\frac{H}{{\mathrm{10}}^{\mathrm{5}}}+f_{\mathrm{0}}\right)}^{-\mathrm{1}},T_{\mathrm{a}}\ge -\mathrm{1}.\end{array}\right.\end{array}$$

The parameter r_lu is the initial resistance for cuticular uptake. Values follow GEOS-Chem Wesely (Table S6).

For Zhang, cuticular resistance (r_cut) is calculated as follows:

$$\begin{array}{}\text{(142)}& r_{\mathrm{cut}}=\left\{\begin{array}{l}\frac{c_{\mathrm{cut},\mathrm{dry}}}{u^{*}{\mathrm{LAI}}^{\mathrm{0.25}}{e}^{\mathrm{3}\mathrm{RH}}},\mathrm{dry}\\ \frac{c_{\mathrm{cut},\mathrm{wet}}}{u^{*}{\mathrm{LAI}}^{\mathrm{0.5}}},\mathrm{wet}.\end{array}\right.\end{array}$$

The parameters c_cut,dry and c_cut,wet [unitless] are empirical coefficients related to dry and wet cuticular uptake, respectively (Table S14). If P is greater than 0.2 mm h⁻¹, cuticles are wet; otherwise, cuticles are dry.

The variable r_cut is adjusted for snow:

$$\begin{array}{}\text{(143)}& r_{\mathrm{cut}}={\left({\displaystyle \frac{\mathrm{1}-f_{\mathrm{snow}}}{r_{\mathrm{cut}}}}+{\displaystyle \frac{\mathrm{2}{f}_{\mathrm{snow}}}{\mathrm{2000}}}\right)}^{-\mathrm{1}}.\end{array}$$

3.8.4 Resistances to the lower canopy and ground (and associated resistances to transport)

For Wesely, the resistance associated with in-canopy convection (r_dc) is calculated as follows:

$$\begin{array}{}\text{(144)}& r_{\mathrm{dc}}=\mathrm{100}\left(\mathrm{1}+{\displaystyle \frac{\mathrm{1000}}{G+\mathrm{10}}}\right).\end{array}$$

The resistance to the lower canopy (r_cl) is expressed as follows:

$$\begin{array}{}\text{(145)}& r_{\mathrm{cl}}={\left({\displaystyle \frac{H}{{\mathrm{10}}^{\mathrm{5}}{r}_{\mathrm{cl},\mathrm{S}}}}+{\displaystyle \frac{f_{\mathrm{0}}}{r_{\mathrm{cl},\mathrm{O}}}}\right)}^{-\mathrm{1}}.\end{array}$$

The parameters r_cl,S and r_cl,O are initial resistances to uptake by the lower canopy and follow GEOS-Chem Wesely (Table S6).

Resistance to the ground (r_g) is calculated as follows:

$$\begin{array}{}\text{(146)}& r_{\mathrm{g}}={\left({\displaystyle \frac{H}{{\mathrm{10}}^{\mathrm{5}}{r}_{\mathrm{g},\mathrm{S}}}}+{\displaystyle \frac{f_{\mathrm{0}}}{r_{\mathrm{g},\mathrm{O}}}}\right)}^{-\mathrm{1}}.\end{array}$$

The parameters r_g,S and r_g,O are initial resistances to the ground and follow GEOS-Chem Wesely (Table S6). The resistance to turbulent transport to the ground (r_ac) follows GEOS-Chem Wesely (Table S6). The changes in resistances when there is snow follow GEOS-Chem Wesely (Table S6).

For Zhang, in-canopy aerodynamic resistance (r_ac) is expressed as follows:

$$\begin{array}{}\text{(147)}& r_{\mathrm{ac}}=r_{\mathrm{ac}\mathrm{0}}{\displaystyle \frac{{\mathrm{LAI}}^{\mathrm{0.25}}}{{\left(u^{*}\right)}^{\mathrm{2}}}}.\end{array}$$

The variable r_ac0 is calculated as follows:

$$\begin{array}{}\text{(148)}& r_{\mathrm{ac}\mathrm{0}}=r_{\mathrm{ac}\mathrm{0},\min }+{\displaystyle \frac{\mathrm{LAI}-{\mathrm{LAI}}_{\min }}{{\mathrm{LAI}}_{\max }-{\mathrm{LAI}}_{\min }}}\left(r_{\mathrm{ac}\mathrm{0},\max }-r_{\mathrm{ac}\mathrm{0},\min }\right).\end{array}$$

Here, the variables LAI_min and LAI_max [m² m⁻²] are the minimum and maximum observed LAI during a specific year, respectively, and r_ac0,min and r_ac0,max are initial resistances (Table S14).

Resistance to the ground (r_g) is expressed as follows:

$$\begin{array}{}\text{(149)}& r_{\mathrm{g}}={\left({\displaystyle \frac{\mathrm{1}-min\left\{\mathrm{1},\mathrm{2}{f}_{\mathrm{snow}}\right\}}{\mathrm{200}}}+{\displaystyle \frac{min\left\{\mathrm{1},\mathrm{2}{f}_{\mathrm{snow}}\right\}}{\mathrm{2000}}}\right)}^{-\mathrm{1}}.\end{array}$$

The variable f_snow is the fraction of the surface covered by snow [unitless]:

$$\begin{array}{}\text{(150)}& f_{\mathrm{snow}}=\min \left\{\mathrm{1},{\displaystyle \frac{\mathrm{SD}}{{\mathrm{SD}}_{\max }}}\right\}.\end{array}$$

3.9 DO₃SE

Deposition of Ozone for Stomatal Exchange (DO₃SE), as described below, is consistent with the parameterization in the European Monitoring and Evaluation Programme (EMEP) model (Simpson et al., 2012). DO₃SE uses two methods to estimate r_st: the multiplicative method based on Jarvis (1976) (“DO₃SE multi”) and the coupled photosynthesis–stomatal conductance method based on Leuning (1995) (“DO₃SE psn”). Unless stated otherwise, the components are the same between DO₃SE multi and DO₃SE psn. Parameters in Table S16 are site-specific.

3.9.1 Surface resistance

Surface resistance (r_c) is calculated as follows:

$$\begin{array}{}\text{(151)}& r_c={\left({\displaystyle \frac{\mathrm{LAI}}{r_{\mathrm{st}}}}+{\displaystyle \frac{\mathrm{StAI}}{r_{\mathrm{cut}}}}+{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{ac}}+r_{\mathrm{g}}}}\right)}^{-\mathrm{1}}.\end{array}$$

The parameter StAI is the stand area index [m² m⁻²].

For forests,

$$\begin{array}{}\text{(152)}& \mathrm{StAI}=\mathrm{LAI}+\mathrm{1}.\end{array}$$

For the other LULC types examined here,

$$\begin{array}{}\text{(153)}& \mathrm{StAI}=\mathrm{LAI}.\end{array}$$

3.9.2 Stomatal resistance

For DO₃SE multi, according to Simpson et al. (2012), stomatal resistance (r_st) is expressed as follows:

$$\begin{array}{}\text{(154)}& r_{\mathrm{st}}={\left(g_{\max }max\left\{f_{\min },f\left(T_{\mathrm{a}}\right)f\left(\mathrm{VPD}\right)f\left(w_{\mathrm{2}}\right)\right\}{a}_{\mathrm{phen}}{a}_{\mathrm{light}}\right)}^{-\mathrm{1}}.\end{array}$$

The parameter g_max is the maximum stomatal conductance [m s⁻¹] (Table S16) and f_min is the minimum factor [unitless] (Table S16). The effects of T_a are expressed as follows:

$$\begin{array}{}\text{(155)}& f\left(T_{\mathrm{a}}\right)=\left\{\begin{array}{l}\frac{T_{\mathrm{a}}-T_{\min }}{T_{\mathrm{opt}}-T_{\min }}{\left(\frac{T_{\max }-T_{\mathrm{a}}}{T_{\max }-T_{\mathrm{opt}}}\right)}^{\frac{T_{\max }-T_{\mathrm{opt}}}{T_{\mathrm{opt}}-T_{\min }}},T_{\min }\le T_{\mathrm{a}}\le T_{\max },\\ \mathrm{0.01},\mathrm{otherwise}.\end{array}\right.\end{array}$$

The parameters T_min, T_max, and T_opt [^∘C] are the minimum, maximum, and optimum temperatures, respectively (Table S16).

The effects of VPD are as follows:

$$\begin{array}{}\text{(156)}& \begin{array}{rl}f\left(\mathrm{VPD}\right)& =\min \mathit{\{}\mathrm{1},\max \mathit{\{}{f}_{\min },f_{\min }\\ & +(\mathrm{1}-f_{\min }){\displaystyle \frac{{\mathrm{VPD}}_{\min }-\mathrm{VPD}}{{\mathrm{VPD}}_{\min }-{\mathrm{VPD}}_{\max }}}\mathit{\}}.\end{array}\end{array}$$

Parameters VPD_min and VPD_max [kPa] are minimum and maximum VPD, respectively (Table S16).

The effects of w₂ are expressed as follows:

$$\begin{array}{}\text{(157)}& \begin{array}{rl}f\left(w_{\mathrm{2}}\right)& =\min \mathit{\{}\mathrm{1},\max \mathit{\{}{f}_{\min },f_{\min }\\ & +(\mathrm{1}-f_{\min }){\displaystyle \frac{w_{\mathrm{wlt}}-w_{\mathrm{2}}}{w_{\max }-\mathrm{0.5}(w_{\mathrm{fc}}-w_{\mathrm{wlt}})}}\mathit{\}}.\end{array}\end{array}$$

The variable a_phen is calculates as follows:

$$\begin{array}{}\text{(158)}& a_{\mathrm{phen}}=\left\{\begin{array}{l}\mathrm{0},d_{\mathrm{y}}\le d_{\mathrm{SGS}}\mathrm{or}{d}_{\mathrm{y}}>d_{\mathrm{EGS}}\\ \begin{array}{c}{\mathrm{\varnothing }}_a+\left(\frac{d_{\mathrm{y}}-d_{\mathrm{SGS}}}{(d_{\mathrm{SGS}}+{\mathrm{\varnothing }}_d)-d_{\mathrm{SGS}}}\right)\left({\mathrm{\varnothing }}_b-{\mathrm{\varnothing }}_a\right),d_{\mathrm{SGS}}\le d_{\mathrm{y}}<d_{\mathrm{SGS}}+{\mathrm{\varnothing }}_d\\ {\mathrm{\varnothing }}_b,d_{\mathrm{SGS}}+{\mathrm{\varnothing }}_d<d_{\mathrm{y}}\le d_{\mathrm{EGS}}-{\mathrm{\varnothing }}_e\end{array}\\ {\mathrm{\varnothing }}_b-\left(\frac{d_{\mathrm{y}}-(d_{EGS-{\mathrm{\varnothing }}_e)}}{d_{\mathrm{EGS}}-{\mathrm{\varnothing }}_e}\right)\left({\mathrm{\varnothing }}_b-{\mathrm{\varnothing }}_c\right),d_{\mathrm{EGS}}-{\mathrm{\varnothing }}_e<d_{\mathrm{y}}\le d_{\mathrm{EGS}}.\end{array}\right.\end{array}$$

Here, d_y is the day of the year, d_SGS is the day of the year that corresponds to the start of the growing season, and d_EGS is the day of the year that corresponds to the end of the growing season. For forests, d_SGS and d_EGS are estimated: d_SGS equals 105 at 50^∘ N and alters by 1.5 d per degree latitude earlier moving south and later moving north, and d_EGS equals 297 at 50^∘ N and alters by 2 d per degree latitude earlier moving north and later moving south. The values of ∅_a, ∅_b, ∅_c, ∅_d, and ∅_e are given in Table S16. For other LULC, we assume a yearlong growing season.

The variable a_light is expressed as follows:

$$\begin{array}{}\text{(159)}& a_{\mathrm{light}}={\displaystyle \frac{{\mathrm{LAI}}_{\mathrm{sun}}}{\mathrm{LAI}}}\left(\mathrm{1}-e^{-\mathit{\alpha }{I}_{\mathrm{PAR}}^{\mathrm{sun}}}\right)+{\displaystyle \frac{{\mathrm{LAI}}_{\mathrm{shd}}}{\mathrm{LAI}}}\left(\mathrm{1}-e^{-\mathit{\alpha }{I}_{\mathrm{PAR}}^{\mathrm{shd}}}\right).\end{array}$$

The parameter α is empirical (Table S16); sunlit and shaded portions of LAI (LAI_sun and LAI_shd, respectively) follow Norman (1979, 1982):

$$\begin{array}{}\text{(160)}& {\displaystyle }{\mathrm{LAI}}_{\mathrm{sun}}=\left(\mathrm{1}-{\mathrm{e}}^{-\mathrm{0.5}\frac{\mathrm{LAI}}{\cos \mathit{\theta }}}\right)\mathrm{2}\cos \mathit{\theta },\text{(161)}& {\displaystyle }{\mathrm{LAI}}_{\mathrm{shd}}=\mathrm{LAI}-{\mathrm{LAI}}_{\mathrm{sun}}.\end{array}$$

The variables $I_{\mathrm{PAR}}^{\mathrm{sun}}$ and $I_{\mathrm{PAR}}^{\mathrm{shade}}$ [W m⁻²] are calculated as follows:

$$\begin{array}{}\text{(162)}& {\displaystyle }{I}_{\mathrm{PAR}}^{\mathrm{shd}}=I_{\mathrm{diff}}{e}^{-\mathrm{0.5}{\mathrm{LAI}}^{\mathrm{0.7}}}+\mathrm{0.07}{I}_{\mathrm{dir}}\left(\mathrm{1.1}-\mathrm{0.1}\mathrm{LAI}\right)e^{-\cos \mathit{\theta }},\text{(163)}& {\displaystyle }{I}_{\mathrm{PAR}}^{\mathrm{sun}}={\displaystyle \frac{I_{\mathrm{dir}}\cos {\mathit{\alpha }}_{\mathrm{1}}}{\cos \mathit{\theta }}}+I_{\mathrm{PAR}}^{\mathrm{shd}}.\end{array}$$

Here, the parameter α₁ is the average inclination of leaves [^∘60], and I_diff and I_dir are the respective diffuse and direct radiation [W m⁻²] estimated as a function of the potential to actual PAR. Potential PAR is estimated using standard solar geometry methods assuming no cloud cover and a sky transmissivity of 0.9.

For DO₃SE psn (Leuning, 1990, 1995), which requires an estimate of net photosynthesis (A_n) [mol CO₂ m⁻² s⁻¹] (Farquhar et al., 1980), stomatal resistance (r_st) is calculated as follows:

$$\begin{array}{}\text{(164)}& \begin{array}{rl}{r}_{\mathrm{st}}& =(g_{\mathrm{0}}+g_{\mathrm{1}}{\displaystyle \frac{A_{\mathrm{n}}}{\left(\right[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{l}}-{\mathrm{\Gamma }}_{*}\left)\right(\mathrm{1}+{\left(\frac{\mathrm{VPD}}{D_{\mathrm{0}}}\right)}^{\mathrm{8}})}}{)}^{-\mathrm{1}}\\ & {\displaystyle \frac{D_{{\mathrm{CO}}_{\mathrm{2}}}}{D_{{\mathrm{O}}_{\mathrm{3}}}}}{\displaystyle \frac{\mathrm{1000.0}\mathit{\rho }}{M_{\mathrm{air}}}}.\end{array}\end{array}$$

Here, the parameter g₀ is minimum conductance [mol air m⁻² s⁻¹] (Leuning, 1990), g₁ is empirical [unitless], D₀ is a parameter related to VPD [kPa] (Leuning et al., 1998) (Table S16), ${\left[{\mathrm{CO}}_{\mathrm{2}}\right]}_{\mathrm{l}}$ is the leaf surface carbon dioxide mixing ratio [mol CO₂ mol air⁻¹], and Γ_* is the carbon dioxide compensation point [mol CO₂ mol air⁻¹]. The ratio of the diffusivities is 0.96. The variable ${\left[{\mathrm{CO}}_{\mathrm{2}}\right]}_{\mathrm{l}}$ is calculated from [CO₂] and leaf boundary layer resistance (r_b,leaf):

$$\begin{array}{}\text{(165)}& r_{\mathrm{b},\mathrm{leaf}}=\mathrm{186}\sqrt{{\displaystyle \frac{u}{l}}}.\end{array}$$

The parameter l is the characteristic dimension of leaves [m].

The variable A_n follows Sharkey et al. (2007):

$$\begin{array}{}\text{(166)}& A_{\mathrm{n}}=min\left\{A_c,A_j,A_{\mathrm{p}}\right\}-R_{\mathrm{d}}.\end{array}$$

The parameter R_d is dark respiration [$\mathrm{0.015}\times {\mathrm{10}}^{-\mathrm{6}}$ mol m⁻² s⁻¹]. The Rubisco-limited rate (A_c) [mol m⁻² s⁻¹] is expressed as follows:

$$\begin{array}{}\text{(167)}& A_c=a_{\mathrm{phen}}f\left(w_{\mathrm{2}}\right)V_{\mathrm{cmax},\mathrm{25}}{\displaystyle \frac{{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}_i-{\mathrm{\Gamma }}_{*}}{{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}_i+K_{\mathrm{c}}\left(\mathrm{1}+\frac{o_{\mathrm{i}}}{K_{\mathrm{o}}}\right)}}.\end{array}$$

Here, the variable ${\left[{\mathrm{CO}}_{\mathrm{2}}\right]}_i$ is the intercellular carbon dioxide partial pressure [Pa]; K_c and K_o are Michaelis–Menten constants for carboxylation and oxygenation [Pa], respectively; o_i is the intercellular oxygen partial pressure [Pa]; Γ_* is the CO₂ compensation point [Pa]; V_cmax,25 is the maximum rate of carboxylation at 25 ^∘C [mol m⁻² s⁻¹] (Table S16); a_phen follows Eq. (158); and f(w₂) follows Eq. (157).

The ribulose-1,5-bisphosphate (RuBP)-limited rate (A_j) [mol m⁻² s⁻¹] is calculated as follows:

$$\begin{array}{}\text{(168)}& A_j=J{\displaystyle \frac{{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}_i-{\mathrm{\Gamma }}_{*}}{a{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}_i+b{\mathrm{\Gamma }}_{*}}}.\end{array}$$

Here, the variable J is the electron transport rate [mol m⁻² s⁻¹], and a and b denote the electron requirements for the formation of nicotinamide adenine dinucleotide phosphate hydrogen (NADPH) and adenosine triphosphate (ATP), respectively. We use a equals 4 and b equals 8 (Sharkey et al., 2007).

The product-limited photosynthetic rate (A_p) [mol m⁻² s⁻¹] is expressed as follows:

$$\begin{array}{}\text{(169)}& A_{\mathrm{p}}=\mathrm{0.5}{V}_{\mathrm{cmax},\mathrm{25}}.\end{array}$$

3.9.3 Cuticular resistance

The resistance to cuticles (r_cut) is prescribed [2500 s m⁻¹].

3.9.4 Resistances to the lower canopy and ground (and associated resistances to transport)

The resistance to in-canopy turbulence (r_ac) follows Erisman et al. (1994):

$$\begin{array}{}\text{(170)}& r_{\mathrm{ac}}=\mathrm{14}{\displaystyle \frac{h\mathrm{StAI}}{u_{*}}}.\end{array}$$

Resistance to the ground (r_g) is calculated as follows:

$$\begin{array}{}\text{(171)}& r_{\mathrm{g}}=\mathrm{200}+\mathrm{1000}{e}^{-T_{\mathrm{a}}-\mathrm{4}}+\mathrm{2000}{\mathit{\delta }}_{\mathrm{snow}}.\end{array}$$

The parameter δ_snow equals 1 when snow is present and 0 when snow is absent.

3.10 MLC-CHEM

The Multi-layer Canopy and Chemistry Exchange Model (MLC-CHEM) has been applied to evaluate the role of in-canopy interactions on atmosphere–biosphere exchanges and atmospheric composition at field sites (e.g., Visser et al., 2021) and the global scale (e.g., Ganzeveld et al., 2010). MLC-CHEM requires a minimum h of 0.5 m, so it has not been configured for all sites. The canopy environment is represented by an understory and crown layer. However, radiation-dependent processes such as biogenic emissions, photolysis, and stomatal conductance are estimated at four canopy layers to consider observed large gradients in in-canopy radiation as a function of the vertical distribution of biomass. For the single-point model, ∼75 % and ∼25 % of the total LAI is present in the crown layer and understory, respectively. These canopy structure settings are used to calculate in-canopy profiles of direct and diffusive radiation as well as the fraction of sunlit leaves from the surface incoming solar radiation (Norman, 1979). Simulated radiation-dependent processes for the four layers are then scaled-up to two layers for in-canopy and canopy-top fluxes and concentrations using the vertical LAI distribution.

MLC-CHEM diagnoses canopy-scale v_d from simulated canopy-top ozone fluxes divided by [O₃], which is the ambient ozone mixing ratio at z_r [ppbv] (Table 1). Turbulent exchanges of ozone between the crown layer (subscript “cl”) and understory (subscript “us”) and between the surface layer (subscript “sl”) and crown layer are calculated from assumed linear [O₃] gradients between heights and from eddy diffusivities. The eddy diffusivity (K_sl→cl) [m² s⁻¹] is expressed as follows (Ganzeveld and Lelieveld, 1995):

$$\begin{array}{}\text{(172)}& K_{\mathrm{sl}\to \mathrm{cl}}=\left(z_{\mathrm{sl}}-z_{\mathrm{cl}}\right)/r_{\mathrm{a}}.\end{array}$$

The eddy diffusivity between the crown layer and understory (K_cl→us) [m² s⁻¹] is calculated as follows:

$$\begin{array}{}\text{(173)}& K_{\mathrm{cl}\to \mathrm{us}}=K_{\mathrm{sl}\to \mathrm{cl}}{u}_{\mathrm{cl}\to \mathrm{us}}/u.\end{array}$$

Here, the variable u_cl→us is wind speed at the crown layer–understory interface [m s⁻¹], calculated as a function of u and canopy structure (Cionco, 1978).

Resistance to leaf-level uptake per layer (r_l,layer) is expressed as follows:

$$\begin{array}{}\text{(174)}& r_{\mathrm{l},\mathrm{layer}}={\displaystyle \frac{r_{\mathrm{b},\mathrm{leaf}}+{\left(\frac{\mathrm{1}}{r_{\mathrm{st}}}+\frac{\mathrm{1}}{r_{\mathrm{cut}}}\right)}^{-\mathrm{1}}}{\max \left\{{\mathrm{LAI}}_{\mathrm{layer}},{\mathrm{10}}^{-\mathrm{5}}\right\}}}.\end{array}$$

Here, the variable r_b,leaf is the resistance to transport through the quasi-laminar boundary layer resistance around leaves (Table S3). Leaf-level stomatal resistance (r_st) is calculated using a photosynthesis–stomatal conductance model (Ronda et al., 2001):

$$\begin{array}{}\text{(175)}& \begin{array}{rl}{r}_{\mathrm{st}}& =f\left(w_{\mathrm{2}}\right)R_{\mathrm{diff},\mathrm{st}}\left({\displaystyle \frac{D_{\mathrm{w}}}{D_{{\mathrm{CO}}_{\mathrm{2}}}}}\right.\\ & {\left.\left(g_{\mathrm{0}}+g_{\mathrm{1}}{\displaystyle \frac{A_n}{\left(\right[{\mathrm{CO}}_{\mathrm{2}}]-{\mathrm{\Gamma }}_{*})(\mathrm{1}+\mathrm{8.09}\frac{\mathrm{VPD}}{D_{\mathrm{0}}})}}{\displaystyle \frac{M_{\mathrm{air}}}{\mathrm{1000}\mathit{\rho }}}\right)\right)}^{-\mathrm{1}}.\end{array}\end{array}$$

The ratio of the diffusivities of water vapor to carbon dioxide is 1.6; g₀ is set to $\mathrm{0.025}\times {\mathrm{10}}^{-\mathrm{3}}$ m s⁻¹ (Leuning, 1990); g₁ is set to 9.09; A_n is net photosynthesis [µmol CO₂ m⁻² s⁻¹], calculated as a function of G, leaf temperature, [CO₂], and soil moisture (Ronda et al., 2001); Γ_* is the CO₂ compensation point [45 ppmv]; and D₀ [kPa] is the VPD at which stomata close (this term is calculated each time step from vegetation-specific constants; Ronda et al., 2001). The soil moisture effect is expressed as follows:

$$\begin{array}{}\text{(176)}& \begin{array}{rl}f\left(w_{\mathrm{2}}\right)& =\mathrm{2}max\mathit{\{}min\mathit{\{}{\mathrm{10}}^{-\mathrm{3}},{\displaystyle \frac{w_s-w_{\mathrm{wlt}}}{\mathrm{0.75}{w}_{\mathrm{fc}}-w_{\mathrm{wlt}}}}\mathit{\}},\mathrm{1}\mathit{\}}\\ & -(max\mathit{\{}min\mathit{\{}{\mathrm{10}}^{-\mathrm{3}},{\displaystyle \frac{w_s-w_{\mathrm{wlt}}}{\mathrm{0.75}{w}_{\mathrm{fc}}-w_{\mathrm{wlt}}}}\mathit{\}},\mathrm{1}\mathit{\}}{)}^{\mathrm{2}}.\end{array}\end{array}$$

Leaf-level cuticular resistance (r_cut) is calculated as follows (Wesely, 1989; Ganzeveld and Lelieveld, 1995; Ganzeveld et al., 1998):

$$\begin{array}{}\text{(177)}& r_{\mathrm{cut}}={\left({\displaystyle \frac{\mathrm{1}-f_{\mathrm{wet}}}{\mathrm{5}\times {\mathrm{10}}^{\mathrm{5}}}}+{\displaystyle \frac{f_{\mathrm{wet}}}{\mathrm{1000}}}\right)}^{-\mathrm{1}}.\end{array}$$

In-canopy aerodynamic resistance (r_ac) considers turbulent transport through the understory to the ground:

$$\begin{array}{}\text{(178)}& r_{\mathrm{ac}}=\mathrm{14}{\displaystyle \frac{\mathrm{0.25}h\mathrm{LAI}}{u^{*}}}.\end{array}$$

To estimate dry deposition to the ground, r_ac is added in series with r_g, which is the resistance to the ground [400 s m⁻¹] (Wesely, 1989; Ganzeveld and Lelieveld, 1995; Ganzeveld et al., 1998). If there is snow, r_g is 2000 s m⁻¹. Resistances are combined with the lowermost understory leaf resistance ($r_{\mathrm{l},\mathrm{layer},\mathrm{1}}$) to create a lowermost understory canopy resistance ($r_{\mathrm{c},\mathrm{layer},\mathrm{1}}$):

$$\begin{array}{}\text{(179)}& r_{\mathrm{c},\mathrm{layer},\mathrm{1}}={\left({\displaystyle \frac{\mathrm{1}}{r_{\mathrm{l},\mathrm{layer},\mathrm{1}}}}+{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{ac}}+r_{\mathrm{g}}}}\right)}^{-\mathrm{1}}.\end{array}$$

In contrast to big-leaf schemes, effective conductances for MLC-CHEM do not add up exactly to v_d because there is an in-canopy [O₃] gradient due to sources and sinks and transport.

4 Measurements for driving and evaluating single-point models

4.1 Turbulent fluxes of ozone

Our best observational constraints on dry deposition are turbulent fluxes, but fluxes integrate the influence of many processes and are not necessarily only reflective of dry deposition. For example, ambient chemical loss of ozone can influence ozone fluxes when the chemistry occurs on the timescale of turbulence. Relevant reactions for ozone fluxes are ozone reacting with highly reactive biogenic volatile organic compounds (BVOCs) or nitrogen oxide (NO). When there are no other sources and sinks aside from dry deposition below the measurement height, dividing the observed turbulent flux by the ambient concentration at the same height can give a measure of the efficiency of dry deposition (“the deposition velocity”). While fluxes provide key constraints on the amount of gas removed by the surface, deposition velocities aid in building the predictive ability of dry deposition given that they indicate how the strength of the removal changes with meteorology and environmental conditions. Turbulent fluxes are mostly measured at individual sites, representing the “ecosystem” scale where the measurement footprint typically extends from the order of 100 m to 1 km. Turbulent fluxes can also be measured from airplanes (e.g., Lenschow et al., 1981; Godowitch, 1990; Mahrt et al., 1995; Wolfe et al., 2015). Turbulent fluxes record changes on hourly or half-hourly timescales, which is important because there is strong sub-daily variability in dry deposition.

Here, we leverage existing long-term and short-term ozone flux datasets over a variety of LULC types to develop a current understanding of model performance and the model spread. Strong observed interannual variability in ozone deposition velocities (Rannik et al., 2012; Clifton et al., 2017; Gerosa et al., 2022), as well as the development of dry deposition schemes based on short-term data (e.g., days to months), motivates our emphasis on multiyear evaluation. Although our evaluation effort would ideally include fluxes of many reactive gases (as well as aerosols), there are not long-term flux measurements of most compounds for which the fluxes primarily represent dry deposition. Such flux observations are oftentimes few and far between and/or challenging to access (Guenther et al., 2011; Fares et al., 2018; Clifton et al., 2020a; Farmer et al., 2021; He et al., 2021). A key reason for this is that obtaining high-frequency concentration measurements of some compounds (e.g., NO₂, SO₂, HNO₃, and H₂O₂) can be challenging due to the detection limits of fast-response sensors, the demands of running research-grade instruments in an eddy covariance configuration (e.g., consumables, dedicated staff, and data storage), and potential flux divergences due to atmospheric chemical consumption or production on the same timescale as deposition processes (Ferrara et al., 2021; Fischer et al., 2021). Nonetheless, recent work further developing or creating new instruments for eddy covariance fluxes of black carbon, ozone, NO₂, ammonia, and a large suite of organic gases (Phillips et al., 2013; Nguyen et al., 2015; Emerson et al., 2018; Fulgham et al., 2019; Novak et al., 2020; Hannun et al., 2020; Ramsay et al., 2018; Schobesberger et al., 2023; Vermeuel et al., 2023) demonstrates the potential for more widespread measurements that would assist in assessing the accuracy of dry deposition schemes more broadly.

Ozone fluxes are the most measured turbulent fluxes of any dry-depositing reactive gas, and they can be measured over seasonal to multiyear timescales. We note that, while the model evaluation component of Activity 2 is only for ozone, the model comparison component of Activity 2 can be performed for other gases.

Ozone turbulent fluxes are measured either via eddy covariance or the gradient method. Eddy covariance is the most fundamental and direct method for measuring turbulent exchange (e.g., Hicks et al., 1989; Dabberdt et al., 1993). Eddy covariance fluxes require concentration analyzers with high measurement frequency to capture the transport of material via turbulent eddies. While fast analyzers are available for ozone, historically they have been resource intensive to operate (note that new techniques like Hannun et al., 2020, are changing this) (Clifton et al., 2020a). However, gradient techniques assume that transport only occurs down the local mean concentration gradient, whereas, in reality, organized turbulent motions can transport material up-gradient (e.g., Raupach, 1979; Gao et al., 1989; Collineau and Brunet, 1993; Thomas and Foken, 2007; Steiner et al., 2011; Patton and Finnigan, 2013). We use some gradient ozone flux datasets, but we caution that they may be particularly uncertain, especially for tall vegetation.

4.2 Site-specific datasets

We simulate ozone deposition velocities by driving single-point models with meteorological and environmental variables measured or inferred from measurements at eight sites. Table 2 summarizes the site locations, LULC types, vegetation composition, and soil types. The set of sites represents a variety of LULC types and climates. The sites include deciduous, evergreen, and mixed forests; shrubs; grasses; and a peat bog. Climate types include Mediterranean, temperate, boreal, and maritime and continental. Dry deposition parameterizations strongly rely on the concept that key processes and parameters are specific to LULC type. While we examine several LULC types here, we emphasize that our measurement test bed is likely insufficient to generalize the results of our study to specific LULC types; thus, we focus our discussion on individual sites. We also cannot discount the fact that differences in ozone flux methods and instrumentation and a lack of coordinated processing protocols across datasets limit meaningful synthesis of our results across sites. Table S17 summarizes details about the ozone flux measurements, the time periods examined, and the post-processing of data. Five of eight sites selected have at least 3 and up to 12 years of ozone flux data (Borden Forest, Easter Bush, Harvard Forest, Hyytiälä, and Ispra). The rest have fewer than 3 years of ozone flux data (Auchencorth Moss, Bugacpuszta, and Ramat Hanadiv) but were included to diversify the climate and LULC types examined.

Table 2 Summary of ozone flux tower sites.

The eddy covariance technique is used for Auchencorth Moss, Bugacpuszta, Harvard Forest, Hyytiälä, Ispra, and Ramat Hanadiv. The gradient technique is used for Borden Forest and Easter Bush. The gradient technique used at Borden Forest is described in Wu et al. (2015, 2016) and was developed for Harvard Forest by comparing gradient and eddy covariance fluxes. Wu et al. (2015) shows that the gradient technique used at Borden Forest strongly overestimates ozone deposition velocities at night and during winter at Harvard Forest, as compared to the ozone deposition velocities calculated from the ozone eddy covariance flux measurements. Wu et al. (2015) also show that parameter choice can strongly influence the deposition velocities inferred from the gradient technique. Thus, seasonal and diel cycle amplitudes as well as the magnitude of observed ozone deposition velocities at Borden Forest are uncertain.

For Activity 2, we selected sites without known influences of highly reactive BVOCs on ozone fluxes. However, there may be unknown influences, especially for coniferous or mixed forests (Kurpius and Goldstein, 2003; Goldstein et al., 2004; Clifton et al., 2019; Vermeuel et al., 2021), and the magnitude of the contribution and how it changes with time are generally uncertain (Wolfe et al., 2011; Vermeuel et al., 2023). Most sites are expected to have very low NO. There may be some influences of NO on ozone fluxes at Ramat Hanadiv (Li et al., 2018) and Ispra, but the magnitude and timing of the contribution are uncertain. Constraining the contributions of highly reactive BVOCs and NO to ozone fluxes is beyond the scope of our work here.

The removal of observed hourly or half-hourly ozone deposition velocity outliers for all sites leverages a univariate adjusted box plot approach following Hubert and Vandervieren (2008), which explicitly accounts for skewness in distributions and identifies the most extreme ozone deposition velocities at each site. Non-Gaussian univariate distributions, or skewness, are present to some degree in each observational dataset used here. This method designates the most extreme 0.7 % of a normal unimodal distribution as outliers, but the exact percentage depends on the degree of skewness. For the datasets used here, which can be highly skewed, we filter 1 %–6 % of ozone deposition velocities across sites. Table S17 describes any other antecedent post-processing of the ozone deposition velocities performed for this effort.

Many dry deposition schemes include adjustments for snow. Table S18 identifies sites with snow depth (SD) measurements. Unless the single-point model directly takes SD input to infer the fractional snow coverage of the surface, we define the presence of snow as a SD greater than 1 cm. Models assume no snow if the SD is less than or equal to 1 cm or is missing.

Canopy wetness is an input to several single-point models. Others do not ingest canopy wetness explicitly as an input variable but rather indicate canopy wetness using a precipitation and/or dew indicator. For the latter type, the fraction of canopy wetness (f_wet) from datasets is not used, and models' indicators are used. Table S18 details canopy wetness measurements at each site. For sites where f_wet data are not available, f_wet values are approximated using an approach used in CMAQ (Table S18).

Soil moisture and soil properties and hydraulic variables are important for stomatal conductance as well as soil deposition processes (Fares et al., 2014; Fumagalli et al., 2016; Stella et al., 2011, 2019). Site-specific details of variables used for near-surface and root-zone volumetric soil water content are described in Table S19. A set of soil hydraulic properties (Table S20) are estimated for each site from soil texture and used across the models employing these parameters. For example, the variable B is an empirical parameter, which is calculated as the slope of the water retention curve in log space (Cosby et al., 1984), that relates volumetric soil water content to soil matric potential and can be referred to as a bulk hydraulic property of the soil (Clapp and Hornberger, 1978; Letts et al., 2000).

Overall, the core description for each site includes the key information needed to drive the single-point models: LULC type, vegetation composition, soil type, and measurement height for ozone fluxes (Table 2 and Table S17). We also describe inputs for snow, canopy wetness, h, and LAI (Table S18). Outside of the core description, other meteorological variables are measured with standard techniques, which are not discussed here. When an input variable is inferred, we detail assumptions involved in the inference because variability in inferred input variables may not be accurately represented and this may need to be accounted for when comparing simulated versus observed ozone deposition velocities.

We note that, in addition to data screening conducted by data providers, driving datasets were visually inspected and clearly erroneous values were set to missing (e.g., in one case T_a was less than −50 ^∘C). Driving datasets are not gap-filled (unless explicitly stated otherwise); therefore, simulated ozone deposition velocities have gaps whenever one or more of a model's input variables is missing. We emphasize that single-point models require different sets of input variables. Thus, output from different models may have different data gaps at a given site. Additionally, because data capture for observed deposition velocities is based on the availability of ozone flux measurements and data gaps in input variables may be different from data gaps in the ozone flux measurements, simulated deposition velocities can have different data gaps from observed deposition velocities. We address data coverage discrepancies across models and observed deposition velocities in two ways: first, we identify the time-averaged observed and simulated deposition velocities with suboptimal coverage in our results (e.g., see Fig. 1); second, we account for diel imbalances in our analysis. Both approaches are described more fully in Sect. 4.3.

Figure 1 Monthly mean ozone deposition velocities (v_d) from the ozone flux observations. The multiyear average is in black, different years are shown using colors, and open symbols indicate months for a given year with low data capture.

4.3 Creation of monthly and seasonal average observed and simulated quantities

We examine averages across 24 h, except for Ramat Hanadiv. For Ramat Hanadiv, many months have missing values during the night and morning; thus, we limit our analysis to 11:00–17:00 (all times are given in local time throughout). Across sites and analyses, we use a weighted averaging approach for daily averages that considers the number of observations for a given hour to avoid the over-representation of any given hour due to sampling imbalances across the diel cycle (e.g., more valid observations during daylight hours).

There are sometimes periods of missing ozone fluxes in the datasets. We indicate year-specific monthly averages with low data capture for observed v_d in Fig. 1. Low data capture is defined as less than or equal to 25 % data capture averaged across 24 h (or 11:00–17:00 for Ramat Hanadiv). In other words, we first compute data capture for each hour of a given month (or season) and then average across hour-specific data capture rates to compare against the 25 % threshold. We indicate multiyear monthly averages with low data capture for observations and models in Figs. 2 and 3. Note that the number of data points used in constructing monthly averages differs between models and observations, and across models. Data capture for each model depends on the availability of the specific measured input data required for driving that model. Data capture for observed v_d is based on the availability of ozone flux measurements.

Figure 2 Multiyear monthly mean ozone deposition velocities (v_d) from ozone flux observations and single-point models. Pink shading denotes the interquartile range across models, red lines denote the minimum and maximum across monthly simulated values, and open symbols on observations indicate months with low data capture.

When we examine multiyear averages, we do not consider sampling biases across years (e.g., more valid observations in 1 year over the other). Thus, more data for 1 year may skew multiyear averages towards values for that year (Fig. 1). However, results are generally similar if we include weighting by years, except when there are only a few years contributing to multiyear averages, and 1 or some of those years have low data coverage. For seasonal averages, months are not given equal weight unless stated otherwise. For example, all non-missing data for a given hour across months of the season are considered equally (e.g., the fact that there may be more data at noon in July than in August is not considered in a summertime average).

5 Results

Figure 1 shows monthly mean observed ozone deposition velocities (v_d) across years, as well as multiyear averages, at all sites. There are a variety of seasonal patterns and magnitudes of observed v_d across sites. Interannual variability is strong in terms of the standard deviation across yearly annual averages normalized by the multiyear average (range of 10 % to 60 % across sites). In some cases, periods with low data coverage contribute to apparent interannual variability and/or seasonality; thus, in these cases, the degree of interannual variability is uncertain. However, more complete ozone flux records also show strong variability from year to year and month to month, suggesting that we can expect strong interannual variability on a monthly basis to be a generally robust feature of the observations. The following discussion focuses on multiyear averages, but we briefly examine summertime (June–August) interannual variability at sites with 3 or more years of data in the individual site subsections below to establish whether models capture the range of interannual variability and/or ranking among different summers.

Figure 2 shows multiyear monthly mean v_d from observations and the spread in multiyear monthly mean v_d across models, whereas Fig. 3 shows multiyear monthly mean values from each individual model and the observations. The minimum and maximum of the monthly averages across the models bracket the observations across most sites and seasons (Fig. 2). The exceptions are Auchencorth Moss (all months except July), Borden Forest (October–November only), and Ispra (October–February only). In some cases, model outliers allow the full set of models to bracket observations (Fig. 3), which suggests limited skill of the model ensemble. If we instead consider the interquartile range across models (hereinafter, “the central models”), there are at least a few months at every site when observations fall out of range. At the same time, at every site except Auchencorth Moss, there are also at least a few months when the observations are within the range, indicating that failure of the central models to capture observations consistently across the seasonal cycle does not suggest a complete lack of skill from the model ensemble that de-emphasizes outliers. Further, the central models are very close to bracketing observations across months at Easter Bush, Hyytiälä, and Harvard Forest.

Figure 3 Multiyear monthly mean ozone deposition velocities (v_d) from ozone flux observations and single-point models. Open symbols indicate months with low data capture.

The model spread in multiyear mean v_d across months and sites is large (Fig. 2). The spread in terms of the model with the highest annual average divided by the model with the lowest ranges from a factor of 1.8 to 2.3 except for Hyytiälä (2.7) and Auchencorth Moss (5). The spread in wintertime (December–February) averages is very high at some sites: Borden (10), Hyytiälä (21), Auchencorth Moss (9.1), and Harvard Forest (6.3). The spread in wintertime averages is a factor of 2 to 3.3 at other sites. The spread is typically lower during summer (June–August) than during winter, on par with annual values. We also use the 75th percentile divided by the 25th percentile as a metric of the spread. This metric for the annual average is a factor of 1.2–1.8. For winter, the metric is also lower for sites with high spreads based on all models (a factor of 3 for Borden Forest, 2.4 for Hyytiälä, 3 for Auchencorth Moss, and 2.7 for Harvard Forest), but it is still higher than the summer and annual spreads (except for Ispra).

Figure 4 Seasonal mean relative biases (simulated minus observed, divided by observed) across models and sites for ozone deposition velocities (v_d), expressed as fractions. Numbers next to model names in the subpanel titles are seasonal mean absolute biases (in cm s⁻¹). DJF is December, January, and February; MAM is March, April, and May; JJA is June, July, and August; and SON is September, October, and November.

Figure 4 shows the relative biases (simulated minus observed, divided by observed) across months, sites, and seasons. When we consider individual model performance, we find that no model is always within 50 % of the observed multiyear averages across sites and seasons (Fig. 4). Models are very low against observations at Auchencorth Moss, but the previous statement holds even excluding this site. In general, a key finding here is that model performance varies strongly by model, season, and site. Below, we first discuss the mean absolute biases across sites and then the drivers of seasonality across models and sites. Following this, in the subsections, we discuss each site, starting with short vegetation and then moving on to forests.

The absolute bias (simulated minus observed) averaged across multiyear seasonal averages and sites is highest for GEM-MACH Wesely (0.22 cm s⁻¹) and lowest for CMAQ M3Dry-psn (0.12 cm s⁻¹) (Fig. 4). GEM-MACH Zhang, WRF-Chem Wesely, GEOS-Chem Wesely, TEMIR Wesely, TEMIR Wesely BB, and TEMIR Wesely Medlyn are on the higher end of the spread with respect to the mean absolute bias across seasons and sites (0.17–0.18 cm s⁻¹), whereas DO₃SE multi, DO₃SE psn, IFS SUMO Wesely (0.13 cm s⁻¹), and CMAQ M3Dry (0.14 cm s⁻¹) are on the lower end, with the rest in between (0.15–0.16 cm s⁻¹). (MLC-CHEM does not simulate three sites, so we exclude it here.)

The absolute biases averaged across seasons may overemphasize model performance when v_d values are high. Given that wintertime v_d tends to be lower in magnitude than during other seasons, we also examine the wintertime mean absolute biases across sites (Fig. 4). Values are highest for GEM-MACH Zhang (0.22 cm s⁻¹), GEM-MACH Wesely (0.20 cm s⁻¹), TEMIR Wesely (0.20 cm s⁻¹), and TEMIR Wesely Medlyn (0.19 cm s⁻¹). Otherwise, model biases are below 0.16 cm s⁻¹.

Figure 5 shows the simulated multiyear wintertime and summertime mean effective conductances as well as the observed multiyear seasonal average v_d (recall that simulated effective conductances sum to simulated v_d). The three main pathways are stomata, cuticles, and soil; even when models simulate lower-canopy uptake, uptake via this pathway tends to be low. Thus, we focus on stomatal, cuticular, and soil pathways. There are three important takeaways from Fig. 5. First, models can disagree in terms of relative contributions from pathways, even when they predict similar v_d; conversely, models can agree in terms of relative contributions of pathways but predict different v_d. Second, stomatal and nonstomatal pathways both have important contributions to v_d across models and are both key drivers of variability across models. Third, models tend to disagree on cuticular versus soil contributions to nonstomatal uptake at some sites while agreeing at others.

Figure 5 Multiyear seasonal mean simulated effective conductances and observed ozone deposition velocities (v_d). Black dots are simulated v_d (black dots should equal the top of the bars). DJF is December, January, and February, and JJA is June, July, and August.

Figure 6 Pathways contributing to variability across simulated multiyear monthly mean ozone deposition velocities. The variance for each effective conductance is a solid color. Twice the covariance between effective conductances is a hatched pattern (the colors of the hatching correspond to the pathways examined). Each value is normalized by the absolute value of the sum of the variances and twice the covariances so that we are comparing the pathways that drive seasonality across models in a relative sense (rather than the seasonal amplitude as well).

Figure 6 shows how the multiyear mean seasonality of effective conductances contributes to the multiyear mean seasonality of simulated v_d across models. Specifically, the variance in each pathway across months is shown, as well as twice the covariance between individual pathways. Negative covariances imply offsetting seasonality between the two pathways (i.e., an anticorrelation in the seasonal cycles of two pathways that acts to dampen the total seasonality). Positive covariances mean that a positive correlation in seasonal cycles of the two pathways acts to amplify total seasonality. Values are normalized by the absolute sum of the variance and twice the covariances so that Fig. 6 does not emphasize differences in the seasonal amplitude, rather what pathways control the seasonality.

The key finding from Fig. 6 is that stomatal uptake is the most important driver of the multiyear mean v_d seasonality for most models and sites. For some models and sites, cuticular uptake also plays a role, albeit mostly just via correlations with stomatal uptake. Correlations between stomatal and cuticular pathways are mostly positive, and thus tend to amplify v_d seasonality. Exceptions are Hyytiälä and Easter Bush, where some models show anticorrelations between stomatal and cuticular uptake seasonal cycles. With a few exceptions (e.g., at Easter Bush and for the GEM-MACH Wesely and DO₃SE models), soil uptake tends to play a more minor role.

In general, the parameters and dependencies driving simulated v_d seasonality are model-dependent. Expected dominant influences include changes in initial resistances with season, cuticular and stomatal dependencies on LAI, stomatal dependencies on soil moisture, temperature response functions (used in Wesely, 1989, to decrease nonstomatal deposition pathways at cold temperatures), and changes with snow.

Figure 7 shows how multiyear monthly mean v_d changes with LAI, for both the models and the observations. Multiyear monthly mean observed and simulated v_d generally increases with LAI across sites during at least some time periods of plant growth (Fig. 7). In general, however, the relationship between v_d and LAI on monthly timescales is nonlinear for both observations and models, distinct between observations versus models, and distinct across models. Many models show a strong sensitivity to the LAI, which has been pointed out in previous work (Cooter and Schwede, 2000; Charusombat et al., 2010; Schwede et al., 2011; Silva and Heald, 2018). Our analysis here, combined with past work, suggests that advancing predictive ability requires a better understanding of observed v_d–LAI relationships in terms of seasonality and site-to-site differences.

Figure 7 Multiyear monthly mean ozone deposition velocities (v_d) versus the leaf area index (LAI).

Figure 8 Multiyear mean ozone deposition velocity (v_d) during all conditions versus when the snow depth is greater than or equal to 1 cm for sites with snow depth records and sufficient time with snow (25 %, averaged across hours per month). Months considered are December–February for Bugacpuszta, December–February for Borden Forest, and November–March for Hyytiälä. Months are given equal weight in the averages.

Figure 8 shows snow's impact on the multiyear mean v_d at sites with snow depth records and sufficient snowy periods. Observations suggest modest reductions with snow at Bugacpuszta and Hyytiälä but not much change at Borden Forest. At Borden Forest, some models show decreases, whereas others show little change. At Hyytiälä and Bugacpuszta, some models capture decreases with snow despite biases, whereas other models understate or exaggerate decreases. Observed reductions with snow are larger at Bugacpuszta than Hyytiälä, and many models capture this. Findings with respect to Borden Forest may reflect that snow is not measured there, rather 15 km away, and thus these measurements do not reflect the exact local conditions. Even though some models do not capture the magnitude of the observed v_d decreases with snow, Fig. 8 shows that models' inability to capture the magnitude of wintertime values (snow or snow-free) at a given site is a much larger problem than models' inability to capture responses to snow, at least at these three sites. The relative model spread (based on the standard deviation across models divided by the average) does not change substantially under snowy versus all conditions, except at Bugacpuszta (27 % versus 70 %), further underscoring the need to better understand wintertime v_d in a more general sense.

The relatively low magnitude of snow-induced observed v_d changes indicates that snow-induced changes are not the main driver of observed v_d seasonality (Fig. 8). For example, observed changes with snow are a small fraction of the observed absolute seasonal amplitude of multiyear monthly averages at these sites, at least for Hyytiälä and Borden Forest. We also note that models simulate v_d reductions with snow at Hyytiälä and Bugacpuszta even when snow is not model input, suggesting that other model dependencies (e.g., temperature response functions) may lead to changes coincident with snow. Recent papers have suggested that better snow cover representation may be key with respect to capturing v_d spatial variability at regional scales and regional average seasonal cycles as well as changes with climate change (Helmig et al., 2007; Andersson and Engardt, 2010; Matichuk et al., 2017; Clifton et al., 2020b). Despite insufficient data to examine spatial variability or responses to climate change, our analysis suggests that it is important to understand drivers of wintertime v_d other than snow.

5.1 Bugacpuszta

Bugacpuszta is a semiarid and seminatural grassland in Hungary that experiences grazing during most of the year. In terms of variability across models, the model spread based on the model with the highest annual average v_d divided by the model with the lowest is a factor of 2.1 (2.8 during summer and 2.2 during winter); however, based on the interquartile range, this value is a factor of 1.3 (1.2 during summer and 1.3 during winter). This model spread at Bugacpuszta is on the lower end of the estimates across the sites examined.

A longer ozone flux data record is needed to assess the interannual variability at Bugacpuszta. Bugacpuszta has only 1 year of data during February–May (from 2013), 2 years of data during August–December (from 2012 and 2013), and 2 years of data during January (from 2013 and 2014) (Fig. 1). Data are always missing during June and July. For time periods with 2 years of data, observed monthly mean v_d values are very close in magnitude between years. The exception is October, for which 2013 values are half of the 2012 values. However, October 2013 has very low data coverage (only ∼2–3 d of coverage), and hourly values exhibit high uncertainty compared with other months (not shown). Agreement between both years for months with sufficient data coverage may suggest that there is low interannual variability at this site; however, more data are really needed to make this assessment, as well as in-depth analyses of shorter-term variability (e.g., diel cycles). In the following, we focus on the “multiyear averages” at this site, acknowledging that there are only 2 years of data during 6 months of the year (and 10 months total with data).

Without June and July observations, we cannot fully assess seasonality at Bugacpuszta. Therefore, we evaluate seasonality across other months. The observed seasonal cycle for the months with data is as follows: v_d maximizes during May, following an increase from March, and minimizes during August, after which v_d increases to November and levels off from December to February (Fig. 1). Seasonal patterns are similar across many models, with midsummer peaks after slow increases from winter and similar values from August to November (Fig. 3). Despite similar seasonal patterns across the models as well as fair agreement in the relative seasonal amplitude across the models (Fig. 9), the models disagree with respect to the pathways dominating the seasonal cycle (Fig. 6). Notably, models disagree the most in terms of the pathway(s) driving seasonality at Bugacpuszta relative to other sites, suggesting that changes in individual pathways on seasonal timescales at this location may be a key uncertainty.

Figure 9 Relative seasonal amplitudes of multiyear monthly mean stomatal uptake (sideways triangles) and ozone deposition velocities (upwards triangles) across models, defined as the maximum across months of multiyear monthly averages minus the minimum, divided by the average. Black triangles denote the relative seasonal amplitude of observations for sites with wintertime minima and summertime maxima. Gray shading denotes the interquartile range across models.

The central models bracket observed v_d at Bugacpuszta during December–May but are too high against the observations during August and September (and only slightly too high during October and November) (Fig. 2). Two clear model outliers during the warm months are the TEMIR Zhang models (Fig. 3), which show relatively low soil and cuticular uptake (Fig. 5). TEMIR psn also shows no stomatal uptake, following very low input root-zone soil moisture (below prescribed wilting point). While the TEMIR Zhang models are clear model outliers during the warm months, they allow the complete set of models to bracket observations during August–November, because the other models are mostly too high (or, in a few cases, just right). Without June and July ozone fluxes, however, it is unclear how the TEMIR Zhang models alter the summertime performance of the model spread.

Only eight models show substantial summertime stomatal uptake at Bugacpuszta (Fig. 5). There is no summertime stomatal uptake simulated by the TEMIR psn, IFS SUMO Wesely, or DO₃SE models, and very little summertime stomatal uptake is simulated by CMAQ STAGE, CMAQ M3Dry, and CMAQ M3Dry-psn. Only these models employ soil moisture dependencies on stomatal conductance (MLC-CHEM does as well but does not simulate values at Bugacpuszta); these models simulate little to no stomatal uptake at Bugacpuszta because input soil moisture is below the prescribed wilting point. We emphasize that the wilting point, which is not a directly measurable quantity, is uncertain across sites. Nonetheless, the magnitude of stomatal uptake at this site is a clear uncertainty. If we instead focus on the models with substantial summertime stomatal uptake, we can see that they show a large spread in the stomatal fraction of v_d – from 12.5 % to 40 % with one model simulating 60 % (Fig. 12) – and produce distinct stomatal uptake seasonal cycles (Fig. 10). On the other hand, many models show similar v_d seasonal cycle shapes (Fig. 3) but dissimilar stomatal uptake seasonal cycle shapes. These results suggest that nonstomatal uptake seasonality plays a role in normalizing differences in v_d seasonal cycles across models, and the models are more distinct than implied by v_d alone.

Figure 10 Multiyear monthly mean effective stomatal conductance (eg_s) from single-point models. Gray shading denotes the multiyear monthly mean LAI (used to emphasize seasonality in this variable; y ranges are not given).

Figure 11 Model spread (standard deviation) across multiyear seasonal mean ozone deposition velocities (v_d) and effective conductances for DJF (stars) and JJA (circles). DJF is December, January, and February, and JJA is June, July, and August.

Bugacpuszta has the most similar summertime model spreads across the top three deposition pathways relative to other sites (except Hyytiälä) (Fig. 11), suggesting a high degree of uncertainty in the magnitude of all pathways during the warm months. Most models show substantial summertime contributions from soil uptake, but the magnitude of soil uptake varies across models (Fig. 5). In contrast, for the summertime cuticular and stomatal pathways, models disagree as to whether contributions are substantial as well as disagreeing on the magnitude of uptake. For example, like how some models show very low stomatal uptake (as discussed above), some models show negligible cuticular uptake. Establishing whether there should be summertime stomatal and/or cuticular uptake at Bugacpuszta would be a first step towards further constraining models.

The multiyear monthly mean LAI at Bugacpuszta shows a sharp summer peak, maximizing during June (∼3.6 m² m⁻²) (Fig. 10). Values are similar during August to November and then decrease from November to March, with a minimum during March. Observed v_d is missing for LAI values greater than 2 m² m⁻² (corresponding to June and July). There is no discernable observed v_d–LAI relationship for LAI values below 1 m² m⁻², and models capture this (Fig. 7). Observations show a strong v_d increase from 1 to 2 m² m⁻². Models show an increase, but most do not capture the large observed slope. This is especially true for models with soil moisture dependencies on stomatal conductance, implying that, during at least some periods of high vegetation density, there should not be soil moisture stress or not as strong soil moisture stress as simulated by some models.

Models simulate that soil uptake dominates wintertime v_d at Bugacpuszta (Fig. 5). The exception is GEM-MACH Wesely, which underestimates wintertime v_d. Wintertime stomatal fractions of v_d can be up to 10 % (due to low v_d overall) but are mostly within 0 %–5 %. Because the central models capture wintertime v_d (Fig. 2) and models also agree that soil uptake dominates, some models may have some skill during cooler months. However, there is variability in soil uptake across models (Fig. 11). Models largely capture observed wintertime v_d decreases with snow, with most slightly overestimating the change but a few (the DO₃SE models, WRF-Chem Wesely, TEMIR Zhang, and GEM-MACH Wesely) underestimating it (Fig. 8). Future attention to the noncentral models should focus on better capturing wintertime nonstomatal uptake generally at this site, rather than changes with snow.

A key outstanding question at Bugacpuszta is as follows: should models simulate low stomatal uptake throughout summer or only during late summer? Most model values are too high compared with observations during August and September. This includes models employing soil moisture dependencies on stomatal conductance (and, thus, simulating very low to no stomatal uptake), implying an overly high simulated nonstomatal uptake. Continuous year-round ozone flux observations, especially during periods of the growing season with and without moisture stress, are needed to better assess model performance at Bugacpuszta. Independent measures of stomatal conductance during periods of missing ozone fluxes would be useful in constraining the absolute stomatal portion of dry deposition, but further constraining nonstomatal uptake, which models indicate is an important fraction of summertime v_d (despite disagreeing on the exact pathway), requires additional ozone flux measurements.

5.2 Auchencorth Moss

Auchencorth Moss is a peat bog in Scotland that is covered with heather, moss, and grass. The model spread in terms of the model with the highest annual average v_d divided by the model with the lowest is a factor of 5 (4.3 during summer and 9.1 during winter); however, based on the interquartile range, this value is a factor of 1.6 (1.5 during summer and 3 during winter). Across sites, for the annual metrics, Auchencorth Moss has the largest spread for the maximum/minimum metric and the second largest for the interquartile range.

There is no clear shape of the observed v_d seasonal cycle at Auchencorth Moss (Fig. 1). Whether this is true on a climatological basis is unclear due to (1) data incompleteness during the 2-year period – observed values during February–May have low data capture mostly because data are missing during 2016 – and (2) strong interannual variability when there are data and (3) the fact that there are only 2 years of data. A longer and more complete ozone flux record is needed to fully assess interannual variability as well as seasonality at Auchencorth Moss. Below, we focus on multiyear averages, acknowledging that only half the months of the year have 2 years of data.

A key finding is that models do not capture the high values of v_d that are observed year-round at Auchencorth Moss (Fig. 2). The exception is TEMIR Zhang Medlyn during July. Auchencorth Moss is the only site examined with negative biases (> 30 % of observed multiyear seasonal averages) across seasons and models (except for TEMIR Zhang Medlyn during July) (Fig. 4). Biases tend to be smallest during summer and largest during winter because many models simulate peak v_d during warm months (Fig. 3). Notably, models differ substantially with respect to their relative seasonal amplitudes, with a very even and wide distribution in relative seasonal amplitude across models (Fig. 9), especially relative to other short-vegetation sites.

Simulated v_d seasonality is mostly due to stomatal uptake (Fig. 6). Some models show that soil uptake plays a role, and all but two models show moderate contributions from correlations between pathways. The seasonality shape of stomatal uptake is very similar across most models, as is the magnitude of stomatal uptake throughout the year (Fig. 10). Major exceptions are the TEMIR Medlyn models, which show peak values of around 0.4 cm s⁻¹ in contrast to the rest that average just under 0.1 cm s⁻¹. For the relative seasonal amplitudes in stomatal uptake, the spread across the central models is low (Fig. 9). The value for GEM-MACH Wesely is very high (> 5), with other models' values spanning a factor of 1.75 to 3. Models deviating from the rest with respect to stomatal uptake's seasonality shape are GEM-MACH Zhang (a strong peak during July and near-zero during August and after) and DO₃SE (low during summer) as well as WRF-Chem Wesely and IFS SUMO Wesely (the latter two are similar and higher than others, especially during spring).

While high summertime stomatal uptake combined with moderately high year-round nonstomatal uptake distinguishes TEMIR Zhang Medlyn from the other models (Fig. 5), we see the best agreement between this model and observations during the warm months. However, TEMIR Zhang Medlyn does not capture observed seasonality (or lack thereof). Thus, TEMIR Zhang Medlyn may have more skill during summer than other models; however, like other models, TEMIR Zhang Medlyn struggles with seasonality. Future work should establish whether there is strong seasonality in stomatal uptake coupled with offsetting seasonality in nonstomatal uptake at Auchencorth Moss or whether stomatal uptake should be higher year-round.

Figure 12 Multiyear seasonal mean stomatal fraction of ozone deposition velocities (v_d) across models during DJF (stars) and JJA (circles). Gray shading denotes the interquartile range across models. DJF is December, January, and February, and JJA is June, July, and August.

For soil uptake, the model spread is large and similar between summer and winter (Fig. 11). During summer, the spread in stomatal uptake is on par with soil uptake; spreads for stomatal and soil uptake are the highest across pathways. During winter, the spread in stomatal uptake is very low, and the spread in soil uptake is the highest. Wintertime stomatal fractions vary from 0 % to 20 % across models (Fig. 12). Models except CMAQ STAGE simulate non-negligible soil uptake (Fig. 5). However, during summer, models disagree on the soil contribution to v_d (0 %–80 %) as well as the magnitude of soil uptake. In contrast, during winter, models agree that soil uptake contributes substantially to v_d (> 60 %) (apart from CMAQ STAGE and GEM-MACH Wesely) but disagree on the magnitude of soil uptake. Snow depth is measured at Auchencorth Moss, but data are missing during half of the ozone flux period, and there is not a substantial amount of time with snow when there are measurements.

Models estimate very low to moderate cuticular uptake at Auchencorth Moss (Fig. 5), which is consistent across low-vegetation sites. Moderate values of cuticular uptake are simulated by the GEM-MACH Zhang and TEMIR Zhang models, and values are similar between summer and winter. Otherwise, models simulate very little cuticular uptake during winter and low cuticular uptake during summer. Nonetheless, the model spread in cuticular uptake is similar between seasons. Summertime stomatal fractions vary across the central models from 25 % to 55 % (Fig. 12). Aside from one model simulating 80 % and two models around 10 %, half are around 20 %–30 % and the other half are around 45 %–60 %. There is a clear division across models in that no model simulates stomatal fractions between 32.5 % and 45 %. The dichotomy seems to be due to variability in both stomatal and soil uptake across models, consistent with high summertime model spreads for these pathways (Fig. 11).

Despite an unclear observed v_d seasonal pattern at Auchencorth Moss, the relationship between the monthly mean LAI and v_d may provide insights into model performance. With strong observed v_d variations at low LAI values (less than 0.6 m² m⁻²), there is no relationship, but there is a positive relationship at moderate LAI values (in the range of 0.6 to 0.9 m² m⁻²) (Fig. 7). Observations then show that v_d decreases with LAI increases above 0.8 m² m⁻², but there is only one data point here. Most models seem to capture the observed relationship at moderate LAI values as well as that there should not be a relationship at low LAI values. However, some models (e.g., the TEMIR models) overestimate the increase's slope at moderate LAI values. Thus, some models may have some skill with respect to simulating seasonality in cuticular and/or stomatal uptake. Nonetheless, strong observed v_d variability at low LAI values and changes with LAI during peak vegetation density require a better understanding. With observational constraints on stomatal uptake, we will be able to understand whether nonstomatal uptake should be higher year-round and/or seasonality in nonstomatal uptake should act to offset seasonality in stomatal uptake.

We close by emphasizing that very high observed v_d values at Auchencorth Moss are uncertain – there is strong interannual and day-to-day variability but a lot of missing data. The peat/bog LULC type does not have many ozone flux measurements at other sites that could be used to provide additional context to Auchencorth Moss measurements. Schaller et al. (2022) show that v_d ranges from 0.05 cm s⁻¹ at night to 0.45 cm s⁻¹ during the day in July 2017 at a peatland in northwestern Germany. El Madany et al. (2017) look at ozone fluxes at the same site during 2014 but do not present v_d values. Fowler et al. (2001) present older measurements at Auchencorth Moss, estimated with the gradient technique (eddy covariance is used for the data examined here), showing much lower observed v_d than examined here (e.g., winter and fall values here are twice what they are during 1995–1998, summer values are almost twice as high, and spring values are higher but not twice as high). It is not clear what drives the higher, more recent v_d measurements at Auchencorth Moss analyzed in this study, and more detailed analysis is needed to figure it out. In general, building an understanding of ozone dry deposition for this LULC type provides a key test of the understanding of soil uptake and of its dependence on its expected drivers (soil organic carbon and water content), given that peat/bog soils are organic-rich and wet.

5.3 Easter Bush

Easter Bush is a managed grassland in Scotland that is used for silage harvest and intensive grazing. In terms of variability across models, the spread based on the model with the highest annual average v_d divided by the model with the lowest is a factor of 1.8 (1.8 during summer and 3.0 during winter); however, based on the interquartile range, this value is a factor of 1.3 (1.3 during summer and 1.4 during winter). Model spreads at Easter Bush are some of the lowest compared with other sites.

Easter Bush has one of the longest ozone flux records (Clifton et al., 2020a), and it has the longest record examined here as well as the strongest interannual variability. For example, the coefficient of variation across years is on average 60 % across months. In contrast, other sites show coefficients of variations across years from 10 % to 30 %. There is also strong interannual variability in the observed seasonal cycle's shape at Easter Bush (Fig. 1). As for other sites with long-term records, we focus on multiyear averages but also touch on summertime interannual variability. Some models capture some low summers, but models do not capture high summers (except GEOS-Chem Wesely, IFS GEOS-Chem Wesely, and TEMIR Wesely, which capture 1 high year) and underestimate interannual spread (Fig. 13). Future work should focus on understanding observed interannual variability and should consider that interannual variability changes strongly by month, both in terms of the spread across years and the ranking of years.

Figure 13 Simulated and observed yearly summertime mean ozone deposition velocities (v_d) for sites with records of at least three summers. Values are normalized by the multiyear average of the respective model or observations to emphasize ranking and spread across years. Colors rank yearly values from low (blue) to high (gold) for the observations. The model year is not shown when the observed year is missing. The highest year for Easter Bush is not shown because it is very high (2 times the multiyear mean observed value).

The central models' spread largely brackets the observed multiyear monthly values across months. Specifically, observed values sit mostly on the lower end of or just below the central models' spread, except during May, November, and December when observed values are on the higher end (Fig. 2). Only CMAQ STAGE consistently shows lower v_d than observed, but the relative bias is low (−18 % to −30 %) (Fig. 4). During winter, GEM-MACH Wesely and TEMIR Wesely psn are too low, and the relative biases are substantial (−51 % to −70 %). With a few exceptions (i.e., winter for GEM-MACH Wesely and TEMIR Wesely psn and summer for WRF-Chem Wesely and TEMIR Wesely Medlyn), models are within ±50 % of observed seasonal averages.

Overall, the following suggests that models may have skill with respect to simulating climatological v_d seasonality at Easter Bush, aside from a clear set of outliers. There is a weak warm-season peak in observed v_d (Fig. 1). Models show weak warm-season maxima (Fig. 3) and relatively similar relative seasonal amplitudes (Fig. 9). However, some models are clear outliers. For example, GEM-MACH Wesely and TEMIR Wesely psn show particularly strong relative seasonal amplitudes (Fig. 9), in part due to low wintertime v_d. The absolute standard deviation across models for v_d is higher during winter than during summer (Fig. 11). This only happens at Easter Bush and Hyytiälä; however, as noted above, the wintertime model spread reduces when considering the full versus interquartile range, suggesting that low outliers may drive the large standard deviation across models.

For most models, the primary driver of v_d seasonality is stomatal uptake (Fig. 6). However, individual contributions from stomatal uptake barely contribute for GEM-MACH Wesely and TEMIR Wesely BB. Several models, including GEM-MACH Wesely, GEM-MACH Zhang, the TEMIR Wesely models, and TEMIR Zhang models, simulate large contributions from soil uptake individually and/or via correlations with other pathways. Only two models, in contrast to seven at the other grassland examined (Bugacpuszta), suggest that individual contributions from cuticular uptake matter for seasonality.

Most models are similar in terms of magnitude and seasonality shape of stomatal uptake (Fig. 10), as well as relative seasonal amplitudes (Fig. 9). Exceptions are GEM-MACH Wesely (a very strong peak during July and almost zero after July, thereby showing an anomalous seasonal amplitude), TEMIR Medlyn (much higher than other models during the warm months), and IFS SUMO Wesely and WRF-Chem Wesely (slightly higher than other models, especially during spring). The DO₃SE models are also an exception: they show very different seasonal cycles from each other, despite both being high and seasonally distinctive relative to other models. DO₃SE psn also shows an anomalous seasonal amplitude.

At Easter Bush, LAI peaks during July, with a broad maximum from May to November and low values during February and March (Fig. 10). With some exceptions, models bound the observed relationship between v_d and LAI, agreeing on a fairly weak but positive dependence (Fig. 7). Outliers with respect to the v_d–LAI relationship (GEM-MACH Wesely and TEMIR Wesely psn) also indicate that stomatal uptake does not strongly influence v_d seasonality, suggesting the latter is incorrect.

During summer, model spreads for v_d and deposition pathways at Easter Bush are highest for soil uptake, then stomatal uptake, and then cuticular uptake (Fig. 11). Most models simulate moderate or substantial stomatal uptake, but there is a division as to whether models simulate very low, low, or moderate cuticular uptake (Fig. 5). Models simulate substantial soil uptake, both in terms of absolute magnitudes and the relative contribution to v_d. Exceptions are the DO₃SE models, which have very low soil uptake. Stomatal fractions range from 10 % to 70 %, with most models around 30 % and only four models above 40 % (Fig. 12). The range across models for stomatal fractions is one of the largest across sites, but the interquartile range is one of the smallest. High agreement regarding the stomatal uptake magnitude, seasonality shape, and relative amplitude as well as the stomatal fractions across most models suggests that an appropriate next step would be to use observation-based estimates of stomatal uptake (e.g., from water vapor fluxes) to evaluate whether models are accurate with respect to this pathway.

During winter, models simulate that v_d is dominated by soil uptake, with some models simulating low to moderate contributions from cuticular uptake (Fig. 5). Only the DO₃SE models and GEM-MACH Wesely show little soil uptake; while soil uptake is still a large fraction of v_d for GEM-MACH Wesely, it is a small fraction for the DO₃SE models. Stomatal uptake is very low except for DO₃SE psn. Stomatal fractions are between 0 % and 10 % except for DO₃SE psn (50 %) (Fig. 12). Because models largely agree that wintertime v_d is dominated by soil uptake and due to the fact that most models overestimate January–April v_d but underestimate November–December values, future work should focus on changes in soil uptake on weekly to monthly timescales. We do not have snow depth measurements at Easter Bush, but we do not expect that accounting for snow would substantially impact simulated values.

5.4 Ramat Hanadiv

Ramat Hanadiv is a shrubland in Israel near the Mediterranean coast. The spread based on the model with the highest annual average v_d divided by the model with the lowest is a factor of 2.2 (2.3 during summer and 2 during winter); however, based on the interquartile range, this value is a factor of 1.4 (1.3 during summer and 1.5 during winter). Metrics are on the lower end of the cross-site range.

There are ozone flux observations at Ramat Hanadiv during January–September only, and only March, August, and September have substantial data coverage (Fig. 1). A total of 3 different years contribute to multiyear averages, with each year only having a few months of data per year. For some months, years have overlapping data coverage. Some months with data for 2 years show interannual variability, whereas others do not. Like Bugacpuszta and Auchencorth Moss, more data are needed to assess interannual variability as well as seasonality at Ramat Hanadiv. Below, we examine multiyear averages, acknowledging that only 6 months of the year have 2 years of data and 3 months have data from 1 year only.

Models show weak relative seasonal amplitudes for v_d (Fig. 9). Values are very similar across models, more so than other sites. Most models also show weak relative seasonal amplitudes for stomatal uptake, but there is a larger spread across the central models and some outliers. The lack of simulated seasonality for most models is likely due to constant LAI. Any simulated v_d seasonality is from stomatal uptake (Fig. 6), more so than (or in contrast to) the other short-vegetation sites. GEM-MACH Wesely and WRF-Chem Wesely, which are two of three models with input initial resistances (i.e., model parameters) varying by season, have very distinct v_d seasonal cycle shapes at this site, compared with the rest of the models (Fig. 3).

The seasonal cycle shape of observed v_d at Ramat Hanadiv is hard to discern, with many months with low or no data coverage (Fig. 1). The current set of observations indicates higher values during early spring and lower values during late summer. Individual models do not capture this, with models simulating near-constant values year-round or increases from winter to early summer (Fig. 3). Exceptions are MLC-CHEM, the DO₃SE models, and GEM-MACH Wesely, which at least somewhat capture that the predominant seasonality feature should be lower late-summer values and higher early-spring values.

Across months with observations, models bracket observed v_d (Fig. 2). In particular, models are within −35 % to +55 % of observed seasonal averages (Fig. 4). Exceptions occur during summer and include GEM-MACH Wesely, IFS GEOS-Chem Wesely, WRF-Chem Wesely, GEOS-Chem Wesely, the TEMIR Wesely models, and the TEMIR Zhang models (biases are higher than +55 %). The central models' spread only brackets observed values during January–April and June, whereas it is too high during May and July–September. The largest deviation happens during August. Thus, like Bugacpuszta, late summer is when the largest model biases occur at Ramat Hanadiv.

The DO₃SE models, MLC-CHEM, and TEMIR psn show weak v_d decreases from spring to fall (Fig. 3). These models and the CMAQ models consider stomatal conductance dependencies on soil moisture. CMAQ models show weaker v_d declines from spring to fall, compared with the DO₃SE models, MLC-CHEM, and TEMIR psn. This behavior is consistent with their soil moisture dependencies. For example, the TEMIR psn and IFS SUMO Wesely models' stomatal conductance is set to zero when input soil moisture is less than the wilting point, but the CMAQ models have more of a taper effect. Future work should aim to understand the role of soil moisture in the observed seasonal variation in v_d and stomatal uptake.

Models with the highest biases during April–September are the TEMIR models, GEM-MACH Wesely, WRF-Chem Wesely, GEOS-Chem Wesely, and IFS GEOS-Chem Wesely (Fig. 3). These models simulate the highest stomatal uptake during this period, apart from a few models with lower-than-average nonstomatal uptake (CMAQ STAGE, the DO₃SE models, and GEM-MACH Zhang) (Fig. 5). Only the CMAQ M3Dry models capture low observed v_d during August. CMAQ M3Dry-psn captures July but CMAQ M3Dry does not, and they do not capture observed values during other months. Notably, the CMAQ M3Dry models show much lower summertime stomatal uptake than other models. CMAQ M3Dry models may have more skill during summer than other models, but, like the other models, they struggle with seasonality.

Lower canopy uptake is the highest for Ramat Hanadiv, during both summer and winter, across sites (Fig. 5). However, the relative and absolute contributions of lower canopy uptake are still low compared with soil and stomatal uptake (and, in some cases, cuticular uptake). Lower canopy uptake is only simulated by the Wesely models. Mostly the Wesely models simulate low cuticular uptake compared with other models, so lower canopy uptake does not necessarily contribute to the very high model biases of the Wesely models.

Uptake by soil and stomata mostly comprises v_d at Ramat Hanadiv during winter and summer (Fig. 5). The model spread is highest for stomatal uptake during winter and summer, compared with other pathways (Fig. 11). The spread for soil uptake is remarkably low given its importance across models (less than 20 % relative spread compared with mostly between 40 % and 75 % of v_d). Ramat Hanadiv is the only site with a large wintertime spread across stomatal uptake estimates and with similar model ranges of stomatal fractions during winter and summer. Models except WRF-Chem Wesely show substantial wintertime stomatal uptake. In general, stomatal uptake is very high compared with other sites during winter, presumably due to the site's Mediterranean climate. Models also show substantial summertime stomatal uptake, except CMAQ M3Dry. Wintertime stomatal fractions range from 20 % to 50 % across models (Fig. 12). The range is only slightly less across central models (25 %–40 %), suggesting that wintertime stomatal uptake is a key uncertainty at this site. The central models simulate a very small range of summertime stomatal fractions (similar to only Easter Bush), centering on 40 %, but the full range spans 12.5 % to 50 %.

At Ramat Hanadiv, most models should simulate lower stomatal and/or nonstomatal uptake during late summer, on par with the CMAQ M3Dry models, which have both lower stomatal and nonstomatal uptake than other models. However, stomatal and/or nonstomatal uptake should be higher than simulated by CMAQ M3Dry during other times of year, and other models bracket observations well at this time so they may provide insight here as to driving processes. Observational constraints on stomatal uptake year-round will help to further narrow uncertainties as to whether and when models need improvement with respect to stomatal versus nonstomatal uptake, including when they capture the absolute magnitude of v_d well.

5.5 Ispra

Ispra is a deciduous broadleaf forest in northern Italy. The model spread in terms of the model with the highest annual average v_d divided by the model with the lowest is a factor of 2.3 (3.1 during summer and 2.9 during winter); however, based on the interquartile range, this value is 1.5 (1.5 during summer and winter). These metrics are towards the higher end of the metrics for other sites.

Observed multiyear monthly mean v_d values are similar year-round except during March and April when values are lower (Fig. 1). This observed climatological seasonal pattern is consistent across years except during October–December. For example, observed v_d is high during October 2013, low during November 2015, and high during December 2014. As discussed below, the causes of high year-round values are uncertain; this, along with strong interannual variability during fall, indicates a need for more years of observations at Ispra, coupled with complementary measurements targeting individual pathways. In the following, we focus on multiyear averages, after briefly evaluating summertime interannual variability.

Summertime observed v_d at Ispra is higher during 2014 than during 2013 and 2015 (Fig. 1). Accordingly, model skill with respect to interannual variability should be determined by whether models capture the much higher summertime average during 2014 versus other years. Some models suggest that v_d should be highest during 2014, but hardly any models capture the large observed relative difference between this year and other years (Fig. 13). The exception is MLC-CHEM, and to a lesser extent GEM-MACH Zhang. Thus, most models have little skill regarding simulating summertime interannual variability at this site.

The v_d seasonality shape is a clear discrepancy between observations and models at Ispra. In contrast to the observations, multiyear monthly mean v_d peaks during the warm months in the central models (Fig. 2). There are similar v_d relative seasonal amplitudes across models, aside from GEM-MACH Wesely (Fig. 9), especially relative to other forests. The central models bracket the observations during April–September, but models show a low bias during October–March. Relative summertime and springtime biases range from −33 % to +32 % except DO₃SE multi, TEMIR Zhang, TEMIR Wesely BB, and GEM-MACH Zhang (lower) as well as GEM-MACH Wesely (higher) (Fig. 4). Relative wintertime and fall biases range from −22 % to −89 % across models. Ispra is the only site besides Auchencorth Moss where models are biased in the same direction for an extended period (i.e., longer than 3 months).

Models show that stomatal uptake largely drives v_d seasonality at Ispra (Fig. 6). Models simulate contributions from cuticular uptake, mostly via positive correlations with the stomatal pathway. Models with nonzero individual contributions from cuticular uptake (GEM-MACH Zhang, the CMAQ models, and the DO₃SE models) are the same as at Harvard Forest and Borden Forest. Models show v_d maxima during the warm months because v_d strongly depends on LAI (Fig. 7), which has a broad maximum during warm months (Fig. 10). Specifically, simulated v_d tends to increase with LAI, which contrasts with observed v_d.

A couple of models deviate from the majority in terms of the v_d seasonal cycles (Fig. 3). For example, GEM-MACH Zhang is low during the warm months and GEM-MACH Wesely is very high during the warm months. WRF-Chem Wesely shows higher wintertime v_d than other models, especially in January–March, due to high soil uptake, as well as high early-springtime uptake due to combined high soil and stomatal uptake (Figs. 5, 10). GEM-MACH Wesely and WRF-Chem Wesely are two of three models with input initial resistances (i.e., model parameters) varying by season, which likely causes these models to produce distinct seasonal cycle shapes. GEM-MACH Zhang has low summertime stomatal and nonstomatal uptake, compared with the rest (Fig. 5).

Even though the central models bracket the observed multiyear monthly mean v_d during April–September at Ispra (Fig. 2) and many individual models capture the increase from April to May, individual models fail to capture the fact that values should be roughly constant from July to September, rather than decreasing (Fig. 3). For example, some models (including DO₃SE psn and MLC-CHEM) simulate the April–July multiyear monthly mean v_d values very well but not August and September when they are low (because they simulate decreases from early to late summer). Models may erroneously simulate decreases from early to late summer because they depend too strongly on LAI, which weakly declines from July to September, or soil moisture.

During summer at Ispra, the model spread is largest for stomatal uptake relative to other pathways (Fig. 11). Models simulate substantial stomatal uptake, with DO₃SE multi and GEM-MACH Zhang simulating the lowest (but non-negligible) values (Fig. 5). The highest stomatal uptake is simulated by GEM-MACH Wesely, GEOS-Chem Wesely, IFS GEOS-Chem Wesely, IFS SUMO Wesely, TEMIR Wesely, and MLC-CHEM. The central models show stomatal fractions of 50 % to 77.5 %, but the full model range is 37.5 % to 87.5 % (Fig. 12). The model spread across pathways is the second largest for cuticular uptake. Soil uptake is very low across models, except for WRF-Chem Wesely, CMAQ STAGE, and GEM-MACH Wesely, where it is higher. The ranking and spread across pathways of pathways' standard deviations at Ispra are very similar to Borden Forest and Harvard Forest but not to Hyytiälä. Given that the central models capture the average magnitude of v_d during the warm season well but disagree mainly on stomatal versus cuticular fractions as well as monthly changes within the warm season (or lack thereof), future work should prioritize using observational constraints on stomatal uptake to further evaluate model performance.

During winter at Ispra, simulated v_d tends not to be dominated by one pathway; instead, there are small contributions from two to four pathways (Fig. 5). Exceptions are WRF-Chem Wesely where soil uptake dominates and a few models where cuticular uptake tends to dominate (e.g., CMAQ STAGE, CMAQ M3Dry, and DO₃SE multi). The model spread in soil uptake is largest across pathways (Fig. 11), and high WRF-Chem Wesely values play a role in this. Otherwise, soil uptake is low or, in a few cases, moderately low (e.g., MLC-CHEM and IFS SUMO Wesely). Cuticular uptake is close behind soil uptake in terms of the spread. Stomatal fractions span 0 % to 47.5 %, with the largest range across the central models (10 %–45 %) across sites (Fig. 12). A total of 11 models show low to moderately low stomatal uptake, but others predict none (GEM-MACH Wesely, GEM-MACH Zhang, CMAQ STAGE, GEOS-Chem Wesely, CMAQ M3Dry, TEMIR Wesely, and DO₃SE multi). More models predict nonzero stomatal uptake at Ispra compared with other sites, apart from Ramat Hanadiv. Whether simulated wintertime stomatal, cuticular, soil, and/or lower canopy uptake should be higher at Ispra is uncertain. There may also be fast ambient losses of ozone. Ispra does not have snow depth observations, but we anticipate that accounting for snow would not substantially change model results. Future attention should be placed elsewhere with respect to developing a better understanding of large wintertime model biases. A key first step is to understand whether there is stomatal uptake during winter and, if so, its magnitude.

5.6 Hyytiälä

Hyytiälä is a boreal evergreen needleleaf forest in Finland. The model spread in terms of the model with the highest annual average v_d divided by the model with the lowest is a factor of 2.7 (1.9 during summer and 21 during winter); however, based on the interquartile range, the value is a factor of 1.6 (1.4 during summer and 2.4 during winter). The metrics of model spread at Hyytiälä are at the higher end of other sites' values, especially for annual and winter values.

Observed multiyear monthly mean v_d maximizes during the warm months, and this is consistent across years (Fig. 1). Most models simulate higher values during warm months relative to cool months (Fig. 3). Outliers with respect to the seasonality are TEMIR Zhang (strong overestimate during cold months leading to near-constant values year-round), GEM-MACH Wesely (strong overestimate during warm months), GEOS-Chem Wesely and TEMIR Wesely (overestimate during summer), and WRF-Chem Wesely (strong overestimate during early spring). Here, we examine observed relative seasonal amplitude for v_d because observed and (most) modeled values have warm-month maxima and cool-month minima as well as full years of observations, allowing meaningful comparisons. The observed relative seasonal amplitude falls within the central models' range, although towards the upper end, and most models predict overly low values (Fig. 9).

In general, the largest relative model v_d biases at Hyytiälä occur during the cool months (Fig. 4) and the wintertime v_d model spread is the highest relative to other sites (Fig. 11), implying that wintertime v_d at this site is a key uncertainty. Wintertime relative biases range from −81 % to +87 % except for a few models that have much higher positive biases: GEM-MACH Zhang (+307 %), the TEMIR Zhang models (+211 % to +245 %), and DO₃SE psn (+104 %). However, most models are biased high, apart from IFS SUMO Wesely (−5 %), IFS GEOS-Chem Wesely (−81 %), GEOS-Chem Wesely (−62 %), and the TEMIR Wesely models (−15 % to −57 %). Models largely simulate that cuticular and soil uptake are dominant contributors (Fig. 5). Most models simulate near-zero wintertime stomatal uptake, despite relatively high LAI (Fig. 10), implying that models have at least rudimentary skill with respect to capturing the seasonality of evergreen vegetation. The central models show stomatal fractions between 0 % and 12.5 %, but a few models show contributions of 17.5 % to 50 % (Fig. 12). The model with the 50 % contribution (TEMIR Wesely BB) in addition to very low stomatal uptake has very low nonstomatal uptake.

During winter, models also show differences in the partitioning and magnitudes of cuticular versus soil uptake (Fig. 5). The model spread in cuticular uptake is larger than soil uptake (Fig. 11); Hyytiälä is the only site where this happens, presumably because LAI remains relatively high year-round and models seem to suggest that cuticular uptake is more important than ground uptake at forests. A total of 10 models show substantial cuticular uptake, whereas only 2 models show low cuticular uptake, and the rest show none. A total of 7 models show substantial soil uptake, whereas 10 show very little to none. Models showing high versus low cuticular and soil uptake are sometimes the same. For example, four simulate substantial cuticular uptake and soil uptake, and five simulate minimal cuticular uptake and soil uptake. In the former case, models overestimate wintertime v_d; in the latter case, models underestimate it. Most models capture small observed decreases in wintertime v_d with snow, but the spread across models during snow and snow-free periods is very large (Fig. 8). Thus, attention should focus on constraining wintertime cuticular versus soil uptake. Establishing whether there is cuticular and/or soil uptake during winter is an important first step towards narrowing model uncertainties.

Within the warm season, whether models show pronounced v_d seasonality varies (Fig. 3). Models also do not capture the fact that observations maximize during August and minimize during March. Specifically, models tend to overestimate late-winter/spring v_d while underestimating fall/early-winter v_d, as indicated by comparing the interquartile range to observations (Fig. 2). The multiyear monthly mean LAI peaks during August (around 3.75 m² m⁻²), after an increase from May (Fig. 10). Then, LAI decreases to November and is constant from November to May (around 2.75 m² m⁻²). Models bound the observed v_d–LAI relationship and largely capture the increase in v_d as LAI increases from 3 to 3.5 m² m⁻² (Fig. 7). However, most models do not capture the v_d change as LAI increases from 3.5 to 3.75 m² m⁻² where observations suggest that the slope should be the same as for 3 to 3.5 m² m⁻² (instead models suggest decreases). Models also overestimate the increase in v_d as LAI increases from 2.75 to 3 m² m⁻². Some effect overrides LAI's influence on seasonality in stomatal uptake in models, given that both the observed LAI and v_d peak during August but simulated stomatal uptake and v_d do not. Simulated declines with soil moisture may play a role here.

Models simulate that stomatal uptake and covariations between pathways are important seasonality drivers (Fig. 6). Only two models suggest that there are not individual contributions by stomatal uptake (GEM-MACH Wesely and GEM-MACH Zhang), but several models suggest that the sum of individual contributions from other pathways and covariations are at least as important as stomatal uptake. There are similarly evenly distributed spreads across models in terms of relative seasonal amplitudes for stomatal uptake and v_d (Fig. 9). Most models' stomatal uptake seasonal cycles show a broad warm-season peak, apart from some models with more pronounced seasonality during the warm months (e.g., GEM-MACH Wesely, GEOS-Chem Wesely, TEMIR Wesely, and the CMAQ M3Dry models) (Fig. 10). IFS SUMO Wesely peaks during May and then declines afterwards. Model outliers in terms of high magnitudes of summertime stomatal uptake include GEOS-Chem Wesely, TEMIR Wesely, MLC-CHEM, and GEM-MACH Wesely.

During summer, relative model biases range from −14 % to +20 % except for GEM-MACH Wesely (+88 %), IFS SUMO Wesely (−25 %), WRF-Chem Wesely (+32 %), TEMIR Wesely (+34 %), and GEOS-Chem Wesely (+40 %) (Fig. 4). Models show substantial stomatal uptake (Fig. 5), with stomatal fractions spanning from 27.5 % to 80 % (Fig. 12). The central models show 42.5 %–65 %. Models that simulate lower canopy uptake show low uptake via this pathway, like other forests. The largest model spread is for soil and stomatal uptake, but this is closely followed by cuticular uptake (Fig. 11), which is distinct from other forests. Soil uptake's high model spread is due to high values from WRF-Chem Wesely and GEM-MACH Wesely and zero values from the DO₃SE models; other models simulate more similar estimates of soil uptake, ranging from low to moderate. Models show non-negligible cuticular uptake but disagree as to whether it is low or moderate. Observational constraints on stomatal uptake will help to further narrow uncertainties as to the magnitude and relative contribution of summertime stomatal uptake, as well as changes on weekly to monthly timescales.

Key findings regarding seasonality at Hyytiälä include the following: models struggle to capture the exact timing of maximum and minimum values, models overestimate wintertime values and thus underestimate the relative seasonal amplitude, and models disagree about seasonality within the warm season while generally capturing that there should higher values during warm months. Silva et al. (2019) use Hyytiälä observations to train a machine learning model and apply the model to predict v_d at Harvard Forest, finding that their model predicts a late summertime peak in v_d, which is observed at Hyytiälä but not at Harvard Forest. Assuming that differences between these two sites are characteristic of sites' broad LULC classifications, both our findings and theirs suggest a need for improved predictive ability of seasonal differences between coniferous and deciduous forests.

Thus far, we have discussed multiyear averages at Hyytiälä. We now turn to summertime interannual variability. Models do not capture the summertime ranking across years (Fig. 13). Several models predict particularly low (high) v_d during some summers, but the observations do not indicate low (high) values for these years. Some models are close to capturing the degree of summertime interannual variability, but these models typically show a more uneven distribution across years than suggested by observations. Notably, models show more variability in their year-to-year rankings at Hyytiälä compared with other sites with longer records. Nonetheless, we conclude that model skill is poor at this site in terms of summertime interannual variability.

5.7 Harvard Forest

Harvard Forest is a temperate mixed forest in the northeastern USA. The model spread in terms of the model with the highest annual average v_d divided by the model with the lowest is a factor of 1.9 (1.8 during summer and 4.8 during winter); however, based on the interquartile range, this value is a factor of 1.2 (1.4 during summer and 2.6 during winter). Like other forests, the wintertime spread is largest. Aside from winter values, the metrics of the spread at Harvard Forest are on the lower end of estimates across sites.

Observed multiyear monthly mean v_d maximizes during May–September (Fig. 1). Observed seasonal cycles vary across years, but values are generally higher during the warmer versus cooler months. We focus on multiyear averages until the subsection end, where we touch on summertime interannual variability. Models capture that v_d peaks during the warm months (Fig. 2). The exception is GEM-MACH Zhang, which has similar monthly averages year-round. Despite capturing the seasonality shape, models overestimate the relative seasonal amplitude (Fig. 9), apart from GEM-MACH Zhang, TEMIR Zhang, and TEMIR Zhang BB (substantial underestimate) as well as DO₃SE psn (slight underestimate). Outliers show high wintertime v_d relative to other models and observations, implying that the full set of models bounding the observed relative seasonal amplitude do not necessarily indicate ensemble skill.

Models are within ±65 % of observed values across seasons (Fig. 4). Exceptions occur during spring and summer for GEM-MACH Wesely, winter and spring for GEM-MACH Zhang, and spring for WRF-CHEM Wesely and TEMIR Zhang Medlyn. The central models bracket observations well (Fig. 2). Specifically, observations fall at the lower end of the spread during the warm months and the upper end during November–January, but they are otherwise are in the middle of the spread. Across models, summertime biases are positive, ranging from +4 to +144 %, except IFS GEOS-CHEM Wesely (−4 %) and TEMIR Zhang (−2 %). Thus, overestimated relative seasonal amplitudes (Fig. 9) are likely due to high summertime v_d. Previous work suggests that GEOS-Chem's overestimate at Harvard Forest is due to overly high model LAI (Silva and Heald, 2018), but there is clearly another issue because models are forced with site-specific LAI here. Most models tend to underestimate v_d at low LAI values and overestimate v_d at high LAI values, overstating v_d increases with LAI (Fig. 7).

During winter, model biases tend to be negative, ranging from −24 % to −71 %, with the exception of GEM-MACH Wesely (+85 %), the TEMIR Zhang models (+25 % to +33 %), and MLC-CHEM (+13 %) as well as two models with very low negative biases (DO₃SE psn and WRF-Chem Wesely) (Fig. 4). The wintertime model spread is highest for soil uptake across pathways, with cuticular uptake close behind (Fig. 11). Soil uptake is always at least 37.5 % (and up to 70 %) of v_d except for GEM-MACH Wesely (20 %) (Fig. 5). Most models show little to no stomatal uptake, but some models show non-negligible values. The central models show stomatal fractions of 5 %–15 % (Fig. 12). Estimates for cuticular uptake vary across models: there are substantial, small, and negligible contributions. Lower canopy uptake is low for models that simulate this pathway but can be an important fraction of v_d. There are no snow depth observations at Harvard Forest. Assuming no snow throughout the time period may influence some models' ability to estimate wintertime v_d well. However, based on our analysis at other sites, we do not anticipate the lack of snow data to be the main driver of model–observation or model–model differences. Establishing whether there should be stomatal or cuticular uptake during winter would be a useful first step in further constraining models. Otherwise, attention should focus on narrowing the uncertainties related to wintertime ground uptake.

Some models capture the broad observed v_d maximum during the warm season, while others show more seasonality within the warm season (Fig. 3). A few models show pronounced declines after July (e.g., MLC-CHEM and TEMIR psn). Pronounced declines after July do not occur in observed multiyear monthly averages but do occur during several individual years (Fig. 1). Simulated pronounced declines may follow these models' soil moisture dependencies. (Note that not all models have soil moisture dependencies, and there are differences among models that do have them.) The fact that models with soil moisture dependencies do not capture the observed multiyear mean seasonality may be due to the soil moisture dependencies themselves and/or owing to uncertainty in soil moisture input. For example, soil moisture was not measured during all years with ozone fluxes at Harvard Forest, and thus we use a climatological average during those years. Future work should examine seasonality during individual years, paying attention to years with climatological average versus year-specific input soil moisture, to determine model strengths and limitations.

Models show that stomatal uptake is an important driver of v_d seasonality at Harvard Forest (Fig. 6). Six models estimate that stomatal uptake largely drives seasonality, with some contributions from covariations between pathways (mainly positive covariations between stomatal and cuticular pathways). The rest estimate moderate contributions from stomatal uptake but at least as much of an influence from individual nonstomatal pathways or covariations (positive or negative). Models show a clear seasonality in stomatal uptake, with a peak during the warm months and zero or near-zero values during winter (Fig. 10). The spread in the relative seasonal amplitude for stomatal uptake across the central models is the smallest across sites (Fig. 9). However, six models deviate from the rest: CMAQ M3Dry, CMAQ STAGE, and GEM-MACH Wesely have high relative seasonal amplitudes for stomatal uptake, whereas GEM-MACH Zhang, IFS SUMO Wesely, and DO₃SE psn have low values. In contrast, the spread in the relative seasonal amplitude for v_d has a more even distribution across models. Thus, while there is a fair amount of agreement across models in terms of seasonality in stomatal uptake, models disagree with respect to nonstomatal uptake seasonality and its role in v_d seasonality. Together with findings that models exaggerate the v_d–LAI relationship and most models overestimate the relative seasonal amplitude for v_d, this result implies that future work should aim to better constrain the nonstomatal influences on seasonality.

During summer, the model spread is highest for stomatal uptake, with cuticular uptake close behind (Fig. 11). Models show substantial contributions from stomatal uptake: the model range spans 30 % to 80 %, but the central models' range spans 50 % to 70 % (Fig. 12). Estimates for cuticular uptake vary across models (Fig. 5): there are substantial, moderate, and low contributions. Soil uptake is low, except for WRF-Chem Wesely and GEM-MACH Wesely. Similar to other forests, lower canopy uptake is low for models that simulate this pathway. Observational constraints on stomatal uptake will help to further narrow model uncertainties as to the magnitude and relative contribution of summertime stomatal uptake.

Interannual variability is strong across months (Fig. 1). A series of papers pointed this out for daytime values and investigated drivers during summer (Clifton et al., 2017, 2019). Models capture neither the large observed spread across years during summer nor the ranking of years (Fig. 13). Most models simulate that some of the summers with the highest observed v_d have low v_d. Previous work points to nonstomatal pathways driving summertime interannual variability (Clifton et al., 2017, 2019), and thus models may be lacking in their ability to simulate the degree to which nonstomatal uptake varies from year to year, and likely key process dependencies.

5.8 Borden Forest

Borden Forest is a mixed forest in the boreal–temperate transition zone in Canada. The model spread in terms of the model with the highest annual average v_d divided by the model with the lowest is a factor of 2.3 (3.4 during summer and 10 during winter); however, based on the interquartile range, this value is a factor of 1.4 (1.8 during summer and 3 during winter). The metrics of model spread are towards the higher end of other sites, except for winter and the summertime interquartile range when they are the highest.

The observed multiyear monthly mean v_d shows a broad maximum during the warm months at Borden Forest (Fig. 1), like Harvard Forest and Hyytiälä. However, uniquely, observations at Borden Forest show particularly large winter versus summer differences and steep changes during spring and fall. Specifically, v_d increases from March to June by 0.5 cm s⁻¹. Then, v_d remains high from June to September (0.6–0.65 cm s⁻¹) and declines steeply from September to November. Models simulate higher v_d during the warmer versus cooler months (Fig. 3), and the observed relative seasonal amplitude lies close to the middle of the central models' spread (Fig. 9). However, there is a clear discrepancy between models and observations in that models do not capture very high v_d across the warm months (Fig. 2). All models except GEM-MACH Wesely have low summertime biases, with a range from −15 % to −74 % (Fig. 4). In general, high observed v_d during the warm months at Borden Forest requires a better understanding, given the uncertainty in ozone flux measurements from the gradient technique (see discussion in Sect. 4.2).

The individual contribution from stomatal uptake is found to be a key driver of v_d seasonality, apart from IFS SUMO Wesely, CMAQ STAGE, and the DO₃SE models (Fig. 6). These four models do, however, show stomatal contributions to seasonality via correlations with other pathways. Notably, there are more individual nonstomatal (e.g., ground, cuticular) contributions to seasonality at Borden Forest than at other forests. There are also a variety of simulated v_d seasonal cycle shapes at Borden Forest, in contrast to Harvard Forest and Ispra. Some models simulate weak changes from cooler to warm months (the DO₃SE models, the TEMIR Zhang models, IFS SUMO Wesely, and GEM-MACH Zhang), whereas others simulate moderate changes (WRF-Chem Wesely, MLC-CHEM, and CMAQ STAGE) or strong changes (GEOS-Chem Wesely, TEMIR Wesely, IFS GEOS-Chem Wesely, GEM-MACH Wesely, the CMAQ M3Dry models, and TEMIR Wesely psn). The TEMIR psn models simulate erratic monthly changes from June to October. Generally, models with the strongest changes from cooler to warm months simulate that stomatal uptake predominately drives v_d seasonality (Fig. 6). Conversely, models with weak changes from cooler to warm months indicate that nonstomatal pathways contribute more predominantly.

With respect to the relationship between the multiyear monthly mean v_d and LAI, observed v_d increases with LAI but the slope varies (Fig. 7). The observed slope is strongest for LAI increases from 0.5 to 1 m² m⁻², and models tend to underestimate the change but do simulate increases. The observed slope weakens but remains positive for LAI increases from 1 to 2 m² m⁻² – most models suggest decreases instead. The observed slope then weakens even further for LAI increases above 2 m² m⁻². Some models capture the slope of LAI increases above 2 m² m⁻², but others exaggerate it (e.g., GEM-MACH Wesely, GEOS-Chem Wesely, TEMIR Wesely, and the CMAQ M3Dry models). The main issue is that individual models tend not to capture that there should be relatively high v_d during May and October (Fig. 3). Specifically, models simulate a later spring onset with respect to the v_d seasonality as well as an earlier fall decline, and thus a shorter season of elevated v_d than observed. Therefore, we suggest that models are too strongly tied to the LAI, which strongly increases from May to June and strongly decreases from September to October (Fig. 10).

Additionally, many models do not capture that the multiyear monthly mean v_d is similar during June–September (Fig. 3). Some models simulate declines from August to September (e.g., CMAQ M3Dry-psn, GEOS-Chem Wesely, TEMIR Wesely, and GEM-MACH Wesely). A weak decline from August to September occurs in the observed multiyear average (the strong decline happens from September to November); some models capture the August–September decline's magnitude, whereas others exaggerate it. Some models show low values during July (e.g., TEMIR psn) in addition to August–September declines. Observations show low values during July, not with respect to multiyear monthly mean seasonal cycles but during 2012 and perhaps 2008 (Fig. 1). Many models show peak v_d during June. Again, this does not happen in the observed multiyear monthly averages, but it occurs in 2010. Thus, models may exaggerate depositional responses (in particular, stomatal responses) to changes in environmental conditions (e.g., soil moisture) on a climatological basis but have some skill in certain years.

During summer, the largest model spread across pathways occurs for stomatal uptake, followed by cuticular uptake and then soil uptake (Fig. 11), similar to Harvard Forest and Ispra. Models show substantial stomatal uptake, apart from two with very low values (IFS SUMO Wesely and DO₃SE multi). Stomatal fractions range from 20 % to 80 % across models but from 40 % to 62.5 % across the central models (Fig. 12). Eight models simulate lower cuticular uptake, whereas the rest simulate higher cuticular uptake (Fig. 5). Models that have the lower canopy uptake pathway show low values of cuticular uptake, with two exceptions: GEM-MACH Wesely, which has high cuticular uptake, and MLC-CHEM, which does not archive the lower canopy uptake diagnostic but has low cuticular uptake. Most models simulate low soil uptake, but a few models simulate moderate to high soil uptake (GEM-MACH Wesely, GEM-MACH Zhang, CMAQ STAGE, WRF-Chem Wesely, and MLC-CHEM). Observational constraints on stomatal uptake will help to further narrow model uncertainties with respect to the magnitude and relative contribution of stomatal uptake.

During winter, models show a mixture of over- and underestimates. Models with overestimates are the TEMIR Zhang models (+68 % to +73 %), GEM-MACH Zhang (+124 %), WRF-Chem Wesely (+13 %), DO₃SE multi (+9 %), and DO₃SE psn (+44 %) (Fig. 4). Otherwise, underestimates span from −20 % to −78 %. Models with high v_d simulate high cuticular uptake, generally high soil uptake, and, in one case (DO₃SE psn), non-negligible stomatal uptake (Fig. 5). Soil and cuticular uptake show the highest spreads across models, with soil uptake the highest, similar to Harvard Forest and Ispra (Fig. 11). The central models show very low stomatal fractions, but outliers span 10 % to 30 % (Fig. 12). Many models largely capture that observations show no v_d change with snow, although some slightly overestimate the change. Thus, the primary issue with wintertime model biases is likely unrelated to responses to snow; rather, it is likely related to mischaracterized magnitudes of pathways or responses to other environmental conditions.

In terms of summertime interannual variability, most models underestimate the relative spread across years (Fig. 13), but some only slightly underestimate it (IFS SUMO Wesely, CMAQ STAGE, TEMIR Zhang, MLC-CHEM, and the DO₃SE models) and a few exaggerate it (TEMIR psn). Models generally struggle to capture the observed relative distribution across summers (i.e., 2 high years, 2 low years, and 1 middle year). No model captures the year-to-year ranking across summers, but many capture 1 of the high years and, in some cases, 1 of the low years. CMAQ STAGE captures a second high year, whereas no other model captures this (or distinguishes it from other years). Given variability within summer in the yearly observations (Fig. 1), future work should examine interannual variability in monthly averages to further establish model skill.

6 Conclusion

We introduce AQMEII4 Activity 2 for the intercomparison and evaluation of 18 dry deposition schemes configured as single-point models driven by the same set of meteorological and environmental conditions at eight sites with ozone flux records. We provide our approach's rationale, document the single-point models, and describe the observational datasets used to drive and evaluate the models. The emphasis on driving models with a consistent set of inputs in Activity 2 allows us to focus on parameter and process uncertainty.

We launch the Activity 2 results by analyzing simulated multiyear mean ozone deposition velocities and effective conductances for plant stomata, cuticles, the lower canopy, and soil, as well as observed multiyear mean ozone deposition velocities. Our focus is monthly and seasonal averages across all hours of the day, apart from one site for which we examine afternoon averages (Ramat Hanadiv). We evaluate the magnitudes and seasonal cycles (e.g., shape and amplitude) of simulated ozone deposition velocities against observations. Moreover, we identify how differences and similarities in the relative and absolute contributions of individual deposition pathways and how some dependencies on environmental conditions influence the model spread and comparison with observations. We encourage future work to examine the roles of parameters, sensitivities, and transport-related processes. For example, previous work has shown that differences in deposition velocities among air quality models under stable conditions may, at least in part, be due to different empirical formulations of Monin–Obukhov similarity theory (Toyota et al., 2016).

There are a variety of observed climatological seasonal patterns and magnitudes of ozone deposition velocities across the sites. We emphasize that our measurement test bed is likely insufficient to generalize results to specific LULC types, so we focus on site-specific results. We also cannot discount the fact that differences in ozone flux methods and instrumentation and a lack of coordinated processing protocols across datasets limit meaningful synthesis of our results across sites. However, given that key processes and parameters are strongly tied to the LULC type in dry deposition parameterizations, a core question is whether the magnitude and dependencies of ozone deposition velocities can be described from a LULC-type perspective. To address this question, future work will need to better understand observed site-to-site differences in ozone deposition velocities, which will likely require new multiscale ozone flux datasets.

We also emphasize the incomplete understanding of observed variations in ozone deposition velocities at several sites. Namely, there are unexpectedly high ozone deposition velocities year-round at Auchencorth Moss, during the cool season at Ispra, and during the warm season at Borden Forest; models do not capture these high values. Further model evaluation at these sites requires a better understanding of these features in the observations as well as whether the models should capture them.

Observed interannual variation in ozone deposition velocities is strong at most sites examined here, demonstrating the importance of long-term ozone flux records for model evaluation. For example, even if a model captures values for a given year, the model may not reproduce the interannual variability nor the multiyear average. Our focus of this first paper is climatological evaluation, with the caveat that three sites (Ramat Hanadiv, Auchencorth Moss, and Bugacpuszta) do not have multiple years of data for several months and two sites are missing some months of data across all years. Of course, full annual records with several years of data are required for confident constraints on climatological seasonality. Nonetheless, sites with short-term records have very similar monthly averages between years when there is good data coverage, with only a few exceptions (October at Auchencorth Moss and fall at Ispra), implying some utility of these datasets towards our aim.

Despite the focus on climatological evaluation, for sites with more than three summers of data, we briefly identify whether models capture the ranking and spread across summers. We find that models do not capture observed summertime interannual variability, a finding that agrees with earlier work with one model at Harvard Forest (Clifton et al., 2017). Our work here shows that the issue is widespread across models and sites. Specifically, we show poor model skill with respect to simulating the degree of the interannual spread as well as the ranking across years.

An important conclusion here is that the individual model performance strongly varies by season and site. Throughout this paper, we examine individual models as well as model ensembles including the full set of models as well as the interquartile range, which helps us to narrow our focus to key common uncertainties across models. The interquartile range across simulated averages of ozone deposition velocities ranges from a factor of 1.2 to 1.9 annually across sites, and it largely, reasonably bounds multiyear monthly mean ozone deposition velocities. Exceptions to the latter finding are times denoted as particularly uncertain at Auchencorth Moss, Ispra, and Borden Forest, in addition to late summer at Bugacpuszta and Ramat Hanadiv. The latter finding, along with our finding that many models that include soil moisture dependencies on stomatal conductance exaggerate late-summer decreases in ozone deposition velocities at forests, suggests a need to focus on refining soil moisture dependencies. Such work should probe interannual variability and seasonality with additional observational constraints on stomatal uptake in the context of uncertainty in soil moisture input data. In general, in some cases, gaps in site-specific measurement data (e.g., soil moisture and characteristics) forced us to make assumptions or derive estimates for key model variables and parameters. This may influence model performance, and it also points to a need for a standard minimum set of observations at future field studies.

Even beyond the differing effects of soil moisture across the ensemble of models, there are differences in the shapes of the simulated seasonal cycles of ozone deposition velocities. Models that rely strongly on seasonally dependent parameters are often identified as outliers; thus, we recommend that related canopy resistance equations should be tied to variables like the leaf area index instead of only seasonally varying parameters. In principle, seasonally varying parameters are not problematic, but a challenge seems to be indicating site-specific phenology accurately. At half the sites, the model spread is highest during the cooler months, implying a need for a better understanding of wintertime deposition processes. Strong wintertime sensitivities of tropospheric ozone abundances in regional to global chemical transport models (Helmig et al., 2007; Matichuk et al., 2017; Clifton et al., 2020b) also point to this requirement. By compositing observed and simulated ozone deposition velocities for all versus snowy conditions during the cool months at sites with snow depth observations, we show that models' inability to capture the magnitude of wintertime values is generally a larger issue than models' inability to capture responses to snow. While our analysis suggests that snow-induced changes are not the main driver of observed seasonality in ozone deposition velocities, we also find models may too strongly rely on the leaf area index to determine seasonality.

Several papers have illustrated challenges with respect to determining which ozone dry deposition parameterization is best given observations compiled from the literature (Wong et al., 2019; Cao et al., 2022; Sun et al., 2022) or comparing seasonal differences for ozone and sulfur dioxide deposition velocities at Borden Forest (Wu et al., 2018). While we agree with these earlier findings based on our more complete and diverse test bed, we take the evaluation a step further by pinpointing how different pathways contribute to the spread. In general, both stomatal and nonstomatal pathways are key drivers of variability in ozone deposition velocities across models. Additionally, in some cases, ozone deposition velocities are similar across models when the partitioning among deposition pathways is very different (i.e., similar results for different reasons).

For the most part, models simulate that stomatal uptake predominately drives seasonality in ozone deposition velocities. Like large model differences in the seasonality of ozone deposition velocities, there are large model differences in the seasonality of stomatal uptake. A few models show that seasonality in nonstomatal uptake terms is also important for seasonality in ozone deposition velocities. Across sites, both stomatal and nonstomatal pathways are important contributors to ozone deposition velocities during the growing season. For example, during summer, the median of the stomatal fraction of the ozone deposition velocity across models ranges from 30 % to 55 % across most sites. Thus, like observationally based estimates of the stomatal fraction over physiologically active vegetation compiled by a recent review (Clifton et al., 2020a), models clearly indicate a codominant role for dry deposition through nonstomatal pathways. Nonetheless, as stated in the previous paragraph, we emphasize large differences in simulated nonstomatal uptake, in addition to stomatal uptake, across models.

In general, we confirm here with our unprecedented full documentation of 18 dry deposition schemes that dry deposition schemes, especially nonstomatal deposition pathways, are highly empirical. While some schemes can capture some of the salient features of observations and could be adjusted to better capture the magnitude of observed ozone deposition velocities at the sites examined here, a better mechanistic understanding of observed variability and a firm grasp on how different deposition pathways change in time and space on different scales are needed to improve the predictive ability of ozone dry deposition. We will continue to chip away at this problem; next for Activity 2 will be to leverage observation-based constraints on stomatal conductance, along with inferred stomatal fractions of ozone deposition velocities, and examine diel, seasonal, and interannual variations to further evaluate single-point models.