Effects of ozone–vegetation interactions on meteorology and air quality in China using a two-way coupled land–atmosphere model

Tropospheric ozone (O₃) is one of the most important air pollutants in China and is projected to continue to increase in the near future. O₃ and vegetation closely interact with each other and such interactions may not only affect plant physiology (e.g., stomatal conductance and photosynthesis) but also influence the overlying meteorology and air quality through modifying leaf stomatal behaviors. Previous studies have highlighted China as a hotspot in terms of O₃ pollution and O₃ damage to vegetation. Yet, few studies have investigated the effects of O₃–vegetation interactions on meteorology and air quality in China, especially in the light of recent severe O₃ pollution. In this study, a two-way coupled land–atmosphere model was applied to simulate O₃ damage to vegetation and the subsequent effects on meteorology and air quality in China. Our results reveal that O₃ causes up to 16 % enhancement in stomatal resistance, whereby large increases are found in the Henan, Hebei, and Shandong provinces. O₃ damage causes more than 0.6 µmol CO₂ m⁻² s⁻¹ reductions in photosynthesis rate and at least 0.4 and 0.8 g C m⁻² d⁻¹ decreases in leaf area index (LAI) and gross primary production (GPP), respectively, and hotspot areas appear in the northeastern and southern China. The associated reduction in transpiration causes a 5–30 W m⁻² decrease (increase) in latent heat (sensible heat) flux, which induces a 3 % reduction in surface relative humidity, 0.2–0.8 K increase in surface air temperature, and 40–120 m increase in boundary-layer height in China. We also found that the meteorological changes further induce a 2–6 ppb increase in O₃ concentration in northern and south-central China mainly due to enhanced isoprene emission following increased air temperature, demonstrating that O₃–vegetation interactions can lead to strong positive feedback that can amplify O₃ pollution in China. Our findings emphasize the importance of considering the effects of O₃ damage and O₃–vegetation interactions in air quality simulations, with ramifications for both air quality and forest management.

1 Introduction

Tropospheric ozone (O₃) is a secondary air pollutant, which is mainly formed from the photochemical oxidation of carbon monoxide (CO), methane (CH₄), and non-methane volatile organic compounds (VOCs) by hydroxyl radicals (OH) in the presence of nitrogen oxides (NO_x= NO + NO₂). O₃ is known as the third most important greenhouse gas, with an estimated radiative forcing of 0.41 W m⁻² for the period of 1750–2010 (IPCC, 2013; Stevenson et al., 2013). As an air pollutant, O₃ is also shown to be harmful to not only human health but also vegetation and crop health (Anenberg et al., 2010; Cohen et al., 2017). Various field experiments and numerical modeling studies have already demonstrated that O₃ can not only reduce gross primary production (GPP) of natural vegetation as well as crop yields (Ainsworth et al., 2012; Lombardozzi et al., 2012; Tai e al., 2014; Feng et al., 2015; Yue et al., 2017; Li et al., 2018) but also decrease transpiration (Arnold et al., 2018), decrease runoff (Li et al., 2016) on larger scales, and therefore affect the global carbon and water cycle (Lombardozzi et al., 2015).

Vegetation can in turn modulate O₃ concentration through influencing the sources and sinks of O₃. Dry deposition of O₃ onto vegetation is a major sink for O₃, mainly via stomatal uptake. Stomata are the pores on plant leaves; they control water exiting and carbon entering the leaf interior and hence influence the water and carbon exchange between the land and atmosphere. When vegetation is exposed to enhanced O₃ levels, cellular and tissue damage can result in a decrease in photosynthesis rate, thus altering CO₂ assimilation. Stomata conductance may decrease subsequently in response to O₃ exposure, thus reducing the dry-depositional sink of O₃ (Sadiq et al., 2017; Zhou et al., 2018), but some studies also suggest that O₃ exposure can cause stomata to respond more sluggishly to changing environmental conditions, such as drought, with complex overall effects on stomatal behaviors and dry deposition (e.g., Huntingford et al., 2018). Moreover, recent studies showed reduced dry deposition velocities of O₃ by drought-stressed vegetation, which affects surface O₃ trends and extremes (Huang et al., 2016; Lin et al., 2019, 2020). Vegetation also affects the sources of O₃; the most abundant biogenic VOC (BVOC) species emitted by vegetation is isoprene (C₅H₈), which is a major precursor for O₃ formation in polluted, high-NO_x environments, but removes O₃ by ozonolysis or by sequestering NO_x in more pristine, low-NO_x regions (Hollaway et al., 2017). Isoprene production is known to be highly coupled with photosynthesis and by extension to stomatal conductance (Arneth et al., 2007). Moreover, transpiration, which is modulated by stomatal behaviors, significantly regulates surface meteorology including water vapor content and air temperature, which further influence the production and loss of O₃. Therefore, through influencing plant ecophysiology (e.g., photosynthesis and stomata behaviors), O₃–vegetation interactions can modulate boundary-layer meteorology and climate, and may further affect O₃ air quality via a series of feedback mechanisms. It is therefore essential to fully understand the O₃–vegetation interactions and the following climatic and biospheric impacts especially in areas with high O₃ concentrations and vegetation density.

In many land surface and biospheric models, such as Noah-Multi Parameterization (Noah-MP) or Community Land Model (CLM), the Farquhar–Ball–Berry model (FBB, Farquhar et al., 1980; Ball et al., 1987) is commonly used to simulate stomatal conductance and photosynthetic rate. In the FBB model, the calculation of stomata conductance is based on the calculation of photosynthesis, which makes them tightly coupled with each other. Therefore, in several land surface models that consider O₃ damage effect on vegetation, the photosynthetic rate is modified first and the stomatal conductance is modified subsequently, which means stomata conductance and photosynthesis will change collinearly under chronic O₃ exposure (Sitch et al., 2007; Yue and Unger, 2014). However, field experiments have shown that, under chronic O₃ exposure, stomata conductance decreases with a smaller magnitude than photosynthetic rate does, which makes the simulations of stomata conductance and photosynthetic rate as well as the following water and carbon cycles in the above models less accurate (Lombardozzi et al., 2012). Modifying stomata conductance and photosynthesis separately in land surface models is therefore more reasonable. Lombardozzi et al. (2012) modified the stomata conductance and photosynthetic rate separately based on the cumulative uptake of O₃ into leaves and has shown a better representation of plant responses to O₃ exposure. Efforts have been made to investigate the effects of O₃ exposure on land biosphere based on the above O₃ damage schemes. For example, based on an offline process-based vegetation model, Yue and Unger (2014) found that O₃ damage decrease GPP by 4 %–8 % on average in the eastern US and leads to significant decreases of 11 %–17 % in east coast hotspots. Using the offline CLM model, Lombardozzi et al. (2015) estimated that the present O₃ exposure reduces GPP and transpiration globally by 8 %–12 % and 2.0 %–2.4 %, respectively.

Several modeling studies conducted so far have demonstrated the importance of considering the interactions and feedbacks between atmosphere and biosphere. By dynamically coupling O₃ and leaf area index (LAI) but without considering the meteorological feedbacks of O₃–vegetation interactions to O₃, Zhou et al. (2018) found that O₃-induced damage on LAI can lead to changes in O₃ concentrations by −1.8 to +3 ppb in boreal summer. By considering the interactions between atmospheric chemistry with biosphere in a two-way coupling model, Lei et al. (2020) quantified the damaging effects of O₃ on vegetation and found a global reduction of annual GPP by 1.5 %–3.6 %, with regional extremes of 10.9 %–14.1 % in the eastern US and eastern China. Based on the Community Earth System Model (CESM) model with fully interactive atmospheric chemistry, biogeochemical, and biogeophysical cycles, Sadiq et al. (2017) estimated that surface O₃ is 4–6 ppb higher in Europe, North America, and China in simulations with O₃–vegetation coupling comparing the surface O₃ concentrations without O₃–vegetation coupling. Based on the modified Weather Research and Forecasting model with chemistry (WRF-Chem), Li et al. (2016, 2018) investigated the effect of O₃ exposure on hydroclimate and crop productivity in the US and highlighted O₃ damage effects on meteorological fields and surface energy balance as well as the crop yields, but the feedbacks of changing meteorology onto surface O₃ were not investigated. Arnold et al. (2018) examined the global climate response to O₃ exposure and found O₃ damage on vegetation can induce widespread surface warming and changes in clouds, which could be critical on regional scales. Although the interactions between O₃ and vegetation are critical to our environment, adequate representation of O₃–vegetation interactions is still missing in most atmospheric models used for climate and atmospheric chemistry simulations, at least in part due to incomplete coupling capacities with land surface or biospheric model components at high resolutions and in part due to limited observations to optimize O₃ damage schemes for wider regional applicability.

With the rapid urbanization and industrialization in the recent decades, China has experienced increasingly severe O₃ pollution, which is expected to continue to worsen in the near future. O₃ concentration in China has been observed to exceed ambient air quality standard by 100 %–200 % (Wang et al., 2017), with the maximum 8 h mean concentration of O₃ (MDA8 O₃) increasing by 4.6 % per year from 2015 to 2017 (Silver et al., 2018). Lu et al. (2019) showed that urban surface O₃ in China during 2013–2017 was significantly higher than that in other regions around the world, and thus vegetation exposure to O₃ is also higher in China. Li et al. (2019) also revealed the increasing trend of O₃ in megacity clusters of China during 2013–2017, which is closely related to meteorology, anthropogenic emissions, and PM_2.5 concentrations. Global-scale studies have highlighted China as a hotspot of O₃ pollution and damage to vegetation compared with other regions (Sadiq et al., 2017; Arnold et al., 2018; Lei et al., 2020). However, a comprehensive study of how O₃ affects meteorology and air quality through O₃–vegetation interactions in China at high spatial resolutions, especially under severe O₃ pollution, is still limited but highly needed. Moreover, there have been limited studies focusing on the feedbacks of O₃–vegetation coupling on O₃ concentration itself, especially in China, which is one of the main scopes of our study.

This study, therefore, first adopted and implemented a semi-mechanistic O₃ damage scheme in a widely used regional atmosphere–land modeling framework and hence used it to simulate and assess the impacts of O₃–vegetation interactions on boundary-layer meteorology and air quality in China at a high spatial resolution. Specifically, O₃-induced damage to vegetation, changes in meteorology in China due to O₃–vegetation coupling, and the subsequent feedback effects onto O₃ concentration itself are examined, which is crucial to fully understand the O₃–vegetation interactions and the following impacts on climate, biosphere, and air quality in areas with both high O₃ concentrations and high vegetation coverage.

2 Methods

2.1 WRF-Chem model setup

The WRF model is a state-of-the-art mesoscale nonhydrostatic meteorological model. An atmospheric chemistry module that includes various gas-phase chemistry and aerosol mechanisms has been implemented into and fully coupled with WRF to create the WRF-Chem model (Grell et al., 2005; Fast et al., 2006). In WRF-Chem, both the air quality and meteorological components use the same transport scheme, model grid, subgrid-scale transport physics, and time step. WRF-Chem has been widely used in previous air quality studies (e.g., Li et al., 2016, 2018; Liu et al., 2018, 2020). In this study, we applied our revised WRF-Chem model based on version 3.8.1 to simulate meteorological fields and O₃ concentration over China. Simulations are conducted from 24 May to 1 September every year from 2014 to 2017 and the days in May were discarded as spin-up. For the land surface component within WRF, we used Noah-MP, which will be described in the next subsection.

The model domain was configured at a horizontal resolution of 27 km on the Lambert conformal projection, centered at 37^∘ N, 108.1^∘ E, and covering all of China. The model has 26 vertical layers, with the lowest layer at 0.17 km and the highest layer at 17.67 km. The meteorological initial and boundary conditions are provided by the 6-hourly Final Operational Global Analysis (FNL) dataset at a horizontal resolution of $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$. The chemical initial and boundary conditions were generated from the Model for Ozone and Related Chemical Tracer version 4 (MOZART-4), which is available at a horizontal resolution of $\mathrm{1.9}{}^{\circ }\times \mathrm{2.5}{}^{\circ }$ with 56 vertical layers (Emmons et al., 2010).

Anthropogenic emissions were from the Multi-resolution Emission Inventory for China (MEIC) compiled at a spatial resolution of 27 km and a 1-hourly temporal resolution suitable for our research domain. Biogenic emissions were calculated online by the Model of Emissions of Gases and Aerosol from Nature (MEGAN) (Guenther et al., 2006). Biomass burning emissions were extracted from the Fire Inventory from NCAR (FINN) version 1.5 datasets (Wiedinmyer et al., 2011). Dust emissions were generated online by the Goddard Global Ozone Chemistry Aerosol Radiation and Transport model (GOCART; Ginoux et al., 2001). Gas-phase chemistry was simulated with second generation Regional Acid Deposition Model (RADM2; Stockwell et al., 1990) mechanism, and the Modal Aerosol Dynamics Model for Europe (MADE; Ackermann et al., 1998), which is coupled with the Secondary Organic Aerosol Model (SORGAM; Schell et al., 2001) for aerosol treatment. Detailed physics schemes used in the simulations are shown in Table S1 in the Supplement.

2.2 Description of the Noah-MP model

Noah-MP is a land surface model that uses multiple options for key land–atmosphere interaction processes (Niu et al., 2011). Noah-MP contains a separate vegetation canopy defined by a canopy top and bottom, crown radius, and leaves with prescribed dimensions, orientation, density, and radiometric properties. The canopy employs a two-stream radiation transfer approach along with shading effects necessary to achieve proper surface energy and water transfer processes (Dickinson, 1983). Noah-MP is capable of distinguishing between C₃ and C₄ photosynthesis pathways and defines vegetation-specific parameters for plant photosynthesis and respiration.

Noah-MP is available for prognostic vegetation growth that combines a Ball–Berry photosynthesis-based stomatal resistance (Farquhar et al., 1980; Ball et al., 1987) that allocates carbon to various parts of vegetation (leaf, stem, wood, and root) and soil carbon pools (fast and slow). GPP, LAI, and canopy height are then predicted downstream from photosynthesis. Noah-MP also considers the photosynthesis of sunlit and shaded leaves separately, whereby sunlit leaves are more limited by CO₂ concentration while shaded leaves are more constrained by insolation, which may thus have different responses to O₃ damage. The dynamic LAI and canopy height calculation will further affect surface energy fluxes, which will then affect the boundary-layer meteorology when coupling with the atmosphere model in WRF-Chem. The land use types and the vegetation parameters are based on the US Geological Survey (USGS) embedded in Noah-MP. Figure 1 shows the spatial distribution of vegetation fraction of dominant vegetation types in China. The distributions of main vegetation groups (broadleaf, needleleaf, crop, and grass) that have different sensitivities to O₃ damage following Lombardozzi et al. (2015) are shown in Fig. 1.

In this study, the O₃ concentration simulated by the chemical module of the WRF-Chem model was also dynamically passed onto the Noah-MP land surface model at every time step to modify the photosynthesis and stomatal conductance due to O₃ damage. The land surface variables simulated by Noah-MP were also dynamically passed back onto the atmospheric components, thus allowing immediate, two-way feedback effects onto meteorological fields, O₃, and other atmospheric chemical constituents. In this way, land surface processes, atmospheric dynamics, and atmospheric chemistry in the WRF-Chem model were fully coupled.

Figure 1 The vegetation fraction of (a) broadleaf, (b) needleleaf, (c) cropland, (d) grass, (e) others, and (f) dominant vegetation types.

2.3 O₃ damage parameterization

In Noah-MP, the Farquhar model (Farquhar et al., 1980) was used to calculate photosynthetic rate, whereas Ball–Berry model was used to calculate stomatal conductance (Ball et al., 1987). The photosynthesis rate, A (µmol CO₂ m⁻² s⁻¹), is calculated separately for sunlit and shaded leaves and is limited by either one of three limiting factors and can be calculated as

$$\begin{array}{}\text{(1)}& A=\min \left(W_c,W_j,W_{\mathrm{e}}\right)I_{\mathrm{gs}},\end{array}$$

where W_c is the RuBisCO-limited photosynthesis rate, W_j is the light-limited photosynthesis rate, and W_e is the export-limited photosynthesis rate. I_gs is the growing season index with values ranging from 0 to 1. Stomatal conductance (g_s) is computed based on the photosynthesis rate from the Farquhar model as

$$\begin{array}{}\text{(2)}& g_{\mathrm{s}}={\displaystyle \frac{\mathrm{1}}{r_{\mathrm{s}}}}=m{\displaystyle \frac{A}{c_{\mathrm{s}}}}{\displaystyle \frac{e_{\mathrm{s}}}{e_{\mathrm{i}}}}{P}_{\mathrm{atm}}+b,\end{array}$$

where g_s is the leaf stomatal conductance (µmol m⁻² s⁻¹); r_s is the leaf stomatal resistance (s m² µmol⁻¹); m is an empirical parameter that relates stomatal conductance and photosynthesis with values ranging from 5 to 9; A is the photosynthesis rate as described above; c_s is the CO₂ partial pressure at the leaf surface (Pa); e_s is the vapor pressure at the leaf surface (Pa); e_i is the saturation vapor pressure inside the leaf (Pa); P_atm is the atmospheric pressure (Pa); and b is the minimum stomatal conductance.

As mentioned above, following Lombardozzi et al. (2015), an O₃ damage scheme was implemented in Noah-MP embedded in WRF-Chem model version 3.8.1. The photosynthesis rate and stomatal conductance are modified independently using two sets of O₃ impact factors, $F_{{\mathrm{pO}}_{\mathrm{3}}}$ and $F_{{\mathrm{cO}}_{\mathrm{3}}}$, respectively, which are then multiplied to the initial A and g_s calculated by the Farquhar–Ball–Berry model, respectively. Lombardozzi et al. (2012) found that independently modifying stomatal conductance and photosynthesis can improve the model prediction of plant response to O₃ damage. The two damage factors are calculated based on the cumulative uptake of O₃ (CUO), which integrates the O₃ flux inside leaves through the stomata throughout the growing season. The CUO (mmol m⁻²) is calculated as

$$\begin{array}{}\text{(3)}& \mathrm{CUO}={\mathrm{10}}^{-\mathrm{6}}\sum {\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{3}}\right]}{k_{{\mathrm{O}}_{\mathrm{3}}}{r}_{\mathrm{s}}+r_{\mathrm{a}}+r_{\mathrm{b}}}}\mathrm{\Delta }t,\end{array}$$

where [O₃] is the surface O₃ concentration (nmol m⁻³); $k_{{\mathrm{O}}_{\mathrm{3}}}=\mathrm{1.61}$ is the ratio of leaf resistance to O₃ to leaf resistance to water (Uddling et al., 2012); r_s is the stomatal resistance, r_a is the aerodynamic resistance and r_b is the boundary-layer resistance (s m⁻¹); Δt is the model time step (s). CUO is only accumulated when LAI is larger than 0.4 and O₃ flux is larger than a threshold value of 0.8 nmol O₃ m⁻² s⁻¹ to consider the detoxification effect of plants to O₃ damage.

The two damage factors have linear relationships with CUO and can be calculated as follows:

$$\begin{array}{}\text{(4)}& {\displaystyle }{F}_{{\mathrm{pO}}_{\mathrm{3}}}=a_p\times \mathrm{CUO}+b_p\text{(5)}& {\displaystyle }{F}_{{\mathrm{cO}}_{\mathrm{3}}}=a_c\times \mathrm{CUO}+b_c,\end{array}$$

where $F_{{\mathrm{pO}}_{\mathrm{3}}}$ is the O₃ damage factor for photosynthesis and $F_{{\mathrm{cO}}_{\mathrm{3}}}$ is the O₃ damage factor for stomatal conductance; a_p, b_p, a_c, and b_c are empirical slopes and intercepts of three different plant groups (broadleaf trees, needleleaf trees, and grasses or crops) from Lombardozzi et al. (2015). The values of these slopes and intercepts are shown in Table 1. The original photosynthesis and stomatal conductance are then multiplied by the two damage factors, respectively, to get the modified photosynthesis and stomatal conductance under O₃ exposure.

Table 1 Slopes (per mmol m⁻²) and intercepts (unitless) used for O₃ damage factors in Eqs. (4) and (5), following Lombardozzi et al. (2015).

2.4 Model experiments and evaluation

Two sets of experiments were conducted in this study. We performed a control simulation (simu_withoutO₃) without O₃ damage on vegetation and a production simulation (simu_withO₃) with O₃ damage on vegetation. Detailed information of the experiments is shown in Table 2. In the simu_withO₃ experiment, the O₃ concentration simulated by the chemical module of the model is dynamically passed onto the land surface model at every time step to modify the photosynthesis and stomatal conductance. The differences between the two sets of experiments including vegetation physiology, meteorological fields, and O₃ concentration can thus be attributed to O₃–vegetation interactions. In this work, each simulation was conducted from 24 May to 1 September every year from 2014 to 2017 and the days in May was discarded as spin-up. For each simulation in the 4 years, anthropogenic emissions were kept at 2014 levels, while meteorological fields were changing every year. The 4-year June–July–August (JJA) averaged results were analyzed and compared. JJA was selected because of the most severe O₃ pollution in this season and because it is within the active growing season of the plants.

Table 2 Description of the two sets of model experiments. ICs are initial conditions; BCs are boundary conditions.

The simulated meteorological variables and air pollutant concentrations were evaluated using available in situ observations in China. The daily meteorological observations including temperature at 2 m (T_2 m), relative humidity at 2 m (RH_2 m), and wind speed at 10 m (WS_10 m) above displacement height were from the National Meteorological Information Center. There are 698 stations in the study domain. The air pollutant observations were provided by the China National Environmental Monitoring Center (CNEMC) network, which offers hourly concentrations of particulate matter with an aerodynamic diameter of less than 2.5 µm (PM_2.5) and 10 µm (PM₁₀), carbon monoxide (CO), O₃, sulfur dioxide (SO₂), and nitrogen dioxide (NO₂). The locations of meteorological stations and the sites of CNEMC network are shown in Fig. 2. The statistical parameters including mean values (mean) of observations and simulated variables, their standard deviations (SDs), indices of agreement (IOAs), mean biases (MBs), and correlation coefficients (CORRs) were computed to evaluate the model performance in this study.

Figure 2 Site locations of air quality monitoring sites (blue dots) and the meteorological monitoring sites (pink dots), and terrain height (m) shown by the color contours.

3 Results

3.1 Model evaluation

Table 3 shows the city-averaged evaluation results of meteorological variables from the modified model. The information of the major cities used for evaluation is shown in Table S4. From Table 3, we can find that T_2 m is underestimated with MB values ranging from −1.00 ^∘C in 2017 to −0.70 ^∘C in 2014. The IOA and CORR are generally higher than 0.8, indicating that the model could reasonably simulate the variations of T_2 m. Unlike temperature, relative humidity is overestimated by the model simulations with MB values ranging from 4.38 in the year 2014 to 7.33 in the year 2016, but the CORR values with observations are still high (CORR>0.7). Wind speed is also overestimated by more than 0.38 m s⁻¹, which might be caused by the underestimation of terrain height as reported in other WRF modeling studies (Brunner et al., 2015; Liu et al., 2020). The detailed evaluation results for each city and for seven major geographic regions of China are shown in Tables S5–S10. The classification of the geographic regions is shown in Fig. S2. As shown in these tables, the model can reasonably capture the spatial distribution of these meteorological variables. For example, the larger values of T_2 m and RH_2 m in cities from southern China compared with the cities in northern China (Table 4) can be reasonably simulated. We also found that the model simulations have better performance in northeastern China, central China, and southern China in terms of IOA and CORR as shown in these tables (Table 4).

Table 3 Evaluation results for the temperature at 2 m (T_2 m), relative humidity at 2 m (RH_2 m), and wind speed at 10 m (WS_10 m) for different years in China. Mean_obs (Mean_simu) is the mean value of observation (model simulation); SD_obs (SD_simu) is the standard deviation of the observation (model simulation); IOA is the index of agreement; CORR is the correlation coefficient; MB is the mean bias.

Table 4 Evaluation results of temperature at 2 m (T_2 m), relative humidity at 2 m (RH_2 m), and wind speed at 10 m (WS_10 m) in seven major geographic regions from the implemented model. NEC is northeast China, NC is north China, CC is central China, EC is east China, SC is south China, SWC indicates southwest China, and NWC is northwest China. Mean_obs (Mean_simu) is the mean value of observations (model simulations); SD_obs (SD_simu) is the standard deviation of the observations (model simulations); IOA is the index of agreement; CORR is the correlation coefficient; MB is the mean bias.

Table 5 shows the city-averaged evaluation results of six air pollutants simulated from the modified model. The information of the major cities used for air pollutant evaluation is shown in Table S11. Form Table 5, positive MB values for O₃, PM_2.5, SO₂, and NO₂, and negative MB values for CO are found. The overestimation of O₃ by WRF-Chem was also reported by Hu et al. (2016) and Gao et al. (2020). For PM₁₀, both positive and negative MB values are found for different years. The results indicate general overestimation by the model of most air pollutants except for CO. The underestimation of CO can be explained by either O₃ chemistry, which points to the problem related to low titration, or in the underestimation of dry deposition by the model, which is also affected by the modification of the model. The IOA of air pollutant concentration ranges from 0.36 (SO₂) to 0.63 (O₃). The correlation coefficient of air pollutants ranges from 0.14 (PM₁₀) to 0.66 (O₃). Detailed evaluation results for each city and major geographic regions of China are shown in Tables S9–S14 and Table 6. In terms of the evaluation for O₃, the model has better performance in northeastern China, eastern and southern China, which may suffer the most severe O₃ damage. Our results are generally consistent with the evaluation results of the Community Multiscale Air Quality Modeling (CMAQ) simulation over China by Liu et al. (2020). MBs of SO₂, NO₂, and CO are consistent in both magnitude and sign with Liu et al. (2020), while the MBs of PM and O₃ are larger than Liu et al. (2020). Correlation coefficients of air pollutants are also of similar magnitude with Liu et al. (2020), showing that our model results can well capture the temporal variations of air pollutants. We also compared the evaluation results between the original model and the modified model, as shown in Tables S2 and S3 in the Supplement and Tables 3 and 5 here. We found no obvious differences in the evaluation results between the original model results and the revised model results. It should be noted that this study might not be able to and was not meant to improve model accuracy, but our modified model is able to capture O₃–vegetation interactions without worsening model performance. Overall, there are systematic biases in simulated variables especially the air pollutant concentrations, but the spatial distributions of both meteorological variables and air pollutant concentrations are reasonably simulated by the model, lending trust to the use of the model for sensitivity studies to examine the effects of O₃–vegetation interactions on the atmospheric environment.

Table 5 Evaluation results for the air pollutants in China. Mean_obs (Mean_simu) is the mean value of observation (model simulation); SD_obs (SD_simu) is the standard deviation of the observation (model simulation); IOA is the index of agreement; CORR is the correlation coefficient; MB is the mean bias.

Table 6 Evaluation results of air pollutants in seven major geographic regions simulated by the implemented model. NEC is northeast China, NC is north China, CC is central China, EC is east China, SC is south China, SWC indicates southwest China, and NWC is northwest China. Mean_obs (Mean_simu) is the mean value of observations (model simulations); SD_obs (SD_simu) is the standard deviation of the observations (model simulations); IOA is the index of agreement; CORR is the correlation coefficient; MB is the mean bias.

3.2 Responses of vegetation to O₃ damage

O₃ can adversely affect photosynthesis rate and stomatal conductance and therefore interfere with vegetation growth, productivity, and transpiration. To understand the O₃-induced damage on vegetation physiology, the spatial distribution, and changes in stomatal resistance (RS), photosynthesis rate (PSN), LAI, GPP, and transpiration rate (TR) during 2014–2017 summer (June–July–August) were analyzed.

Figure 3a and d display the spatial distribution of sunlit stomatal resistance (RSSUN) and shaded stomatal resistance (RSSHA) from the simu_withoutO₃ experiment, respectively. The absolute and relative changes in RSSUN (RSSHA) between simu_withO₃ and simu_withoutO₃ experiments are shown in the middle and the right panel of Fig. 3, separately. In general, simulated stomatal resistance in eastern China is larger than that in western China. Both RSSUN and the RSSHA are enhanced in response to O₃ damage to vegetation. The maximum increases in RSSUN and RSSHA can be up to 1.0×10³ s m⁻¹, which is equivalent to a ∼16 % increase compared to the simu_withoutO₃ simulation. Comparing the changes in RSSUN vs. RSSHA, the changes in RSSHA are larger than that in RSSUN, reflecting the larger sensitivity of shaded leaves to O₃ damage (Kinose et al., 2017). Northern China experiences larger changes in stomatal resistance generally, especially in the Henan, Hebei, and Shandong provinces, where the changes in stomatal resistance are twice as large as the changes in stomatal resistance over other regions.

Figure 3 Spatial distribution of mean stomatal resistance in JJA of 2014–2017 for (a) sunlit leaves (RSSUN) and (d) shaded leaves (RSSHA) from the simu_withoutO₃ experiment. Absolute changes in (b) RSSUN and (e) RSSHA caused by O₃ damage. Relative changes in (c) RSSUN and (f) RSSHA caused by O₃ damage. Absolute changes are the RSSUN (RSSHA) from simu_withO₃ minus RSSUN (RSSHA) from simu_withoutO₃. Relative changes are calculated by absolute changes over the RSSUN (RSSHA) from simu_withoutO₃.

The spatial distribution of 2014–2017 JJA mean PSN, LAI, and GPP from the simu_withoutO₃ simulations and their changes induced by O₃ damage are presented in Fig. 4. From Fig. 4a, we find that the PSN values are generally higher in eastern China compared with western China with the largest values of up to ∼7 µmol CO${}_{\mathrm{2}}^{-\mathrm{1}}$ m⁻² s⁻¹. Similar spatial distribution and hotspot areas can also be observed for LAI (Fig. 4d) and GPP (Fig. 4g), with LAI and GPP values in hotspot areas up to 3.6 and 10 g C m⁻² d⁻¹, respectively. We also find that the Henan, Hebei, Shanxi, and Shandong provinces have smaller values of PSN, LAI, and GPP when compared with other provinces in eastern China.

Figure 4 Spatial distribution of 2014–2017 JJA mean (a) photosynthesis rate (PSN), (d) leaf area index (LAI), and (g) gross primary productivity (GPP) from the simu_withoutO₃ experiment; absolute changes in (b) PSN, (e) LAI, and (h) GPP caused by O₃ damage; and relative changes in (c) PSN, (f) LAI, and (i) GPP caused by O₃ damage. Absolute changes are the results from simu_withO₃ minus results from simu_withoutO₃. Relative changes are calculated from the absolute changes over the results from simu_withoutO₃.

With O₃ damage, PSN decreases in general, with absolute changes in PSN ranging from 0.6 to 3.6 µmol CO${}_{\mathrm{2}}^{-\mathrm{1}}$ m⁻² s⁻¹ (Fig. 4b), representing 20 %–40 % reductions in PSN. For northeastern and southern China, where the original PSN values are large, ∼20 % reductions in PSN are found (Fig. 4c). In western China where the dominant vegetation type is grassland and the original PSN values are small, more than 40 % of PSN is reduced due to O₃ damage (Fig. 4c). In response to the PSN reductions, LAI and GPP also decrease. More than 0.4 reductions in LAI are found in central and northern China (Fig. 4e), corresponding to more than 20 % reductions in LAI; in other regions, 5 %–15 % reductions in LAI are observed. More than 0.8 g C m⁻² d⁻¹ reductions in GPP are found generally in China. Similar to Fig. 3c, we find that GPP decreases by ∼20 % in northeastern and southern China and decreases by more than 40 % in other regions (Fig. 4i). Based on offline models without considering atmosphere–biosphere coupling, O₃ damage was found to decrease GPP at most by 11 %–17 % in the east coast hotspots of the US (Yue and Unger, 2014). Using the offline CLM model, Lombardozzi et al. (2015) estimated that the present O₃ exposure reduces GPP globally by 8 %–12 %. Based on the Regional Climate-Chemistry Model version 4 (RegCM-CHEM4) model coupled with Yale Interactive terrestrial Biosphere (YIBs) model, Xie et al. (2019) revealed that O₃ damage induces a significant reduction (12.1±4.4 %) in the GPP, up to 35 % in summer over China (Table S15). Comparing our results with previous studies, our results are broadly consistent with Xie et al. (2019), but the magnitude is larger than the studies conducted by Yue and Unger (2014) and Lombardozzi et al. (2015). Differences or uncertainties may arise from the different model settings. It appears that offline models as used by Yue and Unger (2014) and Lombardozzi et al. (2015) generally found smaller damage than studies with two-way coupling between the atmosphere and biosphere as used by Xie et al. (2019) and our work; this could be due to the existence of positive biosphere–atmosphere feedbacks that potentially worsen O₃ damage, as will be discussed in subsequent sections. Different O₃ damage schemes employed in the models may also be a source of differences, although we note that both this work and Lombardozzi et al. (2015) used the same scheme, so the differences appear to arise more likely from the effect of coupling and other model settings than from the schemes alone.

The spatial distributions of dominant vegetation types in China are shown in Fig. 1, where we can see that the croplands dominant in eastern China and especially in southern China suffer the greatest GPP reductions, indicating that crop yields in China would also be heavily affected by O₃ damage.

Figure 5 depicts the spatial distribution of transpiration rate (TR) of vegetation and the changes in transpiration rate induced by O₃ damage. TR values are higher in eastern China, where there is larger vegetation coverage (Fig. 5a). As shown in Fig. 5b, TR deceases by 0.2–1.0 mm d⁻¹ generally in eastern China with large reductions in northern China, especially in the Henan, Shandong, Anhui, and Jiangsu provinces. In terms of relative changes, TR decreases by ∼12 % in northeastern and southern China, while more than 24 % reductions are found in other regions. Transpiration is affected by the changes in both RS and LAI. With O₃ damage, both the increases in RS (Fig. 3c and f) and decreases in LAI (Fig. 4f) cause TR to decrease, as shown in Fig. 5b and c. Comparing the changes in RS (Fig. 3c and f), LAI (Fig. 4f) and TR (Fig. 5c), we can find that the distribution of changes in TR is more consistent with that of RS, reflecting the dominance of RS in controlling TR.

Figure 5 Spatial distribution of 2014–2017 JJA mean (a) transpiration rate (TR), and (b) absolute changes and (c) relative changes in TR caused by O₃ damage.

3.3 Changes in meteorology due to O₃–vegetation coupling

Through interacting with vegetation, O₃ has the potential to further affect the meteorological environment in China via modifying, e.g., surface heat fluxes, temperature, humidity, and boundary-layer height. The distribution of meteorological variables from simulations with and without O₃ damage is thus compared and analyzed in this section.

Figure 6 shows the spatial distribution of latent heat (LH) flux and sensible heat (SH) flux, and the changes in LH and SH due to O₃–vegetation coupling. With O₃ included in the model simulations, the LH flux decreases by more than 4 W m⁻² (Fig. 6b) on average following the decreases in transpiration rate. Hotspot areas are found in the Henan, Shandong, Anhui, and Jiangsu provinces, where reductions in LH can be up to 30 W m⁻². Meanwhile, 5–30 W m⁻² increases in SH flux are observed in central and northern China (Fig. 6d). With O₃–vegetation coupling, more than 20 % reductions in LH flux are found in central and northern China (Fig. 6c), 20 % increments in SH flux are found in similar regions (Fig. 6f), indicating that O₃ damage shifts the energy balance toward more net radiation being dissipated by SH flux than LH flux, with ramifications for surface temperature.

Figure 6 Spatial distribution of mean (a) latent heat flux (LH) and (d) sensitive heat flux (SH) from the simu_withoutO₃ experiment; absolute changes in (b) LH flux and (e) SH flux in JJA of 2014–2017 caused by O₃ damage; and relative changes in (c) LH flux and (f) SH flux caused by O₃ damage. Absolute changes are the LH (SH) flux from simu_withO₃ minus LH (SH) flux simu_withoutO₃. Relative changes are calculated by absolute changes over LH (SH) flux from simu_withoutO₃.

Figure 7 shows the distribution and the changes in surface relative humidity, temperature and planetary boundary-layer height (PBLH) in response to O₃ damage. Reductions in transpiration rate can directly cause reductions in relative humidity. As shown in Fig. 7b, relative humidity has at least 3 % absolute reductions. Values of relative humidity decrease more in northern China than in southern China. Similar to the changes in TR (Fig. 5b), larger reductions in relative humidity (3 %–9 %) are found over the Henan, Hebei, Shandong, and Anhui provinces. The decreases in LH flux and increases in SH flux following the changes in transpiration rate drive the increases in temperature and contribute to PBLH growth. As presented in Fig. 7e and h, the distribution and hotspot areas of the changes in temperature and PBLH are similar to those in relative humidity. Generally, northern China has larger increases of temperature and PBLH compared with other regions. Generally, temperature increases by 0.2–0.8 K and PBLH increases by 40–120 m for northern China. The hotspot areas experience at least 0.6 K increases in temperature and 80 m increases in PBLH.

Figure 7 Spatial distribution of mean (a) 2 m relative humidity, (d) 2 m temperature, and (g) planetary boundary-layer height (PBLH) in JJA of 2014–2017 from the simu_withoutO₃ experiment; absolute changes in (b) RH_2 m, (e) T_2 m, and (h) PBLH caused by O₃ damage; and relative changes in (c) RH_2 m, (f) T_2 m, and (i) PBLH caused by O₃ damage. Absolute changes are the results from simu_withO₃ minus results from simu_withoutO₃. Relative changes are calculated by absolute changes over the results from simu_withoutO₃.

As shown in Table S15, our results are comparable with results from a regional simulation conducted by Li et al. (2016), which showed that O₃ damage decreases LH flux by 10–27 W m⁻² and O₃ damage increases temperature by 0.6–2.0 ^∘C in the US. However, in their study, Li et al. (2016) assumed that O₃ damage to plants happens when O₃ concentration is over a threshold of 20 ppb to imitate a weaker detoxifying effect of plants, instead of the 40 ppb threshold that was commonly used in previous studies. Considering the severe O₃ air pollution in China, we resorted to use the more universal O₃ threshold used by previous studies (Lombardozzi et al., 2015; Sadiq et al., 2017; Zhou et al., 2018) to represent a more conventional detoxifying effect, instead of lowering the threshold value that would cause much larger changes in the surface fluxes and meteorological fields. Using a two-way coupling model and the same O₃ damage scheme, Arnold et al. (2018) revealed that O₃ causes less than 8 W m⁻² changes in surface heat fluxes regionally, which is smaller than the changes of surface heat fluxes in our study. One possible reason is that the simulated changes in O₃ and aerosol in Arnold et al. (2018) did not feed back onto radiation and climate simulation or affect LAI.

3.4 O₃–vegetation feedbacks on O₃ concentrations

O₃-induced changes in vegetation, surface fluxes, and the overlying meteorology can also constitute important feedback effects onto O₃ concentration itself. Figure 8 shows the spatial distribution of surface O₃ concentration. The change in surface O₃ concentration during daytime is also shown in Fig. S2. As shown in Fig. 8a (Fig. S2), surface O₃ concentration is higher in central and northern China during summer. In terms of the feedbacks on O₃ concentration, we found generally enhancements in O₃ concentration when O₃–vegetation interactions are accounted for, thus representing a positive feedback that worsens O₃ air quality (Fig. 8b). O₃ concentration increases the most (by up to 6 %) in the Hebei, Shanxi, and Henan provinces, with the maximum increment of 6 ppb. The enhancement in surface O₃ concentration from our study is at a similar magnitude to that from the study conducted by Sadiq et al. (2017), in which both biogeochemical and meteorological feedbacks from O₃–vegetation interactions to O₃ are considered. Without considering the meteorological feedbacks following the changes in transpiration to O₃ concentrations, smaller feedbacks on surface O₃ concentrations are found by the following studies. For instance, by incorporating O₃–LAI coupling in chemical transport model, Zhou et al. (2018) found an O₃ feedback of −1.8 to +3 ppb globally. Another similar work conducted by Gong et al. (2020) showed that O₃-induced inhibition in stomatal conductance increases surface O₃ by 2.1 ppb in eastern China, while considering the addition effects of O₃ on isoprene emission slightly reduces surface O₃ concentrations by influencing the precursors. Soil moisture deficit, which has been shown to reduce stomatal uptake, if considered, will also contribute to the enhancement in O₃ concentration (Rydsaa et al., 2016). Together with previous findings, it is increasingly clear that meteorological feedback could be an important pathway whereby O₃–vegetation interactions can further worsen O₃ air quality, almost doubling the effect of biogeochemical feedback alone (i.e., via changes in O₃-relevant chemical fluxes alone). It should be cautiously noted that in terms of magnitude alone the model biases in O₃ are comparable and sometimes larger than the up to 6 ppb systematic enhancement caused by O₃ damage, which represents be one major source of uncertainties in our study.

Figure 8 Same as Fig. 5 but for surface O₃ concentration.

Reduced dry deposition due to stomatal closure and reduced LAI, as well as increased isoprene emission, are all found to be the drivers for the overall positive O₃ feedback. Reductions in dry deposition velocity, following closely the corresponding reductions in transpiration rate as both processes are modulated by stomatal regulation, contribute in part to the O₃ enhancement. Figure 9 shows the spatial distribution of isoprene emission and its changes due to O₃ damage. We observe general increases in isoprene emission in eastern China, mainly due to increased surface temperature (Fig. 7e and f) that is more than enough to offset reduced isoprene caused by reduced LAI (Fig. 4e and f). All in all, O₃ damage on vegetation can further enhance O₃ levels via an overall positive effect, due to not only the associated reductions in dry deposition velocity, but also the reductions in transpiration, LH flux, and the resulting rise in surface temperature.

Figure 9 Same as Fig. 5 but for isoprene emission.

4 Conclusions

Tropospheric O₃ is one of the most concerning air pollutants due to its global warming effects and its ability to affect human health, vegetation, and crops. O₃ and vegetation closely interact with each other and such interactions may not only affect plant physiology (e.g., stomatal conductance and photosynthesis) but also influence the overlying meteorology and air quality through modifying leaf stomatal behavior, plant structure (e.g., LAI), and subsequently land–atmosphere fluxes. According to previous field experiments and modeling works, China has been recognized as one of the hotspot areas suffering from severe O₃ pollution and the resulting damage on vegetation and crops, but the feedback effects onto air quality and climate have not been fully characterized. Previous studies mainly focused on the global scale with coarse spatial resolutions, which did not fully capture the spatial distribution of O₃ damage on vegetation in China. Based on the results from global studies pointing out that China is a hotspot in terms of O₃ pollution and O₃ damage on vegetation, our model simulations performed at high spatial resolutions were capable of investigating O₃ damage effects on regional and provincial scales in China. In this study, we examined the effects of O₃–vegetation interactions on O₃ air quality and meteorology in China during 2014–2017 based on the two-way coupled WRF-Chem model simulations whereby O₃, meteorology, and vegetation physiology and structure can co-evolve with each other in real time.

We found that in China stomatal resistance is enhanced by up to 16 %, which is the direct response to O₃ damage. Northern China, especially the Henan, Hebei, and Shandong provinces, is identified as a hotspot area. For photosynthesis, more than 20 % reductions are observed in China. Large reductions (>2.4 µmol CO₂ m⁻² s⁻¹) are found in northeastern and southern China. Following reduced photosynthesis, LAI shows relatively small reductions (5 %–15 %), while GPP shows more than 20 % reductions (1.6 g C m⁻² d⁻¹). Changes in transpiration rate are due to both changes in stomatal resistance and changes in LAI. With the increases in stomatal resistance and decreases in LAI, transpiration deceases from 0.2 to 1.0 mm d⁻¹ in eastern China with the largest reductions occur in northern China. We also found that the distribution of changes in transpiration is consistent more with the distribution of stomatal resistance than with those of LAI, indicating the dominance of the former in contributing to the overall transpiration rate.

With O₃ damage, the LH fluxes decrease by more than 4 W m⁻² on average, with hotspot areas appearing in the Shandong, Anhui, and Jiangsu provinces, in which the decreases can be up to 30 W m⁻² following mostly the decreases in transpiration rate. SH fluxes increase in similar areas at comparable magnitudes (10–25 W m⁻²). The decreases in LH and the increases in SH cause the increases in temperature and PBLH. We found that northern China has larger decreases in relative humidity, temperature, and PBLH compared with other regions. Generally, relative humidity shows at least 4 % relative reductions, temperature increases by 0.2–0.8 K, and PBLH increases by 40–120 m for northern China. This indicates that O₃–vegetation interactions will cause a shift in the energy balance toward a state where available net radiation is dissipated more by SH flux than LH flux, with ramifications for surface temperature. This represents an additional pathway whereby anthropogenic O₃ pollution can worsen warming, in addition to O₃ being a greenhouse gas itself and O₃-induced plant damage diminishing the global net carbon sink (e.g., Sitch et al., 2007; Lombardozzi et al., 2015).

O₃ induces changes in vegetation, surface fluxes, and meteorology, and in turn affects its own concentration. In this study, we found that reduced dry deposition in China is mainly due to enhanced stomatal conductance, while enhanced isoprene emission is mainly due to enhanced surface temperature and the corresponding increase in O₃ concentration. O₃ concentration increases the most (up to 6 %) in the Hebei, Shanxi, and Henan provinces, with the maximum value of 6 ppb. Our results demonstrate that O₃–vegetation interactions can lead to strong positive feedback that can amplify O₃ pollution in China, in agreement with the suggestions by previous studies focusing on a global scale (Sadiq et al., 2017; Zhou et al., 2018; Gong et al., 2020). We also found that fully considering the positive O₃–vegetation feedbacks, especially when meteorological changes are also accounted for, generates greater damage on vegetation productivity than found by studies that only considered “offline” O₃ damage on plants without feedbacks (Yue and Unger, 2014; Lombardozzi et al., 2015).

In this study, the summertime simulation period of JJA was selected due to the high O₃ pollution in this season and the overlapping-with-vegetation growing season to capture the severe O₃ damage on vegetation. Nevertheless, uncertainty may still arise from that our simulation period may not cover the growing season of all vegetation types and may not cover all periods that O₃ damage happens, which may represent an underestimation of the full scale of O₃ damage. Future work should be conducted for longer time periods and for all seasons, which will help us better understand O₃–vegetation interactions in China. Uncertainty may also arise from the O₃ scheme employed in this study in terms of the CUO calculation and the consideration of O₃ detoxification mechanism of different vegetation types. The calculation of CUO heavily relies on the O₃ threshold. Considering the sensitivities of different vegetation types to O₃ damage, CUO threshold should be varied with different vegetation types. However, a constant O₃ threshold was employed in our study for the whole simulation domain and for all vegetation types, which may either underestimate or overestimate the actual O₃ damage. Moreover, following the work of Lombardozzi et al. (2015), we classified all the vegetation types into only three groups, which may be too coarse to investigate O₃ damage effects on regional or local scales. For example, Zhou et al. (2018) pointed out that Lombardozzi et al. (2015) treated tropical and temperate plants equivalently, which might lead to possible biases. More studies should be conducted to derive more appropriate O₃ thresholds for CUO calculation and make them available for regional scales or for different vegetation types. Another source of uncertainty may arise from the lack of representation of the direct effect of O₃ on isoprene emission. As pointed out by Gong et al. (2020), including the effect of O₃ damage on isoprene emission may reduce O₃ concentration by influencing precursors, but increase O₃ concentration at the same time through weakening the shortwave radiative forcing of secondary organic aerosols, which would help constitute a more complete feedback mechanism between O₃ and vegetation. Moreover, uncertainties may also come from that the effect of soil moisture deficit was not considered in this study, which may underestimate the reduction in dry deposition sink of O₃. It should also be noted that keeping the anthropogenic emission inventory fixed in 2014 levels may be another limitation because of the nonlinear chemistry involving biogenic and anthropogenic precursors. Despite these uncertainties and limitations, our study provides detailed and comprehensive results whereby O₃–vegetation impacts will adversely affect plant growth and crop production, contribute to global warming, worsen the severe O₃ air pollution in China via feedbacks, and identifies the hotspot areas in the country. Our findings clearly pinpoint the need to consider the O₃ damage effects in both air quality studies and climate change studies.