Co-benefit of forestation on ozone air quality and carbon storage in South China

Introduction

Green vegetation plays a pivotal role in the interaction between the terrestrial biosphere and atmosphere, affecting land-atmosphere energy and water exchanges and modulating the carbon and nitrogen cycle^1,2,3. It also serves as a major source of nonmethane volatile organic compounds (NMVOC) and depositing surfaces to take up air pollutants^4,5,6. Their density is measured by the leaf area index (LAI), defined as the one-sided green leaf area per unit ground area in broadleaf species and half of the total needle surface area per unit ground area in coniferous species^7,8. Recent satellite observations have revealed increasing trends in LAI over the globe, reflecting a generally greening Earth since 2000, despite ongoing deforestation in regions such as South America and Africa^9,10. China accounted for 25% of the global net increase in LAI during the period, mainly driven by forestation^10,11.

Forestation, including afforestation (planting new trees) and reforestation (replanting trees), is an important measure to enhance carbon storage in the biosphere and help achieve carbon neutrality^12,13,14. The recent Intergovernmental Panel on Climate Change (IPCC) report suggests that forestation should remove approximately one-quarter of atmospheric carbon to help limit global warming to 1.5 degrees^12,13. Sustainable forest management also helps achieve the United Nations Sustainable Development Goals (SDG) of food and clean water supply (SDG2 and SDG6) and ecosystem protection (SDG 15)^12,14,15. China has implemented multiple forestation programs since 1978¹⁶, increasing forest stock volume (FSV) from 8.7 to 13.7 billion m³ over 1975–2005¹⁷ (Supplementary Table 1). Afforestation is estimated to contribute over half of the increases in China’s total land carbon storage during the 2000s^18,19. As a necessary solution to protect ecosystems and to meet carbon neutrality, specific forestry development goals have been set in China²⁰, i.e., to expand the FSV by 6 billion m³ in 2030 relative to 2005 (13.7 billion m³) and further increase it to 26.5 billion m³ by 2050 (“Methods”).

Apart from the challenge of carbon neutrality, China still faces the air quality challenge, and greening can profoundly affect surface ozone air quality. Surface ozone is a major air pollutant in China and has continued to grow till recent years^21,22,23,24. Ozone is produced by the photochemical oxidation of carbon monoxide (CO) and NMVOC in the presence of nitrogen oxides (NO_x = NO + NO₂). Previous studies have demonstrated the complexity of biosphere-atmosphere interactions by accounting for counterfactual scenarios assuming changes in the land type and LAI, e.g., greening can change biogenic volatile organic compound (BVOC, biogenic NMVOC) emissions and hence ozone formation^25,26,27. Changes in the land type and LAI can also alter ozone dry deposition velocity, thus affecting ozone pollution and ozone-induced vegetation damage^28,29,30. A recent study found that forests can affect surface ozone concentrations by changing the vertical diffusivities and turbulent exchanges within and below the canopies^31,32. However, this canopy effect is rarely considered in current air quality models, and the net effects of historic actual greening in South China on ozone pollution remain to be determined.

In this work, we quantify the greening effects on carbon storage and ozone air pollution in South China over the past 15 years (2005–2019) and future projections under different shared socioeconomic pathways. We use a regional atmospheric chemistry model (WRF-Chem), driven by the satellite-measured LAI and land cover data, coupled with a forest canopy turbulence parameterization and an improved dry deposition scheme (Fig. 1). The improved model results are evaluated using ozone gradient measurements at a South China forest station. We demonstrate that future forestation-induced greening can benefit both carbon neutrality and ozone mitigation.

Results and discussion

Forestation-induced greening and its impacts on ozone in South China

Figure 2 shows the satellite observations of the growing season (April, July, and October) mean LAI and dense forest cover fraction over China since 2003. Significant increases in LAI have been observed in the past 15 years in South China, with 42% of the net increases attributed to forests and 32% to croplands¹⁰. The regional mean LAI increased from 1.14 m² per m² in 2005 to 1.30 m² per m² in 2019 over China, while it increased from 2.38 m² per m² to 2.74 m² per m² over South China. Increases in annual mean LAI have also been reported in previous studies^9,10,33. Here, we find that the largest LAI increases occur in April, followed by October, which could be caused by an extension of the growing season in response to climate change over this period^34,35. LAI increases in July are small even though LAI peaks in the month (Supplementary Fig. 1). The satellite-observed land fractions of forest cover also increased steadily over this period (Fig. 2). The national mean dense forest cover fraction has increased from 9.8% in 2005 to 11.5% in 2019. The mean dense forest cover fraction in South China has increased from 18.6% in 2005 to 23.9% in 2019, and the increases in dense forests are associated with the decreases in savannas (Fig. 2 and Supplementary Fig. 2). The National Forest Inventory (NFI), which defines forests with the tree cover > 20%, shows a similar trend of 13% increases during this period (20.36% over 2004–2008 vs 22.96% over 2014–2018 in China)¹⁷.

Fig. 2: Greening due to forestation in South China as evidenced by MODIS satellite observations of leaf area index (LAI) and dense forest (tree cover > 60% and canopy height > 2 m) cover fraction.

We apply the Weather Research and Forecasting model coupled with the Chemistry (WRF-Chem) regional air quality model to assess the impacts of greening on ozone air quality over southern China. A number of model improvements are implemented to better quantify the impacts of forestation by using satellite observations of LAI and major land cover fractions (Fig. 1; “Methods”). Recent newly conducted near-surface ozone measurements at the Dinghushan (DHS) Forest Ecosystem Research Station (a forest station in South China; Methods) reveal large ozone vertical gradients from above and below the forest canopy, and we find that capturing the near-surface ozone gradients requires coupling a forest canopy turbulence scheme and improving the dry deposition scheme to the regional air quality model (“Methods”; Supplementary Fig. 3). A series of simulations are conducted using the improved model for the growing season of 2019 (Table 1; “Methods”). As shown in Fig. 3, the baseline simulation well captures the nationwide measured surface ozone concentrations in the metric of maximum daily 8 h average (MDA8) with correlation coefficients of 0.7 and mean biases of within 12% in China.

Table 1 WRF-Chem model simulations conducted in this study

Full size table

**Fig. 3: 2005-2019 greening effects on ozone air quality over South China in the growing season.**

Figure 3 shows that the 2005–2019 greening induces net decreases in surface MDA8 ozone by 1.4 ± 2.3 ppbv (mean ± standard deviation) over South China during the growing season. Large MDA8 ozone reductions of even over 8.0 ppbv occur over areas where the land cover switched from savanna land to forest (Supplementary Fig. 2). Ozone reductions are not confined in forested regions but are widely spread over South China. In addition to the human health metric MDA8, we find much greater effects on vegetation damage based on the ozone exposure metrics (AOT40, W126, and O3RH), gross primary productivity (GPP) reductions, and phytotoxic ozone dose (POD₁; “Methods”)^36,37. The monthly AOT40 and W126 values show high ozone exposure in North China, where ozone concentrations are particularly high, and considerable ozone exposure in South China (AOT40 of 2728 ± 2143 ppbv hour, W126 of 3098 ± 3022 ppbv hour) compared with other regions worldwide^21,37. The metrics of O3RH, GPP reductions, and POD₁ show higher ozone exposure and vegetation damage in South China than in North China, as the low RH conditions in North China limit the stomatal uptake of ozone and the high vegetation abundance in South China³⁷.

The 2005-2019 greening can significantly reduce the exposure of vegetation to ambient ozone. We find that the average AOT40, W126, and O3RH in forests over South China decreased by 39%, 41%, and 15% during the growing season, respectively (Fig. 3 and Supplementary Fig. 4). However, the response of ozone-induced vegetation damage to greening effects differs depending on the metric chosen. The average GPP reductions related to ozone exposure in forests over South China decreased by 24% after considering LAI-induced GPP increases of 5% in 2019 compared to 2005, while the dose-based metric of POD₁ increased by 10% due to increased stomatal uptake (Supplementary Figs. 4–5). This reverse response between the concentration metrics (AOT40, W125, O3RH, and GPP reductions) and the flux metric (POD₁) is also measured by Ronan et al.³⁸.

Greening impacts via biosphere-atmosphere interactions

The impacts of forestation-induced greening on ozone air quality are complicated by the multiple potentially counteracting pathways of biosphere-atmosphere interactions. Here, as shown in Fig. 4, we diagnose the suppression of turbulence transfer below the forest canopy as the “canopy” effect, changes in dry deposition of ozone associated with LAI as the “deposition” effect, and changes in BVOC emissions due to greening as the “emission” effect. These effects can be separately quantified by model sensitivity simulations using fixed land cover data (i.e., using the 2005 level), fixed LAI in the model dry deposition scheme, and fixed LAI in the biogenic emission scheme (Table 1; Methods).

**Fig. 4: Biosphere-atmosphere interaction pathways of greening effects on ozone air pollution.**

Figure 4 summarizes the statistics of the 2005–2019 ozone changes in South China due to greening via different effects averaged over the growing season. Supplementary Fig. 6 shows their spatial and temporal variations. For the regional means, the 2005–2019 greening has reduced the mean surface MDA8 ozone concentrations in South China by 0.7 ppbv (− 0.3–4.5 ppbv; 25–75% percentiles) via the “canopy” effect and by 1.0 ppbv (0.4–1.5 ppbv) via the “deposition” effect, while the “emission” effect could have slightly increased the mean MDA8 ozone by 0.3 ppbv (0.2–0.5 ppbv). The “canopy” and “deposition” effects exceed the “emission” effect, leading to overall net decreases in surface MDA8 ozone levels as induced by the 2005–2019 greening in South China.

Substantial differences in spatial and seasonal patterns are found for the three effects (Supplementary Fig. 6). The “canopy” effects tend to be localized, lowering MDA8 ozone values by up to 10 ppbv in forestation areas, but raising them by up to 10 ppbv in deforestation areas and by ~ 1.3 ppbv in urbanization areas (Supplementary Fig. 2). The ozone changes due to “deposition” and “emission” effects in South China are less localized compared to the “canopy” effect and heavily rely on LAI changes. Different from the “canopy” effect, which has almost no seasonal changes, the “deposition” and “emission” effects show large seasonal variations as driven by seasonal changes in LAI (Supplementary Fig. 6). Greening, as represented by LAI increases, would provide larger depositing surfaces to weaken resistances of stomatal and cuticular removal and thus enhance ozone dry deposition velocities over South China by ~ 2% in July and up to 9% in April (Supplementary Fig. 7). Meanwhile, the 2005–2019 LAI increases further enhanced emissions of BVOC by 4% in July and 20% in April in South China. The resulting monthly mean surface MDA8 ozone changes in South China associated with the 2005–2019 greening are − 0.3 ~ − 1.7 ppbv due to the “deposition” effect and 0.1–0.4 ppbv due to the “emission” effect. A previous study reported up to 0.5 ppbv reduction in surface ozone from 1992 to 2014 due to dry deposition enhancements and small isoprene emission increases over southern China³⁹, consistent with our results.

Our analyses further emphasize the important role of the “canopy” effect. The forest can suppress turbulence and vertical transport below and near its canopy, thus reducing ozone from the free troposphere mixing down to the surface. Such a turbulence suppression effect can be captured by our improved regional air quality model when coupled with the forest canopy turbulence scheme. In the presence of the forest canopy, the regional mean near-surface eddy diffusivity coefficients over South China decreased by 31% in April 2005 and by 41.7% in April 2019 compared with the conditions in the absence of the forest canopy effect (Supplementary Fig. 8). The presence of the canopy turbulence effect leads to a sharp gradient of ozone near the surface, as recorded by the ozone gradient measurements at DHS and simulated by the model, and reduces the mean surface MDA8 ozone concentrations by 0-12.6 ppbv over South China in April 2019 (Supplementary Figs. 3 and 8). Meanwhile, it increases ozone concentrations above the canopy and redistributes the vertical ozone profile in the lower troposphere. The 2005–2019 forestation in South China has led to about 35% stronger turbulence suppression by the forest canopy (Supplementary Fig. 8), which largely drove the “canopy” ozone effect shown in Fig. 4. Forest canopy expansion would also induce a modest increase in ozone dry deposition velocity via increases in forestland cover, while we find that such an effect on ozone pollution is much smaller than the suppression of turbulence (Supplementary Figs. 6, 7).

Co-benefit of future greening on carbon storage and ozone air quality

As a potential solution to facilitate carbon neutrality, future forestation could further green South China and may impact both ozone air quality and carbon stocks. Previous forestation goals have been achieved in China’s five-year plans. China has now announced an increase in FSV to 19.7 billion m³ by 2030 and to 26.5 billion m³ by 2050, based on which we can estimate future increases in forest LAI and carbon stocks (Methods). We find a strong linear relationship between the FSV and the forest LAI using reported historical statistics (Supplementary Fig. 9). Such a relationship can be applied to estimate future LAI given the projected FSV. We estimate that the annual mean forest LAI will increase by 8.1% in 2030 and by 34% in 2050 relative to the 2019 condition in South China (Methods; Fig. 5). The largest increase in LAI is projected to occur in the evergreen broadleaf forests of southwestern China. The estimated maximum LAI is 7.7 in 2030 and 14.9 in 2050, within the range of historically measured high values^40,41.

**Fig. 5: Co-benefits of future forestation on carbon stock and ozone air quality.**

Forest carbon stock (also called carbon storage) is a commonly used metric to represent carbon stored in natural forests. It can be calculated by the product of FSV and the biomass expansion factor (BEF, defined as the ratio of forest biomass to timber volume)^42,43. Based on previous estimates of Chinese forest biomass carbon stocks (Supplementary Table 2), we derive a national mean BEF of 0.96 tons per m³ with a range of 0.85–1.07 tons per m³ (Methods). The forest biomass carbon stocks are thus estimated to be 6.6 Pg C (5.8 to 7.3 Pg C) in the year 2005, 8.4 Pg C (7.4 to 9.4 Pg C) in 2019, increasing to 9.4 Pg C (8.3 to 10.5 Pg C) in 2030, and further to 12.7 Pg C (11.2 to 14.2 Pg C) in 2050. Our estimate of future forest biomass carbon stocks is in line with previous data-driven studies (8.0 to 10.8 Pg C in 2030 and 9.0 to 13.1 Pg C in 2050)^44,45,46.

In addition to carbon uptake, future greening would have a co-benefit on surface ozone air quality as driven by biosphere-atmosphere interactions. To assess such impacts on future ozone air quality, we conduct model sensitivity simulations with projected LAI over China in 2030 and 2050 (Table 1). Future anthropogenic emissions are projected in different shared socioeconomic pathways (SSPs) reported by the IPCC. We use the future (2030) anthropogenic emissions by a low-emission scenario (SSP126) and a high-emission scenario (SSP370) in 2030 (Methods). In the SSP126 scenario, Chinese NO_x emissions would decrease by 31%, and NMVOC emissions would decrease by 42%, while in SSP370, Chinese NO_x and NMVOC emissions would increase by 38% and 24%, respectively, during the period from 2019 to 2030. Future (2030) anthropogenic emissions changes alone would decrease the mean surface MDA8 ozone concentration by 4.3 ppbv under the SSP126 scenario, while they would increase it by 3.7 ppbv under the SSP370 scenario in South China (Fig. 5 and Supplementary Fig. 10)

We find that, for all three anthropogenic emission scenarios (the 2019 condition, SSP126, and SSP370), the 2019–2030 forest greening would lead to decreases in surface ozone levels. As shown in Fig. 5, accounting for the 2019–2030 forest greening would decrease the mean MDA8 ozone in South China by 0.36 ppbv (0.01–0.71 ppbv) under the 2019 emission condition, 0.34 ppbv (0.01–0.67 ppbv) in SSP126, and 0.40 ppbv (0.01–0.79 ppbv) in SSP370. The greening effects are slightly higher in SSP370 due to the high NO_x emission levels in the scenario. 2019–2030 anthropogenic emission increases in the SSP370 scenario would increase the mean ozone by 3.7 ppbv in South China; however, if we consider the impacts of future greening, such increases would be reduced by 11%. As the forest LAI increases extending to 2050, the mean surface MDA8 ozone concentrations would be reduced by 1.4 ppbv (0.2–2.6 ppbv) over South China compared to the 2019 level (Fig. 5). Furthermore, future forest greening would decrease ozone-induced GPP damage by 0.5% ± 2.9% in 2030 and 4.9% ± 6.7% in 2050 though vegetation abundance increases (Supplementary Fig. 5). These ozone benefits occur mainly in South China, where we project the most forest LAI enhancements (Supplementary Figs.9, 10). Such benefits may only present a conservative estimate, as here, the “canopy” effect that could further reduce near-surface ozone is not accounted for by only accounting for future LAI increases. Future forestland cover increases, which are highly uncertain, can further help ozone air pollution mitigation.

Forestation thus offers a win-win strategy for ozone air pollution abatement and carbon neutrality. The complex biosphere-atmosphere interactions associated with forestation, as quantified in this study, can also be important in many other regions over the globe, such as Brazil, India, Mexico, and tropical African countries, where significant forestation has occurred in recent decades (Supplementary Fig. 11). These regions are also facing increasing tropospheric ozone levels as seen by satellite observations⁴⁷, and have committed nearly 300 Mha of degraded land aiming to restore 350 Mha of the forest by 2030 (www.bonnchallenge.org)¹³. In addition to the benefits for carbon storage and ozone air quality, forestation can also contribute to ecosystem services, biodiversity protection⁴⁸, and water management⁴⁹, requiring integral consideration in environmental strategies in future work.

Our analyses include some uncertainties. First, a number of global LAI products based on different satellite instruments and inversion algorithms provide a wide range of leaf area estimates. Most of these LAI products recorded an increasing trend over South China in recent decades, and the MODIS estimate of a 15% increase in LAI over 2005–2019 is in the middle of the various products^10,33,50 (Supplementary Information). Second, land cover changes, such as forest expansion, would affect local surface albedo, energy balance, and meteorological conditions^51,52. Forestation of dryland regions was recently found to have a large warming potential due to the reduced albedo effect⁵³. Forestation in South China is mainly converted from savannas. This may cause small decreases in albedo while increasing latent heat near the surface, leading to slight decreases in surface temperature, and thus the effect on ozone formation can be small (Supplementary Fig. 12). Our model simulations by nudging the meteorological conditions every two days can account for the diurnal meteorological responses due to the land cover changes while excluding long-term and large-scale atmospheric circulation changes. In a model simulation without meteorological nudging, we find that the surface ozone changes due to the 2005–2019 forestation greening would be ~ 6% larger (Supplementary Fig. 13). Furthermore, recent studies have demonstrated strong positive ozone-vegetation feedbacks where ozone-induced vegetation damage and inhibition of stomatal conductance would enhance ozone air pollution^26,54,55. Addressing such feedback requires online coupling of an air quality model with a land vegetation model, which is not accounted for in the study. We expect that surface ozone benefits due to forestation would be enlarged when considering reduced ozone levels to further lower vegetation damage.

Methods

Satellite observations of leaf area index and land cover

Land type and vegetation information are from Moderate Resolution Imaging Spectroradiometer (MODIS) sensors on NASA’s Terra and Aqua satellites. The 8-day composite leaf area index (LAI) product (MCD15A2H Version 6) with a 500 m spatial resolution, the monthly normalized difference vegetation index (NDVI) product (MOD13C2 Version 6.1) with a 0.05° resolution, and the land cover type product (MCD12Q1 Version 6) with a 500 m resolution are used in this study^56,57,58. We filter out LAI and NDVI data under cloudy conditions and replace invalid values with valid ones nearest in time to improve data coverage^59,60. The LAI dataset is averaged monthly, and then a 24-month running average is applied to identify the long-term trend. The LAI and NDVI datasets are then spatially aggregated to a 27 km resolution for model inputs. The monthly mean LAI and NDVI datasets from 2003 to 2020 are calculated and utilized in this study.

The MCD12Q1 dataset contains 17 land cover categories following the International Geosphere and Biosphere Program (IGBP)⁶¹. We reproject the tiled MCD12Q1 dataset from 2003 to 2019 and spatially aggregate it to generate yearly land cover model inputs that include the major land cover types and fractions for the 17 categories at 27 km resolution. The IGBP classification types are further grouped into four plant functional types (broadleaf tree, needleleaf tree, shrub/brush, and herb/crops/grasses) for the calculation of biogenic emissions. We define forest grid cells when the total forestland cover fractions (including evergreen needleleaf, evergreen broadleaf, deciduous needleleaf, deciduous broadleaf, and mixed forest as categorized in MCD12Q1) are higher than 20%. It should be noted that forestland cover types in MCD12Q1 IGBP represent dense forests with tree cover > 60% and canopy height > 2 m. This differs from NFI, which defines forests with tree cover > 20%¹⁷.

Observations of air pollutants and meteorology

Ozone concentration gradients were measured at a flux tower in the Dinghushan (DHS) Forest Ecosystem Research Station (a forest station in South China; 112°32′28′′E, 23°10′25′′N) from 21 August to 6 September 2019⁶². The flux tower is located in the central DHS and is surrounded by coniferous and broad-leaved mixed woodland. The ozone concentrations are concurrently measured on the tower at 7 m and 30 m above the surface using two closed-path eddy covariance ozone flux measurement systems. We calculate ozone concentration differences at the two levels as the 7 m–30 m ozone gradients. More details on the DHS and ozone measurements are given in Liu et al.⁶³. Ozone measurements at five background stations in China have accessed from the Tropospheric Ozone Assessment Report (TOAR) dataset (https://doi.org/10.34730/19cdb0ca3ad046bd88449f38b426542b).

Nationwide hourly measurements of surface ozone concentrations at over 300 cities in China for 2019 are obtained from the China Ministry of Ecology and Environment national network (MEE, http://106.37.208.233:2035/). Observations of the tropospheric formaldehyde (HCHO) column in 2019 retrieved from the TROPOspheric Monitoring Instrument (TROPOMI) onboard the Copernicus Sentinel-5 Precursor (S5P) satellite are used to evaluate BVOC emissions in the model. BVOC emissions are the dominant source of HCHO in South China^64,65 and have proven to be strongly correlated with the satellite HCHO vertical column in the region during the summer months^65,66. We filter out the HCHO column under cloudy conditions and only use the valid data with a QA > 0.5^67,68.

In addition, meteorological observations, including 10-m wind speed (WS10), 2-m air temperature (T2), and 2 m dew-point temperature (Td2) at 366 stations over China in the growing season of 2019 are collected from the National Climatic Data Center (NCDC, https://ncdc.noaa.gov/isd/data-access). Dew-point temperature observations are converted to relative humidity (RH) according to the Tetens empirical formula. Evaluations of model-simulated meteorology with these observations generally show good agreement (Supplementary Table 3).

These available observations are used to assess the performance of the air quality model, employing various statistical metrics such as mean values (OBS and MOD), mean normalized bias (MNB), normalized mean bias (NMB), correlation coefficient (R), and the index of agreement for validation (IOA).

$${{rm{MNB}}}=frac{1}{N}mathop{sum }limits_{i=1}^{N}left(frac{{Mod}left(iright)-{Obs}(i)}{{Obs}(i)}right)$$

(1)

$${{rm{NMB}}}=frac{{sum }_{i=1}^{N}({Mod}left(iright)-{Obs}(i))}{{sum }_{i=1}^{N}({Obs}(i))}$$

(2)

$${{rm{R}}}=frac{{sum }_{i=1}^{N}({Mod}left(iright)-{MOD})times {sum }_{i=1}^{N}({Obs}left(iright)-{OBS})}{sqrt{{sum }_{i=1}^{N}{({Mod}left(iright)-{MOD})}^{2}}times sqrt{{sum }_{i=1}^{N}{({Obs}left(iright)-{OBS})}^{2}}}$$

(3)

$${{rm{IOA}}}=1-frac{{sum }_{i=1}^{N}{({Mod}left(iright)-{Obs}(i))}^{2}}{{sum }_{i=1}^{N}{(left|{Mod}left(iright)-{OBS}right|+left|{Obs}left(iright)-{OBS}right|)}^{2}}$$

(4)

Ozone metrics relevant to air quality and vegetation damage

We calculate the maximum daily average 8 h (MDA8) ozone, three exposure-based vegetation ozone metrics (AOT40, W126, and O3RH), and one dose-based vegetation ozone metric of phytotoxic ozone dose (POD₁). These metrics are commonly used to quantify ozone air pollution and its impacts on human health and vegetation^36,37. The MDA8 ozone is the metric used by the WHO air quality guideline and the air quality standard in China. The mean MDA8 ozone values (in the unit of ppbv) in the growing season are used to estimate long-term ozone exposure in this study. The monthly AOT40 (accumulated hourly ozone concentrations over a threshold of 40 ppbv during 08:00 to 19:59; in the unit of ppbv hour) and the monthly W126 (the sigmoidal weighted sum of all hourly ozone concentrations during 08:00 to 19:59; in the unit of ppbv hour) metrics consider the impacts of cumulative ozone concentrations on vegetation. The metric of O3RH (in the unit of ppbv hours) further accounts for the influence of ambient RH on ozone stomatal uptake and is estimated as MDA8 over 20 ppbv multiplied by RH between 40–-80%³⁷. The metric of POD₁ (in the unit of nmol m⁻² s⁻¹) further considers vegetation damage as ozone uptake fluxes entering stomata rather than ozone exposure and is estimated as accumulated hourly stomatal ozone flux (F_s,O3) that exceeds a threshold of 1 nmol O₃ m⁻² s⁻¹ during 08:00 to 19:59³⁸. We calculate monthly AOT40, monthly W126, daily O3RH, and monthly POD₁ using the equations below.

$${{rm{AOT40}}}=mathop{sum }limits_{{{rm{d}}}=1}^{{{rm{days}}}}mathop{sum }limits_{{{rm{h}}}=8}^{19}({max }({[{O}_{3}]}_{{{rm{d}}},{{rm{h}}}}-40,0))$$

(5)

$${{rm{W}}}126=mathop{sum }limits_{{{rm{d}}}=1}^{days}mathop{sum }limits_{{{rm{h}}}=8}^{19}left({[{O}_{3}]}_{{{rm{d}}},{{rm{h}}}}times left(frac{1}{1+4403times exp (-0.126times {[{O}_{3}]}_{{{rm{d}}},{{rm{h}}}})}right)right)$$

(6)

$${{rm{O3RH}}}=, {max }({{rm{MDA}}}8{[{O}_{3}]}_{d}-20,0)times , {max }({min }({{rm{R}}}{{rm{H}}}_{d}-40 %,40 % ),0)$$

(7)

$${{{rm{POD}}}}_{1}=mathop{sum }limits_{{{rm{d}}}=1}^{{{rm{days}}}}mathop{sum }limits_{{{rm{h}}}=8}^{19}{,max }({[{{{rm{F}}}}_{{{{rm{s}}},{{rm{O}}}}_{3}}]}_{{{rm{d}}},{{rm{h}}}}-1,0)times 3600$$

(8)

WRF-Chem model configuration

The Weather Research and Forecasting model coupled with Chemistry (WRF-Chem) version 3.6.1 is employed to simulate ozone air pollution. The model domain covers China and its adjacent areas with a 27 km spatial resolution. It has 38 layers vertically from the surface to 50 hPa, with the first and second layers centered at 10 m and 30 m above the ground, respectively. The initial and lateral meteorological boundary conditions use the National Center for Environmental Prediction Final (FNL) Analysis data with 1° spatial resolution and 6-hour temporal resolution. To avoid drifting in simulated meteorology, we reinitiate meteorological fields with the FNL Analysis every two days. The chemical initial and boundary conditions are provided by the CAM-Chem simulation output⁶⁹. We use the gas-phase Carbon-Bond Mechanism Z mechanism⁷⁰ to simulate tropospheric ozone chemistry. Dry deposition and vertical mixing processes are calculated by Wesely’s dry deposition scheme⁷¹, the Unified Noah land-surface model⁷², and the Yonsei University planetary boundary layer scheme⁷³. Other physical and chemical parameterization schemes follow Liu et al.⁷⁴.

For emissions in the model, the baseline simulation uses the 2017 Multi-resolution Emission Inventory for China (MEIC) and the 2010 MIX inventory for regions outside China (www.meicmodel.org/), except for Chinese agricultural ammonia (NH₃) emissions from Zhang et al.⁷⁵ with updated statistics for 2017. Biomass burning emissions use the Fire Inventory from the NCAR⁷⁶. Biogenic emissions are calculated online using the Model of Emissions of Gases and Aerosols from Nature⁷⁷ (MEGAN) with improved model inputs of LAI and land type described below. The projections of future anthropogenic emissions are based on two shared socioeconomic pathways (SSPs; SSP1-sustainability and SSP3-regional rivalry) with different 2100 forcing levels (2.6 W·m⁻² and 7.0 W·m⁻²) accessed from Coupled Model Intercomparison Project Phase 6 (CMIP6; https://esgf-node.llnl.gov/projects/input4mips/). The 2015–2030 emission trends in a low-emission scenario (SSP126) and a high-emission scenario (SSP370) are then applied to the 2015 MEIC emissions to produce two sets of hypothetical 2030 Chinese anthropogenic emissions.

Model improvements for the simulation of land-atmospheric interactions

Coupling a forest canopy turbulence parameterization

Forest would suppress vertical turbulence diffusivity within and below its canopy; however, this effect is not accounted for in default WRF-Chem. Here we implement a forest canopy parameterization in the boundary layer mixing scheme of the model. The parameterization modifies the turbulence diffusivity and then the vertical mixing of atmospheric chemical constituents (all gaseous and aerosol species) within the forest canopy, following a localized-near-field approximation proposed by Raupach⁷⁸. It has been recently applied by Makar et al.³², which demonstrated the large effects of the forest canopy on turbulence mixing and surface ozone levels over North America. This Eulerian-based modeling approach with first-order closure K-theory is suitable for WRF-Chem since it is more numerically efficient than the higher-order closure approach and the Lagrangian-based stochastic modeling approach^79,80. Following Makar et al.³², the parameterization can be expressed by the formulas below to estimate the suppression of vertical diffusivity (K) within and below the canopy:

$${{rm{{K}}}}_{{mathrm{mod}}}({{rm{z}}})=frac{{{{rm{K}}}}_{{mathrm{mod}}}({{rm{{z}}}}_{{{rm{base}}}})}{{{{rm{K}}}}_{{rm{{est}}}}({{{rm{z}}}}_{{rm{{base}}}})}times {{{rm{K}}}}_{{rm{{est}}}}({{rm{z}}})$$

(9)

$${{mathrm{{K}}}}_{{rm{{est}}}}({{rm{z}}})={sigma }_{{{rm{w}}}}^{2}({{rm{z}}})times {{{rm{T}}}}_{{rm{{L}}}}({{rm{z}}})$$

(10)

$${sigma }_{{rm{{w}}}}^{2}({{rm{z}}})={{{rm{u}}}}^{ * }times left[{{rm{A}}}+{{rm{B}}}times {f}left(frac{{{rm{z}}}}{{{rm{Hc}}}}right)right]$$

(11)

$${{mathrm{{T}}}}_{{rm{{L}}}}({{rm{z}}})=frac{{{rm{Hc}}}}{{{{rm{u}}}}^{ * }}times left[0.256times left(frac{{{rm{z}}}-0.75times {{rm{Hc}}}}{{{rm{Hc}}}}right)+0.492times exp left(frac{-0.256times {{rm{z}}}times {{rm{C}}}}{0.492times {{rm{Hc}}}}right)right]$$

(12)

Here, K_mod is the modified vertical diffusivity coefficient at height z above the ground below the canopy. z_base is the height of the first model layer above the forest canopy. The vertical diffusivity K_mod(z_base) is not affected by the forest canopy and is calculated by the model boundary layer mixing scheme. We reduce the diffusivity K_mod(z) based on the vertical profile of diffusivity (K_est(z)) determined by the near-field theory. Following refs. ^32,78, K_est profiles are estimated by Eq. (10), where ({{{rm{sigma }}}}_{{{rm{w}}}}^{2}left({{rm{z}}}right)) is the Eulerian vertical velocity variance (Eq. 11) and ({{{rm{T}}}}_{{{rm{L}}}}left({{rm{z}}}right)) is the Lagrangian timescale (Eq. 12), u* is the friction velocity, and Hc is the forest canopy height. The parameters A and B and the function of z/Hc in Eq. (11) are adopted from Makar et al.³², which depend on the stability of the atmosphere as diagnosed in the model by Hc/L with L being the Monin-Obukhov length. The term C in Eq. (12) represents 1-Hc/L for unstable atmospheres and a constant (0.1) for neutral and stable atmospheric conditions. The model diffusivity coefficient below the canopy reduced to 0.4 m² s⁻¹ from 4.3 m² s⁻¹ at the DHS site in August 2019 after coupling the forest canopy turbulence parameterization. This value is slightly lower than the approximately 0.6 m² s^-1 estimated with the concurrent below-canopy ozone eddy covariance flux measurements at the DHS site⁶³. This discrepancy is likely due to the absence of ozone uptake by the upper canopy, leading it to compensate by overly enhancing the suppression of turbulence.

The forest canopy height Hc is computed from the satellite dataset of NDVI by using the equation⁸¹:

$${{rm{Hc}}}={{{rm{Hc}}}_{min}}+frac{{{rm{H}}}{{rm{{c}}}_{max} }-{{{rm{Hc}}}_{min} }}{{{{rm{NDVI}}}_{max} }-{{{rm{NDVI}}}_{min} }} times left({{rm{NDVI}}}-{{rm{NDV}}}{{{rm{I}}}}_{{{min }}}right)$$

(13)

where Hc_max and Hc_min are the prescribed maximum and minimum canopy heights and are set to 35 m and 2 m, respectively, for the entire model domain. NDVI_max and NDVI_min are the maximum and minimum forest NDVI in the dataset, respectively. The spatial distribution of canopy height is determined by the satellite dataset of NDVI. The forest canopy parameterization is applied to the forest grid cells as defined by the MODIS land cover dataset with the forest fraction exceeding 20%. Due to limitations in model vertical resolution, the parameterization is applied to the lowermost two layers in most forest grid cells. We also tested the forest fraction of 15% and 25% and the resulting effect of the forest canopy parameterization is similar. As indicated by the equations above, both the increases in forest coverage fractions and tree height can suppress vertical turbulence diffusivity in model forest grid cells.

Accounting for updated LAI and land cover fractions in land-atmosphere exchanges

The default WRF-Chem model only considers the dominant land cover type in each grid cell in the calculations of land-atmosphere exchange and dry deposition. Here, we implement a mosaic approach into the Unified Noah land-surface scheme following Li et al.⁸². We mosaic each model grid with the land cover fractions of five major land types from the MODIS dataset. The land-surface exchanges in each model grid are calculated for five major land types and then averaged with land type fractions as weights. The total fractions of the five major land cover types are above 99% in most grid cells. The same mosaic approach is applied in the model dry deposition scheme so that the effects of changes in land cover types on dry deposition can be more accurately represented in this study. We further improved the model dry deposition calculation by accounting for the dependence of canopy stomatal and leaf surface resistances on LAI⁸³. In addition, biogenic emissions estimated by the MEGAN algorithm require inputs of LAI and plant functional types (PFTs). These model inputs are now updated using the MODIS satellite datasets as described above. Estimated biogenic isoprene emissions in China are 0.65, 2.07, and 0.58 Tg C in April, July, and October 2019, respectively, in line with previously reported values^65,84,85. The resulting model simulated HCHO column concentrations are in good agreement with satellite observations over South China (Supplementary Fig. 14).

The improved model simulation is supported by surface ozone measurements at the five background sites in China (Supplementary Fig. 15). In particular, the large daily mean near-surface (7 m/30 m) ozone vertical gradient of 6.3 parts per billion by volume (ppbv) measured at DHS and surface ozone concentrations at Shangri-La (XGLL, a forest background station) can now be well captured in the baseline simulation (Supplementary Figs. 3 and 15, and Supplementary Table 4).

WRF-Chem model sensitivity simulations

We conduct WRF-Chem model simulations for three months (April, July, and October) of 2019 each after a 3-day spin-up for initialization due to limited computing resources. We average values of the three months as the growing season mean. Four sets of model simulations are conducted to examine the impacts of greening from 2005 to 2019 and future forestation from 2019 to 2030 and 2050 on ozone air quality in China (Table 1).

First, the Base simulation uses the model configurations described above and applies MODIS LAI, NDVI, land cover datasets, and anthropogenic emissions for 2019. The baseline simulation represents the present-day condition that can be evaluated with observations. Second, a sensitivity simulation replaces LAI, NDVI, and land cover data with 2005 conditions. The differences in ozone concentrations with Base are then used to estimate the overall effects of greening from 2005 to 2019. Third, a set of sensitivity simulations are performed by (1) fixing NDVI and land cover data to the year 2005 conditions, (2) fixing LAI in dry deposition for the year 2005, and (3) fixing LAI and PFT in MEGAN for the year 2005. Their differences with the Base simulation thus quantify, respectively, the “canopy,” “deposition,” and “emissions” effects. The “canopy” effect of greening is attributable to forest canopy height and land cover changes, as determined by NDVI and land cover data. The “deposition” effect of greening is attributable to changes in dry deposition associated with LAI. The “emission” effect of greening is attributable to changes in BVOC emissions associated with LAI and PFT. It should be noted that the “canopy” effect diagnosed here includes not only changes in turbulence mixing due to forestation but also changes in dry deposition and surface energy balance due to changes in land type, while we find the latter to be small for ozone air quality. The overall greening effects slightly differ from the sum of “canopy”, “deposition”, and “emissions” effects due to nonlinear behaviors. Fourth, a group of sensitivity simulations is conducted by (1) replacing LAI with the 2030 condition (denoted LAI2030), (2) using 2030 future anthropogenic emissions (denoted SSP126 and SSP370), (3) replacing both LAI and anthropogenic emissions with 2030 conditions (denoted SSP126LAI2030 and SSP370LAI2030), and (4) replacing LAI with the 2050 condition (denoted LAI2050). Their differences with the Base simulation reflect, respectively, (1) the effects of greening in 2030, (2) the effects of 2030 future anthropogenic emission changes, (3) the effects of 2030 greening and anthropogenic emission changes jointly, and (4) the effects of greening in 2050 over China. In all simulations, we hold other factors such as meteorology and CO₂ levels constant, allowing us to isolate the effects of greening on ozone.

Ozone-induced GPP damage

We apply an offline version the Yale Interactive Terrestrial Biosphere (YIBs) model, a process-based vegetation model, to evaluate the impacts of ozone-induced vegetation damage on gross primary productivity (GPP) in China. The YIBs model simulates photosynthesis with the Michaelis–Menten enzyme-kinetics scheme, the stomatal conductance with the model of Ball and Berry, and the ozone damage ratio with the semi-mechanistic parameterization to generate the GPP and ozone-induced GPP changes⁸⁶. The YIBs model has been verified thoroughly with site-level measurements and satellite retrievals⁸⁶. We performed the YIBs model across China with a 0.25° spatial resolution in 2019 after a 4-year spin-up for initialization. The model was driven by meteorological fields from Version 2 of Modern Era Retrospective-analysis for Research and Application (MERRA2), and WRF-Chem simulated ozone concentrations. Our estimates of GPP after considering ozone-induced damage are, in general, consistent with the MODIS product (MOD17A2HGF Version 6) across China in 2019 (Supplementary Fig. 16).

Forestation-induced GPP increase

We estimate forestation-induced GPP increases using the light-use efficiency approach, which is widely used in satellite-based studies^87,88,89. This approach decomposes GPP into the fraction of absorbed photosynthetic active radiation (fAPAR), photosynthetic active radiation (PAR), and light use efficiency (LUE). fAPAR can be evaluated via Lambert-Beer’s law.

$${{rm{GPP}}}={{rm{fAPAR}}}times {{rm{PAR}}}times {{rm{LUE}}}$$

(14)

$${{rm{fAPAR}}}=1-{e}^{-{{rm{k}}}times {{rm{LAI}}}}$$

(15)

$$Delta {{rm{GPP}}}={{rm{M}}}times (1-{e}^{-{{rm{k}}}times Delta {{rm{LAI}}}})$$

(16)

where, k is the light extinction coefficient, assuming k = 0.6 for all land types⁸⁹. M is assumed to remain constant when estimating LAI-related changes in GPP.

Future LAI in 2030 and 2050

China plans to enhance its forestry development in the future through technological advancements and phased forestry goals. We obtain a conservative estimate of the Chinese LAI in 2030 by forestation according to the goal of increasing the national forest stock volume (FSV) by 6 billion m³ in 2030 relative to 2005 (13.7 billion m³) to peak CO₂ emissions before 2030 and to achieve carbon neutrality before 2060 (http://english.mee.gov.cn/News_service/media_news/202012/t20201214_812718.shtml, proposed at the Climate Ambition Summit). We also estimate an ambitious LAI based on the goal of expanding FSV to 26.5 billion m³ by 2050 for forestry development proposed at the 19th National Congress (https://bioenergyinternational.com/policy/china-sets-ambitious-forestry-development-goals-2050). Using historical FSV statistics from the NFI yearbook (Supplementary Table 1), we find that mean forest LAI is linearly correlated with FSV over China. Thus, we forecast that the mean forest LAI in 2030 (2050) would increase to 2.16 ± 0.02 m² per m² (2.68 ± 0.05 m² per m²) (Supplementary Fig. 9). Here, we assume that the future LAI increments relative to the 2019 LAI follow the spatial distribution of changes in forest cover fractions during 2005–2019 (Supplementary Fig. 9).

Forest biomass carbon stocks

The biomass expansion factor (BEF, defined as the ratio of forest biomass to timber volume) method is applied to estimate forest stand biomass^42,43. On a national scale, total forest biomass can be calculated as:

$${{mathrm{Y}}}={{rm{A}}}times {{rm{TV}}}times {{rm{BEF}}}={{rm{FSV}}}times {{rm{BEF}}}$$

(17)

where Y, A, TV, and BEF are the national total forest biomass, forest area, mean timber volume, and biomass expansion factor, respectively. Using FSV from the forest inventory dataset (Supplementary Table 1) and estimates of Chinese forest biomass carbon stocks in the literature (Supplementary Table 2), we derive a national BEF value of 0.84–1.07 tons per m³ over China. Previous studies reported BEF values ranging from 0.69–1.47 for China’s major forest types^43,90,91. We estimate future Chinese forest biomass with the BEF value and future FSV and then convert it to the forest biomass carbon stock with a conversion ratio of 0.5.