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 GC-YIBs 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 (O₃) by GEOS-Chem affect plant photosynthesis and biophysics in YIBs. The updated stomatal conductance and LAI improve the simulated O₃ 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 O₃ 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.

1 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¹² 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₃) and aerosols (Petroff, 2005; Fowler et al., 2009; Park et al., 2014). Dry deposition accounts for ∼25 % of the total O₃ 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₃ 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₃, 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₂ 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; Yue and Unger, 2017).

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₃-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₃ 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 (Yue et al., 2015; 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₃. We also quantify the detrimental effects of O₃ on gross primary productivity (GPP) using instant pollution concentrations from the chemical module. For the first step, we focus on the coupling between O₃ 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₃ from GC-YIBs with that from the original GC models and explores the causes of differences. Section 4 quantifies the damaging effects of O₃ on global GPP using the GC-YIBs model. The last section summarizes progress and discusses the next steps in optimizing the GC-YIBs model.

2 Methods and data

2.1 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₃ and C₄ grasses, and C₃ and C₄ 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:

$$\begin{array}{}\text{(1)}& g_{\mathrm{s}}={\displaystyle \frac{\mathrm{1}}{r_{\mathrm{s}}}}=m{\displaystyle \frac{A_{\mathrm{net}}}{c_{\mathrm{s}}}}\mathrm{RH}+b,\end{array}$$

where r_s is the leaf stomatal resistance (s m⁻¹); m is the empirical slope of the Ball–Berry stomatal conductance equation and is affected by water stress; c_s is the CO₂ concentration at the leaf surface (µmol m⁻³); RH is the relative humidity of the atmosphere; b (m s⁻¹) represents the minimum leaf stomatal conductance when net leaf photosynthesis (A_net, µmol m⁻² s⁻¹) 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₃ and C₄ plants is computed based on well-established Michaelis–Menten enzyme-kinetics scheme (Farquhar et al., 1980; von Caemmerer and Farquhar, 1981):

$$\begin{array}{}\text{(2)}& A_{\mathrm{net}}=min\left(J_{\mathrm{c}},J_{\mathrm{e}},J_{\mathrm{s}}\right)-R_{\mathrm{d}},\end{array}$$

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₂ 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:

$$\begin{array}{}\text{(3)}& \mathrm{LAI}=f\times {\mathrm{LAI}}_{\max },\end{array}$$

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:

$$\begin{array}{}\text{(4)}& C_{\mathrm{veg}}=C_{\mathrm{l}}+C_{\mathrm{r}}+C_{\mathrm{w}},\end{array}$$

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:

$$\begin{array}{}\text{(5)}& \left\{\begin{array}{l}{C}_{\mathrm{l}}=\mathit{\alpha }\times \mathrm{LAI}\\ C_{\mathrm{r}}=\mathit{\alpha }\times {\mathrm{LAI}}_{\max }\\ C_{\mathrm{r}}=\mathit{\beta }\times {\mathrm{LAI}}_{\max }^{\mathit{\gamma }}\end{array}\right.,\end{array}$$

where α represents the specific leaf carbon density; β and γ represent allometric parameters. The vegetation carbon content C_veg is updated every 10 days:

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{d}{C}_{\mathrm{veg}}}{\mathrm{d}t}}=\left(\mathrm{1}-\mathit{\tau }\right)\times \mathrm{NPP}-\mathit{\phi },\end{array}$$

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:

$$\begin{array}{}\text{(7)}& \mathrm{NPP}=\mathrm{GPP}-{\mathrm{R}}_{\mathrm{a}}.\end{array}$$

In addition, the YIBs model implements the scheme for O₃ damage on vegetation proposed by Sitch et al. (2007). The scheme directly modifies photosynthesis using a semi-mechanistic parameterization, which in turn affects stomatal conductance. The O₃ damage factor is considered as the function of stomatal O₃ flux:

$$\begin{array}{}\text{(8)}& F=\left\{\begin{array}{l}-a\left(F_{{\mathrm{O}}_{\mathrm{3}}}-T_{{\mathrm{O}}_{\mathrm{3}}}\right),F_{{\mathrm{O}}_{\mathrm{3}}}>T_{{\mathrm{O}}_{\mathrm{3}}}\\ \mathrm{0},F_{{\mathrm{O}}_{\mathrm{3}}}\le T_{{\mathrm{O}}_{\mathrm{3}}}\end{array}\right.,\end{array}$$

where a represents the sensitivity to damage and $T_{{\mathrm{O}}_{\mathrm{3}}}$ represents the O₃ flux threshold (µmol m⁻² s⁻¹). For a specific PFT, the values of coefficient a vary from low to high to represent a range of uncertainties for ozone vegetation damage (Table S1 in the Supplement). $T_{{\mathrm{O}}_{\mathrm{3}}}$ is a critical threshold for O₃ damage and varies with PFTs. The F becomes negative only if $F_{{\mathrm{O}}_{\mathrm{3}}}$ is higher than $T_{{\mathrm{O}}_{\mathrm{3}}}$. Stomatal O₃ flux $F_{{\mathrm{O}}_{\mathrm{3}}}$ (µmol m⁻² s⁻¹) is calculated as follows:

$$\begin{array}{}\text{(9)}& F_{{\mathrm{O}}_{\mathrm{3}}}={\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{3}}\right]}{r_{\mathrm{a}}+r_{\mathrm{b}}+k\cdot r_{\mathrm{s}}}},\end{array}$$

where [O₃] represents O₃ concentrations at top of the canopy (µmol m⁻³); r_a is aerodynamic resistance (s m⁻¹); r_b is boundary layer resistance (s m⁻¹); r_s represents stomatal resistance (s m⁻¹). 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, 2017).

Figure 1 Diagram of the GC-YIBs global carbon-chemistry model. Processes with red font are implemented in this study. Processes in the blue dashed box will be developed in the future.

Table 1 Summary of simulations using the GC-YIBs model.

Table 2 List of measurement sites used for dry deposition evaluation.

Figure 2 Annual gross primary productivity (GPP) and leaf area index (LAI) from offline simulations (a, b), observations (c, d), and their differences (e, f) averaged for the period of 2010–2012. Global area-weighted GPP and LAI are shown in parentheses in the panel captions. The correlation coefficients (R) and global normalized mean biases (NMBs) are shown in the bottom panels.

2.2 Descriptions of the GEOS-Chem model

GC is a global 3-D model of atmospheric compositions with fully coupled O₃–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₃ mainly through LAI and canopy stomatal conductance, which affect both the sources and sinks of tropospheric O₃ through changes in BVOC emissions, soil NO_xemissions, 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₂ 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⁻¹) for O₃ is computed based on a resistance-in-series model within GC:

$$\begin{array}{}\text{(10)}& V_{\mathrm{d}}={\displaystyle \frac{\mathrm{1}}{R_{\mathrm{a}}+R_{\mathrm{b}}+R_{\mathrm{c}}}},\end{array}$$

where R_a (m s⁻¹) 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⁻¹) 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₃ over vegetated regions, V_d is mainly driven by R_c (m s⁻¹) 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):

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\mathrm{1}}{R_{\mathrm{c}}}}={\displaystyle \frac{\mathrm{1}}{R_{\mathrm{s}}+R_{\mathrm{m}}}}+{\displaystyle \frac{\mathrm{1}}{R_{\mathrm{lu}}}}+{\displaystyle \frac{\mathrm{1}}{R_{\mathrm{cl}}}}+{\displaystyle \frac{\mathrm{1}}{R_{\mathrm{g}}}},\end{array}$$

where R_s is the stomatal resistance (s m⁻¹), R_m is the leaf mesophyll resistance (R_m=0 s m⁻¹ for O₃), R_lu is the upper canopy or leaf cuticle resistance, and R_cl is the lower canopy resistance (s m⁻¹). R_s is calculated based on minimum stomatal resistance (r_s, s m⁻¹), solar radiation (G, W m⁻²), surface air temperature (T_s, ^∘C), and the molecular diffusivities ($D_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ and D_x) for a specific gas x:

$$\begin{array}{}\text{(12)}& R_{\mathrm{s}}=r_{\mathrm{s}}\left\{\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{{\left[\mathrm{200}\left(G+\mathrm{0.1}\right)\right]}^{\mathrm{2}}}}\right\}\left\{{\displaystyle \frac{\mathrm{400}}{T_{\mathrm{s}}\left(\mathrm{40}-T_{\mathrm{s}}\right)}}\right\}{\displaystyle \frac{D_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}{D_x}}.\end{array}$$

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₂ 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; Lin et al., 2019; Wong et al., 2019). They found that GEOS-Chem can generally capture the diurnal and seasonal cycles except for the amplitude of O₃ dry deposition velocity (Silva and Heald, 2018).

Figure 3 Annual surface O₃ concentrations ([O₃]) from offline simulations (a), observations (b), and their differences (c) averaged for the period of 2010–2012. Global area-weighted surface [O₃] over grids with available observations are shown in parentheses in the panel captions. The correlation coefficient (R) and global normalized mean biases (NMBs) are shown in the bottom panel with indication of grid numbers (N) used for statistics.

Figure 4 Simulated annual surface [O₃] from online GC-YIBs model (a) and its changes (b–d) relative to offline simulations. Changes of [O₃] are caused by (b) jointly coupled LAI and stomatal conductance (Online_ALL – Offline), (c) coupled stomatal conductance alone (Online_ALL – Online_LAI), and (d) coupled LAI alone (Online_ALL – Online_GS). Global area-weighted [O₃] or Δ [O₃] are shown in the panels.

Figure 5 Simulated annual O₃ dry deposition velocity from online GC-YIBs model (a) and its changes caused by coupled LAI and stomatal conductance (b–d) averaged for the period of 2010–2012. The changes of dry deposition velocity are driven by (b) coupled LAI and stomatal conductance (Online_ALL – Offline), (c) coupled stomatal conductance alone (Online_ALL – Online_LAI), and (d) coupled LAI alone (Online_ALL – Online_GS). Global area-weighted annual O₃ dry deposition velocity and changes are shown in the panels.

Figure 6 Comparisons of annual O₃ dry deposition velocity between online GC-YIBs (Online_ALL simulation) and GC (Offline simulation) models for different land types, including (a) deciduous forest (DF), (b) coniferous forest (CF), (c) agricultural land (AL), (d) shrub/grassland (SG), and (e) Amazon forest (AF). The box plots of dry deposition velocity simulated by online GC-YIBs (blue) and GC models (red) for different land types are shown in (f). Each point in (a–e) represents annual O₃ dry deposition velocity at one grid point averaged for the period of 2010–2012. The red lines indicate linear regressions between predictions from GC-YIBs and GC models. The regression fit, correlation coefficient (R), and normalized mean biases (NMBs) are shown in each panel.

2.3 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. In the online GC-YIBs configuration, GC provides the hourly meteorology, aerodynamic resistance, boundary layer resistance, and surface [O₃] to YIBs. Without YIBs implementation, the GC model computes O₃ dry deposition velocity using prescribed LAI and parameterized canopy stomatal resistance (R_s), and as a result ignores feedbacks from ecosystems (details in Sect. 2.2). With YIBs embedded, daily LAI and hourly stomatal conductance are dynamically predicted for the dry deposition scheme within the GC model. The online-simulated surface [O₃] affects carbon assimilation and canopy stomatal conductance; in turn, the online-simulated vegetation variables such as LAI and stomatal conductance affect both the sources and sinks of O₃ by altering 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 $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ 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 parameters. At specific grids ($\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ or $\mathrm{2}{}^{\circ }\times \mathrm{2.5}{}^{\circ }$), dry deposition velocity is calculated as the weighted sum of native resolution ($\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$). Replacing of Olson with YIBs land types induces a global mean difference of −0.59 ppbv on surface [O₃] (Fig. S3). Large discrepancies are found in Africa and the southern Amazon, where the local [O₃] 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).

Figure 7 Comparison between observed and simulated O₃ dry deposition velocity at the observational sites. The different marker types represent different land types. The blue and red markers represent the simulation results from online GC-YIBs (Online_ALL simulation) and GC (Offline simulation) models, respectively. The blue and red lines indicate linear regressions between simulations and observations. The regression fits, root-mean-square errors (RMSEs), normalized mean errors (NMEs), and correlation coefficients (R) for GC-YIBs (blue) and GC (red) models are also shown.

Figure 8 Comparison of monthly O₃ dry deposition velocity at the Harvard (a), Ulborg (b), Hyytiälä (c), and Auchencorth (d) sites. The black lines represent observed O₃ dry deposition velocity. The blue and red lines represent simulated O₃ dry deposition velocity from GC (Offline simulation) and online GC-YIBs (Online_ALL simulation) models, respectively.

2.4 Model simulations

We conduct six simulations to evaluate the performance of GC-YIBs and to quantify global O₃ 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) Online_GS is another sensitive run using YIBs-predicted stomatal conductance but prescribed MODIS LAI. (iv) In Online_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₃ damage to plant photosynthesis with high sensitivities. (vi) Online_ALL_ LS is the same as Online_ALL_HS but with low O₃ 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₃]. Differences between Offline and Online_ALL then represent joint effects of coupled LAI and stomatal conductance. The last three runs are used to quantify the global O₃ damage on ecosystem productivity.

2.5 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₃] 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).

Figure 9 Comparison of multi-site mean diurnal cycle of O₃ dry deposition velocity at the Amazon (a), coniferous (b), and deciduous (c) forests. Error bars represent the range of values from different sites. Black lines represent observed O₃ dry deposition velocity. The blue and red lines represent simulated O₃ 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.

Figure 10 Percentage changes in (a, b) GPP and (c, d) LAI caused by the damaging effects of O₃ with (a, c) low (Online_ALL_LS simulation) and (b, d) high sensitivities (Online_ALL_HS simulation). Both changes of GPP and LAI are averaged for 2010–2012.

3 Results

3.1 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² m⁻², the simulated global area-weighted LAI (1.42 m² m⁻²) is close to observations (1.33 m² m⁻²), 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₃] 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₃] in southern China and predicts lower values in western Europe and western USA, the simulated area-weighted surface [O₃] (45.4 ppbv) is only 6 % higher than observations (42.8 ppbv). Predicted summertime surface [O₃] instead shows positive biases in eastern USA and Europe (Fig. S4), consistent with previous evaluations using the GC model (Travis et al., 2016; Schiferl and Heald, 2018; Yue and Unger, 2018).

3.2 Changes of surface O₃ in the online GC-YIBs model

Surface O₃ is changed by the coupling of LAI and stomatal conductance (Fig. 4). Global [O₃] shows similar patterns between Offline (Fig. 3a) and Online_ALL (Fig. 4a) simulations. However, the online GC-YIBs predicts [O₃] of 0.5–2 ppbv higher in the middle-to-high latitudes of NH, leading to an average [O₃] enhancement of 0.22 ppbv compared to offline simulations (Fig. 4b). Regionally, some negative changes of 1–2 ppbv can be found at the tropical regions. With sensitivity experiments Online_LAI and Online_GS (Table 1), we separate the contributions of LAI and stomatal conductance changes to Δ [O₃]. It is found that Δ [O₃] between Online_ALL and Online_LAI (Fig. 4c) resembles the total Δ [O₃] pattern (Fig. 4b), suggesting that changes in stomatal conductance play the dominant role in regulating surface [O₃]. As a comparison, Δ [O₃] 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₃] 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₃] between online and offline GC-YIBs models. Figure 5 shows simulated annual O₃ dry deposition velocity from the online GC-YIBs model and its changes in different sensitivity experiments. The global average velocity is 0.25 cm s⁻¹ with a regional maximum of 0.5–0.7 cm s⁻¹ 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₃] in the tropics, while smaller dry deposition increases [O₃] in boreal regions. Such spatial patterns are broadly consistent with Δ [O₃] 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₃ 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 Online_ALL simulations, we further evaluate the performance of online GC-YIBs in simulating O₃ dry deposition velocity for specific deposition land types (Fig. 6). For agricultural land and shrub/grassland, the simulated O₃ 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, 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₃ 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⁻¹ is limited to only 0.03 cm s⁻¹. However, the new model does not improve the prediction at the other Amazon forest site. Overall, the simulated daytime O₃ 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-mean-square 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₃ dry deposition velocity (Fig. 8). The GC model well captures the seasonal cycles of O₃ 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₃ 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₃ 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₃ dry deposition velocity at 15 sites (Fig. S6). Observed O₃ 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₃ 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₃ 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.

3.3 Assessment of global O₃ damage to vegetation

An important feature of GC-YIBs is the inclusion of online vegetation damage by surface O₃. Here, we quantify the global O₃ damage to GPP and LAI by conducting Online_ALL_HS and Online_ALL_LS simulations (Fig. 10). Due to O₃ damage, annual GPP declines from −1.5 % (low sensitivity) to −3.6 % (high sensitivity) on the global scale. Regionally, O₃ 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₃] due to anthropogenic emissions and (ii) large stomatal conductance due to active ecosystem productivity in monsoon areas. The O₃ effects are moderate in tropical areas, where stomatal conductance is also high, while [O₃] is very low (Fig. 4a) due to limited anthropogenic emissions. Furthermore, O₃-induced GPP reductions are also small in the western USA and western Asia. Although [O₃] is high over these semi-arid regions (Fig. 4a), the drought stress decreases stomatal conductance and consequently constrains the O₃ 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₃] in the simulations (Yue and Unger, 2014, 2015), likely because GC-YIBs considers O₃–vegetation interactions. The feedback of such an interaction to both chemistry and biosphere will be explored in future studies.

4 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₃-induced damage to vegetation could reduce stomatal conductance and in turn alter ambient O₃ level (Sadiq et al., 2017; Zhou et al., 2018). In this study, we implement YIBs into the GC model with fully interactive surface O₃ and the terrestrial biosphere. The dynamically predicted LAI and stomatal conductance from YIBs are instantly provided to GC, meanwhile the prognostic O₃ simulated by GC simultaneously affects vegetation biophysics in YIBs. With these updates, simulated O₃ 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₃] 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⁻¹ in GC-YIBs but is much larger, by −0.15 to 0.25 cm s⁻¹, 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₃ 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₃ 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.

• 2.

  The validity of Δ [O₃], especially those at high latitudes in NH, cannot be directly evaluated due to a lack of measurements. Although changes of dry deposition show improvements in GC-YIBs, the ultimate effects on surface [O₃] remain unclear within the original GC framework.

• 3.

  [O₃] at the lowest model level is used as an approximation of canopy [O₃]. The current model does not include a sub-grid parameterization of pollution transport within the canopy, leading to biases in estimating O₃ vegetation damage and the consequent feedback. However, development of such parameterization is limited by the availability of simultaneous measurements of microclimate and air pollutants.

• 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 $\mathrm{2}{}^{\circ }\times \mathrm{2.5}{}^{\circ }$ 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₃] in GEOS-Chem. The comparison results of 2007 show that the low resolution of $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ induces a global mean bias of −0.24 ppbv on surface [O₃] compared to the relatively high resolution at $\mathrm{2}{}^{\circ }\times \mathrm{2.5}{}^{\circ }$ (Fig. S7). Compared with surface [O₃], 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₃ 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₃, 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₃–vegetation interactions (Sadiq et al., 2017; Zhou et al., 2018). We can explore the full biosphere–chemistry coupling for both O₃ and aerosols using the GC-YIBs model in the future.