Implementation of trait-based ozone plant sensitivity in the Yale Interactive terrestrial Biosphere model v1.0 to assess global vegetation damage

A major limitation in modeling global ozone (O₃) vegetation damage has long been the reliance on empirical O₃ sensitivity parameters derived from a limited number of species and applied at the level of plant functional types (PFTs), which ignore the large interspecific variations within the same PFT. Here, we present a major advance in large-scale assessments of O₃ plant injury by linking the trait leaf mass per area (LMA) and plant O₃ sensitivity in a broad and global perspective. Application of the new approach and a global LMA map in a dynamic global vegetation model reasonably represents the observed interspecific responses to O₃ with a unified sensitivity parameter for all plant species. Simulations suggest a contemporary global mean reduction of 4.8 % in gross primary productivity by O₃, with a range of 1.1 %–12.6 % for varied PFTs. Hotspots with damage >10 % are found in agricultural areas in the eastern US, western Europe, eastern China, and India, accompanied by moderate to high levels of surface O₃. Furthermore, we simulate the distribution of plant sensitivity to O₃, which is highly linked with the inherent leaf trait trade-off strategies of plants, revealing high risks for fast-growing species with low LMA, such as crops, grasses, and deciduous trees.

1 Introduction

Tropospheric ozone (O₃) has long been recognized as a hazardous pollutant for plants (Richards et al., 1958; Reich and Amundson, 1985). As a strong oxidant, O₃ can cause damage to leaf cells (Feng et al., 2014), impact stomata conductance (Buker et al., 2015), and reduce photosynthesis and biomass (Wittig et al., 2009). These negative impacts dampen global plant productivity (Ainsworth et al., 2012, 2020) and crop yield (Tai et al., 2014; Emberson et al., 2018; Feng et al., 2022), altering multiple ecosystem functions and services across various spatiotemporal scales (Agathokleous et al., 2020; Feng et al., 2021). Thus, it is of crucial importance to quantify O₃ plant damage in global modeling and assess its coupling effects in biosphere–atmosphere systems (Zhou et al., 2018).

To date, O₃ fumigation experiments have been conducted for various plant species. Accordingly, O₃ damage sensitivities, denoted as the dose–response relationships (DRRs), were derived as regressions between O₃ exposure metrics and changes in biotic indicators (Mills et al., 2011). The widely used O₃ metrics include ambient O₃ concentrations for AOT40 (accumulated O₃ concentration above the threshold of 40 ppbv; Fuhrer et al., 1997) or the stomatal O₃ flux for POD_y (phytotoxic O₃ dose above a threshold flux of y; Buker et al., 2015). The biotic indicators include visual leaf states, photosynthetic rate, biomass, or crop yield. Normally, the DRRs were derived for typical tree/grass species in specific regions, for example, Norway spruce, birch, and beech in Europe (Buker et al., 2015) or poplar (Shang et al., 2017) and crops (Peng et al., 2019) in East Asia.

Some assessment studies used DRRs to derive contemporary O₃ plant damage patterns at large scales. Concentration-based DRRs were widely measured and applied on the homogenized land cover, mostly for estimating crop yield loss (Feng et al., 2022; Tai et al., 2021; Hong et al., 2020). However, such DRRs do not include information about biochemical defense and stomatal regulation. In comparison, flux-based DRRs reflect a more detailed consideration in biological processes but are limited by the application scales in both space and time (Mills et al., 2011, 2018b). For example, the estimate of POD_y needs a dry deposition model, “DO₃SE” (Deposition of Ozone for Stomatal Exchange) (LRTAP Convention, 2017), or an equivalent model to account for environmental constraints on plant stomatal uptake during the whole growing season. Furthermore, the application of DRRs might introduce uncertainties due to the omission of complex interactions among biotic and abiotic factors at varied spatiotemporal scales.

Alternatively, more and more mechanistic schemes were developed and implemented in dynamic global vegetation models (DGVMs) to assess the joint effects of environmental factors and O₃ on plants. Felzer et al. (2004) considered both the damaging (through AOT40) and healing (through growth) processes related to O₃ effects within the framework of the Terrestrial Ecosystem Model. They further estimated a reduction of 2.6 %–6.8 % in the net primary productivity by O₃ pollution in the US during 1980–1990. Different from Felzer et al. (2004), Sitch et al. (2007) proposed a flux-based scheme linking the instantaneous POD_y with plant damage through the coupling between stomatal conductance and photosynthetic rate. Implementing this scheme into the Yale Interactive terrestrial Biosphere (YIBs) vegetation model, Yue and Unger (2015) predicted a range of 2 %–5 % reduction in global gross primary productivity (GPP) taking into account the low to high O₃ sensitivities for each vegetation type. Lombardozzi et al. (2015) collected hundreds of measurements and derived the decoupled responses on stomatal conductance and photosynthesis for the same O₃ uptake fluxes. They further implemented the separate response relationships into the Community Land Model and estimated a reduction of 8 %–12 % in GPP by O₃ in the present day. Coupling these schemes with Earth system models, studies have assessed interactive O₃ impacts on the carbon sink (Oliver et al., 2018; Yue and Unger, 2018), global warming (Sitch et al., 2007), and air pollution (Zhou et al., 2018; Gong et al., 2020, 2021; Zhu et al., 2022).

Although different schemes considered varied physical processes (Ollinger et al., 1997; Felzer et al., 2004, 2009; Sitch et al., 2007; Lombardozzi et al., 2015; Oliver et al., 2018), they followed the same principle that different O₃ sensitivities should be applied for varied plant functional types (PFTs), as revealed by many measurements in the past 4 decades (Buker et al., 2015; Mills et al., 2018b) (Table S1 in the Supplement). Generally, needleleaf trees, deciduous woody plants, and crop species show ascending sensitivities to O₃ (Reich and Amundson, 1985; Davison and Barnes, 1998; Buker et al., 2015). But the cause of such variation is not fully understood and thus has not been uniformly described in vegetation models (Massman et al., 2000; Tiwari et al., 2016). As a result, all large-scale assessments of O₃ vegetation damage had to rely on a PFT-based range of sensitivity parameters derived from a limited number of plant species and attempted to envelop the range of O₃ impacts by assuming all species within a PFT have either a “high” or “low” sensitivity to O₃. For example, Felzer et al. (2004) defined empirical sensitivity coefficients for three major plants including deciduous trees, coniferous trees, and crops. In Sitch et al. (2007), the sensitivity coefficients were defined separately for five PFTs with high/low ranges calibrated by DRRs of typical species. These synthesized assumptions cannot resolve the intra-PFT variations in the O₃ sensitivity and thus may cause large uncertainties in regional to global assessments.

Recent observations revealed a uniform plant sensitivity to O₃ if stomatal O₃ flux was expressed based on a leaf mass rather than leaf area basis (Li et al., 2016, 2022; Feng et al., 2018). The trait of leaf mass per area (LMA) is an important metric linking leaf area to mass. In a comparative study with 21 woody species (Li et al., 2016) and a meta-analysis of available experimental data (Feng et al., 2018), the DRR showed convergent O₃ sensitivities for conifer and broadleaf trees if the area-based stomatal uptake was converted to the mass-based flux with LMA. This is likely related to the diluting effect of thicker leaves, which normally have stronger defenses against O₃ in their cross-section. Nowadays, a large number of trait observations have been synthesized by global networks (Gallagher et al., 2020). The TRY initiative (Kattge et al., 2011) was one of the most influential datasets with 2.3 billion trait data by the year 2021. Based on the TRY dataset, global LMA was estimated with upscaling techniques such as Bayesian modeling (Butler et al., 2017) (hereafter B2017) or the random forest model (Moreno-Martinez et al., 2018) (hereafter M2018). These advances in the retrieval of LMA provide the possibility to depict more accurate O₃ vegetation damage at the global scale.

Here, we present a major advance in large-scale assessments of O₃ plant damage using a trait-based approach. We implement LMA into a stomatal flux-based O₃ damage framework aiming at a unified representation of plant O₃ sensitivities over the global grids. We couple this new approach to the Yale Interactive terrestrial Biosphere (YIBs) model (Yue and Unger, 2015) and evaluate the derived O₃ sensitivities against observations. We further assess contemporary O₃ impacts on global GPP in combination with the recently developed LMA datasets (Butler et al., 2017; Moreno-Martinez et al., 2018; Gallagher et al., 2020) (Fig. 1a) and the multi-model ensemble mean surface O₃ concentrations (Fig. 1b). The updated risk map for O₃ vegetation damage is used to identify the regions and vegetation types most at risk of O₃ damage.

Figure 1 Global leaf mass per area (LMA) and surface ozone (O₃) concentrations. The (a) LMA is adopted from Moreno-Martinez et al. (2018) (M2018) and (b) annual mean O₃ is derived from TF-HTAP (Turnock et al., 2018).

2 Scheme development and calibration

2.1 The trait-based O₃ vegetation damage scheme

We develop the new scheme based on the Sitch et al. (2007) (hereafter S2007) framework for transient O₃ damage calculation. In the original S2007 scheme, the undamaged fraction F for net photosynthetic rate is dependent on the excessive area-based stomatal O₃ flux, which is calculated as the difference between $f_{{\mathrm{O}}_{\mathrm{3}}}$ and PFT-specific area-based threshold y, and modulated by the sensitivity parameter a_PFT:

$$\begin{array}{}\text{(1)}& F=\mathrm{1}-a_{\mathrm{PFT}}\times max\left\{f_{{\mathrm{O}}_{\mathrm{3}}}-y,\mathrm{0}\right\},\end{array}$$

where a_PFT is calibrated and varies among PFTs with a typical range from “low” to “high” values indicating uncertainties of plant species within the same PFT. The stomatal O₃ flux $f_{{\mathrm{O}}_{\mathrm{3}}}$ (nmol m⁻² s⁻¹) is calculated as

$$\begin{array}{}\text{(2)}& f_{{\mathrm{O}}_{\mathrm{3}}}={\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{3}}\right]}{r+\left[\frac{k_{{\mathrm{O}}_{\mathrm{3}}}}{g_{\mathrm{p}}\times F}\right]}},\end{array}$$

where [O₃]  is the O₃ concentration at the reference level (nmol m⁻³), and r is the aerodynamic and boundary layer resistance between leaf surface and reference level (s m⁻¹). The setting of $k_{{\mathrm{O}}_{\mathrm{3}}}$ to 1.67 represents the ratio of leaf resistance for O₃ to that for water vapor. g_p represents potential stomata conductance for H₂O (m s⁻¹).

Studies suggested that LMA could be used to unify the area-based plant sensitivities to O₃ (Li et al., 2016; Feng et al., 2018), resulting in a constant mass-based parameter a independent of plant species and PFTs:

$$\begin{array}{}\text{(3)}& a=a_{\mathrm{PFT}}\times \mathrm{LMA}.\end{array}$$

Here, we convert the area-based O₃ stomatal flux expression in Eq. (1) to a mass-based flux as follows:

$$\begin{array}{}\text{(4)}& F=\mathrm{1}-a\times max\left\{{\displaystyle \frac{f_{{\mathrm{O}}_{\mathrm{3}}}}{\mathrm{LMA}}}-x,\mathrm{0}\right\},\end{array}$$

where the new sensitivity parameter a is a cross-species constant (nmol⁻¹ s g), LMA is leaf mass per area (g m⁻²), and the flux threshold is replaced by a mass-based value of x (nmol g⁻¹ s⁻¹) (Feng et al., 2018). Equations (2) and (4) can form a quadratic equation. F can be derived at each time step (i.e., hourly) and applied to net photosynthetic rate and stomatal conductance to calculate the O₃-induced damage. The updated LMA-based framework (YIBs-LMA) reduces the number of O₃ sensitivity parameters from three for each PFT (Sitch et al., 2007) in S2007 to a single parameter a for all PFTs. For the YIBs-LMA framework, the default value of the x threshold in Eq. (4) is set to 0.019 nmol g⁻¹ s⁻¹, as recommended by Feng et al. (2018).

2.2 Dose–response relationship (DRR)

We compare the simulated and observed sensitivities to O₃ so as to calibrate the LMA-based scheme. In field experiments, DRR is used to quantify species-specific damage by O₃ with a generic format as follows:

$$\begin{array}{}\text{(5)}& R=\mathrm{100}+S_{\mathrm{O}}\times {\mathit{\varphi }}_{{\mathrm{O}}_{\mathrm{3}}},\end{array}$$

where R (%) is the relative percentage of a biotic indicator (such as biomass or yield) after and before O₃ damage, ${\mathit{\varphi }}_{{\mathrm{O}}_{\mathrm{3}}}$ is an area-based O₃ metric (e.g., POD_y measured in sunlit leaves at the top of canopy), and S_O (usually negative) is the observed sensitivity derived as the slope of linear relationship between R and ${\mathit{\varphi }}_{{\mathrm{O}}_{\mathrm{3}}}$. We collected S_O from DRRs with conventional criteria (typically POD_y=1 for natural PFTs and POD_y=6 for crops as dose metrics (LRTAP Convention, 2017); the biotic indicators include the relative biomass for natural PFTs and relative yield for crops) among plant species from the International Cooperative Programme on Effects of Air Pollution on Natural Vegetation and Crops (CLRTAP) (LRTAP Convention, 2017) and multiple literature sources (Table S1). Such observations are used to calibrate the LMA-based scheme.

As a comparison with observations, we calculate annual relative GPP percentage (R_GPP, %) and POD_y of sunlit leaves in the first canopy layer (mmol m⁻² yr⁻¹, based on per leaf area) from the vegetation model to derive the slopes (S_S) of simulated DRRs. Here, POD_y is a diagnostic variable calculated as

$$\begin{array}{}\text{(6)}& {\mathrm{POD}}_y=\int max\left\{f_{{\mathrm{O}}_{\mathrm{3}}}-y,\mathrm{0}\right\},\end{array}$$

where $f_{{\mathrm{O}}_{\mathrm{3}}}$ represents the stomatal O₃ flux under instant O₃ stimulus at each time step, which can be calculated following Eq. (2) at the leaf level, and y is the prescribed critical level (1 nmol m⁻² s⁻¹ for natural or 6 nmol m⁻² s⁻¹ for crop species; LRTAP Convention, 2017). Excessive O₃ flux above y is accumulated for the sunlit leaves of the top canopy layer and over the growing season to derive the POD_y. The simulated S_S is calculated as the slope of the regression between simulated R_GPP (%) and POD_y at the PFT level. Only the dominant PFT in each grid is considered for the estimate of S_S at both PFT level and gridded analyses.

Similarly, the mass-based POD_x is derived from O₃-impacted $f_{{\mathrm{O}}_{\mathrm{3}}}$ (nmol m⁻² s⁻¹) in Eq. (2), together with gridded LMA (g m⁻²) and mass-based threshold x (nmol g⁻¹ s⁻¹), as

$$\begin{array}{}\text{(7)}& {\mathrm{POD}}_x=\int \left({\displaystyle \frac{f_{{\mathrm{O}}_{\mathrm{3}}}}{\mathrm{LMA}}}-x\right).\end{array}$$

Table 1 Summary of simulations.

^a Units of thresholds are nmol g⁻¹ s⁻¹ for x and nmol m⁻² s⁻¹ for y. ^b Units of key parameters are nmol⁻¹ s g for a and nmol⁻¹ m² s for a_PFT.

2.3 Simulations and calibrations

We perform two groups of supporting experiments (Table 1). The first group explores modeling uncertainties associated with the mass-based framework: (1) YIBs-LMA_B2017 replaces the default LMA map of M2018 (Moreno-Martinez et al., 2018) with B2017 (Butler et al., 2017). (2) YIBs-LMA_PFT applies PFT-specific LMA values (Table S2) for each PFT without considering global LMA geo-gradient. (3) YIBs-LMA_T replaces the default threshold of x=0.019 nmol g⁻¹ s⁻¹ with x=0.006 nmol g⁻¹ s⁻¹, which is an alternative parameter suggested by observations (Feng et al., 2018). The second group of supporting experiments explores the differences between mass-based and S2007 area-based frameworks. Typically, S2007 has a “low to high” a_PFT range for each PFT. Here, a mean sensitivity parameterization of S2007 (YIBs-S2007_adj) is re-calibrated according to S_O in Table S1.

Table 2 Calibrations of the YIBs-LMA^a experiment with varied a.

^a All runs from the YIBs-LMA experiment use x=0.019 nmol g⁻¹ s⁻¹ and the LMA map from M2018. ^b Slopes of simulated DRRs (S_S) are divided by observations (S_O, Table S1) to derive the model-to-observation ratios (“$S_{\mathrm{S}}/S_{\mathrm{O}}$”). The O₃ dose metric is POD_y=1 for natural PFTs and POD_y=6 for crops. The median and standard deviation (SD) of $S_{\mathrm{S}}/S_{\mathrm{O}}$ ratios of all PFTs are calculated for each set of specific parameter a. The values in parentheses exclude the effect of those numbers marked with ^c that are beyond 1 standard deviation. ^d The slopes (fitting) of linear regressions between S_O and S_S are listed for each a. The optimal a values with 1:1 fitting between S_S and S_O are shown in bold.

For all supporting experiments, the parameter a for YIBs-LMA or the eight mean a_PFT values for YIBs-S2007_adj are derived with the optimal 1:1 fitting between S_S and S_O to minimize the possible biases (Tables 2 and S3–S6). The basic method for calibration is feeding the model with series values of a or a_PFT until the predicted O₃ damage matches observations with the lowest normalized mean bias (NMB). For all LMA-based experiments, S_S values from varied PFTs grouped for the calibration of a, while for a_PFT in YIBs-S2007_adj, each a_PFT is determined individually by matching simulated S_S with S_O. Since S_O values are available only for six out of the eight YIBs PFTs, including evergreen broadleaf forest (EBF), needleleaf forest (NF), deciduous broadleaf forest (DBF), C₃ grass, C₄ grass, and crop (Table S1), S_O values of these PFTs are used for calibration. All runs are summarized in Table 1.

2.4 YIBs model and forcing data

In this study, all O₃ vegetation damage schemes are implemented in the YIBs model (Yue and Unger, 2015), which is a process-based dynamic global vegetation model incorporated with well-established carbon, energy, and water interactive schemes. The model applies the same PFT classifications as the Community Land Model (Bonan et al., 2003) (Fig. S1 in the Supplement). Eight PFTs are employed including evergreen broadleaf forest (EBF), needleleaf forest (NF), deciduous broadleaf forest (DBF), cold shrub (C_SHR), arid shrubland (A_SHR), C₃ grassland (C3_GRA), C₄ grassland (C4_GRA), and cropland (CRO) (Fig. S1). For each PFT, phenology is well evaluated (Yue and Unger, 2015) to generate a reliable growing season, which is crucial for the simulation of stomatal O₃ uptake (Anav et al., 2018). Photosynthesis and stomatal processes are calculated using Farquhar et al. (1980) and Ball et al. (1987) algorithms, respectively. Leaf area index (LAI) and tree height are predicted dynamically based on vegetation carbon allocation. The YIBs model has joined the multi-model ensemble project TRENDY and showed reasonable performance in the simulations of global biomass, GPP, LAI, net ecosystem exchange, and soil carbon relative to observations (Friedlingstein et al., 2022). Key plant biogeochemical parameters of the YIBs model are adjusted for this research (Table S7).

The hourly Modern-Era Retrospective Analysis for Research and Applications version 2 (MERRA2) climate reanalyses (Gelaro et al., 2017) are used to drive the YIBs model. The gridded LMA required for the main mass-based simulation is derived from Moreno-Martinez et al. (2018) (M2018), which shows the highest value of >150 g m⁻² for needleleaf forest at high latitudes while low values of ∼40 g m⁻² for grassland and cropland (Figs. 1a and S1). Grids with missing LMA data are filled with the mean of the corresponding PFT. Contemporary O₃ concentration fields in the year 2010 from the multi-model mean in Task Force on Hemispheric Transport of Air Pollutants (TF-HTAP) experiments (Turnock et al., 2018) (Fig. 1b) are used as forcing data. The original monthly O₃ data are downscaled to an hourly scale using the diurnal cycle predicted by the chemistry–climate–carbon fully coupled model ModelE2-YIBs (Yue and Unger, 2015). Generally, areas of severe O₃ pollution are found in the midlatitudes of the Northern Hemisphere, with the highest annual average O₃ concentration of over 40 ppbv in East Asia. All data are interpolated to a spatial resolution of $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$.

Figure 2 Simulated area-based (a, b, c) or mass-based (d, e, f) DRRs for the YIBs-LMA experiment. Three tests from the YIBs-LMA experiment all adopt x=0.019 nmol g⁻¹ s⁻¹ and gridded LMA from M2018 with parameter a=2.5, 3.5, and 4.5 nmol⁻¹ s g, respectively. Each dot represents estimated POD-R_GPP (POD_y=1 for a–c, POD_x=0.019 for d–e) at a grid with corresponding PFT. The PFT-specific regressions between area- or mass-based POD and R_GPP are displayed with solid lines shown in legend. Regression relationships of all PFTs are displayed by the dashed black line, with coefficients of determination (R²) denoted in each panel. Note the differences of ranges in x axis among PFTs. The YIBs-LMA experiment is summarized in Table 1.

3 Results

3.1 Comparison of simulated sensitivities with observations

Simulated relative GPP percentage (R_GPP) values at global grids were sorted by dominant PFTs (Fig. S1) and plotted against area-based accumulated phytotoxic O₃ dose above a threshold y nmol m⁻² s⁻¹ (POD_y=1) at the corresponding grids (Fig. 2). The DRR shows varied slopes among different PFTs, resulting in a coefficient of determination (R²) around 0.54 for all PFTs (Fig. 2a–c). We further calculated the mass-based accumulated phytotoxic O₃ dose above a threshold of 0.019 nmol g s⁻¹ (POD_x=0.019) and compared it with R_GPP. The updated DRR showed convergent slopes and reached a high R² of 0.77 across all PFTs (Fig. 2d–f), suggesting that the mass-based scheme could better unify O₃ sensitivities among different PFTs.

We then calibrated the single, best-fit a value for the YIBs-LMA framework by minimizing the absolute difference between simulated (S_S) and observed (S_O) slopes of O₃ DRR for all PFTs. With different a parameters, the YIBs-LMA framework yielded considerably high R² values of ∼1.0 but varied biases between simulated and observed O₃ impacts across PFTs (Fig. 3). Both the 1:1 fitting and the lowest bias between S_S and S_O were achieved with an optimal a=3.5 nmol⁻¹ s g (Fig. 3c). Notably, such calibration of a is robust under different O₃ fields (see Fig. S2). Consistent with observations, YIBs-LMA with this optimal a parameter simulated low S_S values of −0.18 % and −0.36 % per mmol m⁻² yr⁻¹ of POD_y=1 for evergreen broadleaf forest and needleleaf forest, respectively (Fig. 4a, b); median S_S values from −0.53 % per mmol m⁻² yr⁻¹ for arid shrubland (Fig. 4e); and high S_S values from −0.64 % to −1.04 % per mmol m⁻² yr⁻¹ for deciduous broadleaf forest, C₃ and C₄ grassland, cropland, and cold shrubland (−3.28 % for crops with POD_y=6; Fig. 4c–d, f–h).

Figure 3 Comparison between S_O (% per mmol m⁻²) and S_S (% per mmol m⁻²) for the YIBs-LMA experiment. Five tests from the YIBs-LMA experiment all adopt x=0.019 nmol g⁻¹ s⁻¹ and gridded LMA from M2018 but with varied parameter a from (a) 2.5 to (e) 4.5 nmol⁻¹ s g. S_O values are from Table S1. S_S values are derived as the slope between R_GPP and POD_y. The linear regression (dashed lines), 1:1 fitting (light-grey lines), normalized mean bias (NMB), and correlation coefficient (r) between S_S and S_O for varied PFTs are shown on each panel. The S_S and S_O of CRO are derived using POD_y=6, while other PFTs use POD_y=1. The YIBs-LMA experiment is described in Table 1.

Figure 4 Comparisons of observed and simulated dose–response relationships. Simulated PFT-specific DRRs are derived from YIBs-LMA with gridded LMA from M2018, x=0.019 nmol g⁻¹ s⁻¹, and a=3.5 nmol⁻¹ s g. Each dot represents results from a grid cell with corresponding PFT. The regressions between relative GPP percentage (R_GPP) and leaf area-based stomatal O₃ uptake fluxes (POD_y=1 for natural PFTs and POD_y=6 for crops) are shown for observations (red; see Table S1) and simulations (blue) with slopes of DRRs denoted as S_O and S_S, respectively. S_O values are missing for (d) cold and (e) arid shrubs. Coefficients of determination (R²) of simulations are displayed in each panel. Note the differences of ranges in x axis among PFTs (PFTs are shown in Fig. S1).

3.2 Global map of O₃ vegetation damage

We estimated contemporary GPP reductions induced by O₃ with the global concentrations of surface O₃ (Fig. 1b) in the year 2010. The YIBs-LMA framework using an increase of a parameter yielded an almost linear enhancement of global GPP reduction (Fig. S3) with consistent spatial distributions (Fig. S4). The simulation with the optimal a=3.5 nmol⁻¹ s g predicted a global GPP reduction of 4.8 % (Fig. 5a), which was similar to the value estimated with the area-based S2007 scheme (YIBs-S2007_adj, Table 1). Large reductions of >10 % were predicted over the eastern US, western Europe, eastern China, and India (Fig. 5a). Hotspots were mainly located in cropland and agricultural areas mixed with deciduous broadleaf forest or grassland, accompanied by moderate to high levels of surface O₃. Few discrepancies between the damage maps of YIBs-LMA and YIBs-S007_adj were found (Figs. 5b and S5), even though the number of parameters was greatly reduced in the YIBs-LMA scheme.

Figure 5 Global O₃ vegetation damage simulated with the LMA-based scheme. Results shown are the (a) GPP reduction percentages by O₃ simulated with the YIBs-LMA framework (gridded LMA from M2018, x=0.019 nmol g⁻¹ s⁻¹, and a=3.5 nmol⁻¹ s g) and (b) their differences compared to the predictions from YIBs-S2007_adj simulation. Blue (red) patches indicate the regions where damage predicted in YIBs-LMA is lower (higher) than in YIBs-S2007_adj.

For YIBs-LMA, PFTs with low LMA such as cropland, grassland, and deciduous broadleaf forest account for 73.3 Pg C yr⁻¹ (50.0 %) of the global GPP (Table S8). However, these PFTs contributed to a total GPP reduction of 5.4 Pg C yr⁻¹ (75.5 % of total GPP loss) by O₃ damage. In contrast, evergreen broadleaf and needleleaf forests with high LMA accounted for 48.8 Pg C yr⁻¹ (33.0 %) of total GPP but yielded only a reduction of 0.75 Pg C yr⁻¹ (10.5 % of total GPP loss). Differences in GPP percentage losses were in part associated with the global pattern of O₃ concentrations, which were usually higher over midlatitudes with populated cities and dense crop plantations (Fig. 1b). However, the differences in LMA and simulated O₃ sensitivities of these PFTs also made important contributions to such discrepancies in GPP damage.

Figure 6 Global O₃-induced GPP reductions simulated in four supporting experiments. All damage maps are based on optimal parameters shown in Table 1. The spatial correlation coefficients between YIBs-LMA and these supporting simulations are shown in each panel as well as the global average damage percentage of each map. Simulations are described in Table 1.

Figure 7 Comparison of $S_{\mathrm{S}}/S_{\mathrm{O}}$ among supporting experiments. Individual ratios for (b) different PFTs are grouped to the box plot in (a). All experiments adopt optimal parameters shown in Table 1.

3.3 Uncertainties of the LMA-based scheme

We quantified the uncertainties in the LMA-based scheme by comparing simulated GPP damage among different experiments (Table 1). The experiment with the alternative LMA map of B2017 (Fig. S6) showed similar spatial patterns but a slightly enhanced GPP reduction of 5.3 % (Fig. 6a) compared to the simulations using LMA map of M2018 (Fig. 5a). The B2017 map contains more LMA data than M2018 (∼40 %), leading to unexpected high O₃ threats over the tundra in the Arctic (Fig. S7). Another experiment using PFT-specific LMA estimated a global GPP reduction of 4.6 % (Fig. 6b) with a consistent spatial pattern similar to the prediction in YIBs-LMA, suggesting that the PFT-level LMA can be used in the case of a lack of regional LMA data. The third experiment with an alternative threshold flux (Feng et al., 2018) of 0.006 nmol g⁻¹ s⁻¹ estimated a high GPP reduction of 6.5 % (Fig. 6c) due to the overestimations of O₃ sensitivities for some tree PFTs (Fig. 7). The fourth run, YIBs-S2007_adj, predicted a similar global GPP damage of 4.79 %, similar to the YIBs-LMA run, with a high spatial correlation coefficient of 0.98 (Fig. 6d). Such good consistency is mainly due to the application of recalibrated PFT-level sensitivities in YIBs-S2007_adj. Finally, we tested a new calibration excluding CRO, the PFT that contributed the most to the calibration biases (shown by dashed orange lines in Fig. S8). The results gave an optimal a of 3.2, with global damage of 4.5 %. All sensitivity experiments achieved results consistent with the YIBs-LMA simulation, with damage ranging from 4.5 % to 6.5 % and spatial correlation coefficients larger than 0.94.

4 Discussion

4.1 Mechanisms behind the LMA-based approach

In recent decades, the plant science community has examined how traits could be used to differentiate and predict the function of plant species (Reich et al., 1997, 1999). LMA, related to leaf density and thickness, is a key trait reflecting many aspects of leaf function (Reich et al., 1998). In the field of O₃ phytotoxicology, experiments have revealed plants with high LMA usually have thick leaves with physical and chemical defenses (Poorter et al., 2009), which can strengthen their resistance to O₃ (Li et al., 2016; Feng et al., 2018). Conversely, plants with low LMA normally have thin leaves which are likely to be less O₃-tolerant (Li et al., 2016; Feng et al., 2018). Moreover, it seems plausible that the oxidative stress caused by a given amount of stomatal O₃ flux per unit leaf area would be distributed over a larger leaf mass, and hence diluted, in a leaf with high LMA. Such an LMA–O₃ sensitivity relationship can be well reproduced by our LMA-based model (Fig. 8a and b). Below we explore the linkage between O₃ plant sensitivities and the mutual adaptation of growth strategies and leaf morphology with plant leaf trade-off theory (Reich et al., 1999; Shipley et al., 2006).

In the natural world, plants often adapt to maximize carbon uptake under prevailing conditions (Reich et al., 1998; Shipley et al., 2006). To make full use of resources in the growing season, leaves under varied living conditions choose either fast photosynthetic rates (fast-growing deciduous types) or long photosynthesis duration (slow-growing evergreen types) with compatible leaf structures (Reich, 2014; Diaz et al., 2016). The former species expand leaf area (low LMA) to maximize light interception, while the latter species produce thick and mechanically strong leaves (high LMA) with ample resistant substances for durable utilization (Poorter et al., 2009) in resource-limited and/or environment-stressed habitats (Wright et al., 2002). As a side effect of such leaf trade-offs, deciduous plants with their high rates of photosynthesis, associated large stomatal conductance (Davison and Barnes, 1998; Henry et al., 2019), and less total defense capacity through the leaf profile (Poorter et al., 2009), are highly O₃-sensitive (Mode1 in Fig. 9). In contrast, moderate photosynthesis, relatively low maximum stomatal conductance (Davison and Barnes, 1998; Henry et al., 2019), and reinforced dense leaves (Poorter et al., 2009) lead to low sensitivity for evergreen plants (Mode2 in Fig. 9). Therefore, in our modeling practice, the mass-based O₃ gas exchange algorithm can be regarded as taking into account several interrelated factors such as growth-driven gas exchange requirements, gas path length, and biochemical reserves, in a unified, simplified, and effective manner via LMA.

Figure 8 Relationships between O₃ sensitivity and LMA. (a) Simulated O₃ sensitivity (S_S) at each grid is compared with LMA for different PFTs. Gridded S_S is derived as GPP change per unit POD_y=1 from the YIBs-LMA simulation. Each point represents the value in a grid cell with a dominant PFT. (b) The PFT-level relationships between LMA and O₃ sensitivity are grouped as box plots, which indicate the median, 25th percentile, and 75th percentile of y-axis variables within the same PFT. The width of box plots represents 1 standard deviation of LMA for a specific PFT.

Figure 9 Illustration of the relationships between leaf trade-off strategy and its sensitivity to O₃.

4.2 Implication of potential risks for fast-growing plants

Our new approach reflected the general experimental findings that deciduous plants are much more vulnerable to O₃ than evergreen species (Li et al., 2017; Feng et al., 2018), and in turn within a PFT, early-successional/pioneers with low LMA are likely more vulnerable than late-successional/canopy trees with high LMA (Fyllas et al., 2012). This law has been neglected in previous modeling studies due to the dependence on the limited observed data used for PFT-specific tuning. Our LMA-based approach bridges this gap through grid-based parameterization, and in addition, our data–model integration specifically emphasizes the broad high risks for fast-growing plants, especially for crops. Among PFTs, crops may endure the largest O₃ threats (Davison and Barnes, 1998; Feng et al., 2021; Mukherjee et al., 2021) because they are artificially bred with high photosynthetic capacities (Richards, 2000), stomatal conductance, generally low LMA (Bertin and Gary, 1998; Wang and Shangguan, 2010; Wu et al., 2018; Li et al., 2018) (roughly 30–60 g m⁻²), and cultivated in populated regions with high ambient O₃ concentrations. Modern technology aims to promote crop yield (Herdt, 2005), but this can potentially elevate crop sensitivities to O₃ (Biswas et al., 2008, 2013). This study estimated the highest annual mean GPP damage for crop, 12.6 %, which is at the high end of the 4.4 %–12.4 % of the O₃-induced yield loss estimated for global modeling of soybean, wheat, rice, and maize (Mills et al., 2018a). Furthermore, human-induced land use activities may also increase O₃ damage risks. The global demand for food and commodities leads to the conversion of natural forests to irrigated croplands, grazing pastures, and economical-tree plantations (Curtis et al., 2018; Zalles et al., 2021). Meanwhile, the urgent actions to combat climate change promote large-scale afforestation and reforestation (Cook-Patton et al., 2020). These land use changes with fast-growing plant species may increase the risks of terrestrial ecosystems to surface O₃.

4.3 Advances in the global O₃ damage assessment

For the first time, we implemented plant trait LMA into a process-based O₃ impact modeling scheme and obtained reasonable interspecific and inter-PFT O₃ responses supported by observations. The similarity between YIBs-S2007 and YIBs-LMA shown in Fig. 5 revealed an advance in the modeling strategy. Simulated O₃ damage in YIBs-S2007 is based on the PFT-level calibrations that tuned sensitivity parameters of each PFT with observed DRRs. Such refinement is a data-driven approach without clear physical reasons. Instead, the YIBs-LMA framework converts the area-based responses to mass-based ones and achieves better unification in O₃ sensitivities among different PFTs. In this algorithm, the O₃ damage efficiency is inversely related to plant LMA, which influences both the O₃ uptake potential and the detoxification capability of the vegetation. The similarity in the global assessment of O₃ vegetation damage between YIBs-S2007 and YIBs-LMA further demonstrated the physical validity of the LMA-based scheme in Earth system modeling because the independent LMA map was applied in the latter approach.

In addition to the advance in physical mechanisms, the LMA-based approach improves global O₃ damage assessments in the following aspects. First, it significantly reduces the number of required key parameters. To account for interspecific sensitivities, many schemes have to define PFT-level parameters to capture the ranges of plant responses (Sitch et al., 2007; Felzer et al., 2009; Lombardozzi et al., 2015). As a result, those schemes rely on dozens of parameters, which increases the uncertainties of modeling and the difficulties for model calibration. The LMA-based approach requires the calibration of one single parameter a, largely facilitating its application across different vegetation models. Second, the new approach accounts for the continuous spectrum of O₃ sensitivities. Previous studies usually categorized species into groups of low or high O₃ sensitivity, depending on very limited data from O₃ exposure experiments. As a result, grid cells for a specific PFT share the same sensitivities regardless of their geographic locations and ecosystem characteristics. In reality, there are hundreds and thousands of plant species in each PFT, and they usually have large variations in biophysical parameters including LMA and O₃ sensitivities. The LMA-based approach takes advantage of the newly revealed unifying concept in O₃ sensitivity (Li et al., 2016, 2022; Feng et al., 2018) and the recent development of a trait-based LMA global map (Fig. 1a). Such configurations present a spectrum of gridded O₃ sensitivities (Fig. 8a) following the variations of LMA distribution.

5 Outlook for future modeling

In nature, all aspects of plant physiochemical processes, such as growth, development, reproduction, and defense, are influenced by abiotic factors like water availability, temperature, CO₂ concentration, and light resources (Kochhar and Gujral, 2020). In our modeling, the cumulative O₃ fluxes are based on dynamic plant simulations with well-established DGVM to calculate the effects of these abiotic factors. LMA is considered a factor representing the vulnerability of each species, by which divergent responses to the same O₃ stomatal dose can be further differentiated. In fact, many other key variables in DGVMs, for example, leaf photosynthetic traits (V_cmax and J_max), nutrient traits (leaf nitrogen and phosphorus), morphological traits (leaf thickness and size), and phenology-related traits (leaf life span) are all more or less interlinked with LMA (Walker et al., 2014). There are some generic regression relationships between them. In addition, efforts are being made to directly predict key traits, including LMA, through environmental factors. As a result, considerable improvements can be made in the direction of trait-flexible modeling within the existing DGVM frameworks. Our study demonstrates the validity of LMA-based approach for the O₃ plant damage modeling.

Although we used the most advanced LMA integrated from available observations, this dataset was developed based on static global grids and revealed the mean state for each pixel. In reality, LMA can vary with biotic/abiotic factors like leaf position in the canopy (Keenan and Niinemets, 2017), phenology, plant health, living environment (Fritz et al., 2018), and climate (Wright et al., 2005; Cui et al., 2020). Even long-term exposure to O₃ can alter leaf morphological characteristics and LMA (Li et al., 2017). In future studies, simulations from local to global scales could implement the spatiotemporal variations in LMA taking into account the demographic information and environmental forcings. We expect a breakthrough in the calculation of reliable LMA to achieve fully dynamic predictions of O₃ plant damage in Earth system modeling, thus facilitating the research of plant response and adaption in changing environments.