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

Abstract

A major limitation in modeling global ozone (O${}_{\mathrm{3}}$) vegetation damage has long been the reliance on empirical O${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ plant injury by linking the trait leaf mass per area (LMA) and plant O${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$, with a range of 1.1 %–12.6 % for varied PFTs. Hotspots with damage $>\mathrm{10}\mathit{\%}$ are found in agricultural areas in the eastern US, western Europe, eastern China, and India, accompanied by moderate to high levels of surface O${}_{\mathrm{3}}$. Furthermore, we simulate the distribution of plant sensitivity to O${}_{\mathrm{3}}$, 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${}_{\mathrm{3}}$) has long been recognized as a hazardous pollutant for plants (Richards et al., 1958; Reich and Amundson, 1985). As a strong oxidant, O${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ plant damage in global modeling and assess its coupling effects in biosphere–atmosphere systems (Zhou et al., 2018).

To date, O${}_{\mathrm{3}}$ fumigation experiments have been conducted for various plant species. Accordingly, O${}_{\mathrm{3}}$ damage sensitivities, denoted as the dose–response relationships (DRRs), were derived as regressions between O${}_{\mathrm{3}}$ exposure metrics and changes in biotic indicators (Mills et al., 2011). The widely used O${}_{\mathrm{3}}$ metrics include ambient O${}_{\mathrm{3}}$ concentrations for AOT40 (accumulated O${}_{\mathrm{3}}$ concentration above the threshold of 40 ppbv; Fuhrer et al., 1997) or the stomatal O${}_{\mathrm{3}}$ flux for POD${}_y$ (phytotoxic O${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$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${}_{\mathrm{3}}$ on plants. Felzer et al. (2004) considered both the damaging (through AOT40) and healing (through growth) processes related to O${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ in the present day. Coupling these schemes with Earth system models, studies have assessed interactive O${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ (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${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ impacts by assuming all species within a PFT have either a “high” or “low” sensitivity to O${}_{\mathrm{3}}$. 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${}_{\mathrm{3}}$ sensitivity and thus may cause large uncertainties in regional to global assessments.

Recent observations revealed a uniform plant sensitivity to O${}_{\mathrm{3}}$ if stomatal O${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ vegetation damage at the global scale.

Here, we present a major advance in large-scale assessments of O${}_{\mathrm{3}}$ plant damage using a trait-based approach. We implement LMA into a stomatal flux-based O${}_{\mathrm{3}}$ damage framework aiming at a unified representation of plant O${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ sensitivities against observations. We further assess contemporary O${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ concentrations (Fig. 1b). The updated risk map for O${}_{\mathrm{3}}$ vegetation damage is used to identify the regions and vegetation types most at risk of O${}_{\mathrm{3}}$ damage.

Figure 1

Global leaf mass per area (LMA) and surface ozone (O${}_{\mathrm{3}}$) concentrations. The (a) LMA is adopted from Moreno-Martinez et al. (2018) (M2018) and (b) annual mean O${}_{\mathrm{3}}$ is derived from TF-HTAP (Turnock et al., 2018).

[Figure omitted. See PDF]

2 Scheme development and calibration

2.1

The trait-based O${}_{\mathrm{3}}$ vegetation damage scheme

We develop the new scheme based on the Sitch et al. (2007) (hereafter S2007) framework for transient O${}_{\mathrm{3}}$ damage calculation. In the original S2007 scheme, the undamaged fraction $F$ for net photosynthetic rate is dependent on the excessive area-based stomatal O${}_{\mathrm{3}}$ 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_{\mathrm{PFT}}$: 1$$F=\mathrm{1}-a_{\mathrm{PFT}}\times max\left\{f_{{\mathrm{O}}_{\mathrm{3}}}-y,\mathrm{0}\right\},$$ where $a_{\mathrm{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${}_{\mathrm{3}}$ flux $f_{{\mathrm{O}}_{\mathrm{3}}}$ (nmol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$) is calculated as 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]}},$$ where [$O_{\mathrm{3}}$]  is the O${}_{\mathrm{3}}$ concentration at the reference level (nmol m${}^{-\mathrm{3}}$), and $r$ is the aerodynamic and boundary layer resistance between leaf surface and reference level (s m${}^{-\mathrm{1}}$). The setting of $k_{{\mathrm{O}}_{\mathrm{3}}}$ to 1.67 represents the ratio of leaf resistance for O${}_{\mathrm{3}}$ to that for water vapor. $g_{\mathrm{p}}$ represents potential stomata conductance for H${}_{\mathrm{2}}$O (m s${}^{-\mathrm{1}}$).

Studies suggested that LMA could be used to unify the area-based plant sensitivities to O${}_{\mathrm{3}}$ (Li et al., 2016; Feng et al., 2018), resulting in a constant mass-based parameter $a$ independent of plant species and PFTs: 3$$a=a_{\mathrm{PFT}}\times \mathrm{LMA}.$$ Here, we convert the area-based O${}_{\mathrm{3}}$ stomatal flux expression in Eq. (1) to a mass-based flux as follows: 4$$F=\mathrm{1}-a\times max\left\{{\displaystyle \frac{f_{{\mathrm{O}}_{\mathrm{3}}}}{\mathrm{LMA}}}-x,\mathrm{0}\right\},$$ where the new sensitivity parameter $a$ is a cross-species constant (nmol${}^{-\mathrm{1}}$ s g), LMA is leaf mass per area (g m${}^{-\mathrm{2}}$), and the flux threshold is replaced by a mass-based value of $x$ (nmol g${}^{-\mathrm{1}}$ s${}^{-\mathrm{1}}$) (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${}_{\mathrm{3}}$-induced damage. The updated LMA-based framework (YIBs-LMA) reduces the number of O${}_{\mathrm{3}}$ 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${}^{-\mathrm{1}}$ s${}^{-\mathrm{1}}$, as recommended by Feng et al. (2018).

2.2 Dose–response relationship (DRR)

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

5$$R=\mathrm{100}+S_{\mathrm{O}}\times {\mathit{\varphi }}_{{\mathrm{O}}_{\mathrm{3}}},$$ where $R$ (%) is the relative percentage of a biotic indicator (such as biomass or yield) after and before O${}_{\mathrm{3}}$ damage, ${\mathit{\varphi }}_{{\mathrm{O}}_{\mathrm{3}}}$ is an area-based O${}_{\mathrm{3}}$ metric (e.g., POD${}_y$ measured in sunlit leaves at the top of canopy), and $S_{\mathrm{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_{\mathrm{O}}$ from DRRs with conventional criteria (typically POD${}_{y=\mathrm{1}}$ for natural PFTs and POD${}_{y=\mathrm{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_{\mathrm{GPP}}$, %) and POD${}_y$ of sunlit leaves in the first canopy layer (mmol m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$, based on per leaf area) from the vegetation model to derive the slopes ($S_{\mathrm{S}}$) of simulated DRRs. Here, POD${}_y$ is a diagnostic variable calculated as 6$${\mathrm{POD}}_y=\int max\left\{f_{{\mathrm{O}}_{\mathrm{3}}}-y,\mathrm{0}\right\},$$ where $f_{{\mathrm{O}}_{\mathrm{3}}}$ represents the stomatal O${}_{\mathrm{3}}$ flux under instant O${}_{\mathrm{3}}$ 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${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ for natural or 6 nmol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ for crop species; LRTAP Convention, 2017). Excessive O${}_{\mathrm{3}}$ 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_{\mathrm{S}}$ is calculated as the slope of the regression between simulated $R_{\mathrm{GPP}}$ (%) and POD${}_y$ at the PFT level. Only the dominant PFT in each grid is considered for the estimate of $S_{\mathrm{S}}$ at both PFT level and gridded analyses.

Similarly, the mass-based POD${}_x$ is derived from O${}_{\mathrm{3}}$-impacted $f_{{\mathrm{O}}_{\mathrm{3}}}$ (nmol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$) in Eq. (2), together with gridded LMA (g m${}^{-\mathrm{2}}$) and mass-based threshold $x$ (nmol g${}^{-\mathrm{1}}$ s${}^{-\mathrm{1}}$), as 7$${\mathrm{POD}}_x=\int \left({\displaystyle \frac{f_{{\mathrm{O}}_{\mathrm{3}}}}{\mathrm{LMA}}}-x\right).$$

Table 1

Summary of simulations.

Experiment${}^{\mathrm{a}}$	Method	Thresholds${}^{\mathrm{a}}$	LMA	LMA	Optimal	Tests
		($x$ or $y$)	format	map	($a$ or $a_{\mathrm{PFT}}$)	($a$ or $a_{\mathrm{PFT}}$)
YIBs-LMA	Mass-based	$x=\mathrm{0.019}$	gridded	M2018	$a=\mathrm{3.5}$	five tests
					(Table 2)	($a=\mathrm{2.5}$, 3, 3.5, 4, 4.5)
YIBs-LMA_PFT		$x=\mathrm{0.019}$	PFT-specific	M2018	$a=\mathrm{2.0}$	five tests
					(Table S3)	($a=\mathrm{2}$, 2.5, 3, 3.5, 4)
YIBs-LMA_T		$x=\mathrm{0.006}$	gridded	M2018	$a=\mathrm{3.0}$	five tests
					(Table S4)	($a=\mathrm{2}$, 2.5, 3, 3.5, 4)
YIBs-LMA_B2017		$x=\mathrm{0.019}$	gridded	B2017	$a=\mathrm{2.8}$	five tests
					(Table S5)	($a=\mathrm{2}$, 2.5, 2.8, 3, 3.5)
YIBs-S2007_adj	Area-based	eight values for $y$	–	–	eight values for $a_{\mathrm{PFT}}$	40 tests
		(Table S6)			(Table S6)	(five each for eight PFTs)

${}^{\mathrm{a}}$ Units of thresholds are nmol g${}^{-\mathrm{1}}$ s${}^{-\mathrm{1}}$ for $x$ and nmol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ for $y$. ${}^{\mathrm{b}}$ Units of key parameters are nmol${}^{-\mathrm{1}}$ s g for $a$ and nmol${}^{-\mathrm{1}}$ m${}^{\mathrm{2}}$ s for $a_{\mathrm{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=\mathrm{0.019}$ nmol g${}^{-\mathrm{1}}$ s${}^{-\mathrm{1}}$ with $x=\mathrm{0.006}$ nmol g${}^{-\mathrm{1}}$ s${}^{-\mathrm{1}}$, 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_{\mathrm{PFT}}$ range for each PFT. Here, a mean sensitivity parameterization of S2007 (YIBs-S2007_adj) is re-calibrated according to $S_{\mathrm{O}}$ in Table S1.

Table 2

Calibrations of the YIBs-LMA${}^{\mathrm{a}}$ experiment with varied $a$.

PFT	$S_{\mathrm{O}}$	$S_{\mathrm{S}}$					$S_{\mathrm{S}}/{S_{\mathrm{O}}}^{\mathrm{b}}$
		$a=\mathrm{2.5}$	$a=\mathrm{3.0}$	$\mathbf{a}=\mathbf{3.5}$	$a=\mathrm{4.0}$	$a=\mathrm{4.5}$	$a=\mathrm{2.5}$	$a=\mathrm{3.0}$	$\mathbf{a}=\mathbf{3.5}$	$a=\mathrm{4.0}$	$a=\mathrm{4.5}$
EBF	$-\mathrm{0.19}$	$-\mathrm{0.13}$	$-\mathrm{0.16}$	$\mathbf{-}$0.18	$-\mathrm{0.21}$	$-\mathrm{0.23}$	0.70	0.83	0.96	1.08	1.20
NF	$-\mathrm{0.23}$	$-\mathrm{0.26}$	$-\mathrm{0.31}$	$\mathbf{-}$0.36	$-\mathrm{0.40}$	$-\mathrm{0.45}$	1.14${}^{\mathrm{c}}$	1.35${}^{\mathrm{c}}$	1.56${}^{\mathrm{c}}$	1.75${}^{\mathrm{c}}$	1.95${}^{\mathrm{c}}$
DBF	$-\mathrm{0.70}$	$-\mathrm{0.51}$	$-\mathrm{0.60}$	$\mathbf{-}$0.69	$-\mathrm{0.78}$	$-\mathrm{0.87}$	0.72	0.86	0.99	1.12	1.24
C_SHR	–	$-\mathrm{0.75}$	$-\mathrm{0.90}$	$\mathbf{-}$1.04	$-\mathrm{1.18}$	$-\mathrm{1.31}$	–	–	–	–	–
A_SHR	–	$-\mathrm{0.38}$	$-\mathrm{0.45}$	$\mathbf{-}$0.53	$-\mathrm{0.60}$	$-\mathrm{0.66}$	–	–	–	–	–
C4_GRA	$-\mathrm{0.85}$	$-\mathrm{0.71}$	$-\mathrm{0.84}$	$\mathbf{-}$0.97	$-\mathrm{1.10}$	$-\mathrm{1.22}$	0.83	0.99	1.14	1.29	1.44
C3_GRA	$-\mathrm{0.62}$	$-\mathrm{0.47}$	$-\mathrm{0.55}$	$\mathbf{-}$0.64	$-\mathrm{0.73}$	$-\mathrm{0.81}$	0.75	0.89	1.03	1.17	1.30
CRO	$-\mathrm{3.35}$	$-\mathrm{1.97}$	$-\mathrm{2.57}$	$\mathbf{-}$3.28	$-\mathrm{4.11}$	$-\mathrm{5.10}$	0.59	0.77	0.98	1.23	1.52
Fitting${}^{\mathrm{d}}$	–	0.61	0.79	0.99	1.23	1.50	–	–	–	–	–
Median	–	–	–	–	–	–	0.74 (0.72)	0.88 (0.86)	1.01 (0.99)	1.20 (1.17)	1.37 (1.30)
SD	–	–	–	–	–	–	0.19 (0.09)	0.21 (0.08)	0.23 (0.07)	0.25 (0.08)	0.28 (0.13)

${}^{\mathrm{a}}$ All runs from the YIBs-LMA experiment use $x=$0.019 nmol g${}^{-\mathrm{1}}$ s${}^{-\mathrm{1}}$ and the LMA map from M2018. ${}^{\mathrm{b}}$ Slopes of simulated DRRs ($S_{\mathrm{S}}$) are divided by observations ($S_{\mathrm{O}}$, Table S1) to derive the model-to-observation ratios (“$S_{\mathrm{S}}/S_{\mathrm{O}}$”). The O${}_{\mathrm{3}}$ dose metric is POD${}_{y=\mathrm{1}}$ for natural PFTs and POD${}_{y=\mathrm{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 ${}^{\mathrm{c}}$ that are beyond 1 standard deviation. ${}^{\mathrm{d}}$ The slopes (fitting) of linear regressions between $S_{\mathrm{O}}$ and $S_{\mathrm{S}}$ are listed for each $a$. The optimal $a$ values with $\mathrm{1}:\mathrm{1}$ fitting between $S_{\mathrm{S}}$ and $S_{\mathrm{O}}$ are shown in bold.

For all supporting experiments, the parameter $a$ for YIBs-LMA or the eight mean $a_{\mathrm{PFT}}$ values for YIBs-S2007_adj are derived with the optimal $\mathrm{1}:\mathrm{1}$ fitting between $S_{\mathrm{S}}$ and $S_{\mathrm{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_{\mathrm{PFT}}$ until the predicted O${}_{\mathrm{3}}$ damage matches observations with the lowest normalized mean bias (NMB). For all LMA-based experiments, $S_{\mathrm{S}}$ values from varied PFTs grouped for the calibration of $a$, while for $a_{\mathrm{PFT}}$ in YIBs-S2007_adj, each $a_{\mathrm{PFT}}$ is determined individually by matching simulated $S_{\mathrm{S}}$ with $S_{\mathrm{O}}$. Since $S_{\mathrm{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${}_{\mathrm{3}}$ grass, C${}_{\mathrm{4}}$ grass, and crop (Table S1), $S_{\mathrm{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${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ grassland (C3_GRA), C${}_{\mathrm{4}}$ 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${}_{\mathrm{3}}$ 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 $>\mathrm{150}$ g m${}^{-\mathrm{2}}$ for needleleaf forest at high latitudes while low values of $\sim \mathrm{40}$ g m${}^{-\mathrm{2}}$ for grassland and cropland (Figs. 1a and S1). Grids with missing LMA data are filled with the mean of the corresponding PFT. Contemporary O${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ pollution are found in the midlatitudes of the Northern Hemisphere, with the highest annual average O${}_{\mathrm{3}}$ 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=\mathrm{0.019}$ nmol g${}^{-\mathrm{1}}$ s${}^{-\mathrm{1}}$ and gridded LMA from M2018 with parameter $a=\mathrm{2.5}$, 3.5, and 4.5 nmol${}^{-\mathrm{1}}$ s g, respectively. Each dot represents estimated POD-$R_{\mathrm{GPP}}$ (POD${}_{y=\mathrm{1}}$ for a–c, POD${}_{x=\mathrm{0.019}}$ for d–e) at a grid with corresponding PFT. The PFT-specific regressions between area- or mass-based POD and $R_{\mathrm{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^{\mathrm{2}}$) denoted in each panel. Note the differences of ranges in $x$ axis among PFTs. The YIBs-LMA experiment is summarized in Table 1.

[Figure omitted. See PDF]

3 Results

3.1 Comparison of simulated sensitivities with observations

Simulated relative GPP percentage ($R_{\mathrm{GPP}}$) values at global grids were sorted by dominant PFTs (Fig. S1) and plotted against area-based accumulated phytotoxic O${}_{\mathrm{3}}$ dose above a threshold $y$ nmol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ (POD${}_{y=\mathrm{1}}$) at the corresponding grids (Fig. 2). The DRR shows varied slopes among different PFTs, resulting in a coefficient of determination ($R^{\mathrm{2}}$) around 0.54 for all PFTs (Fig. 2a–c). We further calculated the mass-based accumulated phytotoxic O${}_{\mathrm{3}}$ dose above a threshold of 0.019 nmol g s${}^{-\mathrm{1}}$ (POD${}_{x=\mathrm{0.019}}$) and compared it with $R_{\mathrm{GPP}}$. The updated DRR showed convergent slopes and reached a high $R^{\mathrm{2}}$ of 0.77 across all PFTs (Fig. 2d–f), suggesting that the mass-based scheme could better unify O${}_{\mathrm{3}}$ 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_{\mathrm{S}}$) and observed ($S_{\mathrm{O}}$) slopes of O${}_{\mathrm{3}}$ DRR for all PFTs. With different $a$ parameters, the YIBs-LMA framework yielded considerably high $R^{\mathrm{2}}$ values of $\sim \mathrm{1.0}$ but varied biases between simulated and observed O${}_{\mathrm{3}}$ impacts across PFTs (Fig. 3). Both the $\mathrm{1}:\mathrm{1}$ fitting and the lowest bias between $S_{\mathrm{S}}$ and $S_{\mathrm{O}}$ were achieved with an optimal $a=\mathrm{3.5}$ nmol${}^{-\mathrm{1}}$ s g (Fig. 3c). Notably, such calibration of $a$ is robust under different O${}_{\mathrm{3}}$ fields (see Fig. S2). Consistent with observations, YIBs-LMA with this optimal $a$ parameter simulated low $S_{\mathrm{S}}$ values of $-\mathrm{0.18}\mathit{\%}$ and $-\mathrm{0.36}\mathit{\%}$ per mmol m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$ of POD${}_{y=\mathrm{1}}$ for evergreen broadleaf forest and needleleaf forest, respectively (Fig. 4a, b); median $S_{\mathrm{S}}$ values from $-\mathrm{0.53}\mathit{\%}$ per mmol m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$ for arid shrubland (Fig. 4e); and high $S_{\mathrm{S}}$ values from $-\mathrm{0.64}\mathit{\%}$ to $-\mathrm{1.04}\mathit{\%}$ per mmol m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$ for deciduous broadleaf forest, C${}_{\mathrm{3}}$ and C${}_{\mathrm{4}}$ grassland, cropland, and cold shrubland ($-\mathrm{3.28}\mathit{\%}$ for crops with POD${}_{y=\mathrm{6}}$; Fig. 4c–d, f–h).

Figure 3

Comparison between $S_{\mathrm{O}}$ (% per mmol m${}^{-\mathrm{2}}$) and $S_{\mathrm{S}}$ (% per mmol m${}^{-\mathrm{2}}$) for the YIBs-LMA experiment. Five tests from the YIBs-LMA experiment all adopt $x=\mathrm{0.019}$ nmol g${}^{-\mathrm{1}}$ s${}^{-\mathrm{1}}$ and gridded LMA from M2018 but with varied parameter $a$ from (a) 2.5 to (e) 4.5 nmol${}^{-\mathrm{1}}$ s g. $S_{\mathrm{O}}$ values are from Table S1. $S_{\mathrm{S}}$ values are derived as the slope between $R_{\mathrm{GPP}}$ and POD${}_y$. The linear regression (dashed lines), $\mathrm{1}:\mathrm{1}$ fitting (light-grey lines), normalized mean bias (NMB), and correlation coefficient ($r$) between $S_{\mathrm{S}}$ and $S_{\mathrm{O}}$ for varied PFTs are shown on each panel. The $S_{\mathrm{S}}$ and $S_{\mathrm{O}}$ of CRO are derived using POD${}_{y=\mathrm{6}}$, while other PFTs use POD${}_{y=\mathrm{1}}$. The YIBs-LMA experiment is described in Table 1.

[Figure omitted. See PDF]

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=\mathrm{0.019}$ nmol g${}^{-\mathrm{1}}$ s${}^{-\mathrm{1}}$, and $a=\mathrm{3.5}$ nmol${}^{-\mathrm{1}}$ s g. Each dot represents results from a grid cell with corresponding PFT. The regressions between relative GPP percentage ($R_{\mathrm{GPP}}$) and leaf area-based stomatal O${}_{\mathrm{3}}$ uptake fluxes (POD${}_{y=\mathrm{1}}$ for natural PFTs and POD${}_{y=\mathrm{6}}$ for crops) are shown for observations (red; see Table S1) and simulations (blue) with slopes of DRRs denoted as $S_{\mathrm{O}}$ and $S_{\mathrm{S}}$, respectively. $S_{\mathrm{O}}$ values are missing for (d) cold and (e) arid shrubs. Coefficients of determination ($R^{\mathrm{2}}$) of simulations are displayed in each panel. Note the differences of ranges in $x$ axis among PFTs (PFTs are shown in Fig. S1).

[Figure omitted. See PDF]

3.2

Global map of O${}_{\mathrm{3}}$ vegetation damage

We estimated contemporary GPP reductions induced by O${}_{\mathrm{3}}$ with the global concentrations of surface O${}_{\mathrm{3}}$ (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=\mathrm{3.5}$ nmol${}^{-\mathrm{1}}$ 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 $>\mathrm{10}\mathit{\%}$ 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${}_{\mathrm{3}}$. 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${}_{\mathrm{3}}$ vegetation damage simulated with the LMA-based scheme. Results shown are the (a) GPP reduction percentages by O${}_{\mathrm{3}}$ simulated with the YIBs-LMA framework (gridded LMA from M2018, $x=\mathrm{0.019}$ nmol g${}^{-\mathrm{1}}$ s${}^{-\mathrm{1}}$, and $a=\mathrm{3.5}$ nmol${}^{-\mathrm{1}}$ 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.

[Figure omitted. See PDF]

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

Figure 6

Global O${}_{\mathrm{3}}$-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 omitted. See PDF]

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.

[Figure omitted. See PDF]

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 ($\sim \mathrm{40}\mathit{\%}$), leading to unexpected high O${}_{\mathrm{3}}$ 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${}^{-\mathrm{1}}$ s${}^{-\mathrm{1}}$ estimated a high GPP reduction of 6.5 % (Fig. 6c) due to the overestimations of O${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ (Li et al., 2016; Feng et al., 2018). Conversely, plants with low LMA normally have thin leaves which are likely to be less O${}_{\mathrm{3}}$-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${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ sensitivity relationship can be well reproduced by our LMA-based model (Fig. 8a and b). Below we explore the linkage between O${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$-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${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ sensitivity and LMA. (a) Simulated O${}_{\mathrm{3}}$ sensitivity ($S_{\mathrm{S}}$) at each grid is compared with LMA for different PFTs. Gridded $S_{\mathrm{S}}$ is derived as GPP change per unit POD${}_{y=\mathrm{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${}_{\mathrm{3}}$ 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 omitted. See PDF]

Figure 9

Illustration of the relationships between leaf trade-off strategy and its sensitivity to O${}_{\mathrm{3}}$.

[Figure omitted. See PDF]

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${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ 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${}^{-\mathrm{2}}$), and cultivated in populated regions with high ambient O${}_{\mathrm{3}}$ concentrations. Modern technology aims to promote crop yield (Herdt, 2005), but this can potentially elevate crop sensitivities to O${}_{\mathrm{3}}$ (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${}_{\mathrm{3}}$-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${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$.

4.3

Advances in the global O${}_{\mathrm{3}}$ damage assessment

For the first time, we implemented plant trait LMA into a process-based O${}_{\mathrm{3}}$ impact modeling scheme and obtained reasonable interspecific and inter-PFT O${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ sensitivities among different PFTs. In this algorithm, the O${}_{\mathrm{3}}$ damage efficiency is inversely related to plant LMA, which influences both the O${}_{\mathrm{3}}$ uptake potential and the detoxification capability of the vegetation. The similarity in the global assessment of O${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ sensitivities. Previous studies usually categorized species into groups of low or high O${}_{\mathrm{3}}$ sensitivity, depending on very limited data from O${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ sensitivities. The LMA-based approach takes advantage of the newly revealed unifying concept in O${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ 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${}_{\mathrm{2}}$ concentration, and light resources (Kochhar and Gujral, 2020). In our modeling, the cumulative O${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ stomatal dose can be further differentiated. In fact, many other key variables in DGVMs, for example, leaf photosynthetic traits ($V_{\mathrm{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${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ 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${}_{\mathrm{3}}$ plant damage in Earth system modeling, thus facilitating the research of plant response and adaption in changing environments.

Code availability

The codes of the YIBs model with the LMA-based O${}_{\mathrm{3}}$ damaging scheme are shared at 10.5281/zenodo.6348731 (Ma and Yue, 2022).

Data availability

Results of all simulations (listed in Table 1) are available upon request. Data for figures in the main article are shared at 10.5281/zenodo.6348731 (Ma and Yue, 2022). The global maps of specific leaf area (SLA) to derive LMA for M2018 and B2017 are from https://www.try-db.org/TryWeb/Data.php#59 (last access: 20 April 2023) with DOI 10.17871/TRY.59 (Moreno-Martínez et al., 2020) and https://github.com/abhirupdatta/global_maps_of_plant_traits (Butler et al., 2017), respectively. Monthly O${}_{\mathrm{3}}$ data is from 10.5194/acp-18-8953-2018 (Turnock et al., 2018). Calibration data are summarized in Table S1.

The supplement related to this article is available online at: https://doi.org/10.5194/gmd-16-2261-2023-supplement.

Author contributions

XY, SS, and NU designed the research, and YM performed modeling, data analyses, and visualization and wrote the draft. JU, LMM, ZF, and AWC advised on concepts and methods. CG helped write the draft. HY and MCDR helped with coding. HZ, CT, YC, YL, and YX helped with data collection. All authors commented on and revised the manuscript.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Acknowledgements

We acknowledge the financial support provided by our funding sources. We also thank the editors and reviewers for their valuable feedback on the revisions.

Financial support

Xu Yue has received funding from the National Natural Science Foundation of China (grant no. 42293323) and the Jiangsu Science Fund for Distinguished Young Scholars (grant no. BK20200040). Stephen Sitch, Nadine Unger, Lina M. Mercado, and Alexander W. Cheesman were supported by NERC funding (grant no. NE/R001812/1). Johan Uddling has been supported by the strategic research area “Biodiversity and Ecosystems in a Changing Climate” (BECC).

Review statement

This paper was edited by Christoph Müller and reviewed by two anonymous referees.