A revised dry deposition scheme for land–atmosphere exchange of trace gases in ECHAM/MESSy v2.54

Dry deposition to vegetation is a major sink of ground-level ozone and is responsible for about 20 % of the total tropospheric ozone loss. Its parameterization in atmospheric chemistry models represents a significant source of uncertainty for the global tropospheric ozone budget and might account for the mismatch with observations. The model used in this study, the Modular Earth Submodel System version 2 (MESSy2) linked to the fifth-generation European Centre Hamburg general circulation model (ECHAM5) as an atmospheric circulation model (EMAC), is no exception. Like many global models, EMAC employs a “resistance in series” scheme with the major surface deposition via plant stomata which is hardly sensitive to meteorology, depending only on solar radiation. Unlike many global models, however, EMAC uses a simplified high resistance for non-stomatal deposition which makes this pathway negligible in the model. However, several studies have shown this process to be comparable in magnitude to the stomatal uptake, especially during the night over moist surfaces. Hence, we present here a revised dry deposition in EMAC including meteorological adjustment factors for stomatal closure and an explicit cuticular pathway. These modifications for the three stomatal stress functions have been included in the newly developed MESSy VERTEX submodel, i.e. a process model describing the vertical exchange in the atmospheric boundary layer, which will be evaluated for the first time here. The scheme is limited by a small number of different surface types and generalized parameters. The MESSy submodel describing the dry deposition of trace gases and aerosols (DDEP) has been revised accordingly. The comparison of the simulation results with measurement data at four sites shows that the new scheme enables a more realistic representation of dry deposition. However, the representation is strongly limited by the local meteorology. In total, the changes increase the dry deposition velocity of ozone up to a factor of 2 globally, whereby the highest impact arises from the inclusion of cuticular uptake, especially over moist surfaces. This corresponds to a 6 % increase of global annual dry deposition loss of ozone resulting globally in a slight decrease of ground-level ozone but a regional decrease of up to 25 %. The change of ozone dry deposition is also reasoned by the altered loss of ozone precursors. Thus, the revision of the process parameterization as documented here has, among others, the potential to significantly reduce the overestimation of tropospheric ozone in global models.

1 Introduction

Ground-level ozone is a secondary air pollutant which is harmful for humans and ecosystems. Besides chemical destruction, a large fraction of it is removed by dry deposition which accounts for about 20 % of the total O₃ loss (Young et al., 2018). The process description of dry deposition considers boundary-layer meteorology (e.g. turbulence), chemical properties of the trace gases and surface types. In most global models, dry deposition of trace gases is parameterized using the “resistance in series” analogy by Wesely (1989). The largest deposition rates of ozone occur over dense vegetation (Hardacre et al., 2015) where it mainly follows two pathways: through leaf openings (stomata) and to leaf waxes (cuticle) (Fares et al., 2012). Thereby, stomatal uptake is commonly parameterized following the empirical multiplicative approach by Jarvis (1976) which uses a predefined minimum resistance and multiple environmental response factors like in Zhang et al. (2003), Simpson et al. (2012) and Emberson et al. (2000). More advanced formulations often used by land surface models (Ran et al., 2017; Val Martin et al., 2014) are based on the CO₂ assimilation by plants during photosynthesis (Ball et al., 1987; Collatz et al., 1992). Both approaches rely on the choice and constraints of ecosystem-dependent parameters and have different advantages (Lu, 2018). A further role in coupling stomata to ecosystems is played by stomatal optimization models, whereas optimal stomatal activity with a maximum amount of carbon gain and a minimum loss of water is calculated based on ecophysiological processes (e.g. Cowan and Farquhar, 1977). Of particular interest are stomatal optimization models which, based on ecophysiological processes, maximize carbon gain while minimizing water loss. According to Wang et al. (2020), these models are promising in representing stomatal behaviour and improving carbon cycle modelling. Non-stomatal deposition has been less investigated by now; therefore, most models use predefined constant resistances or scale it with leaf area index (e.g. Val Martin et al., 2014; Simpson et al., 2012), while some apply an explicit parameterization based on the observational findings of enhance cuticular uptake under leaf surface wetness (Altimir et al., 2006).

The different parameterizations of the (surface) resistances cause main model uncertainties in computing dry deposition fluxes of trace gases, which depend on the response to hydroclimate and land-type-specific properties (Hardacre et al., 2015; Wu et al., 2018; Wesely and Hicks, 2000). Thereby, it has been shown that the original Wesely-based parameterization generally captures well the seasonal and diurnal cycles of dry deposition velocity, whereas model–observation discrepancy at seasonal scales arises from biased land type and leaf area index input data (Silva and Heald, 2018). Wong et al. (2019) stated that discrepancies of up to 8 ppb in ground-level ozone arise from different parameterizations.

The current dry deposition scheme of EMAC uses six surface types, where the parameterized processes represent the forest canopy as a whole (big-leaf approach). Thereby, the uptake over vegetation relies on stomatal deposition as the only pathway determined by the photosynthetically active radiation (Kerkweg et al., 2006). According to Fares et al. (2012) and Rannik et al. (2012), the stomatal uptake in parameterizations often lacks the dependence on meteorological and environmental variables (leaf area index, temperature, vapour pressure deficit). Moreover, several studies (e.g. Hogg et al., 2007; Fares et al., 2012; Clifton et al., 2017) found the contribution of an additional process to dry deposition at the leaf covering of plants. Zhang et al. (2002) firstly derived a parameterization from field studies which establishes the important link of this process to meteorology. In general, findings by Solberg et al. (2008), Andersson and Engardt (2010) and Wong et al. (2019) highlight the importance of considering the dry deposition–meteorology dependence in global models. Such an extension would realistically enhance the sensitivity of dry deposition to climate variability and would result in a more accurate prediction of ground-level ozone.

Given the importance of ozone as a major tropospheric oxidant, air pollutant and greenhouse gas, an accurate representation of dry deposition is desirable (Jacob and Winner, 2009). Additionally, the significance of a realistic representation of land–atmosphere feedbacks rises in light of the changing Earth's climate with the projected increase of extreme events' frequency and intensity (Coumou and Rahmstorf, 2012).

Here, we present a revision of the existing Wesely-based dry deposition scheme in the Modular Earth Submodel System (MESSy), which has a very simplified representation of vegetation and soil. The modifications are done by well-established findings about the controls of stomatal and cuticular uptake of trace gases. The calculation of stomatal deposition fluxes is extended by including the vegetation density, two meteorological adjustment factors and an improved soil moisture availability function for plant stomata following the multiplicative algorithm by Jarvis (1976). For the first time in MESSy, a parameterization for cuticular dry deposition dependent on important meteorological and environmental variables is implemented explicitly (Zhang et al., 2003). In Sect. 2, a description of the model setup and the simulations is provided, whereas especially the transition to the new vertical exchange scheme is described in detail. Subsequently, the new VERTEX scheme is evaluated. In Sect. 4, the impact of the changes on ozone dry deposition is evaluated on daily and seasonal scales by comparison with measurements at four different sites. Here, advantages, uncertainties and missing processes in the revised scheme are identified. Next, the global impact on ground-level ozone is assessed by separating the effect of the different implemented parameterizations. Then, Sect. 6 provides a description of the uncertainties in modelling stomatal conductance and Sect. 7 comprises an investigation of the sensitivity to model resolution. Section 8 summarizes the main findings and the remaining process and model uncertainties which form the basis for the provided recommendations. Section 9 describes planned future developments.

2 Model description

This study uses the ECHAM/MESSy atmospheric chemistry model. MESSy v2.54 (Jöckel et al., 2010) provides a flexible infrastructure for coupling processes to build comprehensive Earth system models (ESMs) and is utilized here with the fifth-generation European Centre Hamburg general circulation model (ECHAM5; Roeckner et al., 2003) as an atmospheric general circulation model. The dry deposition process of gases is calculated within the submodel DDEP (Kerkweg et al., 2006). This is described in Sect. 2.2. It relies on the VERTEX vertical exchange submodel (Sect. 2.1), former E5VDIFF, which contains the calculation of stomatal uptake (Eq. 5) and soil moisture stress (Eq. 12). The stomatal uptake parameterization is the base for the evapotranspiration scheme in VERTEX (Appendix B) which also incorporates the soil moisture stress.

2.1 The new VERTEX vertical exchange submodel

The VERTEX submodel represents land–atmosphere exchange and vertical diffusion as an alternative to the default E5VDIFF submodel in ECHAM5/MESSy. In 2016, Huug Ouwersloot branched VERTEX off from E5VDIFF. He optimized the code and applied bug fixes. This includes changes in calculation of the transfer coefficients for vertical diffusion, the latent heat vaporization, the convective transfer coefficient, the storage of the friction velocity, the roughness length over sea, the kinematic heat and moisture fluxes and the 2 and 10 m friction velocity. A detailed description can be found in the Supplement.

2.2 Dry deposition over vegetation

Dry deposition of trace gases to vegetation is calculated according to the multiple resistance scheme by Wesely (1989) shown in Fig. 1. The scheme, originally designed for a regional model with 11 land types and five seasonal categories, is used here with six generalized land types (Kerkweg et al., 2006). This was adapted by Ganzeveld and Lelieveld (1995) to the surface scheme of the ECHAM climate model (Klimarechenzentrum et al., 1992). The vegetation canopy is represented as one system; i.e. the detailed structure and plant characteristics are neglected (one big-leaf approach). Only one assumption about the canopy structure is made: the leaves are horizontally oriented and the leaf density is uniformly vertically distributed (Sellers, 1985). This is required in the formula for the calculation of stomatal resistance (Eq. 5).

Figure 1 Dry deposition resistance analogy (adapted from Zhang et al., 2003); modified resistors are marked with red boxes.

The resistances (in s m⁻¹) in the big-leaf approach account for mass and energy transfer mainly exerted by the boundary layer turbulence (R_a), molecular diffusion via the quasi-laminar boundary layer (R_qbr) and heterogeneous losses at the surface (R_s) (Kerkweg et al., 2006). With these, the dry deposition velocity v_d of a trace gas X (in s m⁻¹) is defined as follows:

$$\begin{array}{}\text{(1)}& v_{\mathrm{d}}\left(X\right)={\displaystyle \frac{\mathrm{1}}{R_{\mathrm{a}}+R_{\mathrm{qbr}}\left(X\right)+R_{\mathrm{s}}\left(X\right)}}.\end{array}$$

The dry deposition flux f_d(X) (in $\mathrm{molecules}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) is determined by multiplying the dry deposition velocity with the trace gas concentration C(X) (in molecules m⁻³):

$$\begin{array}{}\text{(2)}& f_{\mathrm{d}}\left(X\right)=-v_{\mathrm{d}}\left(X\right)\cdot C\left(X\right).\end{array}$$

The total resistance over land combines the resistances over snow, soil, vegetation (veg) and wet skin (ws) weighted by the respective land-covered fraction of a grid box (Kerkweg et al., 2006). In the following, only the latter two are considered. The resistances R_a and R_qbr are commonly parameterized with standard formulations from micrometeorology (Kerkweg et al., 2006; Wesely and Hicks, 1977). For the surface resistance over vegetation (R_s,veg), the parameterization according to Zhang et al. (2003) is used:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{1}}{R_{\mathrm{s},\mathrm{veg}}\left(X\right)}}={\displaystyle \frac{\mathrm{1}}{R_{\mathrm{can}}+R_{\mathrm{s},\mathrm{soil}}\left(X\right)}}\\ & +\underset{R_{\mathrm{leaf}}\left(X\right)}{\underbrace{{\displaystyle \frac{\mathrm{1}}{R_{\mathrm{cut}}\left(X\right)}}+{\displaystyle \frac{\mathrm{1}}{R_{\mathrm{stom},\mathrm{corr}}\left(X\right)+R_{\mathrm{mes}}\left(X\right)}}}},\end{array}\end{array}$$

which consists of the soil resistance (R_s,soil(X)), the in-canopy aerodynamic resistance (R_can) (as in Kerkweg et al., 2006) and the leaf resistance (R_leaf(X)). The gas uptake by leaves (leaf) can be separated in two parallel pathways: the cuticular (cut) and the stomatal (stom) with its associated mesophilic pathway (mes), where the latter has negligible resistance for ozone and highly soluble species (Wesely, 1989). In contrast to the default formulation in MESSy (Eq. A1), the resistances in the updated scheme are provided at canopy scale in order to avoid linear scaling with the leaf area index (LAI, area of leaves [m²]/surface area [m²]). In fact, the linear scaling of resistances with LAI assumes that the leaves act in parallel and overestimates the uptake for high LAI values (>3–4) (Ganzeveld et al., 1998; Baldocchi et al., 1987). Furthermore, the quasi-laminar boundary resistance of individual leaves is included through the cuticular deposition scheme (see Sect. 2.2.2), whereas R_qbr,veg is a separate term in the old formulation (Eq. A1).

Due to the importance of stomatal and cuticular uptake for ozone dry deposition, their respective parameterizations are modified in this study (see Sect. 2.2.1 and 2.2.2). Also, ozone deposition to soil might be an important pathway (Schwede et al., 2011; Fares et al., 2012) but process understanding remains limited due to scant observational constraints (Clifton et al., 2020b, a). Stella et al. (2011) showed an exponential increase of soil resistance with surface relative humidity in three agricultural data sets which, however, varies much between different sites (Stella et al., 2019) and contradicts previous findings (Altimir et al., 2006; Lamaud et al., 2002; Zhang et al., 2002). Models by, e.g. Mészáros et al. (2009); Lamaud et al. (2009) apply a linear dependence on soil water content for parameterizing soil resistance. These parameterizations rely on input variables like the minimum soil resistance (Stella et al., 2011) which introduce an uncertainty due to measurement constraints. Also, the performance of a mechanistic model as proposed by Clifton et al. (2020b) depends on many input variables and parameters whose estimation is challenging and mostly biome dependent. Due to these uncertainties and limitations, the current parameterization of soil resistance in MESSy (see Kerkweg et al., 2006 for details) was not modified in this study.

2.2.1 Uptake through plant stomata

The stomata are actively regulated openings between the plant cells. They are scattered mostly over the lower (hypostomatous) epidermis of leaves. They control the H₂O and CO₂ exchange by plants which is the essential coupling of vegetation to the atmosphere and therefore to weather and climate. Here, the default parameterization of stomatal resistance (Eq. A2) is extended by adding dependencies on meteorological variables according to the Simple Biosphere Model (SiB) by Sellers et al. (1986) based on previous work by Jarvis (1976) for temperature (T) and vapour pressure deficit (VPD):

$$\begin{array}{}\text{(4)}& R_{\mathrm{stom},\mathrm{corr}}\left(X\right)={\displaystyle \frac{R_{\mathrm{stom}}(\mathrm{PAR},\mathrm{LAI})}{f\left(W_{\mathrm{s}}\right)\cdot f\left(T\right)\cdot f\left(\mathrm{VPD}\right)}}\cdot {\displaystyle \frac{D_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}{D\left(X\right)}}.\end{array}$$

The optimal stomatal resistance for water (R_stom(PAR,LAI)) is corrected with the ratio of the molecular diffusivity of the species (D(X)) and water ($D_{{\mathrm{H}}_{\mathrm{2}}\mathrm{0}}$). The optimal stomatal resistance depends on the photosynthetically active radiation (PAR) and LAI (Ganzeveld and Lelieveld, 1995; Sellers, 1985):

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& R_{\mathrm{stom}}(\mathrm{PAR},\mathrm{LAI})=\\ & {\displaystyle \frac{kc}{\left[\frac{b}{\mathrm{dPAR}}\ln \left(\frac{d\exp \left(k\mathrm{LAI}\right)+\mathrm{1}}{d+\mathrm{1}}\right)-\ln \left(\frac{d+\exp (-k\mathrm{LAI})}{d+\mathrm{1}}\right)\right]}}\end{array},\end{array}$$

where k=0.9 is the extinction coefficient, c=100 s m⁻¹ is the minimum stomatal resistance, and a=5000 J m⁻³, b=10 W m⁻² and $d=\frac{a+b\cdot c}{c\cdot \mathrm{PAR}}$ are fitting parameters (Sellers, 1985). For historical reasons, LAI was set to 1 in order to obtain the stomatal resistance at leaf level (Ganzeveld and Lelieveld, 1995). This has been changed and the seasonal evolution of stomatal resistance now follows the LAI which, in our study, is based on a 5-year climatology of monthly normalized difference vegetation index (NDVI) satellite data (Ganzeveld et al., 2002).

First, the stomatal resistance is corrected by the inverse of the temperature stress factor (1∕f(T)) derived by Jarvis (1976):

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }f\left(T\right)=b_{\mathrm{3}}(T-T_l)(T_h-T{)}^{b_{\mathrm{4}}}\text{(7)}& {\displaystyle }& {\displaystyle }{b}_{\mathrm{3}}=(T_{\mathrm{0}}-T_l)(T_h-T_{\mathrm{0}}{)}^{-b_{\mathrm{4}}}\text{(8)}& {\displaystyle }& {\displaystyle }{b}_{\mathrm{4}}=(T_h-T_{\mathrm{0}})/(T_h-T_l),\end{array}$$

where the empirical parameters are T_h=318.15 K, T_l=268.15 K and T₀=298.15 K.

Secondly, following the analysis by Katul et al. (2009), a stress factor dependent on vapour pressure deficit (1∕f(VPD)) was added to the calculation of stomatal resistance in VERTEX:

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }{p}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O},\mathrm{sat}}\left(T\right)=\mathrm{0.61078}\exp \left({\displaystyle \frac{\mathrm{17.1}\cdot T\left(p_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\right)}{\mathrm{235}+T\left(p_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\right)}}\right)\text{(10)}& {\displaystyle }& {\displaystyle }\mathrm{VPD}=p_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O},\mathrm{sat}}\left(T\right)-p_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}=\left(\mathrm{1}-{\displaystyle \frac{\mathrm{RH}}{\mathrm{100}}}\right)p_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O},\mathrm{sat}}\left(T\right)\text{(11)}& {\displaystyle }& {\displaystyle }f\left(\mathrm{VPD}\right)={\mathrm{VPD}}^{-\frac{\mathrm{1}}{\mathrm{2}}},\end{array}$$

with $T\left(p_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\right)$ (in K) as the surface temperature, $p_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ (in kPa) as the pressure of water vapour and $p_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\left(T\right)$ [kPa] the pressure of saturated air. The vapour pressure deficit is calculated according to Kraus (2007).

While the stomatal resistance at canopy scale is actually calculated within the MESSy VERTEX submodel, the submodel DDEP uses it for the calculation of dry deposition fluxes. Thus, in DDEP, the user can choose between the old scheme based on Ganzeveld and Lelieveld (1995) and the new scheme actually using the stomatal resistance at canopy scale. The latter is activated by setting the DDEP &CTRL namelist parameter l_ganzeori to .FALSE.. How the stomatal resistance is calculated is chosen in VERTEX by the &CTRL namelist parameter irstom.

• irstom=0 activates the original parameterization.

• Separate modifications:

  • irstom=2: variable LAI,

  • irstom=3: T dependency and

  • irstom=4: VPD dependency, respectively.

• irstom=5: all modifications.

• irstom=1: stomatal resistance with variable LAI at leaf scale. Instead of choosing LAI of 1 in Eq. (5) to represent the stomatal resistance at leaf level, as is done by the original code, Eq. (5) is calculated at canopy level using the actual LAI and then multiplied by LAI to obtain the average stomatal resistance at leaf level. For this case, the DDEP namelist parameter l_ganzeori has to be set to .TRUE..

The stomatal activity of plants and the strength of surface–atmosphere coupling strongly depend on the parameterized plant–water stress (Combe et al., 2016). The soil water budget is represented by a “bucket scheme” where the soil water in a single layer is prescribed by a geographically varying predefined field capacity and soil wetness governed by transpiration, precipitation, runoff, snow melt and drainage (Roeckner et al., 2003). This scheme is used by so-called “first-generation” models. However, EMAC controls evapotranspiration through the stomatal resistance (Appendix B), which is the most important feature of biophysical (“second-generation”) land surface models. Thereby, the stomatal resistance is calculated often like the one described here (Eq. 4) including temperature, VPD and soil moisture stress (Seneviratne et al., 2010; Sellers et al., 1997). The originally used plant–water stress function of Jarvis (1976) and Sellers et al. (1986), however, relies on leaf water potential (f(ψ)) for different plant types, which is difficult to estimate. Hence, EMAC uses a plant–water stress function dependent on soil moisture (f(W_s). The default parameterization (Eq. A3, ifws=0 in VERTEX &CTRL), applies the permanent wilting point of plants (W_pwp, 35 % of field capacity¹) as a lower threshold in the calculation of the soil moisture stress factor (f(W_s)). However, soil moisture is significantly underpredicted by the model in some regions and the calculated f(W_s) can be 0 for long periods. This is unrealistic and effectively shuts down dry deposition, e.g. during the dry season in the Amazon region. For this reason, f(W_s) is parameterized here according to the original formulation by Delworth and Manabe (1988) by removing the lower limit:

$$\begin{array}{}\text{(12)}& f\left(W_{\mathrm{s}}\right)=\left\{\begin{array}{ll}\mathrm{1}& W_{\mathrm{s}}\left(t\right)>W_{\mathrm{cr}}\\ \frac{W_{\mathrm{s}}\left(t\right)}{W_{\mathrm{cr}}}& W_{\mathrm{s}}\left(t\right)\le W_{\mathrm{cr}},\end{array}\right.\end{array}$$

where W_s(t) is the surface soil wetness (in m). W_cr (in m) is defined as the critical soil moisture level (75 % of the field capacity) at which the transpiration of plants is reduced. The modified parameterization in Eq. (12) can be applied by setting the &CTRL parameter ifws=  in the VERTEX namelist.

2.2.2 Cuticular deposition

According to several field studies (e.g. Van Pul and Jacobs, 1994; Hogg et al., 2007; Fares et al., 2012), cuticular deposition is an important contributor to ozone uptake and should not be neglected in models. Therefore, an explicit parameterization of cuticular deposition as used in many North American air quality modelling studies (Huang et al., 2016; Kharol et al., 2018) has been implemented. The gas uptake by leaf surfaces is based on two parallel routes, for which an analogy to ozone (highly reactive) and sulfur dioxide (very soluble) is used. The cuticular resistance is calculated as

$$\begin{array}{}\text{(13)}& R_{\mathrm{cut}}\left(X\right)={\displaystyle \frac{R_{\mathrm{cut},\mathrm{d}}\left({\mathrm{O}}_{\mathrm{3}}\right)}{{\mathrm{10}}^{-\mathrm{5}}\cdot H\left(X\right)+s_{\mathrm{reac}}\left(X\right)}},\end{array}$$

where H(X) is the effective Henry's law coefficient as a measure of the solubility. The reactivity of a species is rated by the parameter s_reac. For highly reactive species (s_reac=1), the same property as for ozone is assumed (second term in Eq. 13), while for less reactive species ($s_{\mathrm{reac}}=\mathrm{0.1},\mathrm{0}$) the uptake is effectively reduced (Wesely, 1989). For soluble species, the uptake at wet skin is assumed to be similar to the one of sulfur dioxide and is calculated as

$$\begin{array}{}\text{(14)}& R_{\mathrm{ws}}\left(X\right)={\left[{\displaystyle \frac{\mathrm{1}/\mathrm{3}}{R_{\mathrm{cut},\mathrm{w}}\left({\mathrm{SO}}_{\mathrm{2}}\right)}}+{\mathrm{10}}^{-\mathrm{7}}\cdot H\left(X\right)+{\displaystyle \frac{s_{\mathrm{reac}}\left(X\right)}{R_{\mathrm{cut},\mathrm{w}}\left({\mathrm{O}}_{\mathrm{3}}\right)}}\right]}^{-\mathrm{1}},\end{array}$$

where R_cut,w(SO₂) and R_cut,w(O₃) are the resistances of sulfur dioxide and ozone at wet surfaces, respectively. The constant values of the default formulae (Eqs. A4, A5) are replaced by parameterizations which account for the meteorological dependence of cuticular uptake according to Zhang et al. (2002):

$$\begin{array}{}\text{(15)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{cut},\mathrm{d}}({\mathrm{O}}_{\mathrm{3}}/{\mathrm{SO}}_{\mathrm{2}})={\displaystyle \frac{R_{\mathrm{cut},\mathrm{d}\mathrm{0}}({\mathrm{O}}_{\mathrm{3}}/{\mathrm{SO}}_{\mathrm{2}})}{\exp (\mathrm{0.03}\cdot \mathrm{RH})\cdot {\mathrm{LAI}}^{\mathrm{0.25}}\cdot u_{*}}}\text{(16)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{cut},\mathrm{w}}({\mathrm{O}}_{\mathrm{3}}/{\mathrm{SO}}_{\mathrm{2}})={\displaystyle \frac{R_{\mathrm{cut},\mathrm{w}\mathrm{0}}({\mathrm{O}}_{\mathrm{3}}/{\mathrm{SO}}_{\mathrm{2}})}{{\mathrm{LAI}}^{\mathrm{0.5}}\cdot u_{*}}},\end{array}$$

where the cuticular resistance of O₃ and SO₂, respectively, is distinguished for dry canopies (R_cut,d) and wet canopies (R_cut,w) depending on relative humidity (RH in %), LAI (in m² m⁻²) and friction velocity (u_* in m s⁻¹). The input parameters are R_cut,d0(O₃)=5000 s m⁻¹, $R_{\mathrm{cut},\mathrm{w}\mathrm{0}}\left({\mathrm{O}}_{\mathrm{3}}\right)=\mathrm{300}$ s m⁻¹ and $R_{\mathrm{cut},\mathrm{d}\mathrm{0}}\left({\mathrm{SO}}_{\mathrm{2}}\right)=\mathrm{2000}$ s m⁻¹ (Zhang et al., 2002). For rain and dew conditions, values of 50 s m⁻¹ and 100 s m⁻¹ are prescribed for R_cut,w0(SO₂). In contrast to traditional approaches, these parameterizations also consider the aerodynamic and the quasi-laminar boundary resistances of individual leaves. For the usage in MESSy, this can be switched on via l_ganzeori=.FALSE. in the &CTRL namelist of DDEP.

2.3 Simulations

In order to answer the different research questions of this study, two different types of simulations have been performed (Table 1).

1. The first kind were simulations to investigate dry deposition and the effect of the modifications in VERTEX:
   these simulations are based on the Chemistry-Climate Model Initiative (CCMI) setup (Jöckel et al., 2016). To allow for comparison with measurements, the model dynamics have been nudged towards realistic meteorology by the assimilation of data from the European Centre for Medium-Range Weather Forecasts (ECMWF) (Jöckel et al., 2010). Additionally, the QCTM mode is used; i.e. the chemistry does not feed back to the dynamics, resulting in the same meteorology for all simulations (Deckert et al., 2011). All modifications for the dry deposition scheme are employed in a 7-year simulation (REV, 2009–2015). Additionally, a 1.5-year simulation covering the period 2017 to July 2018 (2017 as spin-up) has been performed to cover the measurement periods (Sect. 4). For the same periods, simulations with the same configuration, except applying the default dry deposition scheme (DEF), have been conducted. The individual effects of the different modifications are investigated by two 2-year simulations employing the different namelist switches (Sect. 2.2). Moreover, a free-running sensitivity simulation with an additional temperature and drought stress factor for evapotranspiration (Appendix B) has been performed aiming at an improved representation of local meteorology especially in the Amazon. The station simulation output and the global output are analysed in Sects. 4 and 5, respectively. In addition, two 2-year simulations are realized for different horizontal resolutions (REST42, REST63) to investigate the resolution dependency of dry deposition (Sect. 7). All these simulations use 31 model layers with the top at 10 hPa and take the first year of simulation as spin-up.

2. The second kind were simulations for the evaluation of VERTEX as the boundary layer scheme:
   two pure dynamical (i.e. without chemistry) 30-year simulations with the old (clim-E5) and the new boundary layer description (clim-VER), respectively, have been performed.

All simulations were performed at the Jülich Supercomputing Centre with the JURECA Cluster (Jülich Supercomputing Centre, 2018).

Table 1 List of EMAC simulations

3 VERTEX evaluation

In order to advise the usage of VERTEX (with the default settings) as the default vertical exchange submodel in MESSy, the dynamics produced by both submodels are compared. Therefore, two dynamical, free-running 30-year simulations have been performed using the E5VDIFF or the VERTEX submodels, respectively. To obtain a comparable radiative imbalance at TOA (top of the atmosphere) with VERTEX, the four cloud parameters have been tuned in advance according to Mauritsen et al. (2012). The tuning factors can be found in Table 2. The radiative imbalance at TOA is slightly positive at present-day conditions (Mauritsen et al., 2012; Stephens et al., 2012); here, E5VDIFF gives a negative value. The difference between the tuned VERTEX and E5VDIFF is small and within the uncertainty range of ±0.4 W m⁻².

Table 2 Overview of tuning parameter settings and global mean properties.

^* Stevens and Schwartz (2012)

Additionally, global mean values of surface temperature, cloud liquid water, relative humidity and planetary boundary layer height of EMAC using E5VDIFF and EMAC using VERTEX with the respective uncertainty range for the period 1979–2008 are represented in Fig. 2. The results for cloud liquid water and planetary boundary height show no significant differences between the VERTEX and E5VDIFF simulations since the annual mean of each falls in the confidence interval of the other. This is not always the case for surface temperature and relative humidity. However, the 30-year means of surface temperature and relative humidity simulated by E5VDIFF and VERTEX are not significantly different.

Figure 2 Global mean properties and the uncertainty range (95 % confidence interval in shaded) of the climatology simulations with E5VDIFF (clim-E5) and with VERTEX (clim-VER) for the period 1979–2008.

4 Evaluation with deposition measurements

To assess the impact of the code revision/modifications on the variability of dry deposition, we compare the sensitivity simulations DEF, REV, REV-fTfVPD, REV-fws and REV-NNTR (see Table 1, all at T106L31 resolution) with dry deposition measurements at four field sites (listed in Table 3). The chosen data sets are the best available of ozone dry deposition (flux data and ozone mixing ratio or velocity data) with the required temporal resolution and coverage of diverse biomes of the world. The analysis is aimed at covering the recent decade, which includes the most extreme drought and heat events (where the stomatal stress factors are aimed). For the reason of uniqueness and importance of atmospheric processes in a remote and pristine forest like the Amazon, we included measurements from, among others, the Amazon Tall Tower Observatory (ATTO). Ozone dry deposition fluxes were measured with the eddy covariance and gradient method (Ontario). From these data, deposition velocities were calculated by the means of ozone concentration data. The eddy covariance technique determines a turbulent flux by the covariance of the measured vertical velocity and the gas concentration. Due to the stochastic nature of turbulence, these measurements have an uncertainty of 10 % to 20 % under typical observation conditions (Rannik et al., 2016). For the gradient method used at the Borden forest research station, the dry deposition flux was estimated from concentration gradients below and above the canopy and the eddy diffusivity according to the Monin–Obukhov similarity theory. The estimated dry deposition velocities (V_d) show an uncertainty of ≈ 20 %, which is due to the assigned canopy, the inherent limitations of the algorithm and the measurement uncertainties in concentrations. However, results are in good agreement with other eddy covariance measurements (Wu et al., 2016).

Keronen et al. (2003)Fares et al. (2012)Fares et al. (2012–2014)Wu et al. (2018)

Table 3 Dry deposition measurements. In the description of vegetation/climate, the reported LAI (in m² m⁻²) is given in brackets; $v_{\mathrm{d}}^{\mathrm{mod}}$ and $v_{\mathrm{d}}^{\mathrm{obs}}$ are the average measured and modelled dry deposition velocity.

^a Meteorological measurement height. ^b Ozone data are not available here.

4.1 Annual cycle of dry deposition

The annual cycle of dry deposition is mainly driven by the evolution of vegetation and is generally represented well in models (Silva and Heald, 2018). Here, we use the long time series measured at Borden and Hyytiälä to identify the impact of the code modifications on the annual cycle of dry deposition velocity. The available micrometeorological data help to distinguish the different effects. From the hourly data, we calculated multiyear (2010–2012) monthly means. To explore the contribution of stomatal and cuticular uptake, the individual velocities are calculated for O₃ according to the model calculations (Kerkweg et al., 2006):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{G}_{\mathrm{cut},\mathrm{d}}={\displaystyle \frac{(\mathrm{1}-\mathrm{ws})\cdot (\mathrm{1}-\mathrm{cvs})\cdot \mathrm{veg}}{R_{\mathrm{cut},\mathrm{d}}\left({\mathrm{O}}_{\mathrm{3}}\right)}}\\ \text{(17)}& {\displaystyle }& {\displaystyle }{G}_{\mathrm{cut},\mathrm{w}}={\displaystyle \frac{\mathrm{ws}\cdot (\mathrm{1}-\mathrm{cvs})}{R_{\mathrm{cut},\mathrm{w}}\left({\mathrm{O}}_{\mathrm{3}}\right)}}\text{(18)}& {\displaystyle }& {\displaystyle }{G}_{\mathrm{ns}}=G_{\mathrm{cut},\mathrm{d}}+G_{\mathrm{cut},\mathrm{w}}\text{(19)}& {\displaystyle }& {\displaystyle }{G}_{\mathrm{stom}}={\displaystyle \frac{(\mathrm{1}-\mathrm{ws})\cdot (\mathrm{1}-\mathrm{cvs})\cdot \mathrm{veg}}{R_{\mathrm{stom},\mathrm{corr}}\left({\mathrm{O}}_{\mathrm{3}}\right)}}\text{(20)}& {\displaystyle }& {\displaystyle }{v}_{\mathrm{p}}={\displaystyle \frac{G_{\mathrm{p}}}{G_{\mathrm{stom}}+G_{\mathrm{ns}}}}\cdot v_{\mathrm{d}},\end{array}$$

where G names the individual conductances (inverse of resistance) of stomata (stom), dry cuticle (cut,d), wet cuticle (cut,w) and non-stomata (ns). Here, veg, ws and cvs give the vegetation fraction, the wet skin fraction and the snow-covered fraction, respectively. G_p and v_p are the individual conductance and the velocity of one pathway. Further terms are described in Sect. 2.2.

The multiyear (2010–2012) annual cycle of the simulated dry deposition velocity at Borden forest (Fig. 3a) captures the observed cycle well until June. The new scheme reproduces the observations better than the old scheme. This is a consequence of the increase in nighttime mean velocities due to the much larger cuticular contribution (Fig. B1a, b). However, due to the overestimated stomatal uptake in the default scheme (see Sect. 2.2.1), only slight deviations from the new dry deposition scheme are visible in the daily mean shown in Fig. 3a. The mismatch of the simulated and measured V_d from August to October is a consequence of the underestimation of relative humidity leading to too-low simulated cuticular deposition (Fig. 3c, e). This effect exceeds the impact of the overestimation of relative humidity (only) in summer, because the LAI is higher in summer. In general, the cuticular uptake parameterization accounting for LAI, friction velocity, RH and surface wetness conditions performs, in our simulations, better than parameterizations without these dependencies as expected from the study of Wu et al. (2018). Unfortunately, the cuticular uptake parameterization also introduces uncertainties to the modelled non-stomatal uptake. Moreover, accounting for biogenic volatile organic compounds (BVOCs) like in Makar et al. (2017) would enhance in-canopy loss of ozone, significantly increase non-stomatal dry deposition and lead to improved simulation results (Wu et al., 2018). The representation of in-canopy air chemistry is outside the scope of the present study but planned within a subsequent study.

In contrast, the amplitude of the annual cycle and the mean of dry deposition fluxes in Hyytiälä are overestimated by both schemes during spring and summer (Fig. 3b). For the default scheme, this is due to the oversimplification of the stomatal uptake that only accounts for a constant LAI of 1 m² m⁻² (see Sect. 2.2.1), which is far from the measured LAI of 3–4 m² m⁻² during this period (Keronen et al., 2003). Enabling the new scheme (REV), increases the dry deposition velocity which reproduces the measured values in autumn better. The contribution of non-stomatal dry deposition of 25 %–45 % during the day reported by Rannik et al. (2012) is represented partly by that. However, the new scheme leads to an even higher overestimation by the model from April to July. The sensitivity simulation REV-fws (default f(W_s)) points to the increase of the soil moisture stress function (see Sect. 2.2.1, Eq. 12) as one reason for the overestimation of V_d in summer (Fig. 3b, f). Moreover, the overestimation in June/July is partly (∼ 10 %) due to the too-high model LAI compared to the measured values of 3–4 (Fig. B2a). The remaining gap (Fig. 3f) can be explained by restricting the analysis to wet conditions (RH > 70 %) only and the analysis of the sensitivity simulation REF-fTfD (no f(T) and f(VPD)). This suggests that the overestimated V_d (Fig. B2c) in summer is due to the stress factors for stomatal uptake since the modelled and measured temperatures are a mismatch. VPD has been identified by Rannik et al. (2012) as a strong driver of daytime total deposition velocity, which confirms the importance of inclusion of VPD dependence for stomatal uptake.

Figure 3 Measured (obs) and modelled (DEF, REV) multiyear mean (2010–2012) and REV-fws (2010) annual cycle. (a, c, e) Borden forest and (b, d, f) Hyytiälä; arrows indicate 1σ.

4.2 Importance of stress factors for the diurnal variation of deposition

The short-term measurements at Lindcove research station and at ATTO are used to assess the impact of the stress factors on the diurnal cycle of dry deposition velocity in spring and summer. Additionally, micrometeorological and additional flux data make possible to consider the stomatal resistance (∼ inverse of the velocity, calculations according to Fares et al., 2012) and the underlying meteorological conditions. Since the respective micrometeorological measurements are not available at ATTO, data extracted from the ERA5 global climate reanalysis at the 1000 hPa pressure level (Copernicus, 2017) are used here.

The diurnal cycle of dry deposition velocity at the Lindcove research station follows the solar variation (Fares et al., 2012) and is generally well reproduced by the model with the best match in spring (Fig. 4). The revised dry deposition scheme reduces the underestimation of measured nighttime V_d due to the inclusion of cuticular uptake, which Fares et al. (2012) identified as an important ozone sink for exactly this measurement site. The measured dry deposition velocity increases at sunrise (around 15:00 UTC) and remains almost constant during the day. This is only reproduced by the revised dry deposition scheme. The comparison of the dry deposition velocity from the revised scheme (red line) and the velocity without stomatal T and VPD stress (gray line) in Fig. 4a illustrates the necessity of accounting for the stress factors. This is consistent with Fares et al. (2012), who reported a high negative correlation of V_d(sto) with VPD and temperature and related it to stomatal stress. The direct comparison of the stomatal resistances calculated from measured and modelled variables (Fig. 4c) shows an improvement of the modelled resistances (comparing DEF and REV). However, the modelled daytime stomatal resistance is still too high compared to the measurements. This points to an underestimation of stomatal uptake by the model during the day. A small fraction can be explained by the direct effect of the stomatal soil moisture stress in the model which does not occur in reality since the citrus orchard was watered during the measurement campaign. Contrastingly, in summer, the model underestimation of V_d is higher than in spring (Fig. 4b). As seen from the comparison of stomatal resistance values (Fig. 4d), the model underestimates the stomatal uptake. This is because the irrigation of the orchard leads to cooling sustained evapotranspiration and keeps f(T) low. Thus, in the model, a too-high temperature stress acts on the stomata. Moreover, neglecting the soil moisture stress on stomata would bring the stomatal resistance values closer since the irrigation at the site ensures a constant and high soil moisture. The irrigation of the citrus orchard during the day also enhances surface wetness and favours deposition at cuticles (Fares et al., 2012; Altimir et al., 2006) which cannot be captured by the model. Fares et al. (2012) estimate the stomatal contribution to only account for 20 %–45 % of the total daytime dry deposition flux during both seasons and point to soil deposition and reactions of ozone with NO and VOCs as major sinks at the citrus orchard, especially during flowering season. The contribution of these pathways is expected to be enhanced by the inclusion of further BVOCs within the chemical mechanism and the explicit parameterization of in-canopy residence and transport.

Figure 4 Diurnal cycles of measured (obs) and modelled (DEF, REV, REV-fTfD) ozone dry deposition velocity and stomatal resistance in spring and summer 2010 at Lindcove research station.

Figure 5 Diurnal cycles of measured (obs) and modelled (DEF, REV, REV-NNTR: free-running f(T) and f(VPD) for evapotranspiration) ozone dry deposition velocities in wet and dry seasons at ATTO (gray: standard deviation).

Figure 6 Multiyear (2010–2015) mean absolute values in boreal summer.

Figure 7 Multiyear (2010–2015) mean absolute changes in boreal summer: i.e. difference between the revised and default scheme (REV – DEF).

Figure 8 Mean changes (2010) of dry deposition velocity in boreal summer. (a) f(W_s) modification; (b) temperature and VPD stress.

Tropical forests are known to be effective O₃ sinks with observed mean midday maximum dry deposition velocity of 2.3 cm s⁻¹ (Rummel et al., 2007) due to much higher LAI compared to other sites (e.g. Lindcove). The measured dry deposition velocity at ATTO shown in Fig. 5a and b is no exception but shows a high variability (standard deviation). The diurnal cycle follows the solar radiation with maximum V_d at 15:00 UTC and highest amplitude during the wet season (April–May 2018). The amplitude of the diurnal cycle is highly underestimated in both EMAC simulations, with the highest mismatch during daytime. This is similar for other models. In fact, Hardacre et al. (2015) report a general and large underestimation of dry deposition velocities by models over tropical forests with highest predicted values of 0.25 cm s⁻¹. Here, the simulation with the revised dry deposition scheme (REV) shows only a minor increase of V_d during the wet season. Since stomatal uptake is known to be an important daytime sink (Freire et al., 2017), the underestimation of the total dry deposition flux is partly attributed to a too-low simulated stomatal uptake caused by the overestimation of temperature and the underestimation of relative humidity (Fig. B3). The increase of dry deposition velocity by the new scheme is mainly due to the lowered soil moisture stress on stomata (f(W_s)) shown in Fig. 5e. Freire et al. (2017) also links stomatal uptake to the efficiency of turbulent mixing in transporting ozone down to the canopy. In general, 10 % of the total ozone sink during daytime and 39 % during night are associated with in-canopy processes (Freire et al., 2017). Freire et al. (2017) and Bourtsoukidis et al. (2018) identified the oxidation of sesquiterpenes as an important contributor to the chemical nighttime sink. Cuticular deposition might also play a role in humid conditions during night (Rummel et al., 2007), which is underestimated by the model due to the biased relative humidity (Fig. 5c).

The uncertainty introduced by the mismatching meteorology becomes even more obvious when comparing measurements and simulations for November 2015. This month was characterized by temperatures of 2 to 3 degreeC above average and unusually little rainfall (compared to usual conditions in this season) due to a strong El Niño event (National Centers for Environmental Information, 2016). The dryness is overestimated by the model with a too-high temperature ($\mathrm{\Delta }=+\mathrm{5}$ to +8 K), too-low relative humidity ($\mathrm{\Delta }=-\mathrm{30}$ % to −40 %)) and too-dry soil. The lack of available soil moisture (f(W_s)=0) effectively shuts down stomatal deposition in the default simulation (DEF), whereas the modification of the soil moisture stress function (neglecting the artificial lower limit; see Eq. 12) in the revised model (REV) allows for an increased deposition (Fig. 5b). The temperature and relative humidity biases result in corresponding mismatching stress factors for the stomata that are double the ones derived from reanalysis data (Fig. 5f). This mismatch leads to an underestimation of stomatal uptake. This result is confirmed by the sensitivity simulation REV-NNTR for which no meteorological nudging has been applied, and the stress factors f(T) and f(VPD) are also used for the calculation of evapotranspiration. The REV-NNTR simulation yields much more realistic results compared to the measurements, capturing at least 50 % of the measured V_d during the day (Fig. 5b). This improvement is partly due to the omission of nudging, as the latter can have a detrimental effect on precipitation and evaporation (Jeuken et al., 1996). The temperature bias of the model is associated with the missing soil moisture buffer simulated by the bucket scheme. Incorporating a 5-layer scheme has been shown to lead to a more realistic soil water storage capacity, especially in the Amazon, and to a removal of this bias (Hagemann and Stacke, 2013). Nevertheless, the REV-NNTR simulation suggests that the stress factors f(T) and f(VPD) significantly contribute to buffer soil moisture and ameliorate the dryness bias.

Figure 9 Relative change [%] and absolute change [Tg yr⁻¹] (numbers on bars) of annual global loss by dry deposition of O₃, SO₂, HNO₃ and HCHO (REV – DEF).

Figure 10 Relative change of multiyear (2010–2015) mean at ground level (DEF – REV).

Figure 11 Relative change of multiyear (2010–2015) zonal mean (DEF – REV).

Figure 12 Relative change of multiyear (2010–2015) boreal summer mean (DEF – REV).

5 Global impact on ground-level ozone

Given the importance of dry deposition for ground-level ozone and the uncertainty of dry deposition parameterizations in models (Young et al., 2018; Hardacre et al., 2015), the global impact of the implemented code changes is assessed in this section. The global (boreal) summer mean distributions of deposition velocity and ground-level mixing ratio for O₃ shown in Fig. 6a–b are generally in the same range as reported for global models (e.g. Val Martin et al., 2014; Hardacre et al., 2015). However, like most global models, EMAC overestimates tropospheric ozone in comparison to satellite observations (Righi et al., 2015). Applying the revised dry deposition scheme increases the mean summer V_d by up to 0.5 cm s⁻¹ (Fig. 7a). The highest fraction of this increase arises from the inclusion of cuticular uptake at wet surfaces (V_cut,w) (Fig. B4b). The effect is large over the most northern continental regions (Fig. 7d) and even more pronounced where LAI is high like in Scandinavia and eastern Canada (for LAI distribution, see Fig. B4a). Additionally, the uptake at dry surfaces (V_cut,d) is enhanced with up to 0.3 cm s⁻¹ higher dry deposition velocity (Fig. 7c). This is because the default scheme applies a very high constant resistance for this process.

Concerning the stomatal deposition, the impacts of three different stress factors are considered. First, over relatively dry soil, i.e. where soil moisture exceeds 35 % of field capacity (wilting point of plants), the soil moisture stress is reduced by the modified parameterization. Neglecting the plants' wilting point as the lower limit for soil moisture stress on stomata weakens the dependency on field capacity. Thus, dry deposition is enhanced by up to 0.32 cm s⁻¹, as illustrated in Fig. 8a. Second, the inclusion of temperature and (third) VPD adjustment factors, indeed, leads to a spatially varying impact of ±0.27 cm s⁻¹ change in V_d (Fig. 8b). In humid and cold temperate regions, like Siberia and Canada, no temperature stress appears and the VPD adjustment factor increases the stomatal uptake. In the eastern US, Kazakhstan and central Amazon during boreal summer, stomata are stressed by temperature and VPD. This effect is overpredicted by the model, as the humidity over the Amazon forest is probably too low in the model (see Fig. B3). The stress factors are shown in Fig. B4c and d.

However, the overall decrease in ozone concentration dampens the impact of the change in dry deposition flux. In total, the changes by the revised dry deposition scheme increase the multiyear mean (2010–2015) loss of ozone by dry deposition from 946 to 1001 Tg yr⁻¹ (Young et al., 2018; Hu et al., 2017). Accordingly, (boreal) summer ground-level ozone over land is reduced by up to 12 ppb (24 %), peaking over Scandinavia, Asia, central Africa and eastern Canada (Fig. 7b). In the Northern Hemisphere, also the zonal mean of the tropospheric ozone mixing ratio show a noticeable reduction far from the ground compared to the default scheme (Fig. 11a). This has the potential to reduce the positive bias of tropospheric ozone on the Northern Hemisphere (20 %) reported by Young et al. (2018). However, besides ozone, also other atmospheric tracer gases are affected by the change in dry deposition. The global annual dry deposition flux of odd oxygen (O_x)², which includes many important tropospheric trace gases, increases from 978 to 1032 Tg yr⁻¹ due to the revision. This is in good agreement with the reported numbers by Hu et al. (2017) and Young et al. (2018). In Fig. 9, we show additionally the absolute and relative change of the multiyear annual average dry deposition loss of SO₂, NO₂, HNO₃ and HCHO. As a very soluble species, the loss of SO₂ is increased by the revised dry deposition scheme, whereas the predefined low cuticular and wet skin resistance of HNO₃ in the old scheme were replaced with the new mechanism, leading to an decrease in dry deposition. The altered loss of NO₂ and HCHO and other ozone precursors at ground level, especially soluble oxygenated VOCs, contributes to the total change in ozone loss. NO₂ is deposited almost 40 % more significantly, contributing to the net reduction in ozone production but is mostly counterbalanced by other processes. The change of HCHO dry deposition flux is small on a global and annual scale and only important regionally, mostly in (boreal) summer, when it decreases HCHO at ground level (Fig. 12b) by up to 25 %. Thereby, the change in wet uptake is highest but is partially counterbalanced by other effects. This leads to lower HO₂ production from HCHO photooxidation and lower NO-to-NO₂ conversion and thus lower ozone production (Seinfeld and Pandis, 2016). These effects also impact the OH mixing ratio (Figs. 10b, 11b) which controls the methane lifetime predicted by the model. However, for a clearer effect, a longer simulated time period would be needed. A detailed analysis of the trace gas budgets is beyond the scope of this paper and will be investigated in a subsequent study.

6 Uncertainties in modelling stomatal conductance

Dry deposition is a highly uncertain term in modelling ozone pollution (Young et al., 2018; Clifton et al., 2020a). Its representation is generally limited by a lack of measurements and process understanding but also to a large extent driven by the quality of land cover information (Hardacre et al., 2015; Clifton et al., 2020b). Although the dry deposition scheme by Wesely (1989) is commonly used in global and regional models (e.g. MOZART, GEOS-Chem), the approach has some constraints (Hardacre et al., 2015). The disadvantage of the big-leaf approach used in MESSy is that a vertical variation of leaf properties, affecting, for instance, the attenuation of solar radiation, is not considered (e.g. Clifton et al., 2020b). Regarding stomatal uptake, we neglect the mesophyll resistance as reactions inside the leaf are commonly assumed to not limit stomatal ozone uptake, whereas, besides mostly supporting laboratory studies (e.g. Sun et al., 2016), a few contradicting findings exist (e.g. Tuzet et al., 2011). The here-used empirical multiplicative algorithm by Jarvis (1976) for stomatal modelling has one general drawback concerning that the environmental responses to stomata are treated clearly in contrast to experimental evidence (Damour et al., 2010). However, Jarvis-type models have been shown to be able to compete with the semi-mechanistic A_net−g_s models which link stomatal uptake to the CO₂ assimilation during plant photosynthesis (Fares et al., 2013; Lu, 2018). The critics in Fares et al. (2013) state that the Jarvis model cannot capture the afternoon depression of ozone dry deposition is due to the original used VPD stress factor which has been replaced here by a mechanistic one based on the optimized exchange of CO₂ and water by plants (Katul et al., 2009). Furthermore, a larger set of land cover types is expected to improve the vegetation-dependent variation of dry deposition. The parameters used to model dry deposition of stomata, cuticle and soil are biome dependent and using generalized ones like for the input cuticular resistance can lead to differences in dry deposition (Hoshika et al., 2018). Exemplary discrepancies for the stomatal conductance calculated with different parameter sets are shown in Fig. 13 as the summer mean of 2010. Thereby, the temperature stress factor has been calculated as in Eq. (6) using the obtained surface temperature by EMAC (Fig. 13a, c) and applied to the model (DEFAULT) stomatal conductance (Eq. 17) with two different parameter sets for coniferous and mixed forest by Simpson et al. (2012)³ and Zhang et al. (2003)⁴. Jarvis (1976) obtained the parameters from a set of measurements in mixed hardwood/coniferous forest in Washington. In general, the parameters are related to measurements where the absolute values are influenced by multiple factors like genotype and local climatic conditions (Sulis et al., 2015; Tuovinen et al., 2009; Hoshika et al., 2018). So, for global modelling, mostly simplified parameters have to be used like in the European Monitoring and Evaluation Programme (EMEP) (Simpson et al., 2012).

Figure 13 Absolute difference of stomatal conductance applied with the temperature stress factor calculated for two different parameter sets by Simpson et al. (2012) (Simp) and Zhang et al. (2003) (Zh) in a comparison with the here-used parameter set by Jarvis (1976) (Jar).

7 Sensitivity to model resolution

The simulation of dry deposition depends on meteorology including boundary layer processes, radiation (cloud distribution and reflectivity) and ozone chemistry as well as on input fields like vegetation density (LAI) (Jones, 1992). Model horizontal resolution inherently affects the amplitude and distribution of (regridded) surface processes and the artificial dilution of ozone precursors that are emitted. This aspect is investigated here by analysing simulations at three different spatial resolutions: 2.8^∘ × 2.8^∘, 1.9^∘ × 1.9^∘ and 1.1^∘ × 1.1^∘ (REST42, REST63, and REV (T106) in Table 1).

In Fig. 14a, the resolution dependency is shown for the annual dry deposition flux of ozone on different continental regions. The annual dry deposition fluxes differ by up to 40 Tg yr⁻¹ globally between the different resolutions, with highest dry deposition at high resolution (T106). For the Northern Hemisphere (and consequently globally), this difference is driven by the higher annual mean ground-level ozone compared to the lower resolutions (Fig. 14c). However, this effect cannot be disentangled from the effect of decreased dry deposition velocity on ground-level ozone. Globally, increasing differences in O₃ are anti-correlated with relative humidity as shown in Fig. 15a ($\mathit{\rho }=-\mathrm{0.8}$). The impact of humidity on ozone chemistry is considered to be relatively weak (Jacob and Winner, 2009), but Kavassalis and Murphy (2017) showed for the US that only dry deposition establishes the observed anti-correlation between ozone and relative humidity. A dominating positive correlation of the dry deposition flux with the velocity only occurs on the Southern Hemisphere extratropics (SH_exT), which is highest between T63 and T106 (Fig. 15c). This can be attributed to discrepancies in stomatal deposition (Fig. 15d) driven by differences in humidity which might be caused by different moisture cycles and transpiration.

Figure 14 Ozone and dry deposition at three different resolutions (T42: 2.8^∘ × 2.8^∘, T63: 1.9^∘ × 1.9^∘, T106: 1.1^∘ × 1.1^∘) and the different regions: Northern Hemisphere extratropics (NH_exT: 90–30^∘ N), tropics (30^∘ N–30^∘ S), Southern Hemisphere extratropics (SH_exT: 90–30^∘ S) and the whole Earth (global).

Figure 15 Correlations of resolution dependent relative differences of ozone, dry deposition and meteorological variables for the whole Earth (global) and the Southern Hemisphere extratropics (SH_exT) for the four boreal seasons: spring (MAM), summer (JJA), autumn (SON) and winter (DJF).

8 Conclusion and recommendations

Dry deposition to the Earth's surface is a key process for the representation of ground-level ozone in global models. Its parameterizations constitutes a relevant part of the model uncertainty (Hardacre et al., 2015; Wu et al., 2018). Revising the dry deposition scheme of EMAC leads to an improved representation of surface ozone in regions with a positive model ozone bias (e.g. Europe). The highest increase in ozone dry deposition is due to the implementation of cuticular uptake whose contribution is important especially during night over moist surfaces. The extension of the stomatal uptake with temperature and VPD adjustment factors accounts for the desired link of plant activity to hydroclimate as recommended by Lin et al. (2019). Especially in drought-stressed regions (e.g. citrus orchards), the dependence on vapour pressure deficit leads to a realistic depression of stomatal uptake at noon. Also the dependence of dry deposition on soil moisture has been modified since the current representation of soil moisture in the model is not satisfactory. Specifically, the model simulates a too-dry soil for the Amazon basin, causing stomatal closure and thus an underestimation of dry deposition (Sect. 4.2). We have indications that the dry bias is a consequence of meteorological nudging in EMAC and also the missing representation of organized convection in the tropics (Mauritsen and Stevens, 2015). The sensitivity of the vegetation to droughts is comparably high in the Amazon region because the model soil cannot hold water in the catchment for a realistic time period and exhibits a memory effect (Hagemann and Stacke, 2013). Deeper root zones or buffering of the soil moisture below the root zone would improve the water holding capacity (Hagemann and Stacke, 2013; Fisher et al., 2007). With an improved representation of soil moisture, the more realistic parameterization of the soil moisture stress on stomatal uptake could be re-enabled. In general, the inclusion of the strong link between dry deposition and meteorology reveals some limitations of the dry deposition scheme associated with the inaccurate representation of local meteorology. The results also indicate that an improved representation of important non-stomatal dry deposition like in-canopy reactions of ozone with volatile organic compounds (e.g. citrus orchards; Sect. 4.2) would lower the positive model–observation discrepancy. This can be achieved with the inclusion of further BVOCs and an explicit parameterization of the transport dynamics in the boundary layer in model simulations (Makar et al., 2017). Explicit field measurements could foster further process understanding, which is required for a detailed process description within the models, especially over tropical rainforests. The seasonal variability of the simulated dry deposition velocity could be further improved by using as model input the time series of vegetation cover from imaging products which also capture land use changes and vegetation trend that are known to impact dry deposition significantly (Wong et al., 2019).

9 Outlook

The representation of gaseous dry deposition in MESSy will be further improved by using the MODIS time series of LAI which captures multi-annual vegetation changes. As the next step of dry deposition modelling in MESSy, a biome-dependent dry deposition model coupled to CO₂ assimilation (White et al., 2004) will be applied. Biome-dependent vegetation cover information, required for this scheme, is then provided by global input data which, however, represent only the annual cycle of vegetation. Coupling MESSy to the recently available dynamic vegetation model LPJ-GUESS, which provides detailed vegetation information with the temporal variability required for a climate model, could be a further improvement. By now, the one-way coupling of LPJ-GUESS as a MESSy submodel is only in the initial evaluation phase of the coupling with the atmospheric model (Forrest et al., 2020).

Appendix A: Default dry deposition scheme

The default dry deposition scheme of MESSy uses the following equations described in Kerkweg et al. (2006).

For surface resistance over vegetation (in s m⁻¹),

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{1}}{R_{\mathrm{s},\mathrm{veg}}\left(X\right)}}={\displaystyle \frac{\mathrm{1}}{R_{\mathrm{can}}+R_{\mathrm{s},\mathrm{soil}}\left(X\right)+R_{\mathrm{qbr},\mathrm{veg}}\left(X\right)}}+\\ & {\displaystyle \frac{\mathrm{LAI}}{r_{\mathrm{cut}}\left(X\right)}}+{\displaystyle \frac{LAI}{r_{\mathrm{stom},\mathrm{corr}}\left(X\right)+r_{\mathrm{mes}}\left(X\right)}},\end{array}\end{array}$$

where R_can(X), R_s,soil(X) and R_qbr,veg(X) are the in-canopy aerodynamic resistance, the soil resistance and the quasi-laminar boundary resistance at canopy scale (in s m⁻¹). r_cut(X), r_stom,corr(X) and r_mes(X) are the cuticular resistance, stomatal resistance and mesophyll resistance at leaf scale scaled with LAI (in m² m⁻²) to canopy scale.
For stomatal resistance,

$$\begin{array}{}\text{(A2)}& r_{\mathrm{stom},\mathrm{corr}}={\displaystyle \frac{r_{\mathrm{stom}}\left(\mathrm{PAR}\right)}{\mathrm{fws}}}\cdot {\displaystyle \frac{D_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}{D\left({\mathrm{O}}_{\mathrm{3}}\right)}}.\end{array}$$

For soil moisture stress function,

$$\begin{array}{}\text{(A3)}& f\left(W_{\mathrm{s}}\right)=\left\{\begin{array}{ll}\mathrm{1}& W_{\mathrm{s}}\left(t\right)\ge W_{\mathrm{cr}}(=\mathrm{75}\mathit{\%})\\ \frac{W_{\mathrm{s}}\left(t\right)-W_{\mathrm{pwp}}}{W_{\mathrm{cr}}-W_{\mathrm{pwp}}}& W_{\mathrm{pwp}}<W_{\mathrm{s}}\left(t\right)<W_{\mathrm{cr}}\\ \mathrm{0}& W_{\mathrm{s}}\left(t\right)\le W_{\mathrm{pwp}}(=\mathrm{35}\mathit{\%}).\end{array}\right.\end{array}$$

For cuticular resistance,

$$\begin{array}{}\text{(A4)}& r_{\mathrm{cut}}\left(X\right)={\displaystyle \frac{r_{\mathrm{cut}}\left({\mathrm{O}}_{\mathrm{3}}\right)}{{\mathrm{10}}^{-\mathrm{5}}\cdot H\left({\mathrm{O}}_{\mathrm{3}}\right)+s_{\mathrm{reac}}\left({\mathrm{O}}_{\mathrm{3}}\right)}},\end{array}$$

where $r_{\mathrm{cut}}{\mathrm{O}}_{\mathrm{3}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ s m⁻¹, H(O₃)=0.01 and s_reac=1.
For wet skin resistance,

$$\begin{array}{}\text{(A5)}& R_{\mathrm{ws}}\left({\mathrm{O}}_{\mathrm{3}}\right)={\left[{\displaystyle \frac{\mathrm{1}/\mathrm{3}}{R_{\mathrm{ws}}\left({\mathrm{SO}}_{\mathrm{2}}\right)}}+{\mathrm{10}}^{-\mathrm{7}}\cdot H\left({\mathrm{O}}_{\mathrm{3}}\right)+{\displaystyle \frac{s_{\mathrm{reac}}\left({\mathrm{O}}_{\mathrm{3}}\right)}{R_{\mathrm{cut},\mathrm{w}}\left({\mathrm{O}}_{\mathrm{3}}\right)}}\right]}^{-\mathrm{1}},\end{array}$$

where R_ws(O₃)=2000 s m⁻¹ and R_ws(SO₂)=100 s m⁻¹.

Appendix B: Evapotranspiration

Plants play a key role in the water and energy cycle and thus contribute to the land–atmosphere coupling, which drives the global climate. In this context, transpiration is an important process, as plants lose water during the necessary CO₂ uptake via their stomata. The amount depends on the aperture behaviour of the respective plant in the respective environmental conditions (Katul et al., 2012). Thus, the latent heat flux incorporates the canopy resistance. The formulation is based on the Monin–Obukhov stability theory:

$$\begin{array}{}\text{(B1)}& \begin{array}{rl}& E=\mathit{\rho }{C}_{\mathrm{h}}\left|\mathit{v}\right|\mathit{\beta }(q_{\mathrm{a}}-h{q}_{\mathrm{s}}(T_{\mathrm{s}},p_{\mathrm{s}}\left)\right)\\ & \mathit{\beta }={\left[\mathrm{1}+{\displaystyle \frac{C_{\mathrm{h}}\left|\mathit{v}\right|R_{\mathrm{stom}}}{\mathrm{fws}}}\right]}^{-\mathrm{1}},\end{array}\end{array}$$

where ρ is the density of air, $\left|\mathit{v}\right|$ is the absolute value of the horizontal wind speed and C_h is the transfer coefficient of heat, whereas $r_a=\mathrm{1}/\left(C_{\mathrm{h}}\right|\mathit{v}\left|\right).$ q_s and q_a are the saturation-specific humidity and the atmospheric specific humidity, whereas the relative humidity h at the surface limits the evapotranspiration from bare soil. β determines the ratio of transpiration between water-stressed plants (β<1) and well-watered plants (β=1) (Giorgetta et al., 2013; Schulz et al., 2001). The formula for the canopy stomatal resistance R_stom is given in Eq. (5). In order to adapt the transpiration to temperature and vapour pressure deficit, the T and VPD adjustment factors can be applied to R_stom inversely like in the new dry deposition scheme via izwet=1 in the VERTEX &CTRL namelist. The modification of the soil moisture stress function f(W_s) (old: Eq. A3; new: Eq. 12) affects evapotranspiration directly.

Figure B1 Measured and modelled (DEF, REV) annual cycles at Borden forest.

Figure B2 Measured (obs) and modelled (DEF, REV) multiyear (2010–2012) and REV-fTfD (2010) annual cycles at Hyytiälä.

Figure B3 Differences of meteorology between EMAC and ERA5 at ATTO.

Figure B4 Boreal summer mean vegetation and meteorological variables predicted by EMAC.