Global assessment of climatic responses to ozone–vegetation interactions

The coupling between surface ozone (O₃) and vegetation significantly influences the regional to global climate. O₃ uptake by plant stomata inhibits the photosynthetic rate and stomatal conductance, impacting evapotranspiration through land surface ecosystems. Using a climate–vegetation–chemistry coupled model (the NASA GISS ModelE2 coupled with the Yale Interactive terrestrial Biosphere, or ModelE2-YIBs), we assess the global climatic responses to O₃–vegetation interactions during the boreal summer of the present day (2005–2014). High O₃ pollution reduces stomatal conductance, resulting in warmer and drier conditions worldwide. The most significant responses are found in the eastern US and eastern China, where the surface air temperature increases by +0.33 ± 0.87 and +0.56 ± 0.38 °C, respectively. These temperature increases are accompanied by decreased latent heat and increased sensible heat in both regions. The O₃–vegetation interaction also affects atmospheric pollutants. The surface maximum daily 8 h average O₃ concentrations increase by +1.46 ± 3.02 ppbv in eastern China and +1.15 ± 1.77 ppbv in the eastern US due to the O₃-induced inhibition of stomatal uptake. With reduced atmospheric stability following a warmer climate, increased cloud cover but decreased relative humidity jointly reduce aerosol optical depth by −0.06 ± 0.01 (−14.67 ± 12.15 %) over eastern China. This study suggests that vegetation feedback should be considered for a more accurate assessment of climatic perturbations caused by tropospheric O₃.

1 Introduction

Tropospheric ozone (O₃), one of the most detrimental air pollutants (Myhre et al., 2013), not only threatens human health (Norval et al., 2011; Nuvolone et al., 2018) but also induces phytotoxic effects on vegetation (Mills et al., 2007; Pleijel et al., 2007). When plants are exposed to certain levels of O₃, plant photosynthesis and stomatal conductance are inhibited because of the O₃ oxidation of cells, enzymes, and chlorophyll (Dizengremel, 2001; Fiscus et al., 2005; Jolivet et al., 2016). Consequently, carbon assimilation in terrestrial ecosystems is limited (Yue and Unger, 2014; Oliver et al., 2018), and the land–air exchange rates of water and heat fluxes are altered (Lombardozzi et al., 2015).

Experimental studies have shown that excessive O₃ exposure reduces both plant photosynthesis and stomatal conductance (Ainsworth et al., 2012; Lombardozzi et al., 2013). The reduction rates are dependent on the O₃ stomatal fluxes and the damage sensitivities, which vary among different vegetation types (Nussbaum and Fuhrer, 2000; Karlsson et al., 2004; Pleijel et al., 2004). Several exposure-based indices, such as accumulated hourly O₃ concentrations over a threshold of 40 ppb (AOT40) and the sum of all hourly average concentrations (SUM00), are used to assess O₃-induced vegetation damage (Fuhrer et al., 1997; Paoletti et al., 2007). In addition, the flux-related POD_y method (phytotoxic O₃ dose above a threshold flux of y) is also widely applied to consider the dynamic adjustment of stomatal conductance (Buker et al., 2015; Sicard et al., 2016). Considering the variability of plant sensitivities, different O₃ damage schemes have been proposed to quantify the impacts of O₃ on land carbon assimilation at regional to global scales (Anav et al., 2011; Lam et al., 2023; Lei et al., 2020). For example, Sitch et al. (2007) calculated simultaneous damage to both photosynthesis and stomatal conductance on the basis of instantaneous O₃ stomatal uptake. In contrast, Lombardozzi et al. (2012) estimated decoupled reductions in plant photosynthesis and stomatal conductance via different response relationships to cumulative O₃ stomatal uptake. The application of different schemes has resulted in a wide range of reductions in gross primary productivity (GPP) of 2 %–12 % globally, with regional hotspots of up to 20 %–30 % (Lombardozzi et al., 2015; Unger et al., 2020; Zhou et al., 2024).

O₃-induced inhibition of stomatal conductance decreases dry deposition and consequently enhances surface O₃ concentrations (Clifton et al., 2020; Wesely and Hicks, 2000; Zhang et al., 2006). Using the Sitch et al. (2007) scheme with high O₃ damage sensitivities in ModelE2-YIBs (NASA GISS ModelE2 coupled with the Yale Interactive terrestrial Biosphere model), Gong et al. (2020) revealed that O₃–vegetation interactions increased regional O₃ concentrations by 1.8 ppbv in the eastern US, 1.3 ppbv in Europe, and 2.1 ppbv in eastern China in 2010. In comparison, Sadiq et al. (2017) reported consistently stronger feedbacks on O₃ concentrations in these polluted regions via the scheme of Lombardozzi et al. (2012) embedded in the Community Earth System Model (CESM). Moreover, the inclusion of online O₃–vegetation interactions in numerical models will also result in a greater loss of simulated land carbon assimilation due to the feedbacks of both ecosystems and surface O₃. This is attributable to several factors. On the one hand, O₃ damage to leaf photosynthesis inhibits plant growth and decreases the leaf area index (LAI), leading to a greater reduction in GPP than in simulations without LAI changes (Yue et al., 2020). On the other hand, O₃ enhancement due to vegetation feedback may cause additional vegetation damage and result in further GPP losses (Lei et al., 2021). As a result, O₃–vegetation interactions should be considered in global estimates of O₃ damage to ecosystem functions.

In addition to affecting surface O₃, O₃–vegetation interactions can also alter water and energy exchanges between the land and atmosphere through the modulation of stomatal conductance. For example, Lombardozzi et al. (2015) used the Community Land Model (CLM) and estimated that the cumulative uptake of O₃ by leaves resulted in a reduction of 2.2 % in transpiration but an increase of 5.4 % in runoff globally. Arnold et al. (2018) used CESM and reported that plant exposure to O₃ could decrease land–air moisture fluxes and atmospheric humidity, which would further reduce shortwave cloud forcing in polluted regions and induce widespread surface warming up to +1.5 K. Two recent studies utilized the WRF-Chem model and revealed considerable warming and associated meteorological perturbations due to O₃–vegetation interactions in China (Zhu et al., 2022; Jin et al., 2023). However, all these modeling studies applied the same O₃ vegetation damage scheme proposed by Lombardozzi et al. (2012). It is necessary to assess the climatic responses to O₃–vegetation interactions via different schemes to explore robust responses and associated uncertainties.

In this study, we quantified the global impacts of O₃–vegetation interactions on climatic conditions and surface air pollutants during the 2010s using ModelE2-YIBs (Yue and Unger, 2015). This fully coupled framework was implemented with the semi-mechanistic O₃ damage scheme proposed by Sitch et al. (2007), which calculates aggregated O₃ damage to photosynthesis on the basis of varied sensitivities to instantaneous stomatal O₃ uptake across eight plant functional types (PFTs). We performed sensitivity experiments to quantify the responses of surface air temperature and precipitation to O₃–vegetation interactions. The feedbacks to aerosols and O₃ concentrations were also examined.

2 Method

2.1 Model descriptions

ModelE2-YIBs is a fully coupled climate–carbon–chemistry model that combines the NASA GISS ModelE2 with the YIBs vegetation model. ModelE2 is a general circulation model with a horizontal resolution of 2° × 2.5° in latitude and longitude and 40 vertical layers up to 0.1 hPa. It dynamically simulates gas-phase chemistry (NO_x, HO_x, O_x, CO, CH₄, and non-methane volatile organic compounds (NMVOCs)), aerosols (sulfate, nitrate, black and organic carbon, dust, and sea salt), and their interactions (Menon and Rotstayn, 2006). Both the physical and chemical processes are calculated every 0.5 h, and the radiation module is called every 2.5 h. The radiation module includes direct and indirect aerosol radiative effects and accounts for the absorption of multiple greenhouse gases (GHGs). For cloud optical parameters, Mie scattering, ray tracing, and matrix theory are used (Schmidt et al., 2006). The model outperforms 20 other IPCC-class climate models in simulating surface solar radiation (Wild et al., 2013) and has been extensively validated for meteorological and hydrological variables against observations and reanalysis data (Schmidt et al., 2014).

The YIBs model employs the well-established Farquhar model for leaf photosynthesis and the Ball–Berry model for stomatal conductance (Farquhar et al., 1980; Ball et al., 1987) as follows:

$$\begin{array}{}\text{(1)}& A_{\text{tot}}=min(J_{\mathrm{c}},J_{\mathrm{e}},J_{\mathrm{s}}).\end{array}$$

Here, the total leaf photosynthesis, denoted as A_tot ($\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}\left[\mathrm{leaf}\right]{\mathrm{s}}^{-\mathrm{1}}$), is calculated considering both C₃ (Collatz et al., 1991) and C₄ plants (Collatz et al., 1992). A_tot is derived from the minimum value of the constraints. The rate of carboxylation limited by ribulose-1,5-bisphosphate carboxylase (Rubisco) is J_c:

$$\begin{array}{}\text{(2)}& J_{\mathrm{c}}=\left\{\begin{array}{ll}{V}_{\text{cmax}}\left(\frac{c_{\mathrm{i}}-{\mathrm{\Gamma }}_{\ast }}{c_{\mathrm{i}}+K_{\mathrm{c}}(\mathrm{1}+O_{\mathrm{i}}/K_{\mathrm{o}})}\right)& \text{ for }{\mathrm{C}}_{\mathrm{3}}\text{ plant},\\ V_{\text{cmax}}& \text{ for }{\mathrm{C}}_{\mathrm{4}}\text{ plant}.\end{array}\right.\end{array}$$

The carboxylation rate restricted by the availability of light is J_e:

$$\begin{array}{}\text{(3)}& J_{\mathrm{e}}=\left\{\begin{array}{ll}{a}_{\text{leaf}}\times \mathrm{PAR}\times \mathit{\alpha }\times \left(\frac{c_{\mathrm{i}}-{\mathrm{\Gamma }}_{\ast }}{c_{\mathrm{i}}+\mathrm{2}{\mathrm{\Gamma }}_{\ast }}\right)& \text{ for }{\mathrm{C}}_{\mathrm{3}}\text{ plant},\\ a_{\text{leaf}}\times \mathrm{PAR}\times \mathit{\alpha }& \text{ for }{\mathrm{C}}_{\mathrm{4}}\text{ plant}.\end{array}\right.\end{array}$$

The export-limited rates for C₃ plants and the phosphoenolpyruvate carboxylase (PEPC)-limited rates of carboxylation for C₄ plants are represented by J_s:

$$\begin{array}{}\text{(4)}& J_{\mathrm{s}}=\left\{\begin{array}{ll}\mathrm{0.5}{V}_{\text{cmax}}& \text{ for }{\mathrm{C}}_{\mathrm{3}}\text{ plant},\\ K_{\mathrm{s}}\times V_{\text{cmax}}\times \frac{c_{\mathrm{i}}}{P_{\text{atm}}}& \text{ for }{\mathrm{C}}_{\mathrm{4}}\text{ plant}.\end{array}\right.\end{array}$$

In these functions, V_cmax ($\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) is the maximum carboxylation capacity. c_i and O_i (Pa) represent the internal leaf CO₂ and oxygen partial pressure, respectively. Γ_∗ (Pa) denotes the CO₂ compensation point, whereas K_c and K_o (Pa) are Michaelis–Menten constants for the carboxylation and oxygenation of Rubisco, respectively. The parameters Γ_∗, K_c, and K_o vary with temperature on the basis of the sensitivity of the vegetation to temperature (Q₁₀ coefficient). PAR ($\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) is the absorbed photosynthetically active radiation, a_leaf is the leaf-specific light absorbance that considers sunlit and shaded leaves, and α is the quantum efficiency. P_atm (Pa) represents the ambient pressure. K_s is set to 4000 as a constant following Oleson et al. (2010) to limit photosynthesis to C₄ plants that become saturated at lower CO₂ concentrations.

The stomatal conductance (g_s, $\mathrm{mol}\left[{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\right]{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) is linked to the variations in A_tot, with several parameters, such as the dark respiration rate (R_d, $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), relative humidity (RH), and CO₂ concentration at the leaf surface (c_s):

$$\begin{array}{}\text{(5)}& g_{\mathrm{s}}=m{\displaystyle \frac{(A_{\text{tot}}-R_{\mathrm{d}})\times \mathrm{RH}}{c_{\mathrm{s}}}}+b.\end{array}$$

The model simulates the biophysical processes of eight PFTs, including tundra, C₃ and/or C₄ grass, shrubland, deciduous broadleaf forest, evergreen broadleaf forest, evergreen needleleaf forest, and cropland. Different values are assigned to parameters m and b for each PFT (Table S1 in the Supplement). The carbon taken up by the leaf then accumulates and is allocated to different organs to support plant development, resulting in dynamic changes in the LAI and tree growth.

2.2 O₃–vegetation damage scheme

The YIBs model employs a semi-mechanistic parameterization proposed by Sitch et al. (2007) to estimate the impact of O₃ on photosynthesis through stomatal uptake. The scheme applies an undamaged factor (F) ($\mathrm{nmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) to calculate both A_tot and g_s as follows:

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }{A}_{\text{totd}}=A_{\text{tot}}\cdot F,\text{(7)}& {\displaystyle }& {\displaystyle }{g}_{\text{sd}}=g_{\mathrm{s}}\cdot F,\end{array}$$

where A_totd and g_sd are the unaffected photosynthesis and stomatal conductance, respectively. The factor F is defined as follows:

$$\begin{array}{}\text{(8)}& F=\mathrm{1}-a_{\mathrm{h}}\cdot max[F_{{\mathrm{O}}_{\mathrm{3}}}-F_{{\mathrm{O}}_{\mathrm{3}},\text{crit}},\mathrm{0.0}],\end{array}$$

where a_h ($\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) is the high O₃ sensitivity coefficient, calibrated by Sitch et al. (2007) on data from field observations by Karlsson et al. (2004) and Pleijel et al. (2004) to represent the “high” sensitivity of the relative species of each PFT. $F_{{\mathrm{O}}_{\mathrm{3}},\text{crit}}$ ($\mathrm{nmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) is the specific threshold for O₃ damage. Both a_h and $F_{{\mathrm{O}}_{\mathrm{3}},\text{crit}}$ vary with vegetation types (Table S1):

$$\begin{array}{}\text{(9)}& F_{{\mathrm{O}}_{\mathrm{3}}}={\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{3}}\right]}{R_{\mathrm{a}}+\left[\frac{k_{{\mathrm{O}}_{\mathrm{3}}}}{g_{\text{sd}}}\right]}},\end{array}$$

where [O₃] represents surface O₃ concentrations, and R_a (s m⁻¹) represents aerodynamic resistance, which expresses the turbulent transport efficiency in transferring sensible heat and water vapor between the land surface and a reference height. The constant $k_{{\mathrm{O}}_{\mathrm{3}}}$ = 1.67 is the ratio of stomatal resistance to O₃, which is estimated on the basis of the theoretical stomatal resistance to water (Laisk et al., 1989). When plants are exposed to [O₃] (Eq. 9), A_tot and g_s decrease (Eqs. 6 and 7) if the excess O₃ enters the leaves (Eq. 8). The increased stomatal resistance protects plants by reducing the O₃ uptake of stomata. Consequently, the damage scheme describes changes in both the photosynthetic rate and stomatal conductance.

2.3 Experiments

To explore the coupled O₃–vegetation effect, we performed two simulations using the ModelE2-YIBs model. The control experiment, O3_offline, was conducted without O₃ damage to vegetation. For comparison, the sensitivity experiment, O3_online, included online O₃–vegetation interactions with high O₃ sensitivity. For both experiments, the anthropogenic emissions from 2010 (the average of 2005–2014) of eight species (BC, OC, CO, NH₃, NO_x, SO₂, alkenes, and paraffin) from eight economic sources (agriculture, energy, industry, transportation, residential, solvent, waste, and international shipping) and biomass burning sources were collected from the Coupled Model Intercomparison Project phase 6 (CMIP6) (van Marle et al., 2017; Hoesly et al., 2018). The ensemble means of the monthly sea surface temperature (SST) and sea ice concentration (SIC) simulated by 21 CMIP6 models during the time period of 2005–2014 were employed as the boundary conditions. The cover fractions of eight PFTs (Fig. S1 in the Supplement) fixed in 2010 were adopted from the land use harmonization (LUH2) dataset (Hurtt et al., 2020). For each time slice simulation, the model was run for 30 years with all the input data fixed, and the first 10 years were used as the spin-up period. We calculated the average of the last 20 years and focused on the boreal summer season (June–July–August, JJA), when the interaction of vegetation and surface O₃ reaches its annual maximum (Fig. S3). To show the uncertainty introduced by the internal variability of the model, all the related global and regional values are denoted as mean/sum ± standard deviation of the last 20 model years. We explored the climatic responses to O₃–vegetation interactions as the differences between O3_online and O3_offline at the global scale, with a focus on hotspot regions, such as the eastern US (30–40° N, 80–90° W) and eastern China (22.5–38° N, 106–122° E).

2.4 Data for model evaluation

We evaluated the simulated air pollutants, carbon fluxes, and meteorological variables from the O3_offline simulation against observational and reanalysis datasets. The worldwide observations of the maximum daily 8 h average O₃ (MDA8 O₃) concentrations were collected from three regional networks: the Air Quality Monitoring Network operated by the Ministry of Ecology and Environment (AQMN-MEE) in China, the Clean Air Status and Trends Network (CASTNET) in the US, and the European Monitoring and Evaluation Programme (EMEP) in Europe. The observations used for validation beyond China, sourced from Sofen et al. (2016), are averaged over the period 2005–2014. This dataset encompasses 7288 station records worldwide and excludes the uncertainty associated with high mountaintop sites. For AQMN-MEE, the mean value of 2014–2018 was used because it was established in 2013. The simulated aerosol optical depth (AOD) and LAI were validated using satellite-based data from the moderate-resolution imaging spectroradiometer (MODIS) retrievals collection 5 (Remer et al., 2005) (http://modis.gsfc.nasa.gov/, last access: 5 August 2022) averaged over the years 2005–2014. The simulated GPP was evaluated against the data product upscaled from the FLUXNET eddy covariance measurements for 2009–2011 (Jung et al., 2011). The daily temperature at 2 m (T_2 m) from 2005 to 2014 was obtained from the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) reanalysis 1 (NCEP1) (Kalnay et al., 1996). For precipitation, we used the monthly data averaged from 2005 to 2014 from the Global Precipitation Climatology Project (GPCP) (Huffman et al., 1997; Adler et al., 2018). All these datasets were interpolated to the same resolution as the ModelE2-YIBs model. The root-mean-square error (RMSE) and normalized mean bias (NMB) were applied to quantify the deviations of the simulations from the observations:

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n(S_{\mathrm{i}}-O_{\mathrm{i}}{)}^{\mathrm{2}}},\text{(11)}& {\displaystyle }& {\displaystyle }\mathrm{NMB}=\sum _{i=\mathrm{1}}^n(S_{\mathrm{i}}-O_{\mathrm{i}})/\sum _{i=\mathrm{1}}^n{O}_{\mathrm{i}}\times \mathrm{100}\mathit{\%}.\end{array}$$

Here, S_i and O_i represent the simulated and observed values, respectively, and n denotes the total number of grid points used in the comparisons.

Figure 1 Evaluation of the present-day boreal summertime (June–August) air pollutants simulated by the ModelE2-YIBs model. The MDA8 O₃ (a–c) and AOD (d–f) from the O3_offline simulation (a, d) and observations (b, e) are compared. The correlation coefficient (r), root-mean-square error (RMSE), normalized mean bias (NMB), and number of grid cells (n) for the comparisons are listed on the mean bias maps (c, f).

3 Results

3.1 Control simulation and model evaluations

We first evaluated the air pollutants simulated by the control simulation (O3_offline) of the ModelE2-YIBs model (Fig. 1). Over a total of 503 grids with site-level O₃ measurements (Fig. 1b), the model replicated both the magnitude and spatial distribution of MDA8 O₃, with a correlation coefficient (r) of 0.59 and an NMB of −2.54 % (Fig. 1c). The simulated summertime surface MDA8 O₃ concentration was high in regions with high anthropogenic emissions, such as western Europe and eastern China (Ohara et al., 2007), as well as in central Africa, which has frequent fire emissions (van der Werf et al., 2017). At the global scale, the model yielded an average MDA8 O₃ concentration of 43.93 ppbv, and observations revealed an average of 44.72 ppbv over the same grids. However, the model overestimated the concentrations over the North China Plain and slightly underestimated them over the US, likely due to biases in the emission inventories and the predicted climate that drive O₃ production. The AOD simulated at 550 nm by O3_offline (Fig. 1d) showed a spatial pattern similar to that of the satellite retrievals (Fig. 1e), with r = 0.75 and an NMB of −7.35 % globally (Fig. 1f). Both simulations and observations revealed AOD hotspots over northern Africa and the Middle East, where dust emissions dominate, and in northern India and eastern China, where anthropogenic emissions are high (Feng et al., 2020).

Figure 2 Same as in Fig. 1 but for gross primary productivity (GPP; a–c) and the leaf area index (LAI; d–f).

We then evaluated the simulated GPP and LAI via a control experiment for the boreal summer period (Fig. 2). The observations revealed GPP hotspots over boreal forests, such as those in the eastern US, Eurasia, and East Asia, and over tropical forests, such as those in the Amazon, central Africa, and Indonesia (Fig. 2b). The seasonal total GPP was estimated to be 41.63 Pg [C], which accounted for 35 % of the annual amount. The simulations captured the observed GPP pattern at the global scale, with r = 0.64 and NMB = −7.81 % across 2581 grids (Fig. 2c), with underestimations in the tundra area and slight overestimations in the tropical-rainforest and evergreen-forest regions. The model simulated a seasonal total GPP of 38.69 Pg [C], equivalent to 34 % of the annual amount. The simulated LAI showed patterns similar to those of GPP (Fig. 2d) and resembled the observed LAI (Fig. 2e), with a spatial correlation of r = 0.79 and a low NMB = −5.43 % across 4435 grids globally (Fig. 2f).

We further validated the simulated meteorology from O3_offline (Fig. S2). For the surface air temperature, the model (Fig. S2a) reproduced the observed (Fig. S2b) pattern, with an RMSE of 3.21 °C and an r of 0.99 compared with the observations (Fig. S2c). For precipitation, the simulation (Fig. S2d) captured the observed spatial pattern (Fig. S2e), with NMB = 17.26 % and r = 0.75 (Fig. S2f). Overall, the model captured the spatial characteristics and magnitudes of air pollutants, biospheric parameters, and meteorological fields, making it a valuable tool for studying O₃–vegetation interactions.

Figure 3 Changes in present-day boreal-summertime biospheric variables induced by O₃–vegetation interactions. The results shown are the changes in (a) GPP, (b) canopy conductance, and (c) LAI between the O3_online and O3_offline simulations. The black dots denote areas with significant changes (p < 0.1).

Figure 4 Changes in present-day boreal-summertime meteorological fields caused by O₃–vegetation interactions. The results shown are changes in (a) surface air temperature, (b) precipitation, (c) surface relative humidity, (d) low-level cloudiness, (e) latent heat flux, and (f) sensible heat flux between the O3_online and O3_offline simulations. For heat fluxes, positive values (shaded in red) indicate that the upward fluxes change. The black dots denote areas with significant changes (p < 0.1).

Figure 5 Vertical profile of the vertical velocity. The results shown are changes in the vertical velocity in (a) eastern China and (b) the eastern US between the O3_online and O3_offline simulations. The solid red line denotes 0. Please note the differences in the scales.

3.2 O₃ damage to terrestrial ecosystems

We assessed the damaging effects of surface O₃ on ecosystems due to online O₃–vegetation interactions (Fig. 3). The impacts of O₃ on biospheric variables were mainly located in regions characterized by abundant vegetation cover and elevated O₃ concentrations. On the global scale, O₃ induced a GPP reduction of −1.80 ± 0.61 Pg C yr⁻¹ (−4.69 ± 1.56 %, Fig. 3a). This deleterious effect was more pronounced in specific regions, notably eastern China and the eastern US, with significant GPP declines of −25.40 ± 1.90 % and −20.14 ± 5.02 %, respectively, under conditions of high O₃ sensitivity (Fig. 3a and Table S2). Moreover, stomatal conductance significantly decreased in the middle latitudes of the Northern Hemisphere (Fig. 3b). The most substantial relative change of −30.62 ± 4.30 % was observed in eastern China, followed by −25.65 ± 9.32 % in the eastern US (Fig. 3b and Table S2). Although there are positive responses in some regions, they are not dominant and are hardly significant. These values were stronger than those for GPP (Fig. 3a), likely because of the climatic feedback to O₃–vegetation interactions. The opening of the plant stoma plays a crucial role in regulating energy and water exchange between the land surface and the atmosphere. The inhibition of stomatal conductance by surface O₃ leads to a warmer (Fig. 4a) and drier (Fig. 4b) climate in those hotspot regions, resulting in even stronger inhibitory effects on stomatal conductance. Following the changes in GPP, the global LAI decreased by, on average, 0.01 ± 0.01 m² m⁻² (−0.62 ± 0.84 %), with regional maxima of −4.53 ± 1.14 % in eastern China and −5.87 ± 3.11 % in the eastern US (Table S2).

3.3 Global climatic responses to O₃–vegetation interactions

In response to the O₃-induced inhibition of stomatal conductance, surface air temperature increased by 0.05 ± 0.20 °C (Fig. 4a), whereas precipitation decreased by −0.01 ± 0.03 mm d⁻¹ (Fig. 4b) at the global scale. The most significant change was the warming of 0.56 ± 0.38 °C and the precipitation reduction of −0.79 ± 1.05 mm d⁻¹ (−16.18 ± 20.38 %) in eastern China (Table S3), followed by the greatest inhibition of stomatal conductance (Fig. 3b). Such warming and rainfall deficits also appeared in the eastern US and western Europe, where the O₃–vegetation interactions were notable. The O₃-induced inhibition of stomatal conductance decreased the latent heat flux (Fig. 4e) and consequent precipitation (Fig. 4b) in those hotspot regions. Moreover, the reduced latent heat flux promoted higher surface air temperatures (Fig. 4a), resulting in an increase in the sensible heat flux (Fig. 4f). Such warming has also been reported in field experiments, where relatively high O₃ exposure resulted in noticeable increases in canopy temperature, along with reductions in transpiration (Bernacchi et al., 2011; VanLoocke et al., 2012). Globally, temperature and precipitation showed patchy responses with both positive and negative anomalies, suggesting that the regional hotspots of O₃-induced meteorological changes propagate to surrounding areas through atmospheric perturbations.

We further examined the changes in air humidity and cloud cover. The surface relative humidity decreased by −0.18 ± 0.53 % globally, with a similar pattern to that of precipitation (Fig. 4c). The most significant reductions occurred over eastern China and the eastern US, where both warming (Fig. 4a) and rainfall deficits (Fig. 4b) contributed to drought. However, in adjacent regions, such as northern China and the central US, both rainfall and surface relative humidity increased. These changes were associated with the regional increase in cloud cover (Fig. 4d). The sensible heat flux increased by 6.3 ± 5.4 W m⁻² (16.54 ± 15.59 %) and 7.12 ± 3.86 W m⁻² (25.46 ± 14.71 %) in the eastern US and eastern China, respectively, suggesting the transfer of thermal energy from the land to the atmosphere via O₃–vegetation interactions (Fig. 4f and Table S3). The warming effect further triggered anomalous updrafts in the lower troposphere, represented by changes in vertical velocity (Fig. 5), leading to enhanced convection; reduced atmospheric stability; and, consequently, an increase in low-level cloud cover (Fig. 4d). However, despite the usual cooling effect associated with increased cloud cover due to reductions in radiation, in regions predominantly influenced by O₃–vegetation interactions, this cooling effect was outweighed by O₃-induced warming through the inhibition of stomatal conductance. Therefore, temperatures exhibited an overall increase of 0.56 ± 0.38 °C in eastern China and 0.33 ± 0.87 °C in the eastern US (Table S3).

Figure 6 Changes in present-day summertime atmospheric pollution caused by O₃–vegetation interactions. The results shown are the changes in (a) O₃, (b) AOD, and (c) PM_2.5 between the O3_online and O3_offline simulations. The black dots denote areas with significant changes (p < 0.1).

3.4 Changes in air pollution caused by O₃–vegetation interactions

Changes in surface water and heat fluxes induced by O₃–vegetation interactions could feed back to affect air pollutants, such as O₃ and aerosols. As shown in Fig. 6a and Table S4, surface MDA8 O₃ concentrations increased by 1.46 ± 3.02 ppbv in eastern China and 1.15 ± 1.77 ppbv in the eastern US due to decreased dry deposition following O₃ inhibition of stomatal conductance. This indicates that high contemporary O₃ pollution may worsen air quality through O₃–vegetation interactions. However, negative O₃ changes were predicted in the central US and western China, where increased rainfall dampened O₃ through chemical reactions and wet deposition. On a global scale, surface MDA8 O₃ concentrations showed a limited increase of 0.03 ± 0.4 ppbv because of the offset between positive and negative feedbacks. The increase in O₃ concentrations in polluted regions may exacerbate the warming effect of O₃ as a greenhouse gas and may cause additional damage to vegetation. For example, the effects of offline O₃ damage on GPP in eastern China and the eastern US were simulated to be −0.52 ± 0.03 Pg [C] (−24.98 ± 0.91 %) and −0.17 ± 0.02 Pg [C] (−16.71 ± 1.16 %), respectively, which are smaller than those induced by O₃–vegetation interactions (Table S2).

Aerosols also exhibited evident changes in response to O₃–vegetation interactions. The AOD significantly decreased over hotspot regions, such as eastern China and the eastern US (Fig. 6b). In the ModelE2-YIBs model, sulfate was especially sensitive to clouds, which could enhance aerosol scavenging through cloud water precipitation (Koch et al., 2006). The large increase in cloud cover removed sulfate more efficiently than the other aerosol species did, leading to an average decrease of −1.94 ± 1.67 µg m⁻³ (−8.52 ± 6.88 %) in the PM_2.5 loading over eastern China (Fig. S4 and Table S4). Moreover, the reduction in surface relative humidity (Fig. 4c) in the regions with strong O₃–vegetation interactions limited the hygroscopic growth of aerosols, leading to a more noticeable decrease in AOD (Petters and Kreidenweis, 2007; Pitchford et al., 2007) of −0.06 ± 0.05 (−14.67 ± 16.75 %) in eastern China (Table S4). Similar aerosol changes were found in the eastern US but with smaller reductions in PM_2.5 and in AOD of −0.27 ± 0.36 µg m⁻³ (−6.01 ± 7.9 %) and −0.01 ± 0.01 (−8.15 ± 9.38 %), respectively (Table S4). In addition to the key O₃–vegetation coupling regions, positive but insignificant changes in AOD were predicted, leading to moderate AOD changes at the global scale (Fig. 6b).

4 Discussion and conclusions

We examined O₃–vegetation feedbacks to climate and air pollution in the 2010s via the fully coupled climate–carbon–chemistry model ModelE2-YIBs. During boreal summer, surface O₃ resulted in strong damage to GPP and inhibited stomatal conductance, with regional hotspots over eastern China and the eastern US. Consequently, surface transpiration was weakened, leading to decreased latent heat fluxes and relative humidity but increased surface air temperature. Moreover, surface warming increased cloud cover by reducing atmospheric stability. However, the increase in cloud cover decreased the surface temperature and promoted precipitation outside the key regions with intense O₃–vegetation interactions. The O₃-induced inhibition of stomatal conductance resulted in a localized increase in O₃ concentration. In contrast, the increased cloud cover and decreased relative humidity jointly reduced the AOD in hotspot regions. At the global scale, the mean changes in both climate and air pollution were moderate because of the offset between the changes with opposite signs.

Our predictions of the changes in water and heat fluxes caused by O₃–vegetation interactions were consistent with those of previous studies (Lombardozzi et al., 2015; Arnold et al., 2018; Gong et al., 2020). For example, simulations by Lombardozzi et al. (2015) revealed that surface O₃ reduces global GPP by 8 %–12 % and reduces transpiration by 2 %–2.4 %, with regional reductions of up to 20 % in GPP and 15 % in transpiration in eastern China and the US. These changes are generally consistent with our results, although we predicted greater reductions in transpiration than in GPP due to O₃–vegetation interactions. Using the same scheme as Lombardozzi et al. (2015), Sadiq et al. (2017) reported that O₃–vegetation coupling induced surface warming of 0.5–1 °C and O₃ enhancement of 4–6 ppbv in eastern China and the eastern US. The magnitude of these responses was much stronger than our predictions, likely because they considered the accumulation effect of O₃. In contrast, regional simulations by Jin et al. (2023) revealed that O₃–vegetation coupling led to increases in temperature and in surface O₃ of up to 0.16 °C and 0.6 ppbv, respectively, in eastern China, both of which were smaller than our predictions. The damage scheme they use, which depends on cumulative O₃ uptake, omits the difference in impact on sunlit or shaded leaves and overestimates the O₃ damage to GPP compared with the scheme we use, which considers transient O₃ flux (Cao et al., 2024). The discrepancies in O₃–vegetation feedbacks when the same O₃ damage schemes were used revealed uncertainties in the climate and chemistry models. Our predictions were within the range of previous estimates for both climatic and O₃ changes.

There were several limitations in our simulated O3–vegetation interactions. First, the semi-mechanistic O₃ damage scheme we used in the present study linked damage to photosynthesis with damage to stomatal conductance (Sitch et al., 2007), leading to a greater percentage of inhibition of stomatal conductance than of photosynthesis considering O₃–vegetation feedbacks. However, some observations have shown that damage to stomatal conductance occurs more slowly and might not be proportional to the decline in photosynthetic rates (Gregg et al., 2006; Lombardozzi et al., 2012). Second, observations have shown large variability in terms of plant sensitivity to O₃ damage. The Sitch et al. (2007) scheme employs low to high ranges of sensitivity to indicate interspecific variabilities. In this study, we employed only high O₃ sensitivity to explore the maximum responses. The possible uncertainties due to varied O₃ damage sensitivities deserve further investigation. Third, large-scale observations were not available to validate the simulated regional to global responses of climate and air pollutants. The O₃ vegetation damage scheme has been extensively validated against site-level measurements of both photosynthesis (Yue and Unger, 2018) and stomatal conductance (Yue et al., 2016). However, we were conservative about the derived global responses given that previous studies have shown large discrepancies when the same O₃ damage scheme is used while being implemented in different climate and/or chemistry models (Lombardozzi et al., 2015; Sadiq et al., 2017; Jin et al., 2023). Furthermore, the 2° × 2.5° resolution of the current version of the ModelE2-YIBs model has limitations due to the high computational demands. However, high-resolution models exhibit improved simulations of extreme events (Chang et al., 2020; Ban et al., 2021), which have certain effects on O₃–vegetation interactions (Mills et al., 2016; Lin et al., 2020). While chemical transport models with relatively coarse resolutions can increase biases in simulated air pollutants, they still capture large-scale patterns similar to those of fine-resolution results and compare reasonably well to observational data (Wang et al., 2013; Li et al., 2016; Lei et al., 2020). Moreover, we omit the slow climatic feedback caused by air–sea interactions in the simulations. Studies have revealed that these interactions may result in different climatic perturbations compared to those in simulations with fast responses by the land surface alone (Yue et al., 2011). A dynamic ocean model would enrich future research. Moreover, this study did not isolate the different impacts of aerosols, even though the radiation module included both direct and indirect radiative effects. We will investigate this further in the future by identifying the main processes involved.

Despite these uncertainties, our simulations revealed considerable changes in both climate and air pollutants in response to O₃–vegetation interactions. The most intense warming, dryness, and O₃ enhancement were predicted in eastern China and the eastern US, affecting the regional climate and threatening public health for these top two economic centers. In contrast, for the first time, we revealed a reduction in aerosol loading in hotspot regions, suggesting both positive and negative effects on air pollutants via O₃–vegetation feedbacks. Such interactions should be considered in Earth system models to better project future changes in climate and air pollutants following anthropogenic interventions against both O₃ precursor emissions and ecosystem functions.