Update and evaluation of the ozone dry deposition in Oslo CTM3 v1.0

High concentrations of ozone in ambient air are hazardous not only to humans but to the ecosystem in general. The impact of ozone damage on vegetation and agricultural plants in combination with advancing climate change may affect food security in the future. While the future scenarios in themselves are uncertain, there are limiting factors constraining the accuracy of surface ozone modeling also at present: the distribution and amount of ozone precursors and ozone-depleting substances, the stratosphere–troposphere exchange, as well as scavenging processes. Removal of any substance through gravitational settling or by uptake by plants and soil is referred to as dry deposition. The process of dry deposition is important for predicting surface ozone concentrations and understanding the observed amount and increase of tropospheric background ozone. The conceptual dry deposition velocities are calculated following a resistance-analogous approach, wherein aerodynamic, quasi-laminar, and canopy resistance are key components, but these are hard to measure explicitly. We present an update of the dry deposition scheme implemented in Oslo CTM3. We change from a purely empirical dry deposition parameterization to a more process-based one which takes the state of the atmosphere and vegetation into account. We examine the sensitivity of the scheme to various parameters, e.g., the stomatal conductance-based description of the canopy resistance and the choice of ozone surface resistance, and evaluate the resulting modeled ozone dry deposition with respect to observations and multi-model studies. Individual dry deposition velocities are now available for each land surface type and agree generally well with observations. We also estimate the impact on the modeled ozone concentrations at the surface. We show that the global annual total ozone dry deposition decreases with respect to the previous model version (−37 %), leading to an increase in surface ozone of more than 100 % in some regions. While high sensitivity to changes in dry deposition to vegetation is found in the tropics and the Northern Hemisphere, the largest impact on global scales is associated with the choice of prescribed ozone surface resistance over the ocean and deserts.

1 Introduction

Ozone is an important trace gas for all life forms on Earth. Depending on the place of its occurrence, it has either a positive or negative connotation. In the stratosphere, ozone absorbs most of the ultraviolet (UV) light from the Sun within the range of 100–315 nm, thus shielding the Earth's surface from the most harmful UV radiation. In addition, ozone is a potent greenhouse gas in both stratosphere and troposphere. With a radiative forcing of 0.40±0.20 W m⁻², it is placed third, only surpassed by CO₂ and CH₄ (IPCC – Intergovernmental Panel on Climate Change, 2013, chap. 8).

In the troposphere, and in particular in ambient air, ozone is considered a highly toxic pollutant. Since the industrial revolution, tropospheric background ozone concentrations have been increasing in the Northern Hemisphere (IPCC – Intergovernmental Panel on Climate Change, 2013, chap. 2). In recent years, the number of episodes of peak concentrations has been, in general, decreasing in North America and Europe due to the implementation of air quality regulations (e.g., Fleming et al., 2018; Mills et al., 2018). At the same time, fast-developing countries, e.g., China or India, saw a significant increase in ozone-related air pollution. Continuously high concentrations of ambient air ozone are hazardous to the whole ecosystem. It is estimated that ozone is cause to an increase in premature deaths (WHO – World Health Organization, 2008), an average global loss of yield in the four major crops (wheat, rice, maize, and soybean) of about 3 %–15 % (Ainsworth, 2017), as well as 7 % loss in primary production in forestry (Wittig et al., 2009; Matyssek et al., 2012). The impact of ozone damage on vegetation and agricultural plants may affect food security in the future especially in Asia (Tang et al., 2013; Tai et al., 2014; Chuwah et al., 2015; Mills et al., 2018) and might be an important additional feedback to climate change (Sitch et al., 2007).

Elevated ozone levels at a site may originate from both the local production of ozone from its precursors, which are transported, and from advection of ozone itself. Long-range ozone transport occurs regularly and might be most important in regions that otherwise lack precursors. Tropospheric ozone is produced in complex photochemical cycles involving precursor gases such as carbon monoxide (CO) or volatile organic substances (VOCs – also known as hydrocarbons) in the presents of nitrogen oxides (NO_x). A typical reaction mechanism for CO is sketched in the following. In a sequence of rapid reactions, a peroxyl radical ${\mathrm{HO}}_{\mathrm{2}}^{\mathrm{•}}$ is formed through an initial reaction of CO with a hydroxyl radical ^•OH. Via a reaction between ${\mathrm{HO}}_{\mathrm{2}}^{\mathrm{•}}$ and NO, NO₂ is formed, which is then photolyzed. The resulting atomic oxygen reacts then with O₂ (and also under the presence of available co-reactants) to form an ozone molecule. Such a cycle leads to a net production via

$$\begin{array}{}\text{(R1)}& \mathrm{CO}+\mathrm{2}{\mathrm{O}}_{\mathrm{2}}+h\mathit{\nu }\to {\mathrm{CO}}_{\mathrm{2}}+{\mathrm{O}}_{\mathrm{3}}.\end{array}$$

Similar cycles involving VOCs exist (Monks et al., 2015). Another source of tropospheric ozone is downward transport from the stratosphere via stratosphere–troposphere exchange (STE) (WMO – Global Ozone Research and Monitoring Project, 2014). Based on observations, STE might only amount to 10 % (550±140 Tg a⁻¹) of the total global ozone budget in the troposphere, while ozone from chemical production is estimated to be 5000 Tg a⁻¹ (Monks et al., 2015). Ozone is removed from the atmosphere by photochemical reactions or scavenging processes. Major sinks are photolysis followed by a reaction with water vapor to form OH, reactions with HO₂, titration reactions, and dry deposition. We will come back to the latter later in this section and cover the implemented scheme in more detail in Sect. 2.1.

Since ozone is highly reactive, its global mean lifetime in the troposphere is roughly 22 d but ranges from a few days in the tropical boundary layer up to 1 year in the upper troposphere (Stevenson et al., 2005; Young et al., 2013). The abundance of tropospheric ozone therefore varies, e.g., with the time of day, season, altitude, location (Schnell et al., 2015), or weather conditions in general (Otero et al., 2018). Typical concentrations of surface ozone range from 10 ppb over the tropical Pacific to 100 ppb in the downwind areas of highly emitting sources (IPCC – Intergovernmental Panel on Climate Change, 2013, chap. 8). This variability poses a challenge for both trend analysis from observations as well as validation and intercomparison of models. From the observational side, the number of long-term observations (starting before the 1950s) is limited and restricted to mainly European sites. Most of these have indicated a doubling of tropospheric ozone since the 1950s (IPCC – Intergovernmental Panel on Climate Change, 2013, chap. 2). But especially the very low pre-industrial ozone abundance cannot be reproduced by the likes of most models. These early observations, however, were subject to interference by other species, e.g., SO₂. Among the participating models in the Atmospheric Chemistry and Climate Model Intercomparison Project (ACCMIP), there is a general tendency to underestimate tropospheric ozone burden (e.g., 10 %–20 % negative bias at 250 hPa in the Southern Hemisphere (SH) tropical region) (IPCC – Intergovernmental Panel on Climate Change, 2013, chap. 8). With respect to surface ozone, Schnell et al. (2015) conclude that all ACCMIP models, which reported hourly surface ozone, tend to overestimate surface ozone values in North America and Europe in comparison with available observations. A key to fathom these slightly contradicting results may lie in the used dry deposition schemes.

Removal of any substance from the atmosphere which is not involving rain, e.g., through gravitational settling or by uptake by plants, soil, and water, is referred to as dry deposition. The process of dry deposition is important for predicting surface ozone concentrations and understanding the observed amount and increase of tropospheric background ozone. It is estimated that 1000±200 Tg a⁻¹ of ozone are removed from the atmosphere by dry deposition processes (Monks et al., 2015). A newer study by Luhar et al. (2018), however, indicates much lower amounts (722.8±87.3 Tg a⁻¹) due to lower dry deposition to the oceans. Conceptually, dry deposition is a product between near-surface ozone concentration [O₃](z₀) (e.g., the lowermost model level) and a dry deposition velocity $v_{\text{DD}}^{{\mathrm{O}}_{\mathrm{3}}}$. Species-dependent dry deposition velocities $v_{\text{DD}}^i$, which are synonymously referred to as conductance Gⁱ, for any gaseous species i, are typically calculated following a resistance-analogous approach:

$$\begin{array}{}\text{(1)}& v_{\text{DD}}^i={\displaystyle \frac{\mathrm{1}}{R_a+R_b^i+R_c^i}},\end{array}$$

wherein aerodynamic R_a, quasi-laminar layer $R_b^i$, and canopy resistance $R_c^i$ are key components (Wesely, 1989; Seinfeld and Pandis, 2006). For all gases, R_a is the same, while $R_b^i$ and $R_c^i$ vary from gas to gas and also depend on land surface types (e.g., ice/snow, water, urban, desert, agricultural land, deciduous forest, coniferous forest). Originally, Wesely (1989) used fixed seasonal average dry deposition resistance for each land surface type. For all three types of resistance in this Wesely-type parameterization, more process-oriented formulations have been developed and validated over the years. Luhar et al. (2017) have validated ozone dry deposition to the ocean with respect to three different formulations of surface resistance. Based on the global atmospheric composition reanalysis performed in the ECMWF project Monitoring Atmospheric Composition and Climate (MACC) (MACC-II Consortium, 2011) and a more realistic process-based oceanic deposition scheme, Luhar et al. (2018) found that the ozone dry deposition to oceans amounts to 98.4±30.0 Tg a⁻¹. In particular, Luhar et al. (2018) found that the average surface resistance of ozone over the ocean (r_c=2200 s m⁻¹) is highly overestimated in most models. An update on the ozone surface resistance over snow- and ice-covered surfaces has been provided from combined model and observation studies (Helmig et al., 2007, $v_{\text{ice/snow}}^{{\mathrm{O}}_{\mathrm{3}}}=\mathrm{1}/\mathrm{10}\mathrm{000}$ m s⁻¹). Canopy conductance is parameterized at the single-leaf level (stomatal conductance) for various plant function types (PFTs) as well as for single plant species based on empirical studies (Jarvis, 1976; Ball et al., 1987; Simpson et al., 2012; Mills et al., 2017). But progress has also been made on process-oriented modeling of stomatal conductance (Anderson et al., 2000; Buckley, 2017). The variety of differing formulations and choices of parameters leads to a wide spread of results in model intercomparisons (Hardacre et al., 2015; Derwent et al., 2018) and about 20 % uncertainty on the resulting total dry deposition (Monks et al., 2015).

In Sect. 2, we will briefly describe Oslo CTM3, give a detailed account of the new dry deposition scheme (Sect. 2.1) as well as present pre-processing of meteorological input data to compute necessary input to the dry deposition scheme such as the beginning and duration of the greening season (GDAY, GLEN) and photosynthetic photon flux density (PPFD) (Sect. 2.2). In Sect. 3, we present sensitivity tests with respect to manifold parameters in the dry deposition scheme (Sect. 3.1) and validate our results with respect to results from the multi-model intercomparison of Hardacre et al. (2015) (Sect. 3.2), the MACC reanalysis (Sect. 3.3), and to surface ozone observations (Sect. 3.4). In Sect. 4, we will summarize and discuss our results and draw conclusions for further development of the model.

2 Model description

Oslo CTM3 is an offline, three-dimensional global chemistry transport model (CTM). The key components of Oslo CTM3 have been described and evaluated by Søvde et al. (2012). A detailed account of the capabilities of Oslo CTM3 in simulating anthropogenic aerosol forcing in the past and recent past using the Community Emission Data System (CEDS) historical emission inventory (Hoesly et al., 2018) is given by Lund et al. (2018). Oslo CTM3 can also be coupled to the Model of Emissions of Gases and Aerosols from Nature (MEGAN v2.10) (Guenther et al., 2006). A publication focusing on this is planned.

While the meteorological data driving Oslo CTM3 are given at a resolution of T159N80L60, with the highest model level at 0.02 hPa, it is very time and memory consuming to run Oslo CTM3 with full chemistry at this resolution. Therefore, we reduced the horizontal resolution to $\mathrm{2.25}{}^{\circ }\times \mathrm{2.25}{}^{\circ }$ in our experiments. In the following, we will give a detailed account of the new dry deposition scheme and the equations that we use.

2.1 Ozone dry deposition scheme

In the original dry deposition scheme, the state of the atmosphere was not taken into account. Dry deposition velocities were rather parameterized following the work of Wesely (1989) with parameter updates from Hough (1991). This means that seasonal day and night average deposition velocities for different land surface types (water, forest, grass, tundra/desert, and ice and snow) were in use. Day was distinguished from night by solar zenith angles below 90^∘. Winter was defined by temperatures below 273.15 K for grid boxes containing land masses. For ocean, winter and summer parameters are equal in this parameterization; therefore, no distinctive treatment was needed for ocean grid boxes. In addition, a reduced uptake due to snow cover above 1 m for forest and 10 cm for grass/tundra, respectively, was taken into account. We will refer to this parameterization as the “Wesely scheme”.

Regarding the new dry deposition scheme, we mainly follow Simpson et al. (2012) in their description of dry deposition used in the European Monitoring and Evaluation Programme (EMEP) MSC-W model (see also Emberson et al., 2000; Simpson et al., 2003; Tuovinen et al., 2004), which is used for air quality modeling implementing the Convention on Long-Range Transboundary Air Pollution (CLRTAP). We will refer to the new scheme as the “mOSaic scheme” throughout the rest of the paper. The mOSaic scheme is a more physical approach compared to the previously used Wesely scheme, because it takes state (e.g., pressure, temperature) of the atmosphere as well as dynamics (e.g., wind stress) of the boundary layer into account. To a certain degree, the global variety of plants and their variability throughout the seasons is also acknowledged. The mOSaic scheme is implemented for the gaseous species O₃, H₂O₂, NO₂, PAN, SO₂, NH₃, HCHO, and CH₃CHO. Since CO has a very small uptake and is not included in Simpson et al. (2003, 2012), the Wesely parameterization is kept. In addition to the gaseous species, Simpson et al. (2012) also modify aerosol deposition velocities, namely black carbon (BC) and organic carbon (OC), sulfuric aerosols (SO₄, MSA), and secondary organic aerosols (SOAs), but we have not updated our model with respect to these.

As displayed in Eq. (1), the dry deposition computation is subdivided into contributions from three different types of resistance. The main idea of a mosaic approach is to calculate these types of resistance separately for each land surface type k in each grid cell: $R_a^k$, $R_b^{i,k}$, and $R_c^{i,k}$. The grid cell average dry deposition velocity ${\stackrel{\mathrm{‾}}{v}}_{\mathrm{DD}}^i$ is then defined by weighting each individual $v_{\mathrm{DD}}^{i,k}$ by the corresponding land fraction factor f_k:

$$\begin{array}{}\text{(2)}& {\stackrel{\mathrm{‾}}{v}}_{\mathrm{DD}}^i=\sum _k{f}_k{v}_{\mathrm{DD}}^{i,k}.\end{array}$$

2.1.1 Aerodynamic resistance

In general, the aerodynamic resistance describes the turbulent transport of any substance down to the surface. To derive $R_a^k$, we follow Simpson et al. (2003, 2012) and compute a local friction velocity at reference height z_ref (Eq. 52, Simpson et al., 2012):

$$\begin{array}{}\text{(3)}& u_{*}^k={\displaystyle \frac{\stackrel{\mathrm{‾}}{u}\left(z_{\mathrm{ref}}\right)\cdot \mathit{\kappa }}{\ln \left(\frac{z_{\mathrm{ref}}-d_k}{z_{\mathrm{0}}^k}\right)-{\mathrm{\Psi }}_{\text{m}}\left(\frac{z_{\mathrm{ref}}-d_k}{L}\right)+{\mathrm{\Psi }}_{\text{m}}\left(\frac{z_{\mathrm{0}}^k}{L}\right)}},\end{array}$$

with the average wind speed $\stackrel{\mathrm{‾}}{u}\left(z_{\mathrm{ref}}\right)$ at reference height, the Kármán constant κ=0.40, the integrated stability equation for momentum Ψ_m (e.g., Garratt, 1992), a grid average Obukhov length L, displacement height d_k, and roughness length $z_{\mathrm{0}}^k$ ($d_k=\mathrm{0.78}\cdot h_k\left(\mathrm{lat}\right)$, $z_{\mathrm{0}}^k=\mathrm{0.07}\cdot h_k\left(\mathrm{lat}\right)$ for forests, $d_k=\mathrm{0.7}\cdot h_k\left(\mathrm{lat}\right)$, $z_{\mathrm{0}}^k=\mathrm{0.1}\cdot h_k\left(\mathrm{lat}\right)$ for vegetation other than forests). Taking the height of vegetation into consideration, we have chosen the model level such that ${\stackrel{\mathrm{‾}}{z}}_{\mathrm{ref}}\approx \mathrm{45}$ m. Using the derived $u_{*}^k$ from Eq. (3), a local Obukhov length L_k can be obtained from (Eq. 8, Simpson et al., 2012):

$$\begin{array}{}\text{(4)}& L_k=-{\displaystyle \frac{\mathit{\rho }{c}_{\text{p}}{T}_{\mathrm{2}\mathrm{m}}{u}_{*}^k}{\mathit{\kappa }gH}}.\end{array}$$

Herein, H is the sensible heat flux, g is the standard gravitational acceleration, c_p the specific heat capacity, and T_2 m the 2 m temperature. With these, we can compute the aerodynamical resistance for each land surface type (Eq. 8.8, Simpson et al., 2003):

$$\begin{array}{}\text{(5)}& R_a^k={\displaystyle \frac{\mathrm{1}}{\mathit{\kappa }{u}_{*}^k}}\left[\ln \left({\displaystyle \frac{z_{\mathrm{ref}}-d_k}{z_{\mathrm{0}}^k}}\right)-{\mathrm{\Psi }}_{\text{h}}\left({\displaystyle \frac{z_{\mathrm{ref}}-d_k}{L_k}}\right)+{\mathrm{\Psi }}_{\text{h}}\left({\displaystyle \frac{z_{\mathrm{0}}^k}{L_k}}\right)\right],\end{array}$$

with the integrated stability equation for heat Ψ_h (e.g., Garratt, 1992). Both integrated stability functions (Ψ_m, Ψ_h) and corresponding parameters are listed in Sect. S1 in the Supplement.

2.1.2 Quasi-laminar layer resistance

The quasi-laminar layer resistance $R_b^{i,k}$ is species specific and differs over land and ocean surfaces. Over land, we use (Eq. 53, Simpson et al., 2012)

$$\begin{array}{}\text{(6)}& R_b^{i,k}={\displaystyle \frac{\mathrm{2}}{\mathit{\kappa }{u}_{*}}}\cdot {\left({\displaystyle \frac{{\text{Sc}}_i}{\text{Pr}}}\right)}^{\frac{\mathrm{2}}{\mathrm{3}}},\end{array}$$

wherein Pr is the Prandtl number (0.72 for air and other gases) and Sc_i is the Schmidt number for a gas i. Equation (6) differs from a similar formulation in Seinfeld and Pandis (2006) by a factor of roughly 1.25. From ${\text{Sc}}_i=\mathit{\nu }/D_i$, with the kinematic viscosity of air ν, we derive a Schmidt number in water equivalent:

$$\begin{array}{}\text{(7)}& {\text{Sc}}_i={\displaystyle \frac{D_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}{D_i}}\cdot {\text{Sc}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}},\end{array}$$

with the molecular diffusivity for any gas D_i, the Schmidt number of water (${\text{Sc}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}=\mathrm{0.6}$) and its molecular diffusivity ($D_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}=\mathrm{0.21}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻¹). The used ratios $D_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}/D_i$ are taken from Simpson et al. (2012, Table S18). Over the ocean, we use (Eq. 54, Simpson et al., 2012)

$$\begin{array}{}\text{(8)}& R_b^i={\displaystyle \frac{\mathrm{1}}{\mathit{\kappa }{u}_{*}}}\cdot \ln \left({\displaystyle \frac{z_{\mathrm{0}}}{D_i}}\cdot \mathit{\kappa }{u}_{*}\right),\end{array}$$

with an imposed lower threshold of 10 s m⁻¹ and an upper limit of 1000 s m⁻¹. The computation of roughness length z₀ over the ocean is divided into a calm and a rough sea case, with a threshold of 3 m s⁻¹. For calm sea, we apply the following upper limit (Hinze, 1975; Garratt, 1992, with a slightly higher coefficient of 0.135):

$$\begin{array}{}\text{(9)}& z_{\mathrm{0}}^{\text{calm}}=\text{min}\left\{\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}},\mathrm{0.135}\cdot {\displaystyle \frac{\mathit{\nu }}{u_{*}}}\right\}.\end{array}$$

The kinematic viscosity of air ν herein can be computed from

$$\begin{array}{}\text{(10)}& \mathit{\nu }={\displaystyle \frac{\mathit{\mu }}{\mathit{\rho }}}={\displaystyle \frac{\mathit{\mu }\left(T\right)}{\frac{P_{\mathrm{0}}}{T\cdot R_{\text{air}}}}}.\end{array}$$

For the temperature-dependent dynamic viscosity of air μ(T), we chose a linear fit to Sutherland's law through the origin within the temperature range $\mathit{\{}T\in \mathbb{R}|(\mathrm{243.15}<T<\mathrm{313.15})\mathrm{K}\mathit{\}}$: $\mathit{\mu }\left(T\right)=\mathrm{6.2}\times {\mathrm{10}}^{-\mathrm{8}}\mathrm{kg}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}\cdot T$. But despite its rough accuracy, we found that the choice of μ(T) has no effect on ${\stackrel{\mathrm{‾}}{R}}_b^{i,k}$ (Sect. S2, Figs. S1–S2). In Eq. (10), ρ is substituted by the air density using the ideal gas law. P₀ is the surface pressure, as T the 2 m temperature is chosen, and R_air is the universal gas constant for air. The rough sea case follows the method of Charnock (1955) and Wu (1980):

$$\begin{array}{}\text{(11)}& z_{\mathrm{0}}^{\text{rough}}=\text{min}\left\{\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}},\mathrm{0.018}\cdot {\displaystyle \frac{u_{*}^{\mathrm{2}}}{g}}\right\},\end{array}$$

with a gravitational acceleration g=9.836 m s⁻². The allowed maximum roughness length in both cases is set to 2 mm. Since the z₀ values computed with this parameterization are rather small ($\mathrm{0}<z_{\mathrm{0}}^{\text{calm}}<\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}\mathrm{m}$, $\mathrm{0}<z_{\mathrm{0}}^{\text{rough}}<\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}$ m), $R_b^i$ is set to its lower limit of 10 s m⁻¹ in about 91 % of all cases (see Sect. S2, Fig. S3).

2.1.3 Surface resistance

The surface resistance consists of both stomatal and non-stomatal resistance.

The stomatal conductance is a measure of the rate of CO₂ exchange and evapotranspiration through the stomata of a leaf. There are several environmental conditions affecting the opening and closing of the stomata and hence the capability of respiration (e.g., light, available water). Stomata sluggishness, a state in which the stomata can no longer fully close, has been reported as ozone-induced damage (Hoshika et al., 2015) but is not taken into account in our formulation. To reflect part of the underlying mechanism, the leaf-level molar stomatal conductance in the mOSaic scheme is computed using a common multiplicative ansatz (Ball et al., 1987; Mills et al., 2017):

$$\begin{array}{}\text{(12)}& g_{\text{sto, m}}^k=g_{\text{max, m}}^k\cdot f_{\text{phen}}^k\cdot f_{\text{light}}^k\cdot max\left\{f_{\text{min}}^k,f_T^k\cdot f_D^k\cdot f_{\text{SW}}^k\right\}.\end{array}$$

The factors herein are normalized and vary within the range 0–1. They account for leaf phenology (f_phen), light (f_light), temperature (f_T), water vapor pressure deficit (f_D), and soil water content (f_SW). All factors differ with land use type k. For clarity reasons, we drop this index in the following, as long as it is not necessary for the equation's completeness. The maximum molar stomatal conductance is given by $g_{\text{max, m}}^k$, which is in units of $\mathrm{mmol}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{2}}$. A unit conversion to m s⁻¹ is necessary in our model:

$$\begin{array}{}\text{(13)}& g_{\text{sto}}^k=g_{\text{sto, m}}^k\cdot R\cdot {\displaystyle \frac{T_{\mathrm{0}}}{P_{\mathrm{0}}}}.\end{array}$$

Herein, R is the universal gas constant. To annotate the differing units, we use the index “m” in Eq. (12). The temperature adjustment f_T is computed from

$$\begin{array}{}\text{(14)}& f_{\text{T}}={\displaystyle \frac{T_{\text{2\hspace{0.17em}m}}-T_{\text{min}}}{T_{\text{opt}}-T_{\text{min}}}}\cdot {\left({\displaystyle \frac{T_{\text{max}}-T_{\text{2\hspace{0.17em}m}}}{T_{\text{max}}-T_{\text{opt}}}}\right)}^{\mathit{\beta }},\end{array}$$

with $\mathit{\beta }=\frac{T_{\text{max}}-T_{\text{opt}}}{T_{\text{opt}}-T_{\text{min}}}$. The parameters T_min, T_max, and T_opt are tabulated for various plant functional types. All parameters are taken from Simpson et al. (2012, Tables S16, S19). Since f_T turns negative outside the range defined by T_min, T_max, we impose a lower limit of 0.01 for numerical reasons.

The water vapor deficit (VPD) is proportional to the saturation partial pressure of water ($P_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}^{\text{s}}$) and relative humidity (RH):

$$\begin{array}{}\text{(15)}& \text{VPD}=P_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}^{\text{s}}\cdot (\mathrm{1}-\text{RH}/\mathrm{100}).\end{array}$$

Using tabulated values of f_min, D_min, and D_max, the water vapor pressure deficit penalty factor f_D can be computed:

$$\begin{array}{}\text{(16)}& f_{\text{D}}=f_{\text{min}}+(\mathrm{1}-f_{\text{min}})\cdot {\displaystyle \frac{D_{\text{min}}-\text{VPD}}{D_{\text{min}}-D_{\text{max}}}}.\end{array}$$

The penalty factor with respect to available soil water (SW) f_SW is defined as

$$\begin{array}{}\text{(17)}& f_{\text{SW}}=\left\{\begin{array}{ll}\mathrm{1}& \text{if}\text{SW}\ge \mathrm{0.5},\\ \mathrm{2}\cdot \text{SW}& \text{if}\text{SW}<\mathrm{0.5}.\end{array}\right.\end{array}$$

SW is evaluated at a soil depths of 0.28–1 m, which corresponds to SWVL3 in OpenIFS.

The phenology of a plant typically describes its life cycle throughout a year; e.g., at midlatitudes and for deciduous species, it starts with the emergence of leaves in spring and ends in fall. In the mOSaic scheme, phenology is parameterized with respect to the start of the greening season (SGS) and its end (EGS). Details about our treatment of these are given in Sect. 2.2.1. In summary, our adaption of the f_phen parameterization reads as follows:

$$\begin{array}{}\text{(18)}& \begin{array}{rl}& f_{\text{phen}}=\\ & \left\{\begin{array}{l}\text{if}\text{GLEN}\ge \mathrm{365}\mathrm{1}\text{(explicitly excluding tropics)}\\ \text{if}\text{GDAY}=\mathrm{0}\mathrm{0}\\ \text{else}\left\{\begin{array}{ll}\text{if}\text{GDAY}\le {\mathit{\varphi }}_{\text{AS}}& {\mathit{\varphi }}_a\\ \text{if}\text{GDAY}\le {\mathit{\varphi }}_{\text{AS}}+{\mathit{\varphi }}_e& {\mathit{\varphi }}_b+({\mathit{\varphi }}_c-{\mathit{\varphi }}_b)\\ & \cdot (\text{GDAY}-{\mathit{\varphi }}_{\text{AS}})/{\mathit{\varphi }}_e\\ \text{if}\text{GDAY}\le \text{GLEN}& {\mathit{\varphi }}_c\\ -{\mathit{\varphi }}_{\text{AE}}-{\mathit{\varphi }}_f& \\ \text{if}\text{GDAY}\le \text{GLEN}-{\mathit{\varphi }}_{\text{AE}}& {\mathit{\varphi }}_d+({\mathit{\varphi }}_c-{\mathit{\varphi }}_d)\cdot (\text{GLEN}\\ & -{\mathit{\varphi }}_{\text{AE}}-\text{GDAY})/{\mathit{\varphi }}_f\\ \text{else}& {\mathit{\varphi }}_d.\end{array}\right.\end{array}\right.\end{array}\end{array}$$

Herein, we use the SGS- and EGS-derived parameters: day of greening season (GDAY), the time elapsed starting at the SGS, and the total length of the greening season (GLEN), the time span between EGS and SGS. The parameters ϕ_a, ϕ_b, ϕ_c, and ϕ_d define start or end points in the five phases of phenology in the mOSaic scheme, while ϕ_e, ϕ_f, ϕ_AS, and ϕ_AE control the temporal timing (Fig. 1). If GLEN is zero, we are, e.g., in Arctic regions, and there is no vegetation anyway; therefore, f_phen=0. Before the start of the greening season (GDAY=0), f_phen=0. Since this phenology is tuned to Northern Hemisphere (NH) midlatitudes, it does not apply to the tropics. We therefore decided to set f_phen=1 if GLEN is greater than or equal to 365, which is the case in the tropics.

Figure 1 Sketch of the five different phases in plant phenology f_phen in accordance with Eq. (18).

Light in the wavelength band 400–700 nm to which the plant chlorophyll is sensitive is called photosynthetic active radiation (PAR). The integral of PAR over these wavelengths is the photosynthetic photon flux density (PPFD). The correction factor f_light in response to varying PPFD is

$$\begin{array}{}\text{(19)}& f_{\text{light}}=\mathrm{1}-\text{exp}(-{\mathit{\alpha }}_{\text{light}}\cdot \text{PPFD}).\end{array}$$

In the mOSaic scheme, non-stomatal conductance is explicitly calculated for O₃, SO₂, HNO₃, and NH₃. For all other species, an interpolation between O₃ and SO₂ values is carried out. The non-stomatal conductance for O₃ consists of two terms: one depending on vegetation type and one depending on the soil/surface. For each land surface type k, we can write

$$\begin{array}{}\text{(20)}& G_{\text{ns}}^{{\mathrm{O}}_{\mathrm{3}},k}={\displaystyle \frac{{\text{SAI}}_k}{r_{\text{ext}}}}+{\displaystyle \frac{\mathrm{1}}{R_{\text{inc}}^k+R_{\text{gs}}^{{\mathrm{O}}_{\mathrm{3}},k}}}.\end{array}$$

SAI_k is the surface area index for vegetation type k, which is leaf area index (LAI) plus a value that represents cuticles and other surfaces. The external leaf resistance is defined by

$$\begin{array}{}\text{(21)}& r_{\text{ext}}=\mathrm{2000}\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}\cdot F_{\text{T}}.\end{array}$$

Herein, F_T is a temperature correction factor for temperatures below $-\mathrm{1}{}^{\circ }\mathrm{C}$ and $\mathit{\{}{F}_{\text{T}}\in \mathbb{R}|(\mathrm{1}\le F_{\text{T}}\le \mathrm{2})\mathit{\}}$:

$$\begin{array}{}\text{(22)}& F_{\text{T}}=\exp (-\mathrm{0.2}\cdot (\mathrm{1}+{\mathit{\theta }}_{\text{2\hspace{0.17em}m}}\left)\right).\end{array}$$

θ_2 m is the 2 m temperature in ^∘C. For most land surface types, SAI≡LAI. Some exceptions are

$$\begin{array}{}\text{(23)}& \begin{array}{rl}& \text{SAI}=\\ & \left\{\begin{array}{ll}\text{LAI}+\mathrm{1}& \text{if}\mathrm{forest}/\mathrm{wetland},\\ \text{LAI}\cdot \mathrm{5}/\mathrm{3.5}& \text{if}\text{cropland, first part of growing season},\\ \text{LAI}+\mathrm{1.5}& \text{if}\text{cropland, second part of growing season},\\ \mathrm{0}& \text{if}\mathrm{cropland},\mathrm{winter}.\end{array}\right.\end{array}\end{array}$$

Extending the mOSaic scheme to the Southern Hemisphere, we use the growing season for crops defined in Table 1.

Table 1 Definition of growing season for crops used in Oslo CTM3 in the NH and SH.

In this way, vegetation affects the conductance also by being there, not only by uptake through the stomata. The in-canopy resistance R_inc (Erisman et al., 1994) is then modified with respect to each (vegetated) land surface type in k:

$$\begin{array}{}\text{(24)}& R_{\text{inc}}=b\cdot {\text{SAI}}_k\cdot {\displaystyle \frac{h_k\left(\text{lat}\right)}{u_{*}}},\end{array}$$

where h_k(lat) is the latitude-dependent vegetation height (see explanation at the end of this section) and b=14 m⁻¹ is an empirical constant. The canopy resistance described in Simpson et al. (2012) does not take temperature and snow into account and is zero for non-vegetated surfaces, but we will adopt the correction previously used in the Oslo CTM3 Wesely scheme.

As initially mentioned, the necessary depth of snow to cover a certain type of vegetation differs. Therefore, we calculate a snow cover fraction f_snow using the snow depth S_D, which is available in units of meters of water equivalent from the meteorological input data, scaled to 10 % of the vegetation height. $R_{\text{gs}}^{{\mathrm{O}}_{\mathrm{3}},k}$ is tabulated. We correct for temperature by F_T and for snow cover fraction:

$$\begin{array}{}\text{(25)}& {\displaystyle \frac{\mathrm{1}}{R_{\text{gs}}^{{\mathrm{O}}_{\mathrm{3}},k}}}={\displaystyle \frac{\mathrm{1}-f_{\text{snow}}^k}{{\widehat{R}}_{\text{gs}}^{{\mathrm{O}}_{\mathrm{3}},k}}}+{\displaystyle \frac{f_{\text{snow}}^k}{R_{\text{snow}}^{{\mathrm{O}}_{\mathrm{3}},k}}}.\end{array}$$

The bulk canopy conductance is then defined as

$$\begin{array}{}\text{(26)}& G_c^k=\text{LAI}\cdot g_{\text{sto, m}}^k+G_{\text{ns}}^k,\end{array}$$

wherein LAI is the one-sided leaf area index taken from ISLSCP2 FASIR, g_sto the leaf-level stomatal conductance, and G_ns the bulk non-stomatal conductance.

2.1.4 Latitude-dependent vegetation height

The vegetation height h_k(lat) as described by Simpson et al. (2012) is linearly decreasing with latitude between 60 and 74^∘ N. To adapt this to a global model, we made a few additional assumptions. The tabulated height for each vegetation type h_k in the mOSaic scheme is regarded as constant at midlatitudes (40–60^∘). Towards the poles, we decrease the height of each vegetation type using the same rate as described in Simpson et al. (2012). At a latitude of 74^∘, a minimum height of $\mathrm{3}/\mathrm{10}\cdot h_k$ is reached and kept constant. Towards the Equator, we increase the height linearly so that at a latitude of 10^∘ a maximum height of 2⋅h_k is reached which is then held constant. We also assume symmetry in both hemispheres. Presuming a typical tree height of 20 m at midlatitudes, this stepwise function yields a height of 8 m at high latitudes and 40 m in the tropics, which is not unrealistic. For four example PFTs, results are shown in the Supplement (Sect. S3, Fig. S4).

2.1.5 Mapping of land surface types

Oslo CTM3 is configured to read land surface types from either the ISLSCP2 (https://daac.ornl.gov/cgi-bin/dataset_lister.pl?p=29, last access: 20 November 2017) product from MODIS or Community Land Model (CLM; http://www.cgd.ucar.edu/tss/clm/, last access: 7 November 2017) 2 categories, which have to be mapped to the land surface types used in the mOSaic scheme (Fig. 2). For both MODIS and CLM 2 land surface categories, snow and ice cover is estimated from input meteorology, while $f_L^{\text{water}}$ is defined as $\mathrm{1}-{\sum }_k{f}_L^k$. From the MODIS category “barren or sparsely vegetated”, everything poleward from 60^∘ is defined as tundra, while everything equatorward is categorized as desert. This mapping differs from the one used in the Wesely scheme.

Figure 2 Mapping of land surface categories. Either land surface categories from ISLSCP2 product of MODIS or the Community Land Model (CLM) 2 can be chosen for mapping to the land surface types we use in the mOSaic scheme. Water bodies of MODIS are actually not mapped. For both MODIS and CLM 2 land surface categories, snow and ice cover is estimated from input meteorology, while water is defined as $\mathrm{1}-{\sum }_k{f}_L^k$. From the MODIS category “barren or sparsely vegetated”, everything poleward from 60^∘ is defined as tundra, while everything equatorward is categorized as desert.

2.2 Pre-processing

As mentioned in the previous section, there are two variables needed for computing the stomatal conductance which are not directly available from the meteorological input data: the greening season, as the time of the year in the mid- and high latitudes when it is most likely for plants to grow, and the photosynthetic photon flux density, as the amount of light that plants need to photosynthesize. In the following, we present the necessary pre-processing of the variables. It is planned to implement an online computation of these variables into Oslo CTM3 later on.

2.2.1 Greening season

In Eqs. (20)–(23), Simpson et al. (2012) use prescribed start of growing season (SGS) and end of growing season (EGS) at 50^∘ N (d_SGS, d_EGS), together with lapse rates (∇d_SGS, ∇d_EGS) to define phenology and dry deposition over agricultural areas. For the growing season of crops in the computation of non-stomatal conductance, we use also prescribed values (Table 1), while for the stomatal conductance, as shown in Eq. (18), we use the SGS- and EGS-derived parameters: GDAY, the time elapsed starting at the SGS, and the total length of the greening season (GLEN), the time span between EGS and SGS. Since the parameterization of SGS and EGS in Simpson et al. (2012) is not applicable in a global model, another latitude-dependent parameterization is needed. First, we used a parameterization which was already implemented in Oslo CTM3 and which had been adopted from the Sparse Matrix Operational Kernel Emissions – Biogenic Emission Inventory System (SMOKE-BEIS; https://www.epa.gov/air-emissions-modeling/biogenic-emission-inventory-system-beis, last access: 24 October 2018). SMOKE-BEIS has fixed values for SGS and EGS for all regions but NH midlatitudes (23^∘ < lat < 65^∘), where it uses lapse rates of ∇d_SGS=4.5 and ∇d_EGS=3.3. As this parameterization is optimized for North America, it does not work well in Europe, e.g., most of northern Scandinavia has no allocated vegetation period. This basically results in a suppression of canopy resistance in northern Scandinavia.

In agriculture, there are different empirical rules to estimate the SGS and EGS. The simplest assumption is that greening starts after 5 consecutive days with a daily average temperature above 5 ^∘C, and vice versa for EGS. Other estimates use growing degree days (Levis and Bonan, 2004; Fu et al., 2014a), include soil moisture (Fu et al., 2014b), or rely on satellite observations. A comprehensive evaluation of different techniques is given by Anav et al. (2017). Another solution would be the usage of a proper land surface model, e.g., LPJ-GUESS, CLM, but the integration of such models into Oslo CTM3 is not planned at the moment.

Based on the empirical rule (5 ^∘C days), we have pre-processed our meteorological input data offline. We added some additional criteria to prevent “false spring”: if, within these 5 d, the average temperature drops below or rises above 5 ^∘C, the counter is reset, respectively. First, we used the 5 ^∘C day criteria for 45^∘ < lat < 85^∘ in the NH but extended them also to 35^∘ < lat < 65^∘ in the SH. In all other cases and where the 5 ^∘C day criteria fail, we still use the SMOKE-BEIS parameterization. The described algorithm written in Python 2.7 has been included in Sect. S4. An example map of the computed GLEN using the 5 ^∘C day criteria in both hemispheres is shown in Fig. 3.

Figure 3 Pre-processing of greening season from meteorological surface temperature fields. Shown is the total length of the greening season (GLEN) for the year 2005. The 5 ^∘C day criteria have been used in both hemispheres' mid–high latitudes. Ocean has been shaded to indicate that greening season will only affect land.

2.2.2 Photosynthetic photon flux density

From OpenIFS an accumulated surface PAR is available. It is integrated both spectrally (presumably 400–700 nm) and temporally. For practical use in Eq. (19), we de-accumulate this field with respect to time and refer to the result as PPFD.

The main obstacle is that PAR has been accumulated since model start, so that the first field kept from the original OpenIFS simulation (00:00 UTC) is 12 h after model start (12:00 UTC on the previous day). In other words, the first time step of each day in Oslo CTM3 has already accumulated PAR from 12:00 UTC on the previous day. De-accumulation of time 03:00 to 21:00 UTC simply means computing the difference:

$$\begin{array}{}\text{(27)}& \text{PPFD}\left(t_i\right)=\text{PAR}\left(t_{i+\mathrm{1}}\right)-\text{PAR}\left(t_i\right).\end{array}$$

For de-accumulation of the remaining time step, the best choice is subtracting the difference between 21:00 and 12:00 UTC of the previous day:

$$\begin{array}{}\text{(28)}& \begin{array}{rl}\text{PPFD}(t=\text{00:00}\mathrm{UTC})& =\text{PAR}(t=\text{00:00}\mathrm{UTC})\\ & -\left[\text{PAR}(t=\text{21:00}\mathrm{UTC}-\mathrm{1}\mathrm{d})\right.\\ & \left.-\text{PAR}(t=\text{12:00}\mathrm{UTC}-\mathrm{1}\mathrm{d})\right]\end{array},\end{array}$$

and limit the result to positive values only. An example PAR de-accumulation for 2 January 2005 is shown in Sect. S5 (Figs. S5–S7). The resulting PPFD fields are still accumulated over a time period of 3 h and should be divided by 3. A known issue (see https://confluence.ecmwf.int/display/CKB/ERA-Interim+known+issues, last access: 29 October 2019) in the OpenIFS (cycles ≤ c41r2) causes surface PAR values to be about 30 % below observations. To counter this, we decided to refrain from the division at this stage but need to bear this in mind for later OpenIFS cycles.

3 Evaluation

In this section, we present results from manifold Oslo CTM3 model integrations testing different parameters of the mOSaic scheme. We focus on changes in ozone total dry deposition $\sum {\mathrm{O}}_{\mathrm{3}}^{\text{DD}}$, dry deposition velocities $v_{\text{DD}}^{{\mathrm{O}}_{\mathrm{3}}}$, concentrations in the lowermost model level [O₃](p₀), and tropospheric burden ∑_tropO₃. We evaluate our results with respect to the multi-model comparison of ozone dry deposition by Hardacre et al. (2015) (Sect. 3.2), the MACC reanalysis (Sect. 3.3), and observations (Sect. 3.4). Oslo CTM3 is driven by meteorological input fields from ECMWF – OpenIFS cy38r1 (https://www.ecmwf.int/en/forecasts/documentation-and-support/evolution-ifs/cycle-38r1-summary-changes, last access: 4 October 2018). The CEDS historical emission inventory (Hoesly et al., 2018) is used for anthropogenic emissions, while biomass burning is covered in daily resolution by NASA's Global Fire Emissions Database version 4 (GFEDv4, Randerson et al., 2018). Biogenic emissions are taken from MEGAN-MACC output (Sindelarova et al., 2014), while emissions from soil and wetlands are computed by MEGAN. Resultant NO_x emissions are upscaled to match the Global Emissions InitiAtive (GEIA) inventory and are estimated to amount to 6.55 Tg (N) a⁻¹. For oceanic emissions of CO, we use predefined global fields from POET (GEIA-ACCENT emission data portal, 2003). Emissions of CH₄ are taken from the EU project (EU GOCE 037048) “Hydrogen, Methane and Nitrous oxide: Trend variability, budgets and interactions with the biosphere (HYMN)” for the year 2003 and scaled to oceanic amounts of CH₄ from NASA. In the following (Sect. 3.1), we will present the various model sensitivity studies.

3.1 Sensitivity studies

Due to significant differences between the mOSaic scheme and the previous Wesely scheme with respect to implementation, it is not possible to fully disentangle and trace back every single difference in results to a respective change. Therefore, we conducted one reference simulation denoted as mOSaic and in total seven sensitivity studies to probe the parameter space for stomatal conductance (mOSaic_offLight, mOSaic_offPhen, and mOSaic_SWVL1), ozone surface resistance $R^{{\mathrm{O}}_{\mathrm{3}}}$ (mOSaic_ice, mOSaic_desert, and mOSaic_hough), and emissions (mOSaic_emis2014). A reference simulation featuring the Oslo CTM3 Wesely scheme has been conducted and will be referred to as Wesely_type, indicating that other implementations of the original work by Wesely (1989) may exist in other models. All model experiments discussed in the following are summarized in Table 2. An “x” therein denotes that the model was run exactly in the configuration and with parameters as has been described in Sect. 2. For all model integrations, the meteorological reference year is 2005. This choice affects the direct comparison with data and studies that either show results based on decadal averages or differing years, because non-linearities in ozone formation and destruction make ozone concentrations sensitive to both differences in local concentration of precursors and meteorological conditions (Jin et al., 2013).

First, we take a closer look at the influence of certain parameters on the stomatal conductance. As indicated by the names, mOSaic_offLight and mOSaic_offPhen are rather extreme scenarios completely switching off the sensitivity to light and phenology in Eq. (12) by setting f_light and f_phen to a fixed value of 1, respectively. Because of the underlying research project's focus on Arctic and alpine ecosystems, where water might only be available from upper soil layers, an experiment was conducted using the uppermost soil water level (SWVL1) in the implementation of f_SW. After this, we want to confirm the importance of choice of $R^{{\mathrm{O}}_{\mathrm{3}}}$ for different land surface types. We conducted three experiments looking at a $R_{\text{ice/snow}}^{{\mathrm{O}}_{\mathrm{3}}}$ update (Helmig et al., 2007) (mOSaic_ice), observed $R_{\text{desert}}^{{\mathrm{O}}_{\mathrm{3}}}$ (Güsten et al., 1996) (mOSaic_desert), and an approximation of $R^{{\mathrm{O}}_{\mathrm{3}}}$ originally used in Wesely_type (Wesely, 1989; Hough, 1991). Finally, we run a simulation with emissions for the year 2014 instead of 2005 (EMEP_emis2014) to characterize the general influence of differing emissions on ozone.

Wesely (1989); Hough (1991)Simpson et al. (2012)Simpson et al. (2012)Simpson et al. (2012)Simpson et al. (2012)Simpson et al. (2012)Helmig et al. (2007)Simpson et al. (2012)Güsten et al. (1996)Simpson et al. (2012)Wesely (1989)Hough (1991)

Table 2 Summary of specifications of all simulations discussed in this section. For simplicity, only the tested parameters are listed. An “x” denotes that the model was run exactly in the configuration as has been described in Sect. 2. Here, “n/a” means that these parameters are not applicable in this experiment for it has been conducted with Wesely scheme.

^a $R_{\text{ice/snow}}^{{\mathrm{O}}_{\mathrm{3}}}=\mathrm{10}\mathrm{000}$ s m⁻¹. ^b $R_{\text{desert}}^{{\mathrm{O}}_{\mathrm{3}}}=\mathrm{800}$ s m⁻¹. ^c For adapted values, see Sect. S6.

Figure 4 Relative difference between reference simulations mOSaic and the Wesely_type with respect to (a) average surface ozone; (b) average ozone dry deposition velocity; (c) total amount of ozone removed from the atmosphere by dry deposition.

In Fig. 4, we show global distributions of the relative difference between mOSaic and Wesely_type for surface ozone, dry deposition velocity, and total ozone dry deposition. The surface ozone increases globally except for some regions covered by tropical forest. Especially in desert regions in Africa, North America, and Asia, the surface ozone increases by more than 100 %. Consistently, dry deposition velocities decrease globally by the same order of magnitude in these regions, while they increase over tropical forest. With respect to total dry deposition, the picture is a bit less clear. We find a decrease of total dry deposition of ozone in desert regions and ocean-covered areas and an increase in regions covered by tropical forest, while at mid- and high latitudes in both hemispheres only small changes are visible. A possible explanation for this divergence especially in desert regions is the difference between the prescribed surface resistance $R^{{\mathrm{O}}_{\mathrm{3}}}$ in the Wesely scheme in comparison to those used in mOSaic. We come back to this in the following sections.

3.2 Comparison with modeling results

In the evaluation of our model, we closely follow suggestions by Hardacre et al. (2015). For the purpose of comparison with the multi-model mean of the Task Force on Hemispheric Transport of Air Pollution (TF HTAP) models, we also have regridded our data to a horizontal resolution of $\mathrm{3}{}^{\circ }\times \mathrm{3}{}^{\circ }$. In Sect. 3.2.1, we look at zonal distributions of [O₃](p₀), $v_{\text{DD}}^{{\mathrm{O}}_{\mathrm{3}}}$, and $\sum {\mathrm{O}}_{\mathrm{3}}^{\text{DD}}$ for all our sensitivity simulations, and study seasonal cycles of hemispheric ozone, as well as for nine land surface types (Sect. 3.2.2). From this, we estimate the total annual ozone dry deposition onto ocean, ice, and land surfaces and compare this also with results from Luhar et al. (2017).

Dry deposition velocities are directly available only for the new model version. For Wesely_type, monthly averaged dry deposition velocities $v_{\text{DD}}^{{\mathrm{O}}_{\mathrm{3}}}$ had to be retrospectively estimated from the ratio between the total ozone dry deposition $\sum {\mathrm{O}}_{\mathrm{3}}^{\text{DD}}\left(p_{\mathrm{0}}\right)$ and monthly averaged ozone amount in the lowermost model level O₃(p₀):

$$\begin{array}{}\text{(29)}& v_{\text{DD}}^{{\mathrm{O}}_{\mathrm{3}}}={\displaystyle \frac{\sum {\mathrm{O}}_{\mathrm{3}}^{\text{DD}}\left(p_{\mathrm{0}}\right)}{{\mathrm{O}}_{\mathrm{3}}\left(p_{\mathrm{0}}\right)}}\cdot c_{\text{month}}.\end{array}$$

Herein, $c_{\text{month}}=\frac{\mathrm{\Delta }{h}_{\text{month}}}{s_{\text{month}}}$, with the monthly average height of the lowermost model level in each grid box Δh_month and the respective number of seconds in a month s_month. In the case of mOSaic, resulting values for $v_{\text{DD}}^{{\mathrm{O}}_{\mathrm{3}}}$ from Eq. (29) are compatible with the values which are directly available from model output.

Figure 5 Comparison of the manifold Oslo CTM3 integrations with respect to (a) ozone concentrations in the lowermost model level, (b) annual average ozone dry deposition velocity, and (c) total annual ozone dry deposition. The different colors indicate sets of simulation with similar baselines. The multi-model mean from the evaluation of TF HTAP models by Hardacre et al. (2015) is shown as a reference (where available).

3.2.1 Zonal distribution

The annual zonal average with respect to surface ozone concentration (Fig. 5a) displays on average, consistent with Fig. 4a, a global increase of surface ozone concentrations by 6 ppb comparing mOSaic to Wesely_type. This increase is largest in the zonal band 25–50^∘ N, which contains the major deserts. In the deep tropics (5^∘ S–5^∘ N), the increase is smallest (𝒪(5 ppb)). We find that the mOSaic scheme further intensifies the strong asymmetry between the Northern Hemisphere and Southern Hemisphere as a consequence of the distribution of the continental land masses and vegetation thereon. Among the sensitivity studies focusing on the stomatal conductance, there is only a low absolute variance. Neglecting the dependence on light in the stomatal conductance formulation (mOSaic_offLight) – or in other words allowing photosynthesis 24/7 – decreases the ozone concentration by 1–2 ppb in the tropics and NH midlatitudes, while choosing soil water at shallower depths (mOSaic_SWVL1) increases [O₃] insignificantly. Rather surprisingly, switching off the phenology completely (mOSaic_offPhen) amounts on average only to a small difference (𝒪(<1 ppb)). Most remarkably, but expected due to the much smaller prescribed dry deposition velocity over ice and snow, mOSaic_ice displays a doubling of surface ozone in the high Arctics compared to Wesely_type (𝒪(20 ppb)) but affects ozone concentrations down to latitudes at about 50^∘ in both hemispheres. Reducing $R_{\text{desert}}^{{\mathrm{O}}_{\mathrm{3}}}$ by 60 % (mOSaic_desert), a reduction on the order of 1 ppb is found mainly limited to the NH. The largest impact on ozone concentrations (𝒪(2–5 ppb)) is found for the experiment mOSaic_hough which is closest to Wesely_type, since we used on average the same $R^{{\mathrm{O}}_{\mathrm{3}}}$ (see Sect. S6). The scenario of differing emissions (2005 in comparison to 2014 or more specifically mOSaic compared to mOSaic_emis2014) yields higher ozone concentrations in the Northern Hemisphere in 2005 in accordance with a reduction in sulfur and NO_x emissions in southeast Asia in later years. An opposite tendency is seen for latitudes south of 30^∘ N, where an increase in ozone precursors is seen in CEDS.

The $v_{\text{DD}}^{{\mathrm{O}}_{\mathrm{3}}}$ are shown in Fig. 5b. The dry deposition velocities in the mOSaic scheme are well below the Wesely scheme and in remarkable agreement with the results shown by Hardacre et al. (2015). In the Arctics, except for mOSaic_ice, all model experiments are slightly above the multi-model mean. This indicates that, with respect to the other models, the Helmig et al. (2007) surface resistance above ice and snow should be considered as the new standard for Oslo CTM3. This may, however, lead to an overcompensation of the current Arctic low bias in surface ozone in Oslo CTM3 and needs further evaluation. The dry deposition velocities are of course independent of the emission scenario but display a strong sensitivity to f_light, f_phen, and especially the choice of $R^{{\mathrm{O}}_{\mathrm{3}}}$. The shape of the normalized zonal average dry deposition velocities of the mOSaic scheme are more similar to the multi-model mean than to Wesely_type (Sect. S7; Fig. S8). The biggest exceptions are the zonal bands 50–70^∘ S (almost entirely covered by ocean), 12–30^∘ S (coinciding with the location of Australia and its desert regions), as well as its counterpart in the Northern Hemisphere (12–30^∘ N).

The annual total ozone dry deposition is shown in Fig. 5c. In accordance with the previously described features, we observe a reduction of the global total ozone dry deposition in all sensitivity studies. In the most extreme case (NH subtropics and midlatitudes), the total ozone dry deposition drops to one-half of the amount given by Wesely_type. The occurrence of this reduction in the zonal bands, where the major deserts are located, points to a substantial difference in $v_{\text{desert}}^{{\mathrm{O}}_{\mathrm{3}}}$. Consulting the parameter file used in the Wesely scheme, we indeed find $v_{\text{desert}}^{{\mathrm{O}}_{\mathrm{3}}}\equiv v_{\text{tundra}}^{{\mathrm{O}}_{\mathrm{3}}}=\mathrm{0.26}$ cm s⁻¹ (Hough, 1991), while in the mOSaic scheme $v_{\text{desert}}^{{\mathrm{O}}_{\mathrm{3}}}=\mathrm{0.05}$ cm s⁻¹ and $v_{\text{tundra}}^{{\mathrm{O}}_{\mathrm{3}}}=\mathrm{0.24}$ cm s⁻¹, respectively. Similarly, dry deposition velocities over ice and snow and ocean have been even higher in the Wesely scheme ($v_{\text{ice/snow}}^{{\mathrm{O}}_{\mathrm{3}}}\equiv v_{\text{water}}^{{\mathrm{O}}_{\mathrm{3}}}=\mathrm{0.07}$ cm s⁻¹) than in the original parameter set ($v_{\text{ice/snow}}^{{\mathrm{O}}_{\mathrm{3}}}\equiv v_{\text{water}}^{{\mathrm{O}}_{\mathrm{3}}}=\mathrm{0.05}$ cm s⁻¹, Simpson et al., 2012). These differences in surface resistance over huge parts of the unvegetated surface of the Earth account for most of the qualitative difference between the Wesely and the mOSaic scheme but do not explain the quantitative difference (compare mOSaic_hough). We further elaborate on this in the following (Sect. 3.2.2).

There seems to be a discrepancy between the Oslo CTM3 response and the multi-model mean, since the Wesely scheme is similar to the multi-model mean with respect to total annual ozone dry deposition, while the $v_{\text{DD}}^{{\mathrm{O}}_{\mathrm{3}}}$ of the mOSaic scheme matches better. This could be a sign of differences in photochemistry and transport (e.g., convective, advective, STE) between Oslo CTM3 and the average TF HTAP model but without comparing to the actual [O₃] of the TF HTAP models that participated in the model intercomparison, we cannot elaborate on this any further. This may also hint at issues in the Oslo CTM3 photochemistry, which may have a too-high ozone production, or the actual removal of ozone from the atmosphere, which might have been adjusted to the less physical dry deposition velocities in the past, but this is subject to further investigations.

In Fig. A1a in the Appendix, the average zonal ozone dry deposition is shown separated by month. Where available, we have added the multi-model mean given by Hardacre et al. (2015) as a reference. As for the global annual comparisons above, the mOSaic scheme matches the multi-model-mean values remarkably well with respect to dry deposition velocities, while it strongly underestimates the total dry deposition. Qualitatively, there are two major phases apparent: NH and SH greening seasons. Spring and summer in the NH are reflected in a pronounced peak of $v_{\text{DD}}^{{\mathrm{O}}_{\mathrm{3}}}$ in the northern midlatitudes, while it is absent in winter (SH summer). Spring and summer in the SH are marked by a southward shift of the tropical peak dry deposition velocity and a slight increase of $v_{\text{DD}}^{{\mathrm{O}}_{\mathrm{3}}}$ in the region 20–40^∘ S. In the Wesely scheme, NH midlatitude peak velocities appear in June compared to July in the mOSaic scheme, indicating that the seasonal cycles differ. The corresponding total monthly ozone dry deposition is shown in Fig. A1b. In general, the seasonal patterns are quite similar in the Wesely scheme and the mOSaic scheme, displaying a strong symmetry around 10^∘ N in January/February and November/December, respectively. What differs most is the molding and intensity of the NH peak dry deposition. Both schemes reach the maximum in June/July but the peak is much more differentiated in March already in the Wesely scheme. Similarly, the SH tropical peak dry deposition is reached in August/September but sustained longer, into October, in the Wesely scheme. Since we have not conducted any simulation with a meteorological year other than 2005, we cannot elaborate on whether this is a special feature of our chosen year or not.

3.2.2 Average seasonal cycles

To further disentangle the contributions of different regions to the global ozone budget, we will look at different projections of seasonal cycles.

In Fig. 6, the total annual ozone dry deposition separated into mid- and high latitudes in the Northern Hemisphere (30–90^∘ N), the tropics and subtropics (30^∘ S–30^∘ N), and the mid- and high latitudes in the Southern Hemisphere (30–90^∘ S) is shown. We have added the multi-model mean by Hardacre et al. (2015) as a reference. While the total ozone dry deposition of Wesely_type agrees well with the multi-model mean in any zonal band, the mOSaic scheme displays a much smaller total ozone dry deposition. This deviation appears to be almost the same for each zonal band (6 %–7 %).

Figure 6 Seasonal cycle of total annual amount of ozone removed from the atmosphere through dry deposition separated into the NH, tropics (TR), and SH. The multi-model mean from the evaluation of HTAP models by Hardacre et al. (2015) is shown as a reference.

As expected, the NH mid- and high latitudes display a strongly pronounced seasonal cycle, while it is less pronounced in the tropics (due to the lack of seasons) and in the SH (due to the small percentage of vegetated surface). The highest ozone dry deposition is found in the tropics and amounts on average to the peak level of dry deposition in the NH for the multi-model mean (Hardacre et al., 2015) and mOSaic scheme. In the Wesely scheme, the average tropical ozone dry deposition diverges by 5 Tg in comparison to its corresponding NH maximum. Compared to the multi-model mean, the seasonal cycle in the Oslo CTM3 NH appears to be shifted towards later in the year. The seasonal cycle in the tropics and subtropics only differs by magnitude; otherwise, the shapes are identical for the mOSaic scheme, the Wesely scheme, and the multi-model mean. The total amount of dry deposition of ozone differs strongly between the different model experiments, with mOSaic_SWVL1 and mOSaic_hough displaying the lowest and highest amounts, respectively. This indicates that surface ozone is much more sensitive to the choice of parameters (𝒪(5 ppb) for mOSaic_hough in the tropics) than to slight changes in precursor emissions (𝒪(1 ppb) for mOSaic_emis2014 in the tropics).

As suggested by Hardacre et al. (2015), we also look at ozone dry deposition velocities with respect to surface types separately. Since dry deposition velocities are not directly available for Wesely_type, we use Eq. (29) to estimate these. Based on a CLM 2 average dynamic land surface map, we generate masks for nine different surface types (Fig. A2a) and use these to select grid boxes with a high percentage of these surface types, ranging from a meager 70 % for cropland in the NH midlatitudes to 100 % for desert, ocean, snow and ice, and tropical forest. Thus, it is not possible to exclusively select grid boxes with 100 % cover for each surface type. Since we have not performed a full unfolding on the data, the results should be treated with slight caution (e.g., over cropland). In the case of the mOSaic scheme, we have preselected the dry deposition velocities in accordance with the land surface type.

In Fig. 7, the seasonal cycles of dry deposition velocities are shown for the nine surface categories. The patterns and absolute numbers differ substantially between the Wesely scheme and the mOSaic scheme and the multi-model mean. The divergence of the average dry deposition velocities between Wesely_type and mOSaic in desert regions ($\mathrm{\Delta }{\stackrel{\mathrm{‾}}{v}}_{\text{desert}}^{{\mathrm{O}}_{\mathrm{3}}}=\mathrm{0.20}$ cm s⁻¹) as well as grassland ($\mathrm{\Delta }{\stackrel{\mathrm{‾}}{v}}_{\text{grassland}}^{{\mathrm{O}}_{\mathrm{3}}}=\mathrm{0.65}$ cm s⁻¹) is quite remarkable. The difference of mOSaic and Wesely_type from the multi-model mean in tropical forest regions is ${\mathrm{\Delta }}_{\text{mOSaic}}{\stackrel{\mathrm{‾}}{v}}_{\text{tropical forest}}^{{\mathrm{O}}_{\mathrm{3}}}=\mathrm{0.61}$ cm s⁻¹ and ${\mathrm{\Delta }}_{\text{Wesely\_type}}{\stackrel{\mathrm{‾}}{v}}_{\text{tropical forest}}^{{\mathrm{O}}_{\mathrm{3}}}=\mathrm{0.49}$ cm s⁻¹, respectively. The multi-model mean displays a rather pronounced seasonal cycle in desert regions ($\mathrm{0.10}\mathrm{cm}{\mathrm{s}}^{-\mathrm{1}}\le v_{\text{desert}}^{{\mathrm{O}}_{\mathrm{3}}}\le \mathrm{0.15}$ cm s⁻¹), which cannot be reproduced with the mOSaic scheme. The dry deposition velocities over desert regions are consistent with the average values from the prescribed ozone surface resistance, which means that in the mOSaic scheme they are 1 order of magnitude lower than in the Wesely scheme. In the mOSaic scheme, dry deposition to deserts is dominated by contribution from R_b. From a limited number of ozone flux measurements in the Sahara, Güsten et al. (1996) deduced $v_{\text{desert, day}}^{{\mathrm{O}}_{\mathrm{3}}}=\mathrm{0.1}$ cm s⁻¹, $v_{\text{desert, night}}^{{\mathrm{O}}_{\mathrm{3}}}=\mathrm{0.04}$ cm s⁻¹, and ${\stackrel{\mathrm{‾}}{v}}_{\text{desert}}^{{\mathrm{O}}_{\mathrm{3}}}=\mathrm{0.065}$ cm s⁻¹. This implies that ozone dry deposition over desert regions is highly overestimated in the Wesely scheme as well as in TF HTAP models, while it may be underestimated in the mOSaic scheme. Similarly, the dry deposition velocities over water differ. From measurements during ship campaigns, a mean value of ${\stackrel{\mathrm{‾}}{v}}_{\text{water}}^{{\mathrm{O}}_{\mathrm{3}}}=\mathrm{0.019}$ cm s⁻¹ over the ocean has been deduced (Helmig et al., 2012). In a model study of different mechanisms of dry deposition to ocean waters by means of prescribed ${\stackrel{\mathrm{‾}}{v}}_{\text{water}}^{{\mathrm{O}}_{\mathrm{3}}}$ and one- and two-layer gas exchange modeling, Luhar et al. (2017) found ${\stackrel{\mathrm{‾}}{v}}_{\text{water}}^{{\mathrm{O}}_{\mathrm{3}}}$ ranging between 0.018 cm s⁻¹ (two-layer scheme) and 0.039 cm s⁻¹ (prescribed). With ${\stackrel{\mathrm{‾}}{v}}_{\text{water}}^{{\mathrm{O}}_{\mathrm{3}}}=(\mathrm{0.046}\pm \mathrm{0.002})$ cm s⁻¹, the mOSaic scheme (Sect. 2.1.2) yields probably a too-strong dry deposition to ocean but is in line with the multi-model mean. This implies that ozone concentrations might even become larger and dry deposition even lower in the model if a more advanced dry deposition scheme to the ocean would be implemented. With respect to vegetation, we might be able to improve the model performance further by allowing more PFTs and phenologies, especially in regions covered by tropical forest (Anav et al., 2017) or in boreal regions.

Figure 7 Average seasonal cycles of ozone dry deposition velocities separated by land use type. Results from Hardacre et al. (2015) are shown as a reference. We refrain from showing the full extent of EMEP_offLight here, since it is an extreme scenario and has been discussed already.

Finally, we take a look at the different global as well as hemispheric dry deposition sinks for ozone (Table 3). Despite its vastness, the ocean amounts only to 35 % of the global ozone sink due to dry deposition in Oslo CTM3 with the operative mOSaic scheme (mOSaic), while permanently ice- and snow-covered regions account for 1.2 %. The remainder is deposited to land surfaces of which deserts might be the most neglected in process modeling. The total annual dry deposition in the mOSaic scheme is one-third below the multi-model-mean result by Hardacre et al. (2015). But also the results of Luhar et al. (2017, 2018) yield a 19 %–27 % lower ozone dry deposition than the models participating in the model intercomparison, with deposition to ocean ranging between 12 % and 21 % of the total annual ozone dry deposition. In particular, Luhar et al. (2018) found that current model estimates of ozone dry deposition to the ocean may be 3 times too high compared to their analysis. This implies that the ozone dry deposition to the ocean in Oslo CTM3 is too high as well.

Table 3 Total ozone dry deposition for the respective model experiment in Tg a⁻¹. The global ozone dry deposition has been weighted by ocean, ice, and land fraction in each grid box, respectively. “Ice” herein refers to regions at high latitudes that are permanently covered by ice and snow.

^* Relative change of global annual total in comparison to mOSaic.

Table 4 displays the average tropospheric ozone burden for all model experiments. Consistent with the previous findings, the mOSaic scheme increases the tropospheric ozone burden by 35 Tg. From various satellite ozone retrieval products, Gaudel et al. (2018, Table 5) deduce a lower limit estimate for global tropospheric ozone burden for the years 2010–2014 of 333–345 Tg but remark that this amount underestimates the actual tropospheric ozone burden, since it is only based on daytime retrievals. Nevertheless, the results of mOSaic lie 17 % above that estimate and also well above the typical modeling range of 302–378 Tg (Young et al., 2013). Despite the strong positive bias in ozone concentrations and accordingly low bias in total dry deposition, the difference between mOSaic and mOSaic_emis2014 (6 Tg) lies well within the range given by satellites for the years 2005 and 2014 (Gaudel et al., 2018, Fig. 26). This indicates that Oslo CTM3 responds well to given changes in global emissions.

Table 4 Annual mean tropospheric ozone burden for all experiments and 1σ standard deviation.

3.3 Comparison with MACC reanalysis

In this section, we conclude the comparison of our results with respect to global ozone by looking at ECMWF's MACC reanalysis (MACC-II Consortium, 2011, data obtained from ECWMF's data center; https://apps.ecmwf.int/datasets/data/macc-reanalysis/levtype=ml/, last access: 8 July 2019). In Fig. 8, we compare mOSaic ozone concentrations in the lowermost model level with surface concentrations deduced from the MACC reanalysis for the year 2005. The MACC reanalysis displays low ozone concentrations above all land masses except for the Greenland ice sheet. The lowest values are found in the deep tropics (e.g., northern South America and central Africa), while the highest values occur within 25–60^∘ N over the oceans. These low values over the oceans are relatively well reproduced by mOSaic (±20 %). Over, e.g., South America, central and southern Africa, the Arabian Peninsula, and the northwestern Indian subcontinent, the ozone concentrations are elevated by up to 200 %. Regarding the global average surface ozone concentration, Wesely_type (Sect. S8) is more consistent with the MACC reanalysis than mOSaic but shows the same tendency of enhanced ozone over the continents as the latter. This enhancement is apparent mostly in the deep tropics and over the northwestern Indian subcontinent, which coincides with regions of high intensity in incoming UV radiation. This may indicate an imbalance in the photochemical production and loss of O₃ in Oslo CTM3.

Figure 8 Mean ozone concentrations for the year 2005. (a) MACC reanalysis (surface); (b) Oslo CTM3 mOSaic (lowermost model level); (c) relative difference.

3.4 Comparison with ground-based observations

In this section, we compare our model results to observations at a selected number of sites which provide ozone flux measurements. For all comparisons, we use the original resolution of Oslo CTM3 ($\mathrm{2.25}{}^{\circ }\times \mathrm{2.25}{}^{\circ }$) instead of the regridded resolution ($\mathrm{3}{}^{\circ }\times \mathrm{3}{}^{\circ }$).

Figure 9 Ozone dry deposition fluxes at different observation sites. (a) Comparison between Wesely scheme and average result from all sensitivity studies with observational averages taken from Hardacre et al. (2015). The model uncertainty of Oslo CTM3 is given as a 1 σ error band. The shaded area around the multi-model mean indicates the broad range of different model results but is not an actual error band. (b) Model divergence from observation and χ²-test results.

In Fig. 9a, seasonal cycles of average ozone dry deposition fluxes for the six selected observation sites are shown. We have computed a model average for all sensitivity studies at the closest grid point and show the 1 σ uncertainty band. The shaded area around the multi-model mean indicates the broad range of model results but is not an actual uncertainty band since such is not given in Hardacre et al. (2015). At four of the six sites, the mOSaic scheme performs better than the Wesely scheme and similar to or better than the multi-model mean. We use a χ² test:

$$\begin{array}{}\text{(30)}& {\mathit{\chi }}^{\mathrm{2}}=\sum _{i=\mathrm{1}}^{\mathrm{12}}{\displaystyle \frac{{\left({\stackrel{\mathrm{‾}}{{\mathrm{O}}_{\mathrm{3}}}}_{\text{DD},i}^{\text{sim}}-{\stackrel{\mathrm{‾}}{{\mathrm{O}}_{\mathrm{3}}}}_{\text{DD},i}^{\text{obs}}\right)}^{\mathrm{2}}}{{\mathit{\sigma }}_i^{\mathrm{2}}}},\end{array}$$

with an estimated standard deviation of observation σ_i=1 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ and divide it by the number of degrees of freedom (NDF) to assess this subjective analysis in a more objective way. The closer to 1 this test scores, the better the simulation represents the observation. A score between 0 and 1 indicates that the estimated σ is too small. The results of the χ² test are shown together with the divergences in Fig. 9b. The χ² test reveals also that in four of the six cases the mOSaic scheme improves the performance of Oslo CTM3 with respect to observed ozone dry deposition fluxes, although a satisfying result is only achieved for two sites (Castel Porziano, Blodgett Forest). With only one full year of simulation, the model uncertainty regarding the seasonal cycle at observational sites cannot be properly quantified. Furthermore, the observational averages comprise at most 9 years worth of data. Statistically, these data may still be subject to interannual variability. Among other aspects, the horizontal as well as vertical resolution play an important role in the model performance. Although we do not explicitly assess the impacts of differing resolutions in our model, we can assume that both high and low biases exist due to dilution of sources and sinks in coarse-resolution models (Schaap et al., 2015). Good matches between observation and model are only to be expected if the station's location is representative of an area similar to the respective model grid box and not substantially affected by differences in modeled and actual topography (e.g., major wind directions).

4 Summary and conclusions

We have presented an update of the dry deposition scheme in Oslo CTM3 from purely prescribed dry deposition velocities (Wesely, 1989; Hough, 1991) to a more process-oriented parameterization taking the state of the atmosphere and vegetation into account. Based on the description of dry deposition in Simpson et al. (2003, 2012), we have implemented a mosaic approach to compute contributions to dry deposition by individual subgrid land surface types. Aerodynamic, quasi-laminar, and surface resistance (with the latter divided into stomatal and non-stomatal contributions) is calculated for each land surface category separately. Based on these, a land fraction-weighted mean is deduced. In addition, the various dry deposition velocities are now directly available from model output for diagnostics and further studies.

The new dry deposition scheme named mOSaic improves the modeled dry deposition velocities which are now compatible with observation and model studies (e.g., Hardacre et al., 2015; Luhar et al., 2018). Dry deposition velocities are reduced by 6 %–60 %. At the same time, the annual amount of ozone dry deposition decreases by more than 100 % over all major desert areas and increases over tropical forest. Compared to results from a multi-model evaluation (Hardacre et al., 2015), the total annual ozone dry deposition in Oslo CTM3 is ${\mathrm{38}}_{-\mathrm{10}}^{+\mathrm{1}}$ % below average. However, there seems to be a tendency that newer TF HTAP models show a lower total annual dry deposition of ozone than older models, indicating that newer developments lead to decreasing ozone dry deposition and increasing tropospheric ozone burden (e.g., Luhar et al., 2017, 2018; Hu et al., 2017).

We found the response of Oslo CTM3 to the changes in dry deposition velocities from the old and the mOSaic scheme to be consistent. A decrease in $v_{\text{DD}}^{{\mathrm{O}}_{\mathrm{3}}}$ leads to a decrease in total ozone dry deposition and an increase in ozone concentration [O₃]. As the new scheme is quantitatively more similar to the multi-model mean (Hardacre et al., 2015) with respect to dry deposition velocities, while the old scheme agrees better in terms of total dry deposition, there is an apparent discrepancy. By means of tropospheric ozone burden (Gaudel et al., 2018) and surface ozone concentrations deduced from the MACC reanalysis (MACC-II Consortium, 2011), Oslo CTM3 with the operational mOSaic scheme shows a pronounced high bias of tropospheric ozone. While the average bias is small or even reversed when using the old scheme, both display elevated ozone in comparison with the MACC reanalysis in continental regions with high average incoming UV radiation (e.g., northern South America, central and southern Africa, the Himalayas). The reason behind this bias has not yet been resolved and may hint at, e.g., issues in photochemistry ([OH]-related ozone production and loss) or previously introduced optimization of ozone removal with respect to the old, less physical dry deposition velocities.

Most of the qualitative change in ozone dry deposition in Oslo CTM3 (−2 %–12 %) can be attributed to changes in dry deposition over the ocean and deserts. This is mainly due to updates of the respective prescribed ozone surface resistance $R^{{\mathrm{O}}_{\mathrm{3}}}$. In the case of desert and grasslands, the difference between the old and new prescribed values is an order of magnitude. Over the ocean, the absolute change in dry deposition is small, but it is accumulated over a large area which is especially amplified in the Southern Hemisphere. Small adjustments to the lower limits in our quasi-laminar layer resistance formulation may help improve Oslo CTM3 performance in this regard. With respect to available measurements of dry deposition velocities of ozone over desert (Güsten et al., 1996) and ocean (Helmig et al., 2012), the new Oslo CTM3 dry deposition scheme slightly underestimates ozone dry deposition velocities over the former and overestimates them over the latter. Regarding the vastness of the ocean and the ongoing desertification, it may be worthwhile to revise the dry deposition scheme for these regimes and add more process-oriented formulations, e.g., two-layer gas exchange with ocean waters (Luhar et al., 2017, 2018), wave breaking and spray (Pozzer et al., 2006).

Although dry deposition to ice and snow amounts to only 1 % of the total global annual ozone dry deposition in mOSaic, a decrease in prescribed dry deposition velocity in accordance with combined measurements and model studies (Helmig et al., 2007) causes almost a doubling in the surface ozone concentrations in the high Arctic and affects surface ozone concentrations down to latitudes at 50^∘ in both hemispheres. Comparing with results from the multi-model evaluation (Hardacre et al., 2015), we conclude that it is important to use this updated ozone dry deposition velocity to counter an Arctic surface ozone low bias in models; however, this currently leads to an overcompensation (high bias) in Oslo CTM3.

We have studied the parameter space of the stomatal conductance parameterization and found that surface ozone in the tropics and the Northern Hemisphere is most sensitive to changes therein. In the most extreme test case, the increase in global total dry deposition amounts to 7.3 %, while the more realistic test cases, e.g., using differing years of emission, amount to changes of the order of ±2 %. This may indicate that future changes in vegetation cover and solar radiation at the surface due to changes in stratospheric ozone, cloud cover, or aerosols could also strongly influence the surface ozone burden in the tropics. Total column ozone in the tropics is predicted to decrease due to changes in the atmospheric circulation (e.g., WMO – Global Ozone Research and Monitoring Project, 2014), while tropospheric and surface ozone increase. The combined effects of increasing emissions of ozone precursors and an increase in UV due to thinning of stratospheric ozone might permit more UV light at ground level and thus increase the ozone production.

An important factor in the global ozone budget is emissions of precursor substances. We cover this by using the same meteorology with different years of CEDS emissions. We chose the years 2005 and 2014 for our comparison. Ozone precursor emissions in 2014 are slightly lower in the NH while enhanced in the tropics and the SH. In 2014, surface ozone burden is higher in the Southern Hemisphere and in the tropics (5 %) compared to 2005, while it is lower in the Northern Hemisphere (2 %).

We also evaluated the model with respect to observed dry deposition fluxes at six sites in the Northern Hemisphere and found that the mOSaic scheme performs better than the old one but is not able to reproduce the measurements at most sites quantitatively. This may be due to several reasons. The model resolution in both horizontal ($\mathrm{2.25}{}^{\circ }\times \mathrm{2.25}{}^{\circ }$) and vertical (L60, P_max=0.02 hPa) does not capture all details in transport, thus affecting the distribution and transport (e.g., long-range, convection, and stratosphere–troposphere exchange) of ozone and its precursors. Depending on the location of the observation site and its respective representativeness for a larger area, ozone dry deposition and ozone concentrations are expected to be over- or underestimated in the model. Because of non-linearities in ozone formation and destruction, ozone concentrations are sensitive to both differences in local concentration of precursors and meteorological conditions (Jin et al., 2013). In addition, a comparison of very few years of measurement to only one specific year of simulation may reflect the year-to-year variability more than the actual model performance.

Future work on Oslo CTM3 should resolve the ozone high bias which may involve revising the photolysis and chemical reaction computation as well as reaction rates. For a better modeling of ozone abundances, ocean emissions of very-short-lived ozone-depleting substances (VSLSs) (Warwick et al., 2006; Ziska et al., 2013) which affect stratospheric ozone (Hossaini et al., 2016; Falk et al., 2017) and a scheme covering Arctic springtime ozone depletion (e.g., Yang et al., 2010; Toyota et al., 2011; Falk and Sinnhuber, 2018) could be worthwhile to implement. The general model performance could also be improved by allowing for more plant functional types and phenologies than currently used or implementing an actual photosynthesis-based modeling of plants. A more efficient parallelization of the code would enable computation at higher horizontal resolutions.

Appendix A: Figures

Figure A1 Comparison of the manifold Oslo CTM3 integrations with respect to (a) zonal average ozone dry deposition velocities and (b) total annual amount of ozone removed from the atmosphere via dry deposition. The multi-model mean from the evaluation of TF HTAP models by Hardacre et al. (2015) is shown as a reference (where available).

Figure A2 Partitioning of land surface types: (a) CLM 2 dynamic land surface types in $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ resolution; (b) zonal distribution of land surface types.