Ozone deposition impact assessments for forest canopies require accurate ozone flux partitioning on diurnal timescales

Dry deposition is an important sink of tropospheric ozone that affects surface concentrations and impacts crop yields, the land carbon sink, and the terrestrial water cycle. Dry deposition pathways include plant uptake via stomata and non-stomatal removal by soils, leaf surfaces, and chemical reactions. Observational studies indicate that ozone deposition exhibits substantial temporal variability that is not reproduced by atmospheric chemistry models due to a simplified representation of vegetation uptake processes in these models. In this study, we explore the importance of stomatal and non-stomatal uptake processes in driving ozone dry deposition variability on diurnal to seasonal timescales. Specifically, we compare two land surface ozone uptake parameterizations – a commonly applied big leaf parameterization (W89; Wesely, 1989) and a multi-layer model (MLC-CHEM) constrained with observations – to multi-year ozone flux observations at two European measurement sites (Ispra, Italy, and Hyytiälä, Finland). We find that W89 cannot reproduce the diurnal cycle in ozone deposition due to a misrepresentation of stomatal and non-stomatal sinks at our two study sites, while MLC-CHEM accurately reproduces the different sink pathways. Evaluation of non-stomatal uptake further corroborates the previously found important roles of wet leaf uptake in the morning under humid conditions and soil uptake during warm conditions. The misrepresentation of stomatal versus non-stomatal uptake in W89 results in an overestimation of growing season cumulative ozone uptake (CUO), a metric for assessments of vegetation ozone damage, by 18 % (Ispra) and 28 % (Hyytiälä), while MLC-CHEM reproduces CUO within 7 % of the observation-inferred values. Our results indicate the need to accurately describe the partitioning of the ozone atmosphere–biosphere flux over the in-canopy stomatal and non-stomatal loss pathways to provide more confidence in atmospheric chemistry model simulations of surface ozone mixing ratios and deposition fluxes for large-scale vegetation ozone impact assessments.

1 Introduction

Ozone (O₃) in the atmospheric surface layer is an air pollutant that is toxic to humans and plants. Ozone is removed by oceans, bare soil, and vegetated areas, which together are called dry deposition and account for ±15 %–20 % of the total tropospheric ozone sink (Hu et al., 2017; Bates and Jacob, 2020). In vegetation canopies, the dominant deposition pathway is stomatal uptake, which typically accounts for 40 %–60 % of the total deposition to vegetation (Fowler et al., 2009). Stomatal ozone uptake reduces carbon assimilation in vegetation (Sitch et al., 2007; Ainsworth et al., 2012), affects the terrestrial water cycle (Lombardozzi et al., 2015; Sadiq et al., 2017; Arnold et al., 2018), and causes economic damage through reduced crop yield (e.g., Tai et al., 2014). Besides stomatal uptake, ozone removal occurs via a range of non-stomatal removal mechanisms, such as uptake by the leaf exterior and soils, and in-canopy chemical removal involving nitrogen oxides (NO_x) or plant-emitted reactive carbon species. The contribution of these ozone removal processes to the total non-stomatal term is uncertain (Fowler et al., 2009) and displays temporal variability on diurnal to interannual timescales that is incompletely understood (Clifton et al., 2020a). Given that these non-stomatal removal processes act in parallel to the stomatal removal of ozone, the characterization and quantification of non-stomatal sinks is important for quantification of total and stomatal ozone uptake.

The contribution of different ozone uptake pathways cannot be routinely measured at the plant canopy level due to the various non-stomatal uptake pathways. Most studies infer stomatal conductance (g_s) from canopy-top micro-meteorological and eddy covariance observations using an inverted form of the Penman–Monteith equation (e.g., Fowler et al., 2001; Clifton et al., 2017, 2019; Ducker et al., 2018), although some studies apply alternative g_s estimation methods based on gross primary production (GPP;  El-Madany et al., 2017; Clifton et al., 2017). In such observation-based studies, the non-stomatal ozone removal component (g_ns) is generally treated as the residual of the total uptake conductance (g_c; inferred based on the ozone dry deposition velocity) and g_s. However, sites with long-term ozone flux measurements are scarce (Clifton et al., 2020a), which limits the characterization of the seasonal to interannual temporal variability in the stomatal and non-stomatal components of ozone removal. Several campaign-based studies partitioned total canopy ozone fluxes by using ozone flux measurements along a vertical gradient to study the in-canopy flux divergence and to relate this to the vertical distribution of ozone sinks in the canopy (Fares et al., 2014; Finco et al., 2018), but these are limited to short timescales. Given the scarce availability of ozone deposition observations that span at least 1 year, and preferentially multiple years, quantifying temporal variability in stomatal and non-stomatal ozone deposition solely based on observations remains challenging.

Studies of ozone deposition (and its impacts) on regional to global scales rely on the application of atmospheric chemistry models and their dry deposition parameterizations. Many models treat deposition in a zero-dimensional manner and do not, or only implicitly, account for the variation in different in-canopy loss pathways as a function of environmental drivers and height within the canopy (the big leaf approach, Clifton et al., 2020a). Recent advances in the description of ozone deposition have been made by improving the simulation of stomatal conductance (Lin et al., 2019; Clifton et al., 2020c) and the representation of various non-stomatal removal terms (e.g., Zhang et al., 2003; Stella et al., 2011, 2019; Potier et al., 2015) and in-canopy turbulence and radiation extinction (Makar et al., 2017). Additionally, some models account for vegetation ozone damage via effects on photosynthesis and stomatal conductance (Lombardozzi et al., 2015; Sadiq et al., 2017; Arnold et al., 2018). Another class of models treats the canopy as a separate exchange regime with different biophysical and chemical conditions compared to the lowermost atmospheric layer and explicitly resolves in-canopy vertical gradients of ozone deposition and its driving variables by using multiple in-canopy layers (e.g., Ganzeveld et al., 2002, 2010; Fares et al., 2014; Otu-Larbi et al., 2020). Despite these advances in the representation of ozone deposition in atmospheric chemistry models, their application for ozone impact assessments remains a challenge. For example, the description of stomatal conductance is an important parameter for understanding year-to-year variability in impact metrics, such as cumulative uptake of ozone (CUO;  Clifton et al., 2020b), but stomatal versus non-stomatal ozone flux partitioning in these models is uncertain. Additionally, spatiotemporal controls of ozone deposition pathways remain incompletely understood (Clifton et al., 2017, 2020a), in part owing to the scarcity of long-term ozone flux observations. Therefore, we here study the temporal controls on stomatal and non-stomatal ozone deposition pathways, and their implications for simulations of CUO, using two multi-year ozone deposition data sets as well as a big leaf and multi-layer parameterization of land surface ozone uptake.

Specifically, we investigate the added value of an explicit multi-layer canopy representation of ozone deposition (MLC-CHEM – the Multi-Layer Canopy–CHemistry Exchange Model;  Ganzeveld et al., 2002) compared to a commonly applied big leaf parameterization (Wesely, 1989) in terms of simulating ozone deposition pathways and ozone impact metrics. We first study long-term (seasonal to annual) and short-term (diurnal) temporal variability in ozone dry deposition to forest canopies at a pristine boreal site (Hyytiälä) and a pre-alpine site that frequently experiences high ozone concentrations (Ispra). We then evaluate the performance of a big leaf and a multi-layer representation of atmosphere–biosphere exchange in simulating ozone dry deposition pathways and their temporal variability. Subsequently, we characterize the relationship of non-stomatal conductance as a function of environmental drivers. Last, we aim to demonstrate how representations of the drivers of long- and short-term variability in ozone stomatal and non-stomatal removal in those different land surface parameterizations affect simulated CUO. To this end, we employ multi-year canopy-top observations of micro-meteorology, ozone mixing ratios, surface energy balance components, and fluxes of ozone to derive the stomatal and non-stomatal components of the total ozone flux, combined with observation-driven ozone dry deposition simulations, using the two aforementioned representations.

2 Data and methods

2.1 Site description

Our study makes use of half-hourly observations of micro-meteorology (net radiation, air pressure, air temperature, relative humidity, precipitation, wind speed, and friction velocity) surface energy balance components and fluxes of CO₂ and ozone from two forested flux observation sites (Ispra and Hyytiälä), which are detailed below.

The Ispra forest flux station is situated in a deciduous forest in northern Italy (45.81^∘ N, 8.63^∘ E) at the European Commission Joint Research Centre (EC-JRC) in a 10 ha almost natural ecosystem mainly consisting of Quercus robur (80 %), Alnus glutinosa (10 %), Populus alba (5 %), and Carpinus betulus (3 %). Leaf area index (LAI) shows an average value of 4.1 m² m⁻² during the growing season (Fumagalli et al., 2016). In our analysis, we rely on continuous LAI measurements unavailable at this site, which we, therefore, take from a remote sensing product derived from MODIS (Xiao et al., 2014). The LAI range at Ispra in this product is 0.7–3.7 m² m⁻², scaled up to a locally measured LAI maximum of 4.5 m² m⁻² in July 2015 (Fumagalli et al., 2016), using a seasonally varying sinusoidal scaling function. The turbulent flux measurements of surface energy balance components and ozone were performed in 2013–2015 at 36 m a.g.l. (above ground level), approximately 10 m above the canopy height of 26 m. More information regarding the measurement setup of this site can be found in Gruening et al. (2012).

The Hyytiälä SMEAR II (Station for Measuring Forest Ecosystem Atmosphere Relations) measurement station is located in a needleleaf forest in southern Finland (61.85^∘ N, 24.28^∘ E), with a forest cover dominated by pine trees. LAI was periodically measured at this site and varies between 2.3 and 4 m² m⁻². Ozone flux measurements are available for 2002–2012, with a 1-year data gap in 2006. Turbulent flux measurements are performed at 23 m a.g.l., which is 5–9 m above the forest top of 14–18 m. Ozone mixing ratios at this altitude are derived by linearly interpolating between observations at 16.8 and 33 m. More information about the measurement setup of this site and eddy covariance flux calculation can be found in Rannik et al. (2012) and Mammarella et al. (2016).

2.2 Observational approach

Our observational analysis, schematically depicted in Fig. 1a, aims to derive bulk canopy stomatal and non-stomatal resistances from canopy-top eddy covariance observations in order to estimate the magnitude of stomatal and non-stomatal ozone removal. We first derive the ozone canopy conductance ($g_{\mathrm{c},{\mathrm{O}}_{\mathrm{3}}}$) from the observed ozone dry deposition velocity (V_d(O₃)), measurement-inferred aerodynamic resistance r_a, and bulk canopy quasi-laminar layer resistance r_b (see the Supplement).

We use the inverted Penman–Monteith equation to derive bulk canopy stomatal conductance (g_s) from canopy-top eddy covariance observations of the latent heat flux complemented with other observed variables as follows (Monteith, 1965; Knauer et al., 2018):

$$\begin{array}{}\text{(1)}& g_{\text{s}}={\displaystyle \frac{\mathit{\lambda }E{g}_{\text{a}}\mathit{\gamma }}{\mathrm{\Delta }(R_{\text{n}}-G)+\mathit{\rho }{c}_{\text{p}}{g}_{\text{a}}\text{VPD}-\mathit{\lambda }E(\mathrm{\Delta }+\mathit{\gamma })}},\end{array}$$

where g_a is the aerodynamic conductance to water vapor (see Sect. S1 in the Supplement), λE is the latent heat flux, γ is the psychrometric constant, which relates the water vapor partial pressure to air temperature, Δ is the slope of the saturation vapor pressure curve, R_n is net radiation, G is the ground heat flux, ρ is the air density, c_p is the specific heat of air, and VPD is the vapor pressure deficit. Note that all components of Eq. (1) are observed or derived from observations. g_s refers to stomatal conductance to H₂O. When we refer to the stomatal conductance for ozone, we scale g_s for the diffusivity (D) ratio of ozone and water vapor as follows: $g_{\mathrm{s},{\mathrm{O}}_{\mathrm{3}}}=\frac{D_{{\mathrm{O}}_{\mathrm{3}}}}{D_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}{g}_{\mathrm{s},{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}=\mathrm{0.61}{g}_{\mathrm{s},{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$. Non-stomatal conductance is derived as the residual of the bulk canopy conductance and the canopy stomatal conductance, assuming that stomatal and (bulk) non-stomatal uptake are two parallel pathways (see Fig. 1a).

Figure 1 Schematic display representing the biophysical controls on surface ozone removal in plant canopies in the three different approaches in this study, together with their input and output variables. The combination of uptake resistances (shown as black rectangles) inside the dashed gray rectangle yields the bulk canopy resistance (r_c). In- and output variables of the mechanisms are shown in blue and green, respectively. Orange rectangles in panel (c) display the derivation of photosynthesis parameters required in MLC-CHEM; this procedure is described in more detail in Appendix A. The shown resistances are the stomatal resistance (r_s), bulk canopy non-stomatal resistance (r_ns), resistance to cuticular uptake (r_cut), the resistance to in-canopy transport (r_a,inc), resistance to soil uptake (r_soil), and resistance to in-canopy transport in the upper canopy layer (r_a,uc), lower canopy layer (r_a,lc), and to the soil (r_a,soil).

2.3 Ozone uptake parameterizations

2.3.1 The big leaf approach

The parameterization of gaseous dry deposition in many atmospheric chemistry models is based on the resistance in the series framework introduced by Wesely (1989), hereafter referred to as W89. The discussion below considers the implementation of the big leaf dry deposition approach in the coupled meteorology–chemistry model WRF-Chem (Grell et al., 2005). Other big leaf parameterizations are available with improved treatments of the stomatal (e.g., Emberson et al., 2000; Val Martin et al., 2014; Lin et al., 2019) and non-stomatal uptake (e.g., Zhang et al., 2003). However, the common use of Wesely's (1989) parameterization in state-of-science 3D atmospheric chemistry and transport models (see, e.g., Galmarini et al., 2021) motivates the choice for this scheme in our experimental setup.

Figure 1b depicts the resistance framework. Note that this dry deposition representation is zero-dimensional, i.e., no explicit in-canopy ozone mixing ratios are calculated. The aerodynamic resistance (r_a) is calculated following Monin–Obukhov similarity theory, and the quasi-laminar layer resistance (r_b) is estimated following Hicks et al. (1987). Stomatal resistance is calculated as follows (Wesely, 1989; Erisman et al., 1994):

$$\begin{array}{}\text{(2)}& r_{\text{s}}=r_{\text{i}}\left(\mathrm{1}+{\left({\displaystyle \frac{\mathrm{200}}{R_{\text{n}}+\mathrm{0.1}}}\right)}^{\mathrm{2}}\right)\left({\displaystyle \frac{\mathrm{400}}{T_{\text{s}}(\mathrm{40}-T_{\text{s}})}}\right),\end{array}$$

where r_i (internal resistance) is a season- and land-use-dependent scaling factor, R_n is net radiation, and T_s is the surface temperature. r_s is corrected for the diffusivity difference between H₂O and ozone, as explained in Sect. 2.2. In this formulation, the resistance to stomatal uptake is lowest during high-radiation conditions and for an optimum temperature of 20 ^∘C, reflecting that stomatal aperture follows a diurnal cycle with a peak around midday. Note that this parameterization does not explicitly account for stomatal closure due to a vapor pressure deficit or soil moisture stress. We use the non-stomatal resistances, following Wesely (1989), which are all constant, except for the resistance to transport to the lower canopy that depends inversely on net radiation. For the soil uptake resistance, we use site-inferred values of 300 s m⁻¹ for Ispra (Fumagalli et al., 2016) and 400 s m⁻¹ for Hyytiälä (Zhou et al., 2017).

2.3.2 The Multi-Layer Canopy–CHemistry Exchange Model (MLC-CHEM)

We also apply the Multi-Layer Canopy–CHemistry Exchange Model (MLC-CHEM) to evaluate simulated long-term canopy-scale ozone deposition at the two sites. This one-dimensional model explicitly simulates the canopy exchange and vertical profiles of ozone concentrations as a function of radiation, turbulent mixing, chemistry (using the Carbon Bond Mechanism, version 4 – CBM-IV), biogenic emissions (following the Model for Emissions of Gases and Aerosols from Nature (MEGAN); Guenther et al., 2006, 2012), soil NO emissions (Yienger and Levy, 1995), and (non-)stomatal uptake and their vertical gradients in the canopy. MLC-CHEM has been applied coupled to single-column and global chemistry–climate modeling studies (Ganzeveld et al., 2002, 2010), as well as in an offline setup for the interpretation of site-scale measurements (e.g., Yanez-Serrano et al., 2018).

In our setup, the model consists of three layers representing the understory and the crown layer, as well as one layer aloft representing a bulk surface layer. In-canopy exchange is represented by two canopy layers whose depth depends on the canopy height (h_c), each with a layer thickness of 0.5 h_c. This two-canopy layer setup allows the simulation of in-canopy concentration and flux profiles using a computationally efficient analytical solution, allowing for the coupling of MLC-CHEM to single-column and global chemistry–climate modeling studies (Ganzeveld et al., 2002, 2010). Given the large gradients in radiation in the canopy, vertical profiles of radiation and radiation-dependent processes (photolysis and biogenic emissions) are calculated considering four canopy layers. The four-layer radiation profiles and biogenic emission rates are subsequently averaged over the two canopy layers for the exchange simulation. The model simulation time step is 30 min, but for processes requiring a higher temporal resolution, a sub-time-step temporal resolution is applied, which depends on the removal rate (Ganzeveld et al., 2002).

Micro-meteorological variables are provided as input to the model, and ozone concentrations in the upper layer are nudged to the observed above-canopy ozone concentrations to represent entrainment and advection. We use a weighting factor of 0.5, which implies that we force simulated above-canopy ozone mixing ratios to observed mixing ratios with a timescale of ±2 h, based on the applied temporal resolution of 0.5 h. The specific procedure to incorporate observations in our model setup is described in Sect. 2.4.

In-canopy aerodynamic resistance (r_a) is calculated as a function of canopy height, LAI, and u_*. Leaf-level stomatal conductance is calculated using the assimilation–stomatal conductance model, A-g_s, as follows (Ronda et al., 2001):

$$\begin{array}{}\text{(3)}& g_{\text{s,c,leaf}}=g_{\text{min,c}}+{\displaystyle \frac{a_{\mathrm{1}}{A}_{\text{g}}}{\left(\right[{\mathrm{CO}}_{\mathrm{2}}]-\mathrm{\Gamma })\left(\mathrm{1}+\frac{D_{\text{s}}(a_{\mathrm{1}}-\mathrm{1})}{D_{\mathrm{0}}}\right)}},\end{array}$$

where g_min,c (cuticular conductance), the constant a₁, and Γ (the CO₂ compensation point) depend on the vegetation type. A_g is gross assimilation, calculated as a function of photosynthetically active radiation (PAR), skin temperature, the internal CO₂ concentration, and the soil water content (SWC). We refer the reader to Appendix A in Ronda et al. (2001) for more details on the calculation of A_g. D_s is the vapor pressure deficit (VPD) leaf level, and D₀ is the VPD at which stomata close. g_s,c,leaf is calculated at the leaf level and subsequently integrated to the specific layer as a function of layer-specific LAI and PAR (Ronda et al., 2001). This stomatal conductance representation accounts for observed increases in g_s for an increase in CO₂ assimilation (which responds to radiation), whereas g_s decreases as the external CO₂ concentration increases (a lower CO₂ uptake rate is needed to maintain the supply of CO₂ to the photosynthesis mechanism). g_s also decreases as the vapor pressure deficit increases in order to minimize plant water loss through transpiration. This is a more mechanistic description of stomatal conductance compared to the big leaf approach (Eq. 2), where g_s is parameterized as a function of radiation and temperature.

The A-g_s model has several degrees of freedom in determining the parameter settings. In order to derive physically appropriate settings, we tested the sensitivity of the MLC-CHEM-simulated canopy stomatal conductance (g_s) and the canopy CO₂ flux to A-g_s parameter settings by comparing with observation-inferred g_s (using Eq. 1) and canopy-top $F_{{\mathrm{CO}}_{\mathrm{2}}}$ observations (see Fig. 1c). This procedure is described in Appendix A, and the final, optimized A-g_s parameters are shown in Table A1. With this approach, we effectively implement a realistic, observation-constrained representation canopy-top CO₂ flux and g_s in MLC-CHEM.

Non-stomatal removal in MLC-CHEM is represented using uptake resistances taken from Wesely (1989), Ganzeveld and Lelieveld (1995), and Ganzeveld et al. (1998). Analogous to W89, we adapt MLC-CHEM's default soil uptake resistance to site-inferred values of 300 s m⁻¹ for Ispra (Fumagalli et al., 2016) and 400 s m⁻¹ for Hyytiälä (Zhou et al., 2017). Experimental evidence suggests increased deposition to dew-wetted leaves (Zhang et al., 2002; Altimir et al., 2006). MLC-CHEM accounts for this by using two distinct uptake resistances for the deposition to leaf cuticles and uptake by water films on leaves of 10⁵ and 2000 s m⁻¹, respectively (Ganzeveld and Lelieveld, 1995). Canopy wetness is represented by inferring the fraction of wet vegetation (f_wet) as a function of relative humidity (RH) as follows (Lammel, 1999):

$$\begin{array}{}\text{(4)}& f_{\text{wet}}=\left\{\begin{array}{ll}\mathrm{1}& \text{RH}\ge \mathrm{0.9}\\ \frac{\text{RH}-\mathrm{0.55}}{\mathrm{0.35}}& \mathrm{0.55}\le \text{RH}<\mathrm{0.9}\\ \mathrm{0}& \text{RH}<\mathrm{0.55}.\end{array}\right.\end{array}$$

2.4 Experimental setup

We apply the W89 big leaf parameterization and the multi-layer ozone atmosphere–biosphere exchange parameterization to simulate total canopy ozone removal and its partitioning into stomatal and non-stomatal removal, at two locations with contrasting climate and pollution regimes, for a total of 12 site years. These simulations are compared against observation-inferred g_s and g_ns. We restrict this analysis to daytime values (08:00–20:00 LT) during April–September, which approximately coincides with the growing season. The observational approach is known to be biased under high canopy wetness conditions due to dew formation or precipitation, and various approaches to correct for this have been reported in the literature (e.g., Rannik et al., 2012; Launiainen et al., 2013; Clifton et al., 2017, 2019). We, therefore, only include data with RH < 90 % and when the accumulated precipitation in the preceding 12 h is less than 0.1 mm. This set of assumptions compromises between data quality and retention of data points.

Figure 2 Time series of April–September monthly average daytime (08:00–20:00 LT) ozone dry deposition velocity for Ispra, for W89 (blue), MLC-CHEM (green), and observations (black). Solid lines and points show monthly daytime medians for simulations and observations, respectively, and shaded areas and whiskers display the interquartile range.

3 Results

3.1 Temporal variability in ozone dry deposition velocity

3.1.1 Monthly and interannual variability

The observed ozone uptake at Ispra is generally highest in June–August, with little interannual variability (Fig. 2. W89 underestimates the observed dry deposition velocity (V_d(O₃)) by ±0.1 cm s⁻¹, while MLC-CHEM reproduces the observed magnitude of V_d(O₃) within 7 % in May–September. On the basis of the statistical model performance metrics in Table 1, there is no parameterization that consistently outperforms the other on monthly timescales. MLC-CHEM systematically overestimates ozone deposition in April. To evaluate this bias further, we performed MLC-CHEM simulations with a deactivated sink to wet leaves, motivated by the considerable uncertainty in this ozone removal pathway (Clifton et al., 2020a). This simulation resulted in the strongest decrease in V_d(O₃) in the relatively humid months of April (Fig. S2), ranging from 0.15 cm s⁻¹ in April 2013 to 0.05 cm s⁻¹ in April 2015. This modification results in an improved representation of seasonality in V_d(O₃), suggesting seasonal variation in the ozone sink to wet leaves that might not be properly captured by the RH-dependent parameterization of wet leaf uptake (Eq. 4).

The observed V_d(O₃) at Hyytiälä is generally lower compared to Ispra, reflecting a lower leaf area and, thus, less stomatal uptake at the Finnish site. W89 and MLC-CHEM both capture the observed magnitude of V_d(O₃) to within the interquartile range of observations (±0.2 cm s⁻¹) in most years, although V_d(O₃) in W89 peaks 1 month early compared to the observations. MLC-CHEM reproduces the seasonal cycle in V_d(O₃) with a Pearson (temporal) correlation coefficient between simulations and observations, which is markedly higher compared to the W89 approach (r²=0.59 for MLC-CHEM; r²=0.11 for W89; Table 1). These results suggest that MLC-CHEM better reproduces stomatal and non-stomatal removal processes, and we will investigate this further below.

The interannual variability in the ozone dry deposition velocity for Hyytiälä is 0.17 cm s⁻¹ and is slightly underestimated in both simulations (0.10–0.11 cm s⁻¹; not shown). We therefore calculated the contributions from stomatal conductance and non-stomatal conductance to the overall deposition velocity, as described in Sect. 2.2. Interannual variability in stomatal conductance is overestimated slightly by W89 and MLC-CHEM, compared to the observation-derived g_s estimates, by 0.02 and 0.05 cm s⁻¹, respectively. Interannual variability in non-stomatal conductance is strongly underestimated in both simulations (0.04–0.07 cm s⁻¹) compared to the observed interannual variability in non-stomatal conductance (0.19 cm s⁻¹). The missing interannual variability in the non-stomatal deposition pathway may be due the chemical, wet leaf, and soil uptake pathways.

Figure 3 As in Fig. 2 but for Hyytiälä.

Table 1 Performance statistics for the monthly averaged simulations of V_d(O₃) with W89 and MLC-CHEM (MLC). The unit is centimeters per second (cm s⁻¹)for mean bias error (MBE), root mean square error (RMSE), and the intercept and unitless for the other metrics. Shown are several conventionally applied performance metrics (MBE, RMSE, slope (s) and intercept (i) of a linear regression fit of simulations against observations and r² from the ordinary least squares regression) and the index of agreement (d).

3.1.2 Diurnal cycles

The observed diurnal cycle of V_d(O₃) at Ispra (Fig. 4a) is characterized by an asymmetrical pattern, with a steep morning increase that plateaus around 0.8 cm s⁻¹ and a decrease in the afternoon that reflects stomatal closure and reduced non-stomatal uptake. W89 underestimates the observed median daytime V_d(O₃) values by ±0.1 cm s⁻¹ (20 %), while MLC-CHEM reproduces the observations within 10 %. The onset of the W89-simulated daytime V_d(O₃) peak shows a 1 h time lag, with an underestimation of around −0.3 cm s⁻¹ (52 %) in the morning (06:00–10:00 LT) and an overestimation of 0.1 cm s⁻¹ (13 %) in the afternoon (12:00–16:00 LT). The contribution of the stomatal and non-stomatal removal in this model–observation mismatch will be discussed in Sect.3.2. MLC-CHEM reproduces the diurnal course of V_d(O₃) within 0.1 cm s⁻¹ throughout the day.

The observed V_d(O₃) diurnal cycle at Hyytiälä (Fig. 4b) increases earlier during the day, compared to Ispra, and decreases later, due to the extended day length during the growing season at the Finnish site. V_d(O₃) peaks at 0.5 cm s⁻¹ between 09:00–12:00 LT and decreases in the early afternoon due to decreasing (non-)stomatal sink ozone removal. W89 overestimates the magnitude of V_d(O₃) by 0.1 cm s⁻¹ (22 %) in the afternoon (12:00–16:00 LT) and underestimates ozone uptake in the morning (03:00–10:00 LT) and evening (after 19:00 LT). Apart from a morning overestimation by up to 0.1 cm s⁻¹, MLC-CHEM reproduces the diurnal evolution of V_d(O₃) well, which is apparently due to a more realistic representation of stomatal and non-stomatal removal processes.

Figure 4 Diurnal cycles of April–September ozone dry deposition velocity at Ispra (a) and Hyytiälä (b) derived from observations (black) and simulations with the W89 parameterization (blue) and MLC-CHEM (green). Lines and points show median values, and shaded areas and whiskers display the interquartile range.

3.2 Diurnal variability in stomatal and non-stomatal uptake

3.2.1 Ispra

Next, we analyze the stomatal and non-stomatal components of ozone deposition to further understand the model–observation agreement on diurnal timescales. Figure 5 shows growing season median diurnal cycles of bulk canopy conductance (g_c), canopy stomatal conductance (g_s), and non-stomatal conductance (g_ns) for Ispra in W89 and MLC-CHEM simulation and observational estimates. At Ispra, the observation-derived daytime median ozone canopy conductance is 0.87 cm s⁻¹ (Fig. 5a). The inferred daytime median stomatal conductance is as small as 0.26 cm s⁻¹ (gray points in Fig. 5b), corresponding to a daytime stomatal uptake fraction of 35 % (Fig. 5d). However, we found a substantial gap (of 56 %) in the energy balance closure (Q_gap), defined as the difference between net incoming radiation (R_n) and the surface energy balance components (Foken, 2008). This indicates underestimations in observed sensible and latent energy fluxes (H and LE), which affects our observation-derived stomatal conductance. The energy balance closure issues remain after filtering the observations based on quality flags and u_* thresholds (Fig. S1).

To resolve these energy balance closure issues, we applied a correction method that partitions Q_gap to H and LE via the evaporative fraction (EF = $LE/(H+LE)$; Twine et al., 2000; Renner et al., 2019). This correction increases LE and H by 156 and 25 W m⁻², respectively, corresponding to an evaporative fraction of 0.86. With these corrected surface energy balance components, we derive a substantially larger daytime median stomatal conductance to ozone of 0.49 cm s⁻¹ (black points in Fig. 5b, d), which is an increase of nearly 90 % with respect to the original observation-derived estimate. The Q_gap correction also leads to a better model–observation agreement for g_s. Ozone fluxes are also affected by the surface energy balance closure gap; additional data filtering based on u_* thresholds leads to increases in observed ozone fluxes, which exceeds 50 % in the morning and evening when absolute fluxes are low, but the effect is smaller (< 15 %) during midday.

The observed diurnal cycle in canopy conductance at Ispra is better captured by MLC-CHEM compared to W89 (Fig. 5a). MLC-CHEM also better captures decreases in g_c observed in the afternoon. MLC-CHEM and W89 simulate a daytime median ozone stomatal conductance of 0.43 and 0.51 cm s⁻¹, respectively, and, thus, agree better with the Q_gap-corrected stomatal conductance estimate derived from observations (Fig. 5b). The observation-derived stomatal fraction during 08:00–20:00 LT (0.62) is overestimated by W89 (0.72) and underestimated by MLC-CHEM (0.52). The observed stomatal uptake fraction increases throughout the day, from ±0.4 at 08:00 LT to ±0.8 at 18:00 LT, and this diurnal course is better reproduced by MLC-CHEM than by W89.

Observation-derived non-stomatal conductance peaks in the morning, levels off at ±0.8 cm s⁻¹ (Fig. 5c; gray points), and decreases in the afternoon before reaching a nighttime value of 0.1 cm s⁻¹. The stomatal conductance increase following Q_gap correction leads to a reduction in the daytime average inferred non-stomatal conductance from 0.57 to 0.35 cm s⁻¹. This correction does, however, not affect the shape of the diurnal cycle in g_ns, which is characterized by a sharp increase in the morning and a more gradual reduction in the afternoon. Daytime non-stomatal conductance is strongly underestimated by W89 and shows little diurnal variability since most in-canopy resistances are constant and, apparently, too high. MLC-CHEM reproduces the observed diurnal evolution in non-stomatal conductance more accurately than W89 (Fig. 5c), apparently due to its representation of diurnal variability in processes involved in non-stomatal removal, wet leaf uptake, and in-canopy turbulence. The contributions of different removal processes to total non-stomatal uptake will be discussed in Sect. 3.3.

Figure 5 April–September median diurnal cycles of ozone bulk canopy conductance (a), canopy stomatal conductance (b), bulk non-stomatal conductance, and (c) the stomatal fraction of total ozone removal ($g_{\text{s}}{g}_{\text{c}}^{-\mathrm{1}}$; d) for Ispra. Observed medians and interquartile ranges after Q_gap correction (OBS^*; see text) are shown as black points and whiskers (the values prior to Q_gap correction, denoted as OBS, are shown in gray). The median and interquartile range of W89 and MLC-CHEM are shown in blue and green, respectively. The shaded area in panel (d) highlights the nighttime period (defined as 08:00–20:00 LT) during which the stomatal flux is calculated.

3.2.2 Hyytiälä

At Hyytiälä, the observation-derived daytime median g_c is 0.53 cm s⁻¹ (Fig. 6a), which is lower compared to Ispra due to lower non-stomatal ozone removal. W89 overestimates canopy conductance by up to 0.2 cm s⁻¹ in the afternoon, while morning and evening g_c are underestimated. Similar to Ispra, MLC-CHEM captures the diurnal evolution in g_c better than W89, with a peak around 09:00 LT, as in the observations, but overestimates morning canopy conductance by 0.1 cm s⁻¹. We did not correct for surface energy balance closure gaps for the Hyytiälä observations, since this gap was considerably smaller (±20 % of r_n, without a distinct diurnal cycle) and in closer agreement to the literature-reported values for tall vegetation (Foken, 2008).

Observed stomatal conductance peaks at ±0.5 cm s⁻¹ at 10:00 LT, followed by a decrease in the afternoon (Fig. 6b). W89 underestimates g_s in the morning (05:00–10:00 LT), and overestimates afternoon values by 20 %–25 %. MLC-CHEM overestimates morning stomatal conductance, but follows the observed diurnal cycle well throughout the rest of the day. The observed stomatal ozone uptake fraction is relatively constant at 0.8 (Fig. 6), comparable to the upper range of stomatal uptake fraction estimates by Rannik et al. (2012). The stomatal fraction is well reproduced by both parameterizations, although for W89 this seems a coincidence given the misrepresented diurnal cycle in g_c and g_s.

Observation-derived non-stomatal conductance at Hyytiälä (Fig. 6c) is relatively constant at 0.1 cm s⁻¹, except for a morning peak of 0.2 cm s⁻¹ around 08:00 LT that likely reflects wet leaf ozone uptake (Altimir et al., 2006; Rannik et al., 2012). W89 reproduces the observed daytime magnitude of g_ns but cannot reproduce its morning peak. MLC-CHEM overestimates the nighttime non-stomatal ozone sink, in line with a study by Zhou et al. (2017) based on a 1-month time series of ozone flux observations (August 2010), indicating that observed nighttime ozone deposition appears to reflect smaller nocturnal soil uptake efficiency than assumed. Except for an overestimation in the morning, MLC-CHEM captures the observation-inferred magnitude of non-stomatal ozone deposition well during daytime.

Figure 6 As in Fig. 5 but for Hyytiälä.

3.3 Dependence of non-stomatal deposition on driving variables

Non-stomatal ozone uptake, and its dependence on micro-meteorological and other environmental drivers, is incompletely understood. Previous studies employed statistical or process-oriented modeling (Rannik et al., 2012; Fares et al., 2014; El-Madany et al., 2017; Clifton et al., 2019) to determine the contribution of driving variables to this ozone sink. In this section, we study observed and simulated relationships between the non-stomatal ozone removal fraction ($g_{\text{ns}}{g}_{\text{c}}^{-\mathrm{1}}$) and two variables (air temperature, T_a, and VPD) that we hypothesize to contribute to temporal variability in non-stomatal ozone removal. This section focuses on non-stomatal ozone removal at Ispra, since Rannik et al. (2012) previously characterized the non-stomatal ozone sink for Hyytiälä, and we compare our findings for Ispra to their results at the end of this section.

We first determine how W89 and MLC-CHEM can reproduce the observed relationship between non-stomatal ozone removal and T_a and VPD. We focus on the average daytime response (08:00–20:00 LT) and, subsequently, on three different periods in the diurnal cycle (06:00–20:00, 10:00–14:00, and 14:00–18:00 LT). In this manner, we can disentangle processes affecting (non-)stomatal uptake that act during different periods of the diurnal cycle (wet leaf uptake in the morning, optimal stomatal functioning during midday, and suppressed stomatal conductance during the afternoon).

The temperature response of the relative contribution of non-stomatal removal to total ozone deposition (expressed by the non-stomatal fraction, $g_{\text{ns}}{g}_{\text{c}}^{-\mathrm{1}}$) during different periods of the diurnal cycle is shown in Fig. 7a–d. Non-stomatal uptake decreases with temperature during the day (Fig. 7a). This decrease is largely driven by the morning temperature sensitivity of $g_{\text{ns}}{g}_{\text{c}}^{-\mathrm{1}}$, which shows less sensitivity to temperature later during the day (Fig. 7b–d). W89 underestimates the observed temperature dependence of the non-stomatal fraction throughout the day by ±0.2, although the morning non-stomatal fraction is higher for the lowest temperature bin (10–15 ^∘C). MLC-CHEM reproduces the daytime response well, as it is characterized by elevated morning non-stomatal uptake under low-temperature conditions. For most temperature bins, W89 strongly underestimates the observed variability in the non-stomatal fraction. The observed variability is also underestimated by MLC-CHEM, although to a smaller extent, and apparently indicates still missing or misrepresented deposition processes.

The observation-derived non-stomatal fraction increases with VPD during daytime (Fig. 7e–h), indicating that non-stomatal ozone removal decreases under dry conditions. This result contradicts an anticipated increase in the contribution by non-stomatal removal to overall canopy removal due to a VPD-induced decrease in stomatal uptake. However, the observed non-stomatal uptake also decreases in the afternoon (Fig. 5c) and, therefore, does not compensate for the decreasing stomatal sink with VPD. The non-stomatal fraction displays the strongest VPD sensitivity in the morning, which mainly reflects simulated wet leaf uptake under humid (i.e., low VPD) conditions. Non-stomatal removal in W89 is insensitive to VPD, and this parameterization particularly underestimates the non-stomatal fraction under humid conditions (Fig. 7a–b). MLC-CHEM reproduces the daytime slope between VPD and the non-stomatal fraction well.

We then perform a number of sensitivity experiments with MLC-CHEM with deactivated non-stomatal sinks to identify the role of each sink in explaining temporal variability in non-stomatal ozone removal and its dependence on T_a and VPD. In these experiments, we exclude the contribution by wet leaf uptake, soil uptake, and in-canopy chemical removal, as well as an experiment with strongly enhanced turbulent exchange between the crown layer and the understory. In the Supplement, Sect. S2 and Fig. S6 display the results from the sensitivity analysis of non-stomatal removal at Ispra. We list the main outcomes of this section and the MLC-CHEM sensitivity analysis as follows:

• Soil deposition accounts for almost 40 % of non-stomatal removal under high-temperature conditions, reflecting a simulated increase in in-canopy turbulent transport with air temperature.

• Non-stomatal uptake is elevated under cold and humid conditions in the morning. This is consistent with MLC-CHEM-simulated wet leaf uptake, which accounts for over 20 % of the morning non-stomatal removal fraction.

• Enhanced turbulent transport from the crown layer to the understory reduces the non-stomatal uptake fraction in MLC-CHEM, as it leads to enhanced stomatal uptake in the understory.

• Chemical removal plays a minor role in the total canopy ozone sink at Ispra.

In their multivariate analysis of environmental drivers of non-stomatal ozone removal at Hyytiälä, Rannik et al. (2012) derived that air temperature and VPD are significantly associated with variations in non-stomatal ozone removal, similar to our findings for Ispra. However, Rannik et al. (2012) also found an explanatory role for monoterpene concentrations at Hyytiälä, while our results suggest a minor role of chemical removal at Ispra.

Figure 7 Non-stomatal ozone removal fraction $g_{\text{ns}}{g}_{\text{c}}^{-\mathrm{1}}$ binned by air temperature (a–d) and vapor pressure deficit (e–h) during June–September 2013–2015 at Ispra. Black dots and whiskers show $g_{\text{ns}}{g}_{\text{c}}^{-\mathrm{1}}$ from observations and simulations by W89 (blue points and whiskers) and MLC-CHEM (green points and whiskers). Dots and whiskers display the median and interquartile range per bin, respectively, and the number of observations in the bin is displayed at the bottom of the panels. Each column corresponds to a different time period in the diurnal cycle, namely all day (08:00–20:00 LT; a, e), morning (06:00–10:00 LT; b, f), midday (10:00–14:00 LT; c, g), and afternoon (14:00–18:00 LT; d, h). Dots and whiskers are only shown if the number of samples in the bin exceeds 10.

3.4 Cumulative uptake of ozone (CUO)

In the previous sections, we have shown that the seasonal evolution of ozone deposition in W89 and MLC-CHEM is relatively similar. However, there are consistent differences in daytime ozone stomatal and non-stomatal sinks between the deposition representations. In this section, we evaluate the implications of these differences in the representation of (non-)stomatal removal for determining the cumulative stomatal uptake of ozone (CUO) over the growing season, which is often used for ozone impact assessments (Musselman et al., 2006; Mills et al., 2011). We here use the term CUO_st to refer to cumulative stomatal uptake and to distinguish this from cumulative non-stomatal ozone removal (CUO_ns). In Fig. 8, we compare growing-season-integrated stomatal and non-stomatal ozone fluxes from W89 and MLC-CHEM to observation-derived estimates of total seasonal ozone uptake for both sites (Fig. 8a, c). Our observation-based derivation of stomatal conductance requires dry conditions (RH < 90 %, and no precipitation in the preceding 12 h) to avoid overestimations in the observation-inferred stomatal conductance which lead to overestimations in CUO. However, the application of these data selection criteria also leads to a reduction in data points that hinders the calculation of CUO_st based on observations. In order to derive a first-order CUO_st estimate, we divide the cumulative stomatal uptake inferred from valid observations by the fraction of valid observations. This method serves mainly to perform a site-to-site comparison of inferred CUO_st. Inferred CUO_st at Ispra varies between 61 and 72 mmol m⁻² (Fig. 8b). The inferred CUO_st in 2014 was lower compared to 2013 and 2015 due to comparatively low ozone mixing ratios, while stomatal conductance displayed less year-to-year variability. At Hyytiälä, inferred CUO_st varies between 39 and 41 mmol m⁻² (Fig. 8d), where the lower value in 2005 (29 mmol m⁻²) is caused by missing data during June–August, when stomatal ozone uptake peaks. The higher inferred CUO_st values at Ispra compared to Hyytiälä reflect both higher stomatal conductance and ozone mixing ratios at the Italian site.

The differences between W89- and MLC-CHEM-simulated conductances are also manifested in the simulated growing season cumulative (stomatal) uptake (Fig. 8a–c). The cumulative total ozone flux for Ispra is underestimated by W89 (−10 %), while this parameterization overestimates cumulative stomatal uptake by 14 %–22 %. MLC-CHEM accurately reproduces observation-derived CUO_st (within 7 %) but overestimates the cumulative total flux by 15 % (Fig. 8a). Therefore, the model–observation agreement of the two parameterizations for the simulated cumulative total ozone removal largely reflects non-stomatal uptake differences, which deviates from observation-inferred values by −64 % and 51 %, respectively. At Hyytiälä, the observed cumulative total ozone flux is 15.1 mmol m⁻² and is overestimated by 34 % by W89, reflecting overestimated stomatal uptake (Fig. 8c). The observation-derived CUO_st is 12.6 mmol m⁻², compared to 15.8 mmol m⁻² in W89 (+28 %) and 12.3 mmol m⁻² in MLC-CHEM (−2.4 %). We conclude that the better representation of canopy stomatal conductance in MLC-CHEM compared to W89, particularly during the afternoon peak in ozone mixing ratios, may lead to a substantially reduced bias in the simulated growing-season-integrated (stomatal) ozone flux.

Figure 8 (a, c) Growing-season-integrated daytime (08:00–20:00 LT) and stomatal (CUO_st; dark colors) and non-stomatal (CUO_ns; light colors) ozone fluxes for the different years in the study period for Ispra (a) and Hyytiälä (c). Results from the W89 parameterization are shown in blue, MLC-CHEM in green, and observations in gray. Only the data points with valid observation-inferred stomatal conductance estimates are selected for this comparison; the fraction of valid data points per growing season that remains is shown at the bottom of panels (a, c). (b, d) Inferred cumulative stomatal ozone uptake (CUO) estimate at both sites (hatched bars) after dividing the season-integrated daytime stomatal ozone flux (dark gray bars) by the fraction of valid data points. Note the different y axis ranges across the four panels.

4 Discussion

This study evaluates the potential added value of a multi-layer representation of vegetation canopies with respect to a commonly applied big leaf approach (W89;  Wesely, 1989) for simulating ozone deposition and ozone impact metrics for forest canopies. We focus on short- to long-term temporal variability in V_d(O₃) and its partitioning into stomatal and non-stomatal components, as well as the simulation of ozone impact metrics. We find that both parameterizations reasonably reproduce the observed seasonal cycle in V_d(O₃), which is in agreement with previous chemistry transport model evaluations (e.g., Hardacre et al., 2015). Despite their comparable performance on seasonal timescales, the parameterizations deviate in their simulation of the diurnal cycle; the W89 parameterization particularly underestimates morning ozone removal by 52 % (Ispra) and 37 % (Hyytiälä), due to a combination of underestimated stomatal removal and a missing non-stomatal sink, which is likely wet leaf uptake. In the afternoon, W89 deviates less from observations at both sites (−13 % at Ispra; +22 % at Hyytiälä). Consequently, cumulative stomatal ozone uptake is overestimated by, on average, 18 % (Ispra) and 28 % (Hyytiälä) in W89 simulations, while cumulative total ozone removal deviates by −10 % (Ispra) and 20 % (Hyytiälä). Ozone mixing ratios typically peak in the afternoon and, thus, occur simultaneously with stomatal conductance misrepresentations, which may lead to simulated ozone fluxes overestimates using this mechanism. The multi-layer mechanism, constrained with latent energy and net ecosystem exchange (NEE) observations to optimally represent stomatal exchange, displays a better agreement with the observed ozone deposition velocity (within 10 %) and inferred cumulative stomatal and total uptake (within 15 % and 9 % for Ispra and Hyytiälä, respectively). Therefore, an accurate representation of diurnal variability in ozone uptake partitioned to stomatal and non-stomatal sinks is essential for reproducing cumulative (stomatal) ozone uptake at the land surface.

We applied a big leaf parameterization that is commonly used in (regional) atmospheric chemistry models, for example in WRF-Chem (Grell et al., 2005; Galmarini et al., 2021). Big leaf parameterizations advantageously depend on a limited number of routinely available meteorological variables and a simplified description of land use characteristics and can be readily applied at any location without location-specific parameter derivations (Clifton et al., 2020a). However, the empirical nature of these schemes leads to an oversimplification of in-canopy physical and chemical processes that affect atmosphere–biosphere exchange of ozone, e.g., by not accounting for stomatal closure based on the vapor pressure deficit (VPD) and soil moisture or in-canopy chemical reactions. There are big leaf versions available with a more process-based description of ozone deposition processes, particularly for stomatal conductance (e.g., Lin et al., 2019; Clifton et al., 2020c; Emberson et al., 2001; Büker et al., 2012) and non-stomatal ozone removal (Zhang et al., 2003).

To further explore the effect of model assumptions in big leaf parameterizations, we performed a comparison between W89 and another commonly used big leaf dry deposition scheme by Zhang et al. (2003, referred to as Z03) in Appendix B. This parameterization includes a separate treatment of sunlit versus shaded leaves and the explicit treatment of water stress in the stomatal conductance calculation and includes variations in non-stomatal resistances as a function of LAI and u_*. We find that both parameterizations overestimate afternoon stomatal conductance compared to observations, while Z03 better reproduces morning g_s (Fig. B1). The differences between these parameterizations are, therefore, largely driven by differences in non-stomatal ozone removal (Fig. B1). The agreement with observation-inferred non-stomatal removal depends on site-specific conditions, particularly friction velocity. Our analyses highlight potential areas of improvement in process representation that can be considered in future larger-scale modeling studies to improve the simulations of ozone deposition pathways and their temporal variability. This is particularly important for season-integrated (stomatal) ozone fluxes with big leaf parameterizations.

Our results suggest that A_net-g_s parameterizations, as applied in MLC-CHEM, simulate stomatal conductance in good agreement with observation-inferred values throughout the diurnal cycle. Such models are sensitive to parameters typically derived at leaf level that display spatiotemporal variability. Further observational constraints on these parameters, e.g., from leaf-level ecophysiological measurements, improve the representation of stomatal conductance and biosphere–atmosphere exchange (Vilà-Guerau De Arellano et al., 2020), benefitting simulations of CO₂ and ozone exchange as simulated by A_net-g_s within MLC-CHEM. Determining these parameters from canopy-top observations is an underdetermined problem in a mathematical sense, which we circumvented by deriving a realistic set of model parameters based on a comparison with canopy-top observed NEE and observation-derived stomatal conductance, while remaining as close as possible to the original parameter set in Ronda et al. (2001). Choosing A_net-g_s parameters could be formalized by applying mathematical techniques such as data assimilation (Raoult et al., 2016).

MLC-CHEM can be driven by diagnostic variables available from the chemical transport model (CTM) output (or their driving meteorological models), favoring its implementation to represent atmosphere–biosphere fluxes of reactive compounds (Ganzeveld et al., 2002, 2010). In such a coupled setup, MLC-CHEM would use simulated stomatal conductance from the driving model to represent an atmosphere–biosphere exchange consistent with the model's representation of (micro-)meteorology. An implementation of A-g_s with CO₂ mixing ratios, calculated online or offline, can be tested if simulated stomatal conductance estimates are unavailable.

Our analysis did not include soil moisture as a predictor of stomatal conductance. Sensitivity simulations in MLC-CHEM with observation-constrained soil water content (SWC) at different depths resulted in strong reductions in simulated NEE and g_s during summer compared to observations, which suggests that these SWC observations are not indicative of root zone soil moisture. Nonetheless, simulations of ozone deposition and mixing ratios at various spatial scales suggest a higher predictive skill when accounting for SWC (Anav et al., 2018; Lin et al., 2019, 2020; Clifton et al., 2020c; Otu-Larbi et al., 2020). Including this stress term is especially important in the context of projected drought risk and intensity increases in future climate scenarios (Cook et al., 2018) that may aggravate ozone smog episodes due to a decreased stomatal sink (Lin et al., 2020).

Our analysis of non-stomatal ozone removal as a function of micro-meteorological drivers (air temperature and VPD) for Ispra reveals that the non-stomatal sink is elevated under low-VPD (i.e., high-RH) morning conditions, likely indicating uptake at the leaf surface in water films formed by dew (Zhang et al., 2002; Potier et al., 2015). This sink is reproduced by MLC-CHEM by applying a wet canopy fraction dependent on RH and a constant wet skin uptake resistance. Observations suggest that this non-stomatal ozone sink is less important at Hyytiälä, which could be due to a lower RH threshold for the development of wet canopy conditions in MLC-CHEM compared to previous work (Altimir et al., 2006; Zhou et al., 2017). Since wet leaf uptake may affect simulated diurnal cycles of ozone in chemistry transport models (Travis and Jacob, 2019), uptake parameterizations would benefit from better observation-based constraints on this removal process, both in terms of canopy wetness and wet leaf uptake efficiency.

Our sensitivity analysis also reveals an important role of soil deposition during the afternoon due to more active in-canopy transport. We applied a constant soil resistance to ozone uptake in our simulations, despite various environmental controls that have been identified, including air temperature, soil water content, near-surface air humidity, and soil clay content (Fares et al., 2014; Fumagalli et al., 2016; Stella et al., 2011, 2019). Our results suggest a minor importance of chemical ozone removal at the two considered sites. However, we did not investigate the role of ozone scavenging by reactive sesquiterpenes (Zhou et al., 2017; Hellén et al., 2018; Vermeuel et al., 2021) or soil-emitted nitric oxide (Finco et al., 2018). Since most (big leaf) parameterizations work with a poorly constrained resistance to transport from the canopy-top to the soil (e.g., Makar et al., 2017), the importance of the chemical and soil ozone sinks for total canopy ozone removal can be best explored with better resolved in-canopy turbulent exchange in model simulations.

We have shown that stomatal and non-stomatal sinks are not accurately reproduced using the W89 big leaf parameterization compared to observations at two forested ozone flux sites, leading to structurally biased instantaneous and growing season cumulated (stomatal) ozone flux simulations. Improved methods (e.g., the DO3SE mechanism,  Emberson et al., 2001; Büker et al., 2012) do correct for soil moisture and VPD in the stomatal conductance calculation. Overestimated stomatal ozone fluxes also likely have implications for simulated ozone mixing ratios. Many models underestimate midday ozone mixing ratios in Europe (Solazzo et al., 2012; Im et al., 2015; Visser et al., 2019), and a misrepresentation of land surface uptake may contribute to this bias. Therefore, an overestimated ozone deposition flux may also affect the simulation of concentration-based vegetation ozone impact metrics, such as AOT40, in the opposite direction compared to flux-based metrics. An improved model representation of the ozone deposition process will provide more confidence in the application of atmospheric chemistry models for surface air quality and vegetation ozone damage assessments.

To stimulate improvement of big leaf and multi-layer parameterizations, modelers may benefit from evaluations against existing long-term dry deposition observations in various ecosystems (e.g., forests and grassland) and for contrasting environmental conditions (e.g., during dry vs. wet seasons). Such an assessment is currently underway in stage 4 of the Air Quality Model Evaluation International Initiative (AQMEII4; Galmarini et al., 2021). Additionally, evaluation against in- and above-canopy ozone flux measurements (Fares et al., 2014; Finco et al., 2018) can reveal information about non-stomatal sinks in these parameterizations, such as soil deposition and in-canopy chemical removal. Last, the application of proposed parameterizations for non-stomatal ozone sinks, such as for wet leaf uptake (Potier et al., 2015) and soil uptake (Stella et al., 2019), should be tested in 3D and in single-point models of ozone deposition.

5 Conclusions

We compare ozone deposition simulations to multi-year observations at two European forested flux sites, with a focus on temporal variability, contributions from stomatal and non-stomatal sinks, and metrics for the damage incurred by ozone on vegetation. The widely used big leaf parameterization Wesely (W89; 1989) and the in-canopy process-resolving MLC-CHEM model both reproduce the seasonal cycle of daytime ozone deposition velocity reasonably well, but there are important differences in the skill of the two approaches to capture the diurnal changes in ozone deposition. Specifically, W89 consistently underestimates ozone deposition velocities in the morning (by 37 %–52 %), while the afternoon model observation is somewhat smaller (−13 %–22 %). MLC-CHEM captures the diurnal cycle much better, with relatively small biases, in the morning (−9 % at Ispra, +17 % at Hyytiälä) and has good agreement (within 10 %) in the afternoon. Accounting for stomatal closure, wet leaf removal and in-canopy turbulent transport followed by soil uptake turns out to be important for accurately simulating ozone deposition on diurnal timescales.

The structural errors in W89 are explained by a misrepresentation of the diurnal cycle in stomatal and non-stomatal conductance. Simulations with a more recent big leaf parameterization result in similar biases regarding stomatal and non-stomatal uptake. The MLC-CHEM model, constrained by local observations of diurnal CO₂ and latent energy fluxes, captures stomatal and non-stomatal ozone conductance better. As a result, W89 systematically overestimates cumulative ozone uptake by 20 %–30 % in the growing season at Ispra and Hyytiälä, whereas MLC-CHEM reproduces cumulative ozone uptake within 3 % at both sites. We conclude that MLC-CHEM, nudged with observation-inferred stomatal conductance, accurately describes non-stomatal uptake processes and vegetation ozone impact metrics.

Sensitivity tests with MLC-CHEM for Ispra point out that, in relatively cold and humid conditions, ozone deposition on wet leaves appears to explain up to 20 % of the non-stomatal ozone sink. During high-temperature conditions characterized by efficient in-canopy transport, enhanced uptake by soils accounts for up to 40 % of non-stomatal ozone deposition. The tests suggest a minor role for chemical destruction of ozone at Ispra.

Our results indicate that current model representations of stomatal and non-stomatal ozone uptake by vegetation, often based on W89, should be thoroughly evaluated. This study provides a strategy for such evaluations and shows how a more detailed, canopy-resolving model driven by ancillary measurements of CO₂ and energy fluxes can provide more realistic estimates of ozone deposition and vegetation ozone impact metrics.

Appendix A: A-g_s optimization

Prior to applying MLC-CHEM to analyze ozone fluxes at our study sites, we first paid attention to simulations of the canopy CO₂ flux ($F_{{\mathrm{CO}}_{\mathrm{2}}}$) and canopy stomatal conductance (g_s) to ensure that the photosynthesis parameterization (A-g_s) functions satisfactorily. An initial simulation with the default settings for the C₃ vegetation class resulted in a strongly overestimated $F_{{\mathrm{CO}}_{\mathrm{2}}}$ compared to observations at both sites (see Table A2). This is accompanied by strong overestimation of the canopy stomatal conductance at Ispra, while MLC-CHEM slightly underestimates stomatal conductance at Hyytiälä.

The default A-g_s settings were derived for low vegetation such as grassland and crops (Ronda et al., 2001) and are therefore not necessarily representative for forest canopies. We performed a sensitivity analysis of simulated $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and g_s to A-g_s model parameters in order to determine optimized parameter sets for our simulations. These settings are given in Table A1. We found that strongly overestimated $F_{{\mathrm{CO}}_{\mathrm{2}}}$ is largely caused by a presumed high reference mesophyll conductance (g_m,298), leading to overestimated transport of CO₂ in the plant's chloroplast. Our reductions in g_m,298 are in better correspondence with previously reported estimates of 0.8–2.0 mm s⁻¹ for different forest plant functional types (Steeneveld, 2002; Voogt et al., 2006; ECMWF, 2020). At Ispra, we additionally modified the mesophyll conductance temperature response curve, which differs between plant species (Calvet et al., 1998; von Caemmerer and Evans, 2015), to improve the amplitude of the seasonal cycle in simulated $F_{{\mathrm{CO}}_{\mathrm{2}}}$. At Hyytiälä, the maximum internal CO₂ concentration (f₀, given as a fraction of the external CO₂ concentration) was increased to improve the correspondence with observation-derived g_s.

Our observational constraints to A-g_s lead to improved simulations of g_s and $F_{{\mathrm{CO}}_{\mathrm{2}}}$ (Table A2). The parameter changes additionally affect the simulation of the ozone dry deposition velocity (V_d(O₃)), as shown in Table A3. At Ispra, the strong reduction in stomatal conductance leads to an underestimation in V_d(O₃) (MBE = −0.12 cm s⁻¹), while the other statistical metrics indicate a modest model improvement. At Hyytiälä, the growing season model overestimation is slightly reduced from 0.04 to 0.02 cm s⁻¹. Our approach results in a reduced model bias at the two study sites, particularly for $F_{{\mathrm{CO}}_{\mathrm{2}}}$, while taking care to stay as close as possible to the original parameter set.

Table A1 A-g_s parameter settings used in MLC-CHEM simulations. The first column indicates the default C₃ settings from Ronda et al. (2001), and the other two columns show the optimal settings from our analysis. A dash (–) indicates that a parameter is unchanged with respect to the default C₃ value.

Table A2 Model performance statistics of MLC-CHEM before and after A-g_s optimization for canopy stomatal conductance and CO₂ flux. Shown are several conventionally applied performance metrics (MBE, RMSE, slope (s), and intercept (i) of a linear regression fit of simulations against observations and r² from ordinary least squares regression), as well as the index of agreement (d). The units are centimeters per second (cm s⁻¹) and micromoles per square meter per second (µmol m⁻² s⁻¹), respectively, unless indicated otherwise.

Table A3 As in Table A2 but for V_d(O₃) (in centimeters per second; cm s⁻¹).

Appendix B: Comparison between two big leaf parameterizations

In order to derive more generic conclusions about big leaf parameterizations, we considered another commonly applied parameterization (Zhang et al., 2003) that has recently been extended to different gases by Wu et al. (2018). This big leaf formulation (hereafter Z03) differs compared to the Wesely (1989) parameterization (hereafter W89) in several aspects, as follows: (1) Z03 calculates stomatal conductance for sunlit and shaded leaves differently, (2) stomatal conductance is affected by VPD and soil moisture stress, and (3) non-stomatal resistances contain seasonal and diurnal variability due to dependencies on leaf area index and friction velocity (u_*). This model version was derived from Zhang and Wu (2021) with two modifications. First, we adapted the soil resistance to locally derived values of 400 s m⁻¹ (Hyytiälä) and 300 s m⁻¹ (Ispra), similar to W89 and MLC-CHEM (see Sect. 2). The implementation by Zhang and Wu (2021) relies on observed canopy wetness, which is not available for our two study sites. We therefore parameterize canopy wetness as a function of relative humidity, analogous to MLC-CHEM (Eq. 4). In this section, we compare simulations by W89 and Z03 to observations of the ozone dry deposition velocity and observation-inferred stomatal and non-stomatal conductance.

Figure B1 shows multi-year growing season median diurnal cycles of V_d(O₃), g_s and g_ns for Ispra and Hyytiälä. From this analysis, we conclude that W89 and Z03 perform similarly for Ispra compared against observed V_d(O₃) (Fig. B1a). Z03 better captures the early morning onset of V_d(O₃) for Hyytiälä than W89, but more strongly overestimates midday and afternoon V_d(O₃) compared to observations (Fig. B1b). Both parameterizations overestimate midday and afternoon g_s, while Z03 better captures the observed morning and afternoon g_s values than W89 (Fig. B1c, d). For g_ns, there is no parameterization that performs best for the two sites. Both parameterizations underestimate observation-inferred g_ns at Ispra (corrected for energy balance closure gaps; see Sect. 3.2), while W89 better captures the magnitude of observation-inferred g_ns (although Z03 better reproduces the shape of the diurnal cycle). This suggests that the g_ns dependence on u_* is less strong in the observations than is suggested in the Z03 parameterization; i.e., a sensitivity experiment with doubled u_* values for Ispra results in daytime g_ns values of 0.2–0.35 cm s⁻¹, which is an increase by a factor 2.3–2.8. Based on our findings, we conclude that the different representation of non-stomatal ozone removal drives the differences between W89 and Z03, but the magnitude of these differences depends on site-specific conditions.

Figure B1 Comparison of the dry deposition parameterizations W89 (Wesely, 1989) and Z03 (Zhang et al., 2003) against the observed dry deposition velocity (a, b) and observation-inferred stomatal conductance (c, d) and non-stomatal conductance (e, f) for Ispra and Hyytiälä (left and right panels, respectively). Lines and shaded areas (points and whiskers) show April–September median and interquartile range of the simulations (observations).