Importance of dry deposition parameterization choice in global simulations of surface ozone

. Dry deposition is a major sink of tropospheric ozone. Increasing evidence has shown that ozone dry deposition actively links meteorology and hydrology with ozone air quality. However, there is little systematic investigation on the performance of different ozone dry deposition parameter-izations at the global scale and how parameterization choice can impact surface ozone simulations. Here, we present the results of the ﬁrst global, multidecadal modelling and evaluation of ozone dry deposition velocity ( v d ) using multiple ozone dry deposition parameterizations. We model ozone dry deposition velocities over 1982–2011 using four ozone dry deposition parameterizations that are representative of current approaches in global ozone dry deposition modelling. We use consistent assimilated meteorology, land cover, and satellite-derived leaf area index (LAI) across all four, such that the differences in simulated v d are entirely due to differences in deposition model structures or assumptions about how land types are treated in each. In addition, we use the surface ozone sensitivity to v d predicted by a chemical transport model to estimate the impact of mean

tropical rainforests and seasonally dry tropical forests in Indochina in December.Parameterization-specific biases based on individual land cover type and hydroclimate are found to be the two main drivers of such discrepancies.We find statistically significant trends in the multiannual time series of simulated July daytime v d in all parameterizations, driven by warming and drying (southern Amazonia, southern African savannah, and Mongolia) or greening (high latitudes).The trend in July daytime v d is estimated to be 1 % yr −1 and leads to up to 3 ppbv of surface ozone changes over 1982-2011.The interannual coefficient of variation (CV) of July daytime mean v d in NH is found to be 5 %-15 %, with spatial distribution that varies with the dry deposition parameterization.Our sensitivity simulations suggest this can contribute between 0.5 to 2 ppbv to interannual variability (IAV) in surface ozone, but all models tend to underestimate interannual CV when compared to long-term ozone flux observations.We also find that IAV in some dry deposition parameterizations is more sensitive to LAI, while in others it is more sensitive to climate.Comparisons with other published estimates of the IAV of background ozone confirm that ozone dry deposition can be an important part of natural surface ozone variability.Our results demonstrate the importance of ozone dry deposition parameterization choice on surface ozone modelling and the impact of IAV of v d on surface ozone, thus making a strong case for further measurement, evaluation, and model-data integration of ozone dry deposition on different spatiotemporal scales.

Introduction
Surface ozone (O 3 ) is one of the major air pollutants that poses serious threats to human health (Jerrett et al., 2009) and plant productivity (Ainsworth et al., 2012;Reich, 1987;Wittig et al., 2007).Ozone exerts additional pressure on global food security and public health by damaging agricultural ecosystems and reducing crop yields (Avnery et al., 2011;McGrath et al., 2015;Tai et al., 2014).Dry deposition, by which atmospheric constituents are removed from the atmosphere and transferred to the Earth's surface through turbulent transport or gravitational settling, is the second-largest and terminal sink of tropospheric O 3 (Wild, 2007).Terrestrial ecosystems are particularly efficient at removing O 3 via dry deposition through stomatal uptake and other non-stomatal pathways (Wesely and Hicks, 2000) (e.g.cuticle, soil, reaction with biogenic volatile organic compounds (BVOCs); Fares et al., 2010;Wolfe et al., 2011).Meanwhile, stomatal uptake of O 3 inflicts damage on plants by initiating reactions that impair their photosynthetic and stomatal regulatory capacity (Hoshika et al., 2014;Lombardozzi et al., 2012;Reich, 1987).Widespread plant damage has the potential to alter the global water cycle (Lombardozzi et al., 2015) and suppress the land carbon sink (Sitch et al., 2007), as well as generate a cascade of feedbacks that affect atmospheric composition including ozone itself (Sadiq et al., 2017;Zhou et al., 2018).Ozone dry deposition is therefore key in understanding how meteorology (Kavassalis and Murphy, 2017), climate, and land cover change (Fu and Tai, 2015;Ganzeveld et al., 2010;Geddes et al., 2016;Heald and Geddes, 2016;Sadiq et al., 2017;Sanderson et al., 2007;Young et al., 2013) can affect air quality and atmospheric chemistry at large.
Analogous to other surface-atmosphere exchange processes (e.g.sensible and latent heat flux), O 3 dry deposition flux (F O 3 ) is often expressed as the product of ambient O 3 concentrations at the surface ([O 3 ]) and a transfer coefficient (dry deposition velocity, v d ) that describes the efficiency of transport (and removal) to the surface from the measurement height: (1) Also analogous to other surface fluxes, F O 3 , [O 3 ], and hence v d , can be directly measured by the eddy covariance (EC) method (e.g.Fares et al., 2014;Gerosa et al., 2005;Lamaud et al., 2002;Munger et al., 1996;Rannik et al., 2012) with random uncertainty of about 20 % (Keronen et al., 2003;Muller et al., 2010).Apart from EC, F O 3 and v d can also be estimated from the vertical profile of O 3 by exploiting flux-gradient relationship (Foken, 2006) (termed the gradient method, GM) (e.g.Gerosa et al., 2017;Wu et al., 2015Wu et al., , 2016)).A recent study (Silva and Heald, 2018) complied 75 sets of ozone deposition measurement from the EC and GM across different seasons and land cover types over the past 30 years.
At the site level, ozone dry deposition over various terrestrial ecosystems can be simulated comprehensively by 1-D chemical transport models (Ashworth et al., 2015;Wolfe et al., 2011;Zhou et al., 2017), which are able to simulate the effects of vertical gradients inside the canopy environment, and gas-phase reaction with BVOCs in addition to surface sinks.Regional and global models, which lack the fine-scale information (e.g.vertical structure of canopy, incanopy BVOC emissions) and horizontal resolution for resolving the plant canopy in such detail, instead represent plant canopy foliage as one to two big leaves, and v d is parameterized as a network of resistance, which accounts for the effects of turbulent mixing via aerodynamic (R a ), molecular diffusion via quasi-laminar sublayer resistance (R b ), and surface sinks via surface resistance (R c ): A diverse set of parameterizations of ozone dry deposition is available and used in different models and monitoring networks.Examples include the Wesely parameterization (1989) and modified versions of it (e.g.Wang et al., 1998), the Zhang et al. (2003) parameterization (Zhang et al., 2003), the Deposition of O 3 for Stomatal Exchange model (Emberson et al., 2000;Simpson et al., 2012), and the Clean Air Status and Trends Network (CASTNET) deposition estimates (Meyers et al., 1998).The calculation of R a (mostly based on Monin-Obukhov similarity theory) and R b across these parameterizations often follow a standard formulation from micrometeorology (Foken, 2006;Wesely andHicks, 1977, 2000;Wu et al., 2011) and thus does not vary significantly.
The main difference between the ozone dry deposition parameterizations lies on the surface resistance, R c .This resistance includes stomatal resistance (R s ), which can be computed by a Jarvis-type multiplicative algorithm (Jarvis, 1976) where R s is the product of its minimum value and a series of response functions to individual environmental conditions.Such conditions typically include air temperature (T ), photosynthetically available radiation (PAR), vapour pressure deficit (VPD), and soil moisture (θ), with varying complexity and functional forms.Such formalism is empirical in nature and does not adequately represent the underlying ecophysiological processes affecting R s (e.g.temperature acclimation).An advancement of these efforts includes harmonizing R s with that computed by land surface models (Ran et al., 2017a;Val Martin et al., 2014), which calculate R s by coupled photosynthesisstomatal conductance (A n -g s ) models (Ball et al., 1987;Collatz et al., 1991Collatz et al., , 1992)).Such coupling should theoretically give a more realistic account of ecophysiological controls on R s .Indeed, it has been shown that the above approach may better simulate v d than the multiplicative algorithms that only consider the effects of T and PAR (Val Martin et al., 2014;Wu et al., 2011).
Furthermore, terrestrial atmosphere-biosphere exchange is also directly affected by CO 2 , as CO 2 can drive increases in LAI (Zhu et al., 2016) while inhibiting g s (Ainsworth and Rogers, 2007).These can have important implications on v d , as shown by Sanderson et al. (2007), where doubling current CO 2 level reduces g s by 0.5-2.0mm s −1 , and by Wu et al. (2012) where v d increases substantially due to CO 2 fertilization at 2100.Observations from the Free Air CO 2 Enrichment (FACE) experiments also confirm CO 2 fertilization and inhibition of g s effects, but the impacts are variable and species specific such that extrapolation of these effects to global forest cover is cautioned (Norby and Zak, 2011).
Various efforts have been made to evaluate and assess the uncertainty in modelling ozone dry deposition using field measurements.Hardacre et al. (2015) evaluate the performance of simulated monthly mean v d and F O 3 by 15 chemical transport models (CTMs) from the Task Force on Hemispheric Transport of Air Pollutant (TF HTAP) against seven long-term site measurements, 15 short-term site measurements, and modelled v d from 96 CASTNET sites.This work suggests that the difference in land cover classification is the main source of discrepancy between models.In this case, most of the models in TF HTAP use the same class of dry deposition parameterization (Wang et al., 1998;Wesely, 1989), so a global evaluation of different deposition parameterizations was not possible.Also, the focus in this intercomparison study was on seasonal but not other (e.g.diurnal, daily, interannual) timescales.Using an extended set of measurements, Silva and Heald (2018) evaluate the v d output from the Wang et al. (1998) parameterization used by the GEOS-Chem chemical transport model.They show that diurnal and seasonal cycles are generally well captured, while the daily variability is not well simulated.They find that differences in land type and LAI, rather than meteorology, are the main reason behind model-observation discrepancy at the seasonal scale, and eliminating this model bias results in up to 15 % change in surface O 3 .This study is also limited to a single parameterization.Using parameterizations that are explicitly sensitive to other environmental variables (e.g. Simpson et al., 2012;Zhang et al., 2003) could conceivably lead to different conclusions.
Other efforts have been made to compare the performance of different parameterizations.Centoni (2017) find that two different dry deposition parameterizations, Wesely (1989) versus Zhang et al. (2003), implemented in the same chemistry-aerosol model (United Kingdom Chemistry Aerosol, UKCA, model), result in up to a 20 % difference in simulated surface O 3 concentration.This study demonstrates that uncertainty in v d can have a large potential effect on surface O 3 simulation.Wu et al. (2018) compare v d simulated by five North American dry deposition parameterizations to a long-term observational record at a single mixed forest in southern Canada and find a large spread between the simulated v d , with no single parameterization uniformly outperforming others.They further acknowledge that as each parameterization is developed with its own set of limited observations, it is natural that their performance can vary considerably under different environments, and advocate for an "ensemble" approach to dry deposition modelling.This highlights the importance of parameterization choice as a key source of uncertainty in modelling ozone dry deposition.Meanwhile, in another evaluation at a single site, Clifton et al. (2017) show that the GEOS-Chem parameterization largely underestimates the interannual variability (IAV) of v d in Harvard Forest based on the measurement from 1990 to 2000, although they do not show how the IAV of v d may contribute to the IAV of O 3 .
These developments have made a substantial contribution to our understanding of the importance of O 3 dry deposition in atmospheric chemistry models.Still, pertinent questions remain about the impact of a dry deposition model on simulations of the global distribution of ozone and its long-term variability.Here, we build on previous works by posing and answering the following questions: The answers to such question could have important consequences on our ability to predict long-term changes in atmospheric O 3 concentrations as a function of changing climate and land cover characteristics.In general, there is a high computational cost to thorough and large-scale evaluations of different dry deposition parameterizations embedded in CTMs.In this study, we explore these questions using a strategy that combines an offline dry deposition modelling framework incorporating long-term assimilated meteorological and land surface remote sensing data, in combination with a set of CTM sensitivity simulations.Here, we consider several "big-leaf" models commonly used by global chemical transport models.More complex multilayer models require the vertical profiles of leaf area density for different biomes which are generally not available for regional and global models.From the wide range of literature on dry deposition studies, we observe that R s is commonly modelled through one of the following approaches: 1. a multiplicative algorithm that considers the effects of LAI, temperature, and radiation (Wang et al., 1998); 2. a multiplicative algorithm that considers the effects of LAI, temperature, radiation, and water stress (e.g.Meyers et al., 1998;Pleim and Ran, 2011;Simpson et al., 2012;Zhang et al., 2003); or 3. a coupled A n -g s model, which exploits the strong empirical relationship between photosynthesis (A n ) and stomatal conductance (g s ) (e.g.Ball et al., 1987;Lin et al., 2015) and simulates A n and g s = 1/R s simultaneously (e.g.Ran et al., 2017b;Val Martin et al., 2014).
Similarly, their functional dependence of non-stomatal surface resistance can be classified into two classes: 1. mainly scaling with LAI, with in-canopy aerodynamics parameterized as function of friction velocity (u * ) or radiation (Meyers et al., 1998;Simpson et al., 2012;Wang et al., 1998); and 2. additional dependence of cuticular resistance on relative humidity (Pleim and Ran, 2011;Zhang et al., 2003) With these considerations, we identify four common parameterizations that are representative of the types of approaches described above: 1. the version of Wesely (1989) with the modification from Wang et al. (1998) (hereafter referred to as W98), which is used extensively in global CTMs (Hardacre et al., 2015) and comprehensively discussed by Silva and Heald (2018) -this represents Type 1 in both stomatal and non-stomatal parameterizations; 2. the Zhang et al. (2003) parameterization (hereafter referred to as Z03), which is used in many North American air quality modelling studies (e.g.Huang et al., 2016;Kharol et al., 2018) and Canadian Air and Precipitation Monitoring Network (CAPMoN) (e.g.Zhang et al., 2009) -this represents Type 2 in both stomatal and non-stomatal parameterizations; 3. W98 with R s calculated from a widely used coupled A n -g s model, the Ball-Berry model (hereafter referred to as W98_BB) (Ball et al., 1987;Collatz et al., 1991Collatz et al., , 1992)) Another important consideration in choosing Z03 and W98 is that they both have parameters for all major land types over the globe, making them widely applicable in global modelling.We extract the source code (Wang et al., 1998) and parameters (Baldocchi et al., 1987;Jacob et al., 1992;Jacob and Wofsy, 1990;Wesely, 1989) of W98 from GEOS-Chem CTM (http://wiki.seas.harvard.edu/geos-chem/index. php/Dry_deposition, last access: 24 January 2019).The source code of Z03 was obtained through personal communication with Zhiyong Wu and Leiming Zhang, which follows the series of papers that described the development and formalism of the parameterization (Brook et al., 1999;Zhang et al., 2001Zhang et al., , 2002Zhang et al., , 2003)).The Ball-Berry A n -g s model (Ball et al., 1987;Collatz et al., 1991Collatz et al., , 1992;;Farquhar et al., 1980) and its solver are largely based on the algorithm of CLM (Community Land Model) version 4.5 (Oleson et al., 2013), which is numerically stable (Sun et al., 2012).We use identical formulae of R a and R b (Paulson, 1970;Wesely and Hicks, 1977) for each individual parameterization, allowing us to focus our analysis on differences in parameterizations of R c alone.Table S1 in the Supplement gives a brief description on the formalism of each of the dry deposition parameterizations.

Dry deposition model configuration, inputs, and simulation
The above parameterizations are re-implemented in R language (R core team, 2017) in the modelling framework of the Terrestrial Ecosystem Model in R (https://github.com/amospktai/TEMIR, last access: 15 November 2019) and driven by gridded surface meteorology and land surface datasets.The meteorological forcing chosen for this study is the Modern-Era Retrospective analysis for Research and Applications version 2 (MERRA-2) (Gelaro et al., 2017), an assimilated meteorological product at hourly time resolution spanning from 1980 to the present day.MERRA-2 contains all the required surface meteorological fields except VPD and RH, which can be readily computed from T , specific humidity (q), and surface air pressure (P ).We use the CLM land surface dataset (Lawrence and Chase, 2007), which contains information for land cover, per-grid cell coverage of each plant functional type (PFT), and PFT-specific LAI, which are required to drive the dry deposition parameterizations, and soil property, which is required to drive the A n -g s model in addition to PFT and PFT-specific LAI.CLM land types are mapped to the land type of W98 following Geddes et al. (2016).The mapping between CLM and Z03 land types is given in Table S2.Other relevant vegetation and soil parameters are also imported from CLM 4.5 (Oleson et al., 2013), while land-cover-specific roughness length (z 0 ) values follow Geddes et al. (2016).Leaf is set to be wet when either latent heat flux < 0 W m −2 or precipitation > 0.2 mm h −1 .Fractional coverage of snow for Z03 is parameterized as a landtype-specific function of snow depth following the original paper of Z03, while W98 flags grid cells with albedo > 0.4 or permanently glaciated as snow covered.
As the IAV of LAI could be an important factor in simulating v d , the widely used third-generation Global Inventory Modelling and Mapping Studies leaf area index product (GIMMS LAI3g, abbreviated as LAI3g in this paper) (Zhu et al., 2013), which is a global time series of LAI with 15 d temporal frequency and 1/12 • spatial resolution spanning from late 1981 to 2011, is incorporated in this study.We derive the interannual scaling factors that can be applied to scale the baseline CLM-derived LAI (Lawrence and Chase, 2007) for each month over 1982 to 2011.All the input data are aggregated into horizontal resolution of 2 • × 2.5 • to align with the CTM sensitivity simulation described in the next subsection.To represent subgrid land cover heterogeneity, grid-cell-level v d is calculated as the sum of v d over all subgrid land types weighted by their percentage coverage in the grid cell (a.k.a tiling or mosaic approach, e.g.Li et al., 2013).This reduces the information loss when land surface data are aggregated to coarser spatial resolution and allows us to retain PFT-specific results for each grid box in the offline dry deposition simulations.
We run three sets of 30-year (1982-2011) simulations with the deposition parameterizations to investigate how v d simulated by different parameterizations responds to different environmental factors over multiple decades.The settings of the simulations are summarized in Table 1.The first set, [Clim], focuses on meteorological variability alone, driven by MERRA-2 meteorology and a multiyear (constant) mean annual cycle of LAI derived from LAI3g.The second set, [Clim + LAI], combines the effects of meteorology and IAV in LAI, driven by the same MERRA-2 meteorology plus the LAI time series from LAI3g.As the increase in atmospheric CO 2 level over multidecadal timescales may lead to significant reduction in g s as plants tend to conserve water (e.g.Franks et al., 2013;Rigden and Salvucci, 2017), we introduce the third set of simulations, [Clim + LAI + CO 2 ], which is driven by varying meteorology and LAI, plus the annual mean atmospheric CO 2 level measured in Mauna Loa (Keeling et al., 2001) (for the first two sets of simulations, atmospheric CO 2 concentration held constant at 390 ppm).Since W98 and Z03 do not respond to changes in CO 2 level, only W98_BB and Z03_BB are run with [Clim + LAI + CO 2 ] to evaluate this impact.We focus on the daytime (solar elevation angle > 20 • ) v d , as both v d and surface O 3 concentration typically peak around this time.We calculate monthly means, filtering out the grid cells with monthly total daytime < 100 h.
In summary, we present for the first time a unique set of global dry deposition velocity predictions over the last 30 years driven by identical meteorology and land cover, so that discrepancies (in space and time) among the predicted v d are a result specifically of dry deposition parameterization choice or assumptions about how land cover is treated in each.

Chemical transport model sensitivity experiments
We quantify the sensitivity of surface O 3 to variations in v d using a global 3-D CTM, GEOS-Chem version 11.01 (http: //acmg.seas.harvard.edu/geos/,last access: 20 June 2019), which includes comprehensive HO x -NO x -VOC-O 3 -BrO x chemical mechanisms (Mao et al., 2013) and is widely used to study tropospheric ozone (e.g.Hu et al., 2017;Travis et al., 2016;Zhang et al., 2010).The model is driven by the assimilated meteorological data from the GEOS-FP (forward processing) Atmospheric Data Assimilation System (GEOS-5 ADAS) (Rienecker et al., 2008), which is jointly developed by National Centers for Environmental Prediction (NCEP) of National Oceanic and Atmospheric Administration (NOAA) and the Global Modelling and Assimilation Office (GMAO).The model is run with a horizontal resolution of 2 • × 2.5 • and 47 vertical layers.The dry deposition module, which has been discussed above (W98), is driven by the monthly mean LAI retrieved from Moderate Resolution Imaging Spectroradiometer (MODIS) (Myneni et al., 2002) and the 2001 version of Olson land cover map (Olson et al., 2001).Both of the maps are remapped from their native resolutions to 0.25 • × 0.25 • .We propose to estimate the sensitivity of surface O 3 concentrations to uncertainty/changes in v d by the following equation: where O 3 is the response of monthly mean daytime surface O 3 to fractional change in v d ( v d /v d ), and β accounts for the sensitivity of surface O 3 concentration in a grid box to the perturbation in v d within that grid box.To estimate β, we run two simulations for the year 2013: one with the default setting and another where we perturb v d by +30 %.Thus, this approach could represent a conservative estimate of O 3 sensitivity to v d if the impacts on other species result in additional effects on O 3 .We use this sensitivity to identify areas where local uncertainty and variability in v d is expected to affect local surface O 3 concentration, and we use the assumption of linearity to estimate those impacts to a first order (e.g.Wong et al., 2018).In the Supplement, we justify this first-order assumption mathematically, as well as demon- LAI3g monthly time series Mauna Loa time series strate the impact of using a second-order approximation, and estimate the uncertainty using an assumption of linearity to be within 30 %.However, we note this first-order assumption may not be able to capture the effects of chemical transport, changes in background ozone, and non-linearity in chemistry, which can contribute to response of O 3 concentration to v d .Our experiment could help identify regions where more rigorous modelling efforts could be targeted in future work.We limit our analysis to grid cells where the monthly average v d is greater than 0.25 cm s −1 in the unperturbed GEOS-Chem simulation, since changes in surface O 3 elsewhere are expected to be attributed more to change in background O 3 rather than the local perturbation of v d (Wong et al., 2018).

Evaluation of dry deposition parameterizations
We first compare our offline simulations of seasonal mean daytime average v d that result from the four parameterizations in the [Clim] and [Clim + LAI] scenarios with an observational database largely based on the evaluation presented in Silva and Heald (2018).We do not include the evaluation of v d from the [Clim + LAI + CO 2 ] scenario, as we find that the impact of CO 2 concentration on v d is negligible over the period of concern, as we will show in subsequent sections.We use two unbiased and symmetrical statistical metrics, normalized mean bias factor (NMBF) and normalized mean absolute error factor (NMAEF), to evaluate our parameterizations.Positive NMBF indicates that the parameterization overestimates the observations by a factor of 1 + NMBF and the absolute gross error is NMAEF times the mean observation, while negative NMBF implies that the parameterization underestimates the observations by a factor of 1−NMBF and the absolute gross error is NMAEF times the mean model prediction (Yu et al., 2006).We use the simulated subgrid land-type-specific predictions of v d that correctly match the land type and the averaging window indicated by the observations.We exclude instances where the observed land type does not have a match within the model grid box.While this removes one-third of the original datasets used in Silva and Heald (2018), this means that mismatched land cover types can be ignored as a factor in model bias.
Figure 1 shows the fractional coverage within each grid cell and the geographic locations of O 3 flux observation sites for each major land type.Nearly all the observations are clustered in Europe and North America, except three sites in the tropical rainforest and one site in a tropical decidu-ous forest in Thailand.For most major land types, there are significant mismatches between the locations of flux measurements and the dominant land cover fraction, which may hinder the spatial representativeness of our evaluation.The resulting NMBF and NMAEF for five major land-type categories are shown in Table 2, and the list of sites and their descriptions are given in Table S3.In general, the numerical ranges of both NMBF and NMAEF are similar to that of Silva and Heald (2018), and no single parameterization of the four parameterizations outperforms the others across all five major land types.
The performance metrics of each parameterization at each land type are summarized in Table 2. Comparing the two multiplicative parameterizations (W98 and Z03), we find that W98 performs satisfactorily over deciduous forests and tropical rainforests, while strongly underestimating daytime v d over coniferous forests.In contrast, Z03 performs better in coniferous forests but worse in tropical rainforests and deciduous forests.The severe underestimation of daytime v d by Z03 over tropical rainforests has previously been attributed to persistent canopy wetness and hence stomatal blocking imposed by the parameterization (Centoni, 2017).We also note that even for the same location, v d can vary significantly between seasons (Rummel et al., 2007) and management practices (Fowler et al., 2011), which models may fail to capture due to limited representations of land cover.Given the small sample size (N = 5), diverse environments, and large anthropogenic intervention in the tropics, the disparity in performance metrics may not fully reflect the relative model performance.Baseline cuticular resistance in Z03 under dry and wet canopy is 1.5 and 2 times that of coniferous forests, respectively (Zhang et al., 2003), such that the enhancement of cuticular uptake by wetness may not compensate the reduced g s over tropical rainforests, and, to a lesser extent, deciduous forests.
Over grasslands, W98 has higher positive biases, while Z03 has higher absolute errors.This is because for datasets at high latitudes, the dominant grass PFT is arctic grass, which is mapped to "tundra" land type (Geddes et al., 2016).While tundra is parameterized similarly to grasslands in W98, this is not the case in Z03.Combined with the general high biases at other sites for these parameterizations, the large low biases for "tundra" sites in Z03 lower the overall high biases but lead to higher absolute errors.
Over croplands, the positive biases and absolute errors are relatively large for both W98 and Z03 (with Z03 performing  Substituting the native g s in W98 and Z03 by that simulated by the Ball-Berry model (the W98_BB and Z03_BB runs) generally, though not universally, leads to improvement in model performance against the observations.W98_BB has considerably smaller biases and absolute errors than W98 over grassland.While having little effect on the absolute error, W98_BB improves the biases over coniferous forest and cropland compared to W98 but worsens the biases over rainforests and deciduous forests.In contrast, Z03_BB is able to improve the model-observation agreement over all five land types when compared to Z03.This finding echoes that from Wu et al. (2011), who explicitly show the advantage of replacing the g s of Wesely (1989)  tential of the Ball-Berry model in improving the spatial distribution of mean v d .The different responses to substituting native g s with that from the Ball-Berry model highlight the significant differences in parameterizing non-stomatal uptake between W98 and Z03, which further suggests that the uncertainty in non-stomatal deposition should not be overlooked.
The minimal impact that results from using LAI that matches the time of observation is not unexpected, since the meteorological and land cover information from a 2 • × 2.5 • grid cell may not be representative of the typical footprint of a site measurement (on the order of 10 −3 to 10 1 km 2 , e.g.Chen et al., 2009Chen et al., , 2012)).The mismatch between model resolution and the footprint of site-level measurements has also been highlighted in previous evaluation efforts in globalscale CTMs (Hardacre et al., 2015;Silva and Heald, 2018).Furthermore, the sample sizes for all land types are small (N ≤ 16) and the evaluation may be further compromised by inherent sampling biases.
In addition to the evaluation against field observation, we find good correlation (R 2 = 0.94) between the annual mean v d from GEOS-Chem at 2013 and the 30-year mean v d of W98 run with static LAI, providing further evidence that our implementation of W98 is reliable.Overall, our evaluation shows that the quality of our offline simulation of dry deposition across the four parameterizations in this work is largely consistent with previous global modelling evaluation efforts.

Impact of dry deposition parameterization choice on long-term averages
Here, we summarize the impact that the different dry deposition parameterizations may have on simulations of the spatial distribution of v d and on the inferred surface O 3 concentrations.We begin by comparing the simulated long-term mean v d across parameterizations, then use a chemical transport model sensitivity experiment to estimate the O 3 impacts.
Figure 2 shows the 30-year July daytime average v d simulated by W98 over vegetated surfaces (defined as the grid cells with > 50 % plant cover), and Fig. 3 shows the difference between the W98 and the W98_BB, Z03, and Z03_BB predictions, respectively.We first focus on results from July because of the coincidence of high surface O 3 level, biospheric activity, and v d in the Northern Hemisphere (NH) and will subsequently discuss the result for December, when such a condition holds for the Southern Hemisphere (SH).W98 simulates the highest July mean daytime v d in Amazonia (1.2 to 1.4 cm s −1 ), followed by other major tropical rainforests, and temperate forests in the northeastern US.July mean daytime v d in other temperate regions in North America and Eurasia typically ranges from 0.5 to 0.8 cm s −1 , while in South American and the African savannah, and most parts of China, daytime v d is around 0.4 to 0.6 cm s −1 .In India, Australia, the western US, polar tundra, and Mediterranean region, July mean daytime v d is low (0.2-0.5 cm s −1 ).
The other three parameterizations (W98_BB, Z03, Z03_BB) simulate substantially different spatial distributions of daytime v d .In North America, we find W98_BB, Z03, and Z03_BB produce lower v d (by −0.1 to −0.4 cm s −1 ) compared to W98 in the deciduous-forest-dominated northeastern US and slightly higher v d in boreal-forest-dominated regions of Canada.Z03 and Z03_BB produce noticeably lower v d (by up to −0.2 cm s −1 ) in arctic tundra and grasslands in the western US.In the southeastern US, W98_BB and Z03_BB simulate a slightly higher v d (by up to +0.1 cm s −1 ), while Z03 suggests a slightly lower v d (by up to −0.1 cm s −1 ).W98_BB simulates a lower (−0.1 to −0.4 cm s −1 ) v d in tropical rainforests, with larger reductions concentrated in southern Amazonia, where July is within the dry season, while the northern Amazonia is not (Malhi et al., 2008).Z03 and Z03_BB simulate much smaller (−0.4 to −0.6 cm s −1 ) v d in all tropical rainforests.
Over the midlatitudes in Eurasia, Australia, and South America except Amazonia, W98_BB, Z03, and Z03_BB generally simulate a lower daytime v d by up to 0.25 cm s −1 , possibly due to the dominance of grasslands and deciduous forests, where W98 tends to be more high biased than other parameterizations when compared to the observations of v d .In the southern African savannah, W98_BB and Z03_BB suggest a much lower daytime v d (by −0.1 to −0.4 cm s −1 ) because of explicit consideration of soil moisture limitation to A n and g s (demonstrated by the spatial overlap with soil moisture stress factors shown in Fig. S2).Z03_BB simulates a particularly high daytime v d over the high-latitude coniferous forests (+0.1 to +0.3 cm s −1 ).W98_BB and Z03_BB produce higher daytime v d (up to +0.15 cm s −1 ) in India and south China due to temperature acclimation (Kattge and Knorr, 2007), which allows more stomatal opening under the high temperature that would largely shut down the stomatal deposition in W98 and Z03, as long as the soil does not become too dry to support stomatal opening.This is guaranteed by the rainfall from summer monsoon in both regions.Low  Our results suggest that the global distribution of simulated mean v d depends substantially on the choice of dry deposition parameterization, driven primarily by the response to hydroclimate-related parameters such as soil moisture, VPD, and leaf wetness, in addition to land-type-specific parameters, which could impact the spatial distribution of surface ozone predicted by chemical transport models.To estimate the impact on surface ozone of an individual parameterization "i" compared to the W98 predictions (which we use as a baseline), we apply the following equation: where O 3,i is the estimated impact on simulated O 3 concentrations in a grid box, v d,i is the difference between parameterization i and W98 simulated mean daytime v d in that grid box, v d W98 is W98 output mean daytime v d for that grid box, and β is the sensitivity of surface ozone to v d calculated by the method outlined in Sect.2.3.Figure 4 shows the resulting estimates of O 3 globally.We find O 3 is the largest in tropical rainforests for all the parameterizations (up to 5 to 8 ppbv).Other hotspots of substantial differences are boreal coniferous forests, the eastern US, continental Europe, Eurasian steppe, and the grassland in southwestern China, where O 3 is either relatively large or the signs disagree among parameterizations.In India, Indochina, and south China, O 3 is relatively small but still reaches up to up to −2 ppbv.We find that O 3 is not negligible (1-4 ppbv) in many regions with relatively high population density, which suggests that the choice of dry deposition parameterization can be relevant to the uncertainty in the study of air quality and its implication on public health.We note that we have not estimated O 3 for some regions with low GEOS-Chem-predicted v d (< 0.25 cm s −1 , as described in Sect.2.3) but where the disagreement in v d between parameterizations can be large (e.g. the southern African sa-vannah; see Fig. 3).Given this limitation, the impacts on O 3 we have summarized may therefore be spatially conservative.
To explore the impact of different prediction of v d on surface O 3 in different seasons, we repeat the above analyses for December.Figure 5 shows the 1982-2011 mean December daytime v d predicted by W98, while Fig. 6 shows the difference between W98 and the Z03, W98_BB, and Z03_BB, respectively.High latitudes in the NH are excluded due to the small number of daytime hours.Z03 and Z03_BB simulate substantially lower in daytime v d at NH midlatitudes because Z03 and Z03_BB allow partial snow cover but W98 and W98_BB only allow total or no snow cover.At midlatitudes, the snow cover is not high enough to trigger the threshold of converting vegetated to snow-covered ground in W98 and W98_BB, resulting in lower surface resistance, and hence higher daytime v d compared to Z03 and Z03_BB; in Amazonia, the hotspot of difference in daytime v d shifts from the south to the north relative to July, which is in the dry season (Malhi et al., 2008).These results for December, together with our findings from July, suggest that the discrepancy in simulated daytime v d between W98 and other parameterizations is due to the explicit response to hydroclimate in the former compared to the latter.Given that field observations indicate a large reduction of v d in dry season in Amazonia (Rummel et al., 2007), the lack of dependence of hydroclimate can be a drawback of W98 in simulating v d in Amazonia.
Figure 7 shows the resulting estimates of O 3 globally for December using Eq. ( 3).In all major rainforests, O 3 is smaller in December due to generally lower sensitivity compared to July.A surprising hotspot of both daytime v d and O 3 is the rainforest/tropical deciduous forest in Myanmar and its eastern bordering region, which also has distinct wet and dry seasons.The proximity of December to the dry season, which starts in January (e.g.Matsuda et al., 2005), indicates that the consistent v d between W98 and other parameterizations is driven by hydroclimate as in Amazonia.
Comparison with field measurements (Matsuda et al., 2005) suggests that the W98_BB and Z03_BB capture daytime v d better than W98, while Z03 may overemphasize the effect of such dryness.The above reasoning also explains some of  the v d in India and south China across the three parameterizations.These findings identify hydroclimate as a key driver of process uncertainty of v d over tropics and subtropics, and therefore its impact on the spatial distribution of surface ozone concentrations, independent of land-type-based biases, in these regions.
Overall, these results demonstrate that the discrepancy in the spatial distribution of simulated mean daytime v d resulting from choice of dry deposition parameterization can have an important impact on the global distribution of surface O 3 predicted by chemical transport models.We find that the response to hydroclimate by individual parameterization not only affects the mean of predicted surface ozone but also has different impacts in different seasons, which is complementary to the findings of Kavassalis and Murphy (2017) that mainly focus on how shorter-term hydrometeorological variability may modulate surface O 3 through dry deposition.

Impact of dry deposition parameterization choice on trends and interannual variability
Here, we explore the impact that different dry deposition parameterizations may have on predictions of IAV and trends in v d and on the inferred surface O 3 concentrations.We use the Theil-Sen method (Sen, 1968), which is less susceptible to outliers than least-square methods, to estimate trends in July daytime v d (and any underlying meteorological variables) and use p value < 0.05 to estimate significance.
Figure 8 shows the trend in July mean daytime v d from 1982 to 2011 predicted by each of the parameterizations and scenarios ([Clim], [Clim + LAI], and [Clim + LAI + CO 2 ]). Figure 9 shows the potential impact of these trends in v d on July daytime surface ozone, which we estimate to a first order using the following equation: where O 330y,i and m v d,i are the absolute change in ozone inferred to a first order as a result of the trend of v d and the normalized Theil-Sen slope (% yr −1 ) of v d for parameterization i over the 30 years .
In [Clim] simulations (where LAI is held constant), significant decreasing trends in July daytime v d are simulated by the Z03, W98_BB, and Z03_BB in Mongolia, where significant increasing trend in T (warming) and decreasing trend in RH (drying) detected in the MERRA-2 surface meteorological field in July daytime.This trend is not present in the W98 parameterization as this formulation does not respond to the long-term drying.We find some decreasing trends in v d across parts of central Europe and the Mediterranean to varying degrees across the parameterizations.In the SH, we find consistent decreasing trends across all four parameterizations in southern Amazonia and southern African savannah due to warming and drying, which we estimate could produce a concomitant increase in July mean surface ozone of between 1 to 3 ppbv (Fig. 9).
In [Clim + LAI] scenario, all four parameterizations simulate a significant increasing trend of v d over high latitudes, which is consistent with the observed greening trend over the region (Zhu et al., 2016).We estimate this could pro-  duce a concomitant decrease in July mean surface ozone of between 1 and 3 ppbv.The parameterizations generally agree in terms of the spatial distribution of these trends in O 3 .Exceptions include a steeper decreasing trend in most of Siberia predicted by W98, while the trend is more confined in eastern and western Siberia in the other three parameterizations.Including the effect of CO 2 -induced stomatal closure ([Clim + LAI + CO 2 ] runs) partially offsets the increase of v d in high latitudes but does not lead to large changes in both the magnitudes and spatial patterns of v d trend.We find negligible trends in daytime v d for December in all cases.These results show that across all dry deposition model parameterizations, LAI and climate, more than increasing CO 2 , can potentially drive significant long-term changes in v d and should not be neglected when analysing the long-term change in air quality over 1982-2011.We note that the importance of the CO 2 effect could grow as the period of study extends to allow a larger range of atmospheric CO 2 concentrations (Hollaway et al., 2017;Sanderson et al., 2007).
We go on to explore the impact of parameterization choice in calculations of IAV in v d .Figure 10 shows the coefficient of variation of linearly detrended July daytime v d (CV v d ).
Figure 11 shows the potential impact this has on IAV in surface ozone, which we estimate to a first order by the following equation: where σ O 3 ,i is the estimated interannual standard deviation in surface ozone resulting from IAV in v d given predicted by dry deposition parameterization i.In both cases, we show only the [Clim] and [Clim + LAI] runs, since IAV in CO 2 has negligible impact on interannual variability in v d .
Using the W98 parameterization, IAV in predicted v d and O 3 is considerably smaller in the [Clim] run than that for the [Clim + LAI] run, since both the stomatal and non-stomatal conductance in W98 are assumed to be strong functions of LAI rather than meteorological conditions.This implies that long-term simulations with W98 and constant LAI can potentially underestimate the IAV of v d and surface ozone.In contrast, IAV in v d calculated by the Z03 parameterization is nearly the same for the [Clim] and [Clim + LAI] runs.In Z03, g s is also directly influenced by VPD in addition to temperature and radiation, and non-stomatal conductance in Z03 is much more dependent on meteorology than W98, leading to high sensitivity to climate.Though the Ball-Berry model also responds to meteorological conditions, it considers relatively complex A n -g s regulation and includes temperature acclimation, which could dampen its sensitivity to meteorological variability compared to the direct functional dependence on meteorology in the Z03 multiplicative algorithm.Thus, the climate sensitivity of W98_BB and Z03_BB is in between Z03 and W98, as is indicated by a more moderate difference between σ O 3 ,i from [Clim] and [Clim + LAI] runs in Fig. 11.All parameterizations produce larger CV v d (and therefore larger σ O 3 ) in southern Amazonia compared to northern and central Amazonia, but we find substantial discrepancies across parameterizations.The estimated impact on IAV in O 3 in southern Amazonia ranges from less than 1 ppbv, predicted by the W98 and W98_BB parameterizations, to exceeding 1.5-2.5 ppbv predicted by the Z03 parameterization.IAV is also relatively large in central Africa.We find that the parameterizations which include a Ball-Berry formulation (W98_BB and Z03_BB) estimate higher IAV in this region (with σ O 3 varying between 1 and 4 ppbv), compared to the W98 and Z03 parameterizations (σ O 3 up to 2 ppbv).We also note that the Ball-Berry formulations show more spatial heterogeneity compared to W98 and Z03.In our implementation of the Ball-Berry model, impact of soil moisture on g s is parameterized as a function of root-zone soil matric potential, which makes g s very sensitive to variation in soil  wetness when its climatology is near the point that triggers limitation on A n and g s .Given the large uncertainty in global soil property map (Dai et al., 2019), such sensitivity could be potentially artificial, which should be taken into consideration when implementing Ball-Berry parameterizations in large-scale models despite their relatively good performance in site-level evaluation (Wu et al., 2011).
Across Europe, the magnitudes of IAV predicted by all four parameterizations show relatively good spatial consistency.Simulated CV v d is relatively low in western and northern Europe (< 10 %), which we estimate translates to less than 1 ppbv of σ O 3 .We find larger CV v d (and therefore large σ O 3 ) over parts of southern Russia and Siberia (σ O 3 up to 2.5 ppbv) from all parameterizations except W98.The local geographic distribution of CV v d and σ O 3 also significantly differs among the parameterizations there.Z03 and Z03_BB simulate larger CV v d in eastern Siberia than W98_BB, while W98 BB and Z03_BB predict larger CV v d over the southern Russian steppe then Z03.Finally, all four parameterizations estimate relatively low CV v d and σ O 3 in India, China, and southeast Asia.
We compare the simulated IAV July CV v d from all four deposition parameterizations with those recorded by publicly available long-term observations.Hourly v d is calculated using Eq. ( 1) from raw data.We filter out the data points with extreme (> 2 cm s −1 ) or negative v d , and without enough turbulence (u * < 0.25 m s −1 ).As v d values in each daytime hour are not uniformly sampled in the observational datasets, we calculate the mean diurnal cycle and then calculate the daytime average July v d for each year from the mean diurnal cycle, from which CV v d can be calculated.
The IAV predicted by all four parameterizations at Harvard Forest is between 3 % and 7.9 %, which is 2 to 6 times lower than that presented in the observations (18 %).We find similar underestimations by all four parameterizations compared to the long-term observation from Hyytiälä (Junninen et al., 2009;Keronen et al., 2003; https://avaa.tdata.fi/web/smart/smear/download, last access: 25 July 2019), where observed CV v d (16 %) is significantly higher than that predicted by the deposition parameterizations (3.5 %-7.1 %).In Blodgett Forest, we find that the models underestimate the observed annual CV v d more seriously (∼ 1 %-3 % compared to 18 % in the observations).This suggests that the IAV of v d may be underestimated across all deposition parameterizations we investigated (and routinely used in simulations of chemical transport).Clifton et al. (2019) attribute this to the IAV in deposition to wet soil and dew-wet leaves, and incanopy chemistry under stressed conditions for forests over the northeastern US.Some of these processes (e.g.in-canopy chemistry, wetness slowing soil ozone uptake) are not represented by existing parameterizations, contributing to their difficulty in reproducing the observed IAV.The scarcity of long-term ozone flux measurements (Fares et al., 2010(Fares et al., , 2017;;Munger et al., 1996;Rannik et al., 2012) limits our ability to benchmark the IAV in our model simulations with observational datasets.
In summary, when both the variability in LAI and climate are considered, the IAV in simulated v d translates to IAV in surface O 3 of 0.5-2 ppbv in July for most regions.Such variability is predicted to be particularly strong in southern Amazonian and central African rainforest, where the predicted IAV in July surface O 3 due to dry deposition can be as high as 4 ppbv.This suggests that IAV of v d can be an important part of the natural variability of surface O 3 .The estimated magnitude of IAV is also dependent of the choice of v d parameterization, which highlights the importance of v d parameterization choice on modelling IAV of surface O 3 .

Discussion and conclusion
We present the results of multidecadal global modelling of ozone dry deposition using four different ozone deposition parameterizations that are representative of the major types of approaches of gaseous dry deposition modelling used in global chemical transport models.The parameterizations are driven by the same assimilated meteorology and satellitederived LAI, which minimizes the uncertainty of model input across parameterization and simplifies interpretation of intermodel differences.The output is evaluated against field observations and shows satisfactory performance.One of our main goals was to investigate the impact of dry deposition parameterization choice on long-term averages, trends, and IAV in v d over a multidecadal timescale and estimate the potential concomitant impact on surface ozone concentrations to a first order using a sensitivity simulation approach driven by the GEOS-Chem chemical transport model.
We find that the performance of the four dry deposition parameterizations against field observations varies considerably over land types, and these results are consistent with other evaluations, reflecting the potential issue that dry deposition parameterizations can often be overfit to a particular set of available observations, requiring caution in their application at global scales.We also find that using more ecophysiologically realistic output g s predicted by the Ball-Berry model can generally improve model performance but at the cost of high sensitivity to relatively unreliable soil data.
However, the number of available datasets of ozone dry deposition observation is still small and concentrated in North America and Europe.We know of only one multi-season direct observational record in Asia (Matsuda et al., 2005) and none in Africa, where air quality can be an important issue.To better constrain regional O 3 dry deposition, effort must be made in making new observations of gaseous dry deposition (Fares et al., 2017) especially in the undersampled regions.Evaluation and development of ozone dry deposition parameterizations will continue to benefit from publicly available ozone flux measurements and related micrometeorological variables that allow for partitioning measured flux into individual deposition pathways (e.g.Clifton et al., 2017Clifton et al., , 2019;;Fares et al., 2010;Wu et al., 2011Wu et al., , 2018)).
We find substantial disagreement in the spatial distribution between the mean daytime v d predicted by the different parameterizations we tested.We find that these discrepancies are in general a function of both location and season.In NH summer, v d values simulated by the four parameterizations are considerably different in many regions over the world.We estimate that this could lead to around 2 to 5 ppbv in uncertainty of surface ozone concentration simulations over a vast majority of land in the NH.In tropical rainforests, where leaf wetness is prevalent and the dry-wet season dynamics can have a large impact on v d (Rummel et al., 2007), we estimate the uncertainty due to dry deposition model choice could even lead to an uncertainty in surface ozone of up to 8 ppbv.We also find noticeable impacts in parameterization choice during SH summer, but we note that due to the unreliability of β at low v d , we have not assessed its impact on surface ozone in many high-latitude regions of the NH.In general, we find hydroclimate to be an important driver of the uncertainty.This demonstrates that the potential impact of parameterization choice (or process uncertainty) of v d is neither spatiotemporally uniform nor negligible in many regions over the world.More multi-seasonal observations are especially needed over seasonally dry ecosystems where the role of hydroclimate in deposition parameterizations needs to be evaluated.Recently, standard micrometeorological measurements have been used to derive g s and stomatal deposition of O 3 over North America and Europe (Ducker et al., 2018), highlighting the potential of using global networks of micrometeorological observation (e.g.FLUXNET; Baldocchi et al., 2001) to benchmark and calibrate g s of dry deposition parameterizations, which could at least increase the spatiotemporal representativeness, if not the absolute accuracy, of dry deposition parameterizations, since it would be difficult to constrain non-stomatal sinks with this method.Further research is required to more directly verify whether betterconstrained g s leads to improved v d simulation.
Over the majority of vegetated regions in the NH, we estimate the IAV of mean daytime v d is generally on the order of 5 % to 15 % and may contribute between 0.5 and 2 ppbv of IAV in July surface O 3 over the 30-year period considered here, with each parameterization simulating different geo-graphic distribution of where IAV is highest.The predicted IAV from all four models is smaller than what long-term observations suggest, but its potential contribution to IAV in O 3 is still comparable to the long-term variability of background ozone over similar timescales in US summer (Brown-Steiner et al., 2018;Fiore et al., 2014).This would seem to confirm that v d may be a substantial contributor to natural IAV of O 3 in summer, at least in the US.In the Southern Hemisphere, the IAV mainly concentrates in the drier part of tropical rainforests.The Ball-Berry parameterizations simulate large and spatially discontinuous CV v d and σ O 3 due to their sensitivity to soil wetness.Globally, we find that IAV of v d in W98 is mostly driven by LAI, while in other parameterizations climate generally plays a more important role.We therefore emphasize that temporal matching of LAI is important for consistency when W98 is used in long-term simulations.While our results show notable impacts across the globe, in many regions there are no available long-term observations to evaluate the model predictions over interannual timescales.This information is helpful in designing and identifying sources of error in model experiments that involve variability of v d .
We are also able to detect statistically significant trends in July daytime v d over several regions.The magnitudes of trends are up to 1 % per year and both climate and LAI contribute to the trend.All four deposition parameterizations identify three main hotspots of decreasing July daytime v d (southern Amazonia, southern African savannah, Mongolia), which we link mainly to increasing surface air temperature and decreasing relative humidity.Meanwhile, extensive areas at high latitudes experience LAI-driven increasing July daytime v d , consistent with the greening trend in the region (Zhu et al., 2016).We do not find a strong influence of CO 2induced stomatal closure in the trend over this time period.Over the 30 years, we estimate the trend in July daytime v d could translate approximately to 1 to 3 ppbv of ozone changes in the areas of impact, indicating the potential effect of long-term changes in v d on surface ozone.This estimate should be considered conservative, since we are unable to reliably test the sensitivity of ozone to regions with low v d with our approach.
While the approach we have presented here allows us to explore the role of dry deposition parameterization choice on simulations of long-term means, trends, and IAV in ozone dry deposition velocity, there remain some limitations and opportunities for development.First, we only used one LAI and assimilated meteorological product.The geographic distribution of trend and IAV of v d may vary considerably as the LAI and meteorological products used due to their inherent uncertainty (e.g.Jiang et al., 2017).While we expect the qualitative conclusions about how LAI and climate controls the modelled trend and IAV of v d to be robust to the choice of dataset, the magnitude and spatial variability could be affected.Second, the estimated effects on surface O 3 are a first-order inference based on a linear approximation of the impact that v d has directly on O 3 .We have not applied our analysis to regions with low GEOS-Chem v d , where other components of parameterization (e.g.definition and treatment of snow cover, difference in ground resistance) may have major impact on v d prediction (Silva and Heald, 2018), nor accounted for the role that v d variability can have on other chemical species which would have feedbacks on O 3 .Moreover, the sensitivity of surface ozone to v d may be dependent on the choice of chemical transport model (here, the GEOS-Chem model has been used) and possibly the choice of simulation year for the sensitivity simulation.Finally, we have neglected the effect of land use and land cover change on global PFT composition at this stage, which can be another source of variability for v d and even long-term LAI retrieval (Fang et al., 2013).Nevertheless, the relatively high NMAEF of simulated v d and the inherent uncertainty in input data (land cover, soil property, assimilated meteorology, and LAI) are considered as the major sources of uncertainty in our predictions of v d .
The impact of dry deposition parameterization choice may also have impacts which we have not explored in this study on other trace gases with deposition velocity controlled by surface resistance, and for which stomatal resistance is an important control of surface resistance (e.g.NO 2 ).As v d has already been recognized as a major source of uncertainty in deriving global dry deposition flux of NO 2 and SO 2 (Nowlan et al., 2014), systematic investigation on the variability and uncertainty of v d for other relevant chemical species does not only contribute to understanding the role of gaseous dry deposition on air quality but also to biogeochemical cycling.Particularly, gaseous dry deposition has been shown to be a major component in nitrogen deposition (Geddes and Martin, 2017;Zhang et al., 2012), highlighting the potential importance of understanding the role of v d parameterization in modelling regional and global nitrogen cycles.
Here, we have built on the recent investigations of modelled global mean (Hardacre et al., 2015;Silva and Heald, 2018) and observed long-term variability (Clifton et al., 2017) of O 3 v d .We are able to demonstrate the substantial impact of v d parameterization on modelling the global mean and IAV of v d , and their non-trivial potential impact on simulated seasonal mean and IAV of surface ozone.We demonstrate that the parameterizations with explicit dependence on hydroclimatic variables have higher sensitivity to climate variability than those without.Lin et al. (2019) likewise recently demonstrated the importance of accounting for water availability in O 3 dry deposition modelling.Difficulties in evaluating predictions of v d for many regions of the world (e.g.most of Asia and Africa) persist due to the scarcity of measurements.This makes a strong case for additional measurement and model studies of ozone dry deposition across different timescales, which would be greatly facilitated by an open data-sharing infrastructure (e.g.Baldocchi et al., 2001;Junninen et al., 2009).

1.
How does the global distribution of mean v d vary with different dry deposition parameterizations, and what drives the discrepancies among them?How much might the choice of deposition parameterization affect spatial distribution of surface ozone concentration simulated by a chemical transport model? 2. How are the IAV and long-term trends of v d different across deposition parameterizations, and what drives the discrepancies among them?Do they potentially contribute different predictions of the long-term temporal variability in surface ozone?

Figure 1 .
Figure 1.Fractional coverage of each major land type at each grid cell.Blue dots indicate the locations of the observational sites.

v
d is simulated by Z03 and Z03_BB in the grasslands near the Tibetan Plateau because the grasslands are mainly mapped to the tundra land type, which typically has low v d as discussed in Sect.3.

Figure 4 .
Figure 4.Estimated difference in July mean surface ozone ( O 3 ) due to the discrepancy of simulated July mean daytime v d among the parameterizations.

Figure 5 .
Figure 5.The 1982-2011 December mean daytime v d (solar elevation angle > 20 • ) over vegetated land surface simulated by W98.The data over high latitudes over Northern Hemisphere are invalid due to insufficient daytime hours over the month (< 100 h month −1 ).

Figure 7 .
Figure 7.Estimated difference in December mean surface ozone ( O 3 ) due to the discrepancy of simulated December mean daytime v d among the parameterizations.

Figure 8 .
Figure 8. Trends of July mean daytime v d during 1982-2011 over vegetated land surface.Black dots indicate statistically significant trends (p < 0.05).

Figure 9 .
Figure 9.Estimated impact of trends of July mean daytime v d on July mean surface ozone ( O 330y ) during 1982-2011 over vegetated land surface.Only grid points with statistically significant trends (p < 0.05) in July mean daytime v d are considered.

Figure 10 .
Figure 10.Interannual coefficient of variation of linearly detrended July mean daytime v d (CV v d ) during 1982-2011 over vegetated land surface.

Figure 11 .
Figure 11.Estimated contribution of IAV in July mean daytime v d to IAV of July mean surface ozone (σ O 3 ) during 1982-2011 over vegetated land surface.

Table 1 .
List of v d simulations with input data.