Implementation of Yale Interactive terrestrial Biosphere model v1.0 into GEOS-Chem v12.0.0: a tool for biosphere–chemistry interactions

The terrestrial biosphere and atmospheric chemistry interact through multiple feedbacks, but the models of vegetation and chemistry are developed separately. In this study, the Yale Interactive terrestrial Biosphere (YIBs) model, a dynamic vegetation model with biogeochemical processes, is implemented into the Chemical Transport Model GEOS-Chem (GC) version 12.0.0. Within this GCYIBs framework, leaf area index (LAI) and canopy stomatal conductance dynamically predicted by YIBs are used for dry deposition calculation in GEOS-Chem. In turn, the simulated surface ozone (O3) by GEOS-Chem affect plant photosynthesis and biophysics in YIBs. The updated stomatal conductance and LAI improve the simulated O3 dry deposition velocity and its temporal variability for major tree species. For daytime dry deposition velocities, the model-to-observation correlation increases from 0.69 to 0.76, while the normalized mean error (NME) decreases from 30.5 % to 26.9 % using the GC-YIBs model. For the diurnal cycle, the NMEs decrease by 9.1 % for Amazon forests, 6.8 % for coniferous forests, and 7.9 % for deciduous forests using the GC-YIBs model. Furthermore, we quantify the damaging effects of O3 on vegetation and find a global reduction of annual gross primary productivity by 1.5 %–3.6 %, with regional extremes of 10.9 %–14.1 % in the eastern USA and eastern China. The online GC-YIBs model provides a useful tool for discerning the complex feedbacks between atmospheric chemistry and the terrestrial biosphere under global change.


Introduction
The terrestrial biosphere interacts with atmospheric chemistry through the exchanges of trace gases, water, and energy (Hungate and Koch, 2015;Green et al., 2017). Emissions from the terrestrial biosphere, such as biogenic volatile organic compounds (BVOCs) and nitrogen oxides (NO x ) affect the formation of air pollutants and chemical radicals in the atmosphere (Kleinman, 1994;Li et al., 2019). Globally, the terrestrial biosphere emits ∼ 1100 Tg (1 Tg = 10 12 g) BVOC annually, which is approximately 10 times more than the total amount of VOC emitted worldwide from anthropogenic sources including fossil fuel combustion and industrial activities (Carslaw et al., 2010). Meanwhile, the biosphere acts as a major sink through the dry deposition of air pollutants, such as surface ozone (O 3 ) and aerosols (Petroff, 2005;Fowler et al., 2009;Park et al., 2014). Dry deposition ac-Y. Lei et al.: A tool for biosphere-chemistry interactions counts for ∼ 25 % of the total O 3 removed from the troposphere (Lelieveld and Dentener, 2000).
In turn, atmospheric chemistry can also affect the terrestrial biosphere (McGrath et al., 2015;Schiferl and Heald, 2018;Yue and Unger, 2018). Surface O 3 has a negative impact on plant photosynthesis and crop yields by reducing gas exchange and inducing phytotoxic damage to plant tissues (Van Dingenen et al., 2009;Wilkinson et al., 2012;Yue and Unger, 2014). Unlike O 3 , the effect of aerosols on vegetation is dependent on the aerosol concentrations. Moderate increase of aerosols in the atmosphere is beneficial to vegetation (Mahowald, 2011;Schiferl and Heald, 2018). The aerosol-induced enhancement in diffuse light results in more radiation reaching surface from all directions than solely from above. As a result, leaves in the shade or at the bottom of the canopy can receive more radiation and are able to assimilate more CO 2 through photosynthesis, leading to an increase of canopy productivity (Mercado et al., 2009;Yue and Unger, 2018). However, excessive aerosol loadings reduce canopy productivity because the total radiation is largely weakened (Alton, 2008;. Models are essential tools to understand and quantify the interactions between the terrestrial biosphere and atmospheric chemistry at the global and/or regional scales. Many studies have performed multiple global simulations with climate-chemistry-biosphere models to quantify the effects of air pollutants on the terrestrial biosphere (Mercado et al., 2009;Yue and Unger, 2015;Oliver et al., 2018;Schiferl and Heald, 2018). In contrast, very few studies have quantified the O 3 -induced biogeochemical and meteorological feedbacks to air pollution concentrations (Sadiq et al., 2017;Zhou et al., 2018). Although considerable efforts have been made, uncertainties in biosphere-chemistry interactions remain large because their two-way coupling is not adequately represented in the current generation of terrestrial biosphere models or global chemistry models. Global terrestrial biosphere models usually use prescribed O 3 and aerosol concentrations (Sitch et al., 2007;Mercado et al., 2009;Lombardozzi et al., 2012), and global chemistry models often apply fixed offline vegetation variables (Lamarque et al., 2013). For example, stomatal conductance, which plays a crucial role in regulating the water cycle and altering pollution deposition, responds dynamically to vegetation biophysics and environmental stressors at various spatiotemporal scales (Hetherington and Woodward, 2003;Franks et al., 2017). However, these processes are either missing or lack temporal variations in most current chemical transport models (Verbeke et al., 2015). The fully two-way coupling between biosphere and chemistry is necessary to better quantify the responses of ecosystems and pollution to global changes.
In this study, we develop the GC-YIBs model by implementing the Yale Interactive terrestrial Biosphere (YIBs) model version 1.0 (Yue and Unger, 2015) into the chemical transport model (CTM) GEOS-Chem version 12.0.0 (http://wiki.seas.harvard.edu/geos-chem/index. php/GEOS-Chem_12#12.0.0, last access: 10 January 2020). The GEOS-Chem (GC) model has been widely used in episode prediction (Cui et al., 2016), source attribution (D'Andrea et al., 2016;Dunker et al., 2017;Ni et al., 2018;Lu et al., 2019), future pollution projection Ramnarine et al., 2019), health risk assessment (Xie et al., 2019), and so on. The standard GC model uses prescribed vegetation parameters and as a result cannot depict the changes in chemical components due to biosphere-pollution interactions. The updated GC-YIBs model links atmospheric chemistry with biosphere in a two-way coupling such that changes in chemical components or vegetation will simultaneously feed back to influence the other systems. Here, we evaluate the dynamically simulated dry deposition and leaf area index (LAI) from GC-YIBs and examine the consequent impacts on surface O 3 . We also quantify the detrimental effects of O 3 on gross primary productivity (GPP) using instant pollution concentrations from the chemical module. For the first step, we focus on the coupling between O 3 and vegetation. The interactions between aerosols and vegetation will be developed and evaluated in the future. The next section describes the GC-YIBs model and the evaluation data. Section 3 compares simulated O 3 from GC-YIBs with that from the original GC models and explores the causes of differences. Section 4 quantifies the damaging effects of O 3 on global GPP using the GC-YIBs model. The last section summarizes progress and discusses the next steps in optimizing the GC-YIBs model.

Descriptions of the YIBs model
YIBs is a terrestrial vegetation model designed to simulate the land carbon cycle with dynamical prediction of LAI and tree height (Yue and Unger, 2015). The YIBs model applies the Farquhar et al. (1980) scheme to calculate leaf level photosynthesis, which is further upscaled to the canopy level by the separation of sunlit and shaded leaves (Spitters, 1986). The canopy is divided into an adaptive number of layers (typically 2-16) for light stratification. Sunlight is attenuated and becomes more diffusive when penetrating the canopy. The sunlit leaves can receive both direct and diffuse radiation, while the shaded leaves receive only diffuse radiation. The leaf-level photosynthesis, calculated as the sum of sunlit and shaded leaves, is then integrated over all canopy layers to derive the GPP of ecosystems.
The model considers nine plant functional types (PFTs), including evergreen needleleaf forest, deciduous broadleaf forest, evergreen broadleaf forest, shrubland, tundra, C 3 and C 4 grasses, and C 3 and C 4 crops. The satellite-based land types and cover fraction are aggregated into these nine PFTs and used as input (Fig. S1 in the Supplement). The initial soil carbon pool and tree height used in YIBs are from the 140-year spin-up processes (Yue and Unger, 2015). The YIBs is driven with hourly 2-D meteorology and 3-D soil variables (six layers) from the Modern-Era Retrospective analysis for Research and Applications, version 2 (MERRA2).
The YIBs uses the model of Ball and Berry (Baldocchi et al., 1987) to compute leaf stomatal conductance: where r s is the leaf stomatal resistance (s m −1 ); m is the empirical slope of the Ball-Berry stomatal conductance equation and is affected by water stress; c s is the CO 2 concentration at the leaf surface (µmol m −3 ); RH is the relative humidity of the atmosphere; b (m s −1 ) represents the minimum leaf stomatal conductance when net leaf photosynthesis (A net , µmol m −2 s −1 ) is 0. For different PFTs, appropriate photosynthetic parameters are derived from the Community Land Model (CLM; Bonan et al., 2011). The net leaf photosynthesis for C 3 and C 4 plants is computed based on well-established Michaelis-Menten enzymekinetics scheme (Farquhar et al., 1980;von Caemmerer and Farquhar, 1981): where J c , J e , and J s represent Rubisco-limited photosynthesis, RuBP-limited photosynthesis, and product-limited photosynthesis, respectively. R d is the rate of dark respiration. They are all parameterized as functions of the maximum carboxylation capacity (Collatz et al., 1991) and meteorological variables (e.g., temperature, radiation, and CO 2 concentrations). The YIBs model applies the LAI and carbon allocation schemes from the TRIFFID model (Cox, 2001;Clark et al., 2011). On the daily scale, canopy LAI is calculated as follows: where f represents the phenological factor controlled by meteorological variables (e.g., temperature, water availability, and photoperiod); LAI max represents the available maximum LAI related to tree height, which is dependent on the vegetation carbon content (C veg ). The C veg is calculated as follows: where C l , C r , and C w represent leaf, root, and stem carbon contents, respectively. And all carbon components are parameterized as the function of LAI max : where α represents the specific leaf carbon density; β and γ represent allometric parameters. The vegetation carbon con-tent C veg is updated every 10 days: where τ and ϕ represent partitioning parameter and litter fall rate, respectively; their calculation methods have been documented in Yue and Unger (2015). Net primary productivity (NPP) is calculated as the residue of subtracting autotrophic respiration (R a ) from GPP: In addition, the YIBs model implements the scheme for O 3 damage on vegetation proposed by Sitch et al. (2007). The scheme directly modifies photosynthesis using a semimechanistic parameterization, which in turn affects stomatal conductance. The O 3 damage factor is considered as the function of stomatal O 3 flux: where a represents the sensitivity to damage and T O 3 represents the O 3 flux threshold (µmol m −2 s −1 ). For a specific PFT, the values of coefficient a vary from low to high to represent a range of uncertainties for ozone vegetation damage ( where [O 3 ] represents O 3 concentrations at top of the canopy (µmol m −3 ); r a is aerodynamic resistance (s m −1 ); r b is boundary layer resistance (s m −1 ); r s represents stomatal resistance (s m −1 ). The Sitch et al. (2007) scheme within the YIBs framework has been well evaluated against hundreds of observations globally (Yue and Unger, 2018) and regionally Yue et al., 2016.

Descriptions of the GEOS-Chem model
GC is a global 3-D model of atmospheric compositions with fully coupled O 3 -NO x -hydrocarbon-aerosol chemical mechanisms (Gantt et al., 2015;Lee et al., 2017;Ni et al., 2018). In this study, we use GC version 12.0.0 driven by assimilated meteorology from MERRA2 with a horizontal resolution of 4 • latitude by 5 • longitude and 47 vertical layers from the surface to 0.01 hPa.
In GC, terrestrial vegetation modulates tropospheric O 3 mainly through LAI and canopy stomatal conductance, which affect both the sources and sinks of tropospheric O 3 through changes in BVOC emissions, soil NO x emissions, and dry deposition (Zhou et al., 2018). BVOC emissions are  calculated based on a baseline emission factor parameterized as the function of light, temperature, leaf age, soil moisture, LAI, and CO 2 inhibition within the Model of Emissions of Gases and Aerosols from Nature (MEGAN v2.1; Guenther et al., 2006). Soil NO x emission is computed based on the scheme of Hudman et al. (2012) and further modulated by a reduction factor to account for within-canopy NO x deposition (Rogers and Whitman, 1991). The dry deposition velocity (V d , m s −1 ) for O 3 is computed based on a resistance-inseries model within GC: where R a (m s −1 ) is the aerodynamic resistance representing the ability of the airflow to bring gases or particles close to the surface and is dependent mainly on the atmospheric turbulence structure and the height considered. R b (m s −1 ) is the boundary resistance driven by the characteristics of the surface (surface roughness) and gas/particle (molecular diffusivity). R a and R b are calculated from the meteorological variables of global climate models (GCMs; Jacob et al., 1992). The surface resistance R c is determined by the affinity of the surface for the chemical compound. For O 3 over vegetated regions, V d is mainly driven by R c (m s −1 ) during the daytime because the effects of R a and R b are generally small. Surface resistances R c are computed using the Wesely (1989) canopy model with some improvements, including explicit dependence of canopy stomatal resistances on LAI (Gao and Wesely, 1995) and direct/diffuse PAR within the canopy (Baldocchi et al., 1987): where R s is the stomatal resistance (s m −1 ), R m is the leaf mesophyll resistance (R m = 0 s m −1 for O 3 ), R lu is the upper canopy or leaf cuticle resistance, and R cl is the lower canopy resistance (s m −1 ). R s is calculated based on minimum stomatal resistance (r s , s m −1 ), solar radiation (G, W m −2 ), surface air temperature (T s , • C), and the molecular diffusivities (D H 2 O and D x ) for a specific gas x: In GC, the above parameters related to R c have prescribed values for 11 deposition land types: snow/ice, deciduous forest, coniferous forest, agricultural land, shrub/grassland, Amazon forest, tundra, desert, wetland, urban and water (Wesely, 1989;Jacob et al., 1992). Although the stomatal conductance scheme of Wesely (1989) has been widely used in chemical transport and climate models, considerable limits still exist because this scheme does not consider the response of stomatal conductance to phenology, CO 2 concentrations, or soil water availability (Rydsaa et al., 2016;Lin et al., 2017). Previous studies have well evaluated the dry deposition scheme used in the GEOS-Chem model against observations globally and regionally (Hardacre et al., 2015;Silva and Heald, 2018;Wong et al., 2019). They found that GEOS-Chem can generally capture the diurnal and seasonal cycles except for the amplitude of O 3 dry deposition velocity (Silva and Heald, 2018).

Implementation of YIBs into GEOS-Chem (GC-YIBs)
In this study, GC model time steps are set to 30 min for transport and convection and 60 min for emissions and chemistry.  precursor emissions and dry deposition at the 1 h integration time step. The above processes are summarized in Fig. 1.
To retain the corresponding relationship between vegetation parameters and land cover map in the GC-YIBs model, we replace the Olson 2001 land cover map in GC with satellite-retrieved land cover dataset used by YIBs (Defries et al., 2000;Hanninen and Kramer, 2007). The conversion relationships between YIBs land types and GC deposition land types are summarized in Table S2. The global spatial pattern of deposition land types converted from YIBs land types is shown in Fig. S2. The Olson 2001 land cover map used in GC version 12.0.0 has a native resolution of 0.25 • × 0.25 • and 74 land types (Olson et al., 2001). Each of the Olson land types is associated with a corresponding deposition land type with prescribed parameters. There are 74 Olson land types but only 11 deposition land types, suggesting that many of the Olson land types share the same deposition parame-ters. At specific grids (4 • × 5 • or 2 • × 2.5 • ), dry deposition velocity is calculated as the weighted sum of native resolution (0.25 • × 0.25 • ). Replacing of Olson with YIBs land types induces a global mean difference of −0.59 ppbv on surface [O 3 ] (Fig. S3). Large discrepancies are found in Africa and the southern Amazon, where the local [O 3 ] decreases by more than 2 ppbv with the new land types. However, limited differences are shown in the middle-to-high latitudes of the Northern Hemisphere (NH, Fig. S3).

Model simulations
We conduct six simulations to evaluate the performance of GC-YIBs and to quantify global O 3 damage to vegetation (Table 1): (i) Offline, a control run using the offline GC-YIBs model. The YIBs module shares the same meteorological forcing as the GC module and predicts both GPP and     LAI. However, predicted vegetation variables are not fed into GC, which is instead driven by prescribed LAI from Moderate Resolution Imaging Spectroradiometer (MODIS) product and parameterized canopy stomatal conductance proposed by Gao and Wesely (1995). (ii) Online_LAI is a sensitive run using online GC-YIBs with dynamically predicted daily LAI from YIBs but prescribed stomatal conductance. (iii) On-line_GS is another sensitive run using YIBs-predicted stomatal conductance but prescribed MODIS LAI. (iv) In On-line_ALL, both YIBs-predicted LAI and stomatal conductance are used for GC. (v) Online_ALL_HS is the same as Online_ALL except it predicts surface O 3 damage to plant photosynthesis with high sensitivities. (vi) Online_ALL_ LS is the same as Online_ALL_HS but with low O 3 damage sensitivities. Each simulation is run from 2006 to 2012 with the first 4 years for spin up, and the results from 2010 to 2012 are used to evaluate the online GC-YIBs model. The differences between Online_ALL and Online_GS (Online_LAI) represent the effects of coupled LAI (stomatal conductance) on simulated [O 3 ]. Differences between Offline and On-line_ALL then represent joint effects of coupled LAI and stomatal conductance. The last three runs are used to quantify the global O 3 damage on ecosystem productivity.

Evaluation data
We use observed LAI data for 2010-2012 from the MODIS product. Benchmark GPP product of 2010-2012 is estimated by upscaling ground-based FLUXNET eddy covariance data using a model tree ensemble approach, a type of machine learning technique (Jung et al., 2009). Although these products may have certain biases, they have been widely used to evaluate land surface models because direct observations of GPP and LAI are not available on the global scale (Yue and Unger, 2015;Slevin et al., 2017;Swart et al., 2019). Measurements of surface [O 3 ] over North America and Europe are provided by the global gridded surface ozone data set of Sofen et al. (2016), and those over China are interpolated from data at ∼ 1500 sites operated by China's Ministry of Ecology and Environment (http://www.cnemc.cn/en/, last access: 10 January 2020). We perform literature research to collect data of dry deposition velocity from eight deciduous forest, two Amazon forest, and nine coniferous forest sites (Table 2).

Evaluation of the offline GC-YIBs model
With the offline simulation, the simulated GPP and LAI are compared with observed LAI and benchmark GPP for the period of 2010-2012 (Fig. 2). Observed LAI and benchmark GPP both show high values in the tropics and medium values in the northern middle-to-high latitudes. Compared to observations, the GC-YIBs model forced with MERRA2 meteorology depicts similar spatial distributions, with spatial correlation coefficients of 0.83 (p < 0.01) for GPP and 0.86 (p < 0.01) for LAI. Although the model overestimates LAI in the tropics and northern high latitudes by 1-2 m 2 m −2 , the simulated global area-weighted LAI (1.42 m 2 m −2 ) is close to observations (1.33 m 2 m −2 ), with a normalized mean bias (NMB) of 6.7 %. Similar to LAI, the global NMB for GPP is only 7.1 %, though there are substantial regional biases, especially in the Amazon and central Africa. Such differences are in part attributed to the underestimation of GPP for tropical rainforests in the benchmark product, because the recent simulations at eight rainforest sites with YIBs model reproduced ground-based observations well (Yue and Unger, 2018).
We then evaluate simulated annual mean surface [O 3 ] during 2010-2012 based on an offline simulation (Fig. 3). The simulated high values are mainly located in the mid-latitudes of NH (Fig. 3a). Compared to observations, simulations show reasonable spatial distribution with a correlation coefficient of 0.63 (p < 0.01). Although the offline GC-YIBs model overestimates annual [O 3 ] in southern China and predicts lower values in western Europe and western USA, the sim-  (Fig. S4), consistent with previous evaluations using the GC model (Travis et al., 2016;Schiferl and Heald, 2018;Yue and Unger, 2018).

Changes of surface O 3 in the online GC-YIBs model
Surface O 3 is changed by the coupling of LAI and stomatal conductance (Fig. 4). Global [O 3 ] shows similar patterns between Offline (Fig. 3a) and Online_ALL (Fig. 4a) (Fig. 4b), suggesting that changes in stomatal conductance play the dominant role in regulating surface [O 3 ]. As a comparison, [O 3 ] values between Online_ALL and Online_GS show limited changes globally (by 0.05 ppbv) and moderate changes in tropical regions (Fig. 4d), mainly because the LAI predicted by YIBs is close to MODIS LAI used in GC (Fig. 2). It is noticed that the average [O 3 ] in Fig. 4b is not equal to the sum of Fig. 4c and d, because of the non-linear effects.
We further explore the possible causes of differences in simulated [O 3 ] between online and offline GC-YIBs models. Figure 5 shows simulated annual O 3 dry deposition velocity from the online GC-YIBs model and its changes in dif-ferent sensitivity experiments. The global average velocity is 0.25 cm s −1 with a regional maximum of 0.5-0.7 cm s −1 in tropical rainforests (Fig. 5a), especially over the Amazon and central Africa, where high ecosystem productivity is observed (Fig. 2). With implementation of YIBs into GC, simulated dry deposition velocity increases over tropical regions but decreases in the middle-to-high latitudes of NH (Fig. 5b). Larger dry deposition results in lower [O 3 ] in the tropics, while smaller dry deposition increases [O 3 ] in boreal regions. Such spatial patterns are broadly consistent with [O 3 ] in online GC-YIBs (Fig. 4b). In a comparison, updated LAI induces limited changes in the isoprene and NO x emissions (Fig. S5), suggesting that changes of dry deposition velocity are the dominant drivers of O 3 changes. Both the updated LAI and stomatal conductance influence dry deposition. Sensitivity experiments further show that changes in dry deposition are mainly driven by coupled canopy stomatal conductance (Fig. 5c) instead of LAI (Fig. 5d), though the latter contributes to the enhanced dry deposition in the tropics.
The original GC dry deposition scheme applies fixed parameters for stomatal conductance of a specific land type. The updated GC-YIBs model instead calculates stomatal conductance as a function of photosynthesis and environmental forcings (Eq. 1). As a result, predicted dry deposition exhibits discrepancies among biomes. With Offline and On-line_ALL simulations, we further evaluate the performance of online GC-YIBs in simulating O 3 dry deposition velocity for specific deposition land types (Fig. 6). For agricultural land and shrub/grassland, the simulated O 3 dry deposition velocity for online GC-YIBs model is close to the GC model, with NMBs of 3 % and −2 % and correlation coefficients of 0.96 and 0.97, respectively. However, the simulated dry deposition velocity in online GC-YIBs is lower than GC by 18 % for deciduous forests and 14 % for coniferous forests, The blue and red lines represent simulated O 3 dry deposition velocity by GC (Offline simulation) and online GC-YIBs (Online_ALL simulation) models, respectively. The site number (N), R, and NME are shown for each panel.
but larger by 17 % for Amazon forests. Such changes match the spatial pattern of dry deposition shown in Fig. 5b.
Since the changes of O 3 dry deposition velocity are mainly found in deciduous forests, coniferous forests, and Amazon forests, we collect 27 samples across these three biomes to evaluate the online GC-YIBs model (Table 2). For the 11 samples in deciduous forests, the normalized mean error (NME) decreases from 29 % in the GC model to 24 % in GC-YIBs with lower relative errors at eight sites (Fig. 7). Predictions with the GC-YIBs also show large improvements over coniferous forests, where 8 out of 14 samples show lower (decreases from 27 % in GC to 25 % in GC-YIBs) errors. For Amazon forests, the GC-YIBs model significantly improves the prediction at one site (117.9 • E, 4.9 • N), where the original error of −0.17 cm s −1 is limited to only 0.03 cm s −1 . However, the new model does not improve the prediction at the other Amazon forest site. Overall, the simulated daytime O 3 dry deposition velocities in online GC-YIBs model are closer to observations than those in the GC model with smaller NME (26.9 % vs. 30.5 %), root-meansquare errors (RMSEs, 0.2 vs. 0.23) and higher correlation coefficients (0.76 vs. 0.69). Such improvements consolidate our strategies in updating GC model to the fully coupled GC-YIBs model.
We collect long-term measurements from four sites across North America and western Europe to evaluate the model performance in simulating seasonal cycle of O 3 dry deposition velocity (Fig. 8). The GC model well captures the seasonal cycles of O 3 dry deposition velocity in all sites with the correlation coefficients of 0.95 at Harvard, 0.8 at Hyytiälä, 0.68 at Ulborg, and 0.71 at Auchencorth. However, the magnitude of O 3 dry deposition velocity is overestimated at the Harvard and Hyytiälä sites (NME of 60 % and 42 %, respectively), but underestimated at the Ulborg and Auchencorth sites (NME of 48.7 % and 58.9 %, respectively) at growing seasons. Compared to the GC model, simulated O 3 dry deposition velocity with the GC-YIBs model shows large improvements over Harvard (Hyytiälä), where the model-to-observation NME decreases from 60 % (42 %) to 32 % (28 %).
Additionally, we investigate the diurnal cycle of O 3 dry deposition velocity at 15 sites (Fig. S6). Observed O 3 dry deposition velocities show a single diurnal peak with the maximum from 08:00 to 16:00 local time (Fig. 9). Compared to observations, the GC model has good performance in simulating the diurnal cycle with correlation coefficients of 0.94 for Amazon forests, 0.96 for coniferous forests, and 0.95 for deciduous forests. The GC model underestimates daytime O 3 dry deposition velocity for Amazon forests (NME of 29.8 %) but overestimates it for coniferous and deciduous forests (NME of 21.9 % and 22.9 %, respectively). Compared to the GC model, the simulated daytime O 3 dry deposition velocities using the GC-YIBs model are closer to observations in all three biomes. The NMEs decrease by 9.1 % for Amazon forests, 6.8 % for coniferous forests, and 7.9 % for deciduous forests.

Assessment of global O 3 damage to vegetation
An important feature of GC-YIBs is the inclusion of online vegetation damage by surface O 3 . Here, we quantify the global O 3 damage to GPP and LAI by conducting On-line_ALL_HS and Online_ALL_LS simulations (Fig. 10). Due to O 3 damage, annual GPP declines from −1.5 % (low sensitivity) to −3.6 % (high sensitivity) on the global scale. Regionally, O 3 decreases GPP by as much as 10.9 % in the eastern USA and up to 14.1 % in eastern China at high sensitivity (Fig. 10a, b). Such strong damage is related to (i) high ambient [O 3 ] due to anthropogenic emissions and (ii) large stomatal conductance due to active ecosystem productivity in monsoon areas. The O 3 effects are moderate in tropical areas, where stomatal conductance is also high, while [O 3 ] is very low (Fig. 4a) due to limited anthropogenic emissions. Furthermore, O 3 -induced GPP reductions are also small in the western USA and western Asia. Although [O 3 ] is high over these semi-arid regions (Fig. 4a), the drought stress decreases stomatal conductance and consequently constrains the O 3 uptake. The damage to LAI (Fig. 10c, d) generally follows the pattern of GPP reductions (Fig. 10a, b) but with lower magnitude. The reductions of GPP are slightly higher than our previous estimates using prescribed LAI and/or surface [O 3 ] in the simulations (Yue andUnger, 2014, 2015), likely because GC-YIBs considers O 3 -vegetation interactions. The feedback of such an interaction to both chemistry and biosphere will be explored in future studies.

Conclusions and discussion
The terrestrial biosphere and atmospheric chemistry interact through a series of feedbacks (Green et al., 2017). Among biosphere-chemistry interactions, dry deposition plays a key role in the exchange of compounds and acts as an important sink for several air pollutants (Verbeke et al., 2015). However, dry deposition is simply parameterized in most current CTMs (Hardacre et al., 2015). For all chemical species considered in the GC model, stomatal resistance R c is simply calculated as the function of minimum stomatal resistance and meteorological forcings. Such parameterization not only induces biases, but it also ignores the feedbacks from biosphere-chemistry interactions. For example, recent studies revealed that O 3 -induced damage to vegetation could reduce stomatal conductance and in turn alter ambient O 3 level (Sadiq et al., 2017;Zhou et al., 2018). In this study, we implement YIBs into the GC model with fully interactive surface O 3 and the terrestrial biosphere. The dynamically predicted LAI and stomatal conductance from YIBs are instantly provided to GC, meanwhile the prognostic O 3 simulated by GC simultaneously affects vegetation biophysics in YIBs. With these updates, simulated O 3 dry deposition velocities and their temporal variability (seasonal and diurnal cycles) in GC-YIBs are closer to observations than those in the original GC model.
An earlier study updated the dry deposition scheme in the Community Earth System Model (CESM) by implementing the leaf and stomatal resistances (Val Martin et al., 2014). Compared to that work, the magnitudes of [O 3 ] in our simulations are smaller in North America, eastern Europe, and southern China. This might be because the original dry deposition scheme in the GC model (see validation in Fig. 7) is better than that in CESM, leaving limited potential for improvements. In GC, the leaf cuticular resistance (R lu ) is dependent on LAI (Gao and Wesely, 1995), while the original calculation of R lu in CESM does not include LAI (Wesely, 1989). In addition, differences in the canopy schemes for stomatal conductance between YIBs and the Community Land Model (CLM) may cause different responses in dry deposition, which is changed by −0.12 to 0.16 cm s −1 in GC-YIBs but is much larger, by −0.15 to 0.25 cm s −1 , in CESM (Val Martin et al., 2014). Moreover, the GC-YIBs is driven with prescribed reanalysis, while CESM dynamically predicts climatic variables. Perturbations of meteorology in response to terrestrial properties may further magnify the variations in the atmospheric components of CESM.
Although we implement YIBs into GC with fully a interactive surface O 3 and terrestrial biosphere, it should be noted that considerable limits still exist, and further developments are required for GC-YIBs.
1. Atmospheric nitrogen alters plant growth and further influences both the sources and sinks of surface O 3 through surface-atmosphere exchange processes (Zhao et al., 2017). However, the YIBs model currently utilizes a fixed nitrogen level and does not include an interactive nitrogen cycle, which may induce uncertainties in simulating carbon fluxes. 4. The current GC-YIBs is limited to a low resolution due to slow computational speed and high computational costs for long-term integrations. The GC model, even at the 2 • × 2.5 • resolution, takes days to simulate 1 model year due to comprehensive parameterizations of physical and chemical processes. Such a low speed constrains the long-term spin up required by dynamical vegetation models. The low resolution will affect local emissions (e.g., NO x and VOC) and transport, leading to changes in surface [O 3 ] in GEOS-Chem. The comparison results of 2007 show that the low resolution of 4 • × 5 • induces a global mean bias of −0.24 ppbv on surface [O 3 ] compared to the relatively high resolution at 2 • × 2.5 • (Fig. S7). Compared with surface [O 3 ], low resolution causes limited differences in vegetation variables (e.g., GPP and LAI, not shown).
Despite these deficits, the development of GC-YIBs provides a unique tool for studying biosphere-chemistry interactions. In the future, we will extend our applications via the following steps.
1. Air pollution impacts on biosphere, including both O 3 and aerosol effects will be included. The GC-YIBs model can predict atmospheric aerosols, which affect both direct and diffuse radiation through the Rapid Radiative Transfer Model for GCMs (RRTMG) in the GC module (Schiferl and Heald, 2018). The diffuse fertilization effects in the YIBs model have been fully evaluated (Yue and Unger, 2018), and as a result we can quantify the impacts of aerosols on terrestrial ecosystems.
2. Multiple schemes for BVOC emissions will be added. The YIBs model incorporates both MEGAN (Guenther et al., 2006) and photosynthesis-dependent (Unger, 2013) isoprene emission schemes (Yue and Unger, 2015). The two schemes within the GC-YIBs framework can be used and compared for simulations of BVOC and consequent air pollution (e.g., O 3 , secondary organic aerosols).
3. We will include biosphere-chemistry feedbacks to air pollution. The effects of air pollution on the biosphere include changes in stomatal conductance, LAI, and BVOC emissions, which in turn modify the sources and sinks of atmospheric components. Only a few studies have quantified these feedbacks for O 3 -vegetation interactions (Sadiq et al., 2017;Zhou et al., 2018). We can explore the full biosphere-chemistry coupling for both O 3 and aerosols using the GC-YIBs model in the future.
Code availability. The YIBs model was developed by Xu Yue and Nadine Unger with code sharing at https://github.com/YIBS01/ YIBS_site (last access: 10 January 2020; Yue, 2015). The GEOS-Chem model was developed by the Atmospheric Chemistry Modeling Group at Harvard University led by Daniel Jacob and improved by a global community of atmospheric chemists. The source code for the GEOS-Chem model is publicly available at https: //github.com/geoschem/geos-chem (last access: 10 January 2020; GEOS-Chem, 2020). The source codes for the GC-YIBs model is archived at https://doi.org/10.5281/zenodo.3659346 (Lei and Yue, 2020).
Author contributions. XY conceived the study. YL and XY were responsible for model coupling, simulations, results analysis, and paper writing. All co-authors improved and prepared the manuscript.
Competing interests. The authors declare that they have no conflict of interest.