Tracking surface ozone responses to clean air actions under a warming climate in China using machine learning

Surface ozone, a major air pollutant with important implications for air quality, ecosystems, and climate, shows long-term trends shaped by both anthropogenic and climatic drivers. Here, we developed a machine learning-based approach, namely the fixed emission approximation (FEA), to decouple the effects of meteorological variability and anthropogenic emissions on summertime ozone trends in China under the clean air actions. Anthropogenic emissions drove an approximately +23.2 ± 1.1 µg m⁻³ increase in summer maximum daily 8 h average ozone during 2013–2017, followed by an approximately −4.6 ± 1.5 µg m⁻³ decrease between 2017 and 2020 in response to strengthened emission controls. In contrast, meteorological anomalies, including heatwaves and rainfall conditions, emerged as substantial drivers of ozone variability during 2020–2023. Satellite-derived formaldehyde-to-nitrogen dioxide ratios revealed widespread urban volatile organic compounds-limited regimes for ozone production, with a shift toward nitrogen oxides-limited sensitivity under influence of heatwaves. Extending the FEA framework to assess long-term climate influences from 1970 to 2023, we find that sustained climate warming has driven a substantial increase in urban summertime ozone in China. These results demonstrate that climate change was increasingly offsetting the benefits of emission reductions and highlight the need for integrated ozone mitigation strategies that jointly address emission controls and climate adaptation in a warming world.

1 Introduction

Surface ozone (O₃) is a critical air pollutant that poses significant threats to human health (Knowlton et al., 2004), ecosystems (Agathokleous et al., 2020), and climate (Fishman et al., 1979; Hauglustaine et al., 1994). It forms through complex photochemical reactions involving nitrogen oxides (NO_x) and volatile organic compounds (VOCs) in the presence of sunlight (Jacob, 2000; Wang et al., 2017), exhibiting a nonlinear response to its precursors (Guo et al., 2023; Liu and Shi, 2021; Wang et al., 2023a). Controlling ozone pollution remains a global environmental challenge. In recent years, China has implemented a series of national clean air actions, most notably the Air Pollution Prevention and Control Action Plan (2013–2017) and the Three-Year Action Plan for Winning the Blue-Sky War (2018–2020) (Geng et al., 2024; Zhang et al., 2019; Zheng et al., 2018), that have markedly improved air quality, particularly by reducing fine particulate matter (PM_2.5) (Geng et al., 2024; Zhang et al., 2019). However, surface ozone levels have continued to rise in many regions, raising concerns over the complex drivers of ozone trends and highlighting the need for scientific attribution to guide effective mitigation strategies (Li et al., 2019a; Liu et al., 2023; Wang et al., 2023a; Weng et al., 2022).

Long-term ozone variability is jointly influenced by anthropogenic emissions and weather conditions as well as regional climate (Hallquist et al., 2016; Li et al., 2019b; Wang et al., 2022a). While emission controls directly regulate precursor abundance, climate change modulates ozone through chemical feedbacks, meteorological dynamics, and biosphere–atmosphere interactions (Ma and Yin, 2021; Xue et al., 2020). Over the past century, global surface temperatures have substantially increased relative to the pre-industrial baseline (1850–1900), driven largely by human activities (IPCC, 2021). In such a warming world, extreme climate anomalies – such as heatwaves and persistent rainfall shifts – were expected to be intensified (Li et al., 2025a, b). These events were increasingly recognized as critical modulators of ozone variability through their impacts on photochemistry and precursor emissions (Gao et al., 2023; Pu et al., 2017; Wang et al., 2022a).

Quantifying the respective roles of anthropogenic emissions and meteorological variability in driving ozone trends is therefore essential for evaluating the effectiveness of clean air policies (Li et al., 2019a; Liu et al., 2023). Previous studies have reported rapid increases in surface ozone concentrations in urban cluster regions in China – such as the Beijing–Tianjin–Hebei (BTH) and Yangtze River Delta (YRD) – during the Phase I (2013–2017), with increases of approximately 28 % and 18 %, respectively (Chen et al., 2020; Li et al., 2019a; Liu et al., 2023). In contrast, a modest decline in ozone levels was observed during 2018–2020, largely attributed to emission reductions (Li et al., 2021; Liu and Wang, 2020b; Wang et al., 2024b, 2023a). However, since 2021, observations indicate a renewed increase in ozone concentrations (Fig. S1 in the Supplement). These fluctuations suggest oscillating trends over the past decade, the drivers of which remain poorly constrained.

Two main approaches have been applied to attribute air pollution trends: chemical transport models (CTMs) (Li et al., 2021; Liu et al., 2023; Liu and Wang, 2020a) and data-driven statistical frameworks (Li et al., 2019a, b, 2020). The CTMs simulate atmospheric composition based on emission inventories, meteorological fields, and chemical mechanisms (Ivatt et al., 2022; Liu and Shi, 2021; Liu et al., 2023; Ye et al., 2024). They allow attribution of trend components to emissions or meteorology, and can also resolve sector-specific impacts. However, these models face challenges, including uncertainties and temporal lags in emission inventories. Statistical models, on the other hand, rely on observational datasets and predictor-response relationships without requiring explicit emissions or chemical schemes (Li et al., 2019a, b, 2020; Zhai et al., 2019). With the growing availability of atmospheric big data, machine learning models have emerged as useful tools for trend attribution (Dai et al., 2023; Grange et al., 2018; Vu et al., 2019; Zhang et al., 2025; Zheng et al., 2023). For instance, Grange et al. (2018) developed a random forest-based framework to isolate meteorological influences on particulate matter. Similarly, Wang et al. (2023a) used an enhanced extreme gradient boosting (XGBoost) model to analyze spatial and temporal ozone patterns in China from 2010 to 2021, confirming that emission reductions played a key role in recent declines. Other recent efforts have extended statistical models to long-term assessments of air pollution drivers under climate change (Wang et al., 2022b).

Here, we developed a machine learning-based model framework – fixed emission approximation (FEA) – to quantify the relative contribution of anthropogenic emissions and meteorological conditions in shaping summertime surface ozone trends in China. Applying the FEA to nationwide observational data from 2013 to 2023, we identified three phases of ozone evolution corresponding to the major clean air actions. We further analyzed short-term ozone anomalies associated with extreme weather events, such as the 2022 heatwave and seasonal rainfall. To characterize ozone production regimes, we integrated satellite-derived formaldehyde-to-nitrogen dioxide (HCHO $/$ NO₂, FNR) ratios from the tropospheric monitoring instrument (TROPOMI), revealing spatiotemporal shifts in ozone formation sensitivity across China. Finally, we extend our FEA analysis to evaluate climate-driven ozone trends from 1970 to 2023, using historical meteorological reanalysis data. Collectively, these analyses provide an integrated understanding of how anthropogenic and climatic factors jointly shape surface ozone dynamics under a warming climate.

2 Data and methods

2.1 Data sources and methodology overview

Figure 1 provides an overview of the data analysis and methodological framework employed in this study. We first integrated multi-dimensional datasets, including hourly surface air pollutant concentrations, meteorological reanalysis fields, and satellite remote sensing data. Hourly surface observations of ozone, NO₂, carbon monoxide (CO), and PM_2.5 were accessed from the National Environmental Monitoring Center of China through the open website https://air.cnemc.cn:18007/ (last access: 20 May 2024). Hourly meteorological data with a spatial resolution of 0.25° × 0.25° were sourced from the ERA5 reanalysis dataset provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) and are available for download at https://cds.climate.copernicus.eu (last access: 20 March 2025). Detailed variables are listed in Table S1. The time variables – hour (hour of day) and month (month of year) – are used as emission surrogates to capture regular diurnal and seasonal variations in anthropogenic activity. A similar strategy is widely applied in previous studies about long-term trends in air pollutants (e.g., Grange et al., 2018; Vu et al., 2019) to separate short-term cyclical emission variability from long-term trends. For 2013–2014, the surface MDA8 ozone data were obtained from the Tracking Air Pollution in China (TAP) dataset (Geng et al., 2021), which can be downloaded from http://tapdata.org (last access: 20 May 2024). The TROPOMI on the Sentinel-5P satellite provides global continuous observation data for two indicators of ozone precursors: NO₂ and HCHO column concentrations (Lamsal et al., 2014; Shen et al., 2019). The FNR was used as a proxy for VOC $/$ NO_x reactivity and as a diagnostic indicator of ozone formation sensitivity (Sillman, 1995) to explain and verify the impact of extreme weather and anthropogenic emissions on ozone. Details of the ozone sensitivity diagnostic method are provided in Sect. S1.

Figure 1 Schematic framework of data analysis and methodology. This study integrates multi-dimensional datasets, including ground-based observations, meteorological reanalysis, and satellite remote sensing. A fixed emission approximation (FEA) approach, developed based on the random forest (RF) model, is employed to quantitatively disentangle the contributions of meteorological conditions (MET) and anthropogenic emissions (ANT) to ozone trend variations, and its performance is compared with the conventional meteorological normalization method. The SHAP technique is further applied to assess the influence of extreme weather events, such as heatwaves (HW) and prolonged rainfall (PR). The satellite-derived formaldehyde-to-nitrogen dioxide ratio (FNR) is used to diagnose ozone production sensitivity, to explain and verify the impact of extreme weather and anthropogenic emissions on ozone. Finally, the FEA framework is extended to evaluate the long-term impacts of climate change on ozone trends since 1970.

2.2 Machine learning-based FEA approach

In this study, we develop a machine learning-based FEA approach to assess the impacts of meteorological factors and anthropogenic emissions on the year-round variations in ozone concentrations (Fig. 1). First, a regression model is constructed using the random forest algorithm to relate ozone concentrations to temporal emission surrogates and to meteorological parameters at multiple atmospheric levels. These temporal emission surrogates, including month and hour, represent short-term regular emission patterns (e.g., diurnal cycles), thereby enabling the model to isolate the long-term emission-driven component of ozone changes (Grange et al., 2018; Meng et al., 2025; Shi et al., 2021; Vu et al., 2019). The meteorological parameters include 18 distinct variables at different altitudes (see Table S1). It should be noteworthy that surface air pollutant observations for each city represent multi-site averages rather than data from a single monitoring station, which reduces the influence of local representativeness uncertainty. The meteorological data are obtained from the nearest grid cell corresponding to each city, ensuring spatial consistency between the pollutant and meteorological datasets. This approach was similar to the methodologies widely adopted in previous studies (Shi et al., 2021; Wang et al., 2025; Yao et al., 2024; Zheng et al., 2023). Our modeling strategy involves building and predicting models for individual cities and for each year from 2015 to 2023, which helps in minimizing the uncertainty caused by surface heterogeneity. Due to the lack of available observational data for many cities in 2013 and 2014, we did not develop models for these two years. In our approach, 80 % of the dataset is used for model training, while the remaining 20 % is reserved for testing. We perform ten-fold cross-validation and assess model performance using seven statistical metrics, as listed in Table S2.

Following the construction of the machine learning models for individual cities and years, we introduce the FEA approach. The key principle of FEA is the assumption that the total emissions of ozone precursors remain unchanged from the baseline year. Specifically, we establish hourly-resolution models for the baseline year (i) during the summer season (June to August) as a reference for anthropogenic emissions, represented by the pink solid line in Fig. 1. These models are then applied to predict ozone concentrations under the meteorological conditions of the prediction year (j), while holding the emission levels constant at those of the baseline year (i), as shown by the blue solid line in Fig. 1. The difference between the predicted values (Pred_i) and the observed values (OBS_i) for the baseline year (i) represents the model residuals (RES_i), as shown in Eq. (1). The difference in observed MDA8 ozone concentrations between baseline year i and prediction year j is driven by the differences in meteorological conditions MET_i(j) and anthropogenic emission controls ANT_i(j) (Eq. 2).The prediction result Pred_i(j) obtained by applying the model trained with data from year i to the meteorological conditions of year j, the difference between Pred_i(j) and Baseline (Pred_i) is driven by MET_i(j), while the difference between Pred_i(j) and the observed levels in year j (OBS_j), minus the RES_i, yields the ozone variation driven by (ANT_i(j)). Therefore, MET_i(j) and ANT_i(j) can be quantified and calculated using Eqs. (3)–(4).

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathrm{OBS}}_i={\mathrm{Pred}}_i+{\mathrm{RES}}_i,\text{(2)}& {\displaystyle }{\mathrm{OBS}}_{i\left(j\right)}={\mathrm{MET}}_{i\left(j\right)}+{\mathrm{ANT}}_{i\left(j\right)},\text{(3)}& {\displaystyle }{\mathrm{MET}}_{i\left(j\right)}={\mathrm{Pred}}_{i\left(j\right)}-{\mathrm{Pred}}_i,\text{(4)}& {\displaystyle }{\mathrm{ANT}}_{i\left(j\right)}={\mathrm{OBS}}_j-{\mathrm{Pred}}_{i\left(j\right)}-{\mathrm{RES}}_i,\end{array}$$

The difference in observed MDA8 ozone concentrations between two different prediction years (j₁, j₂) is driven by the differences in meteorological conditions (ΔMET_i(j1,j2)) and anthropogenic emission controls (ΔANT_i(j1,j2)) (Eq. 5). The term ΔMET_i(j1,j2) represents the changes in meteorological conditions and can be calculated by the difference between the predicted values, Pred_i(j1) and Pred_i(j2), for the corresponding years (Eq. 6). Similarly, the value of ΔANT_i(j1,j2), representing the change in anthropogenic emissions between the two years j₁ and j₂, can be therefore calculated using Eq. (7). By performing these calculations, we can isolate and quantify the contributions of meteorological conditions and anthropogenic emission controls to the observed ozone trends. We used a cross-matrix research method to assess the uncertainty of FEA, with specific formulas available in Sect. S2 in the Supplement.

$$\begin{array}{}\text{(5)}& {\displaystyle }\mathrm{\Delta }{\mathrm{OBS}}_{(j\mathrm{1},j\mathrm{2})}=\mathrm{\Delta }{\mathrm{MET}}_{i(j\mathrm{1},j\mathrm{2})}+\mathrm{\Delta }{\mathrm{ANT}}_{i(j\mathrm{1},j\mathrm{2})},\text{(6)}& {\displaystyle }\mathrm{\Delta }{\mathrm{MET}}_{i(j\mathrm{1},j\mathrm{2})}={\mathrm{Pred}}_{i\left(j\mathrm{2}\right)}-{\mathrm{Pred}}_{i\left(j\mathrm{1}\right)},\end{array}$$

$$\begin{array}{}\text{(7)}& \begin{array}{rl}\mathrm{\Delta }{\mathrm{ANT}}_{i(j\mathrm{1},j\mathrm{2})}& ={\mathrm{ANT}}_{i\left(j\mathrm{2}\right)}-{\mathrm{ANT}}_{i\left(j\mathrm{1}\right)}\\ & =\left({\mathrm{OBS}}_{j\mathrm{2}}-{\mathrm{Pred}}_{i\left(j\mathrm{2}\right)}-{\mathrm{RES}}_i\right)\\ & -\left({\mathrm{OBS}}_{j\mathrm{1}}-{\mathrm{Pred}}_{i\left(j\mathrm{1}\right)}-{\mathrm{RES}}_i\right)\\ & =\left({\mathrm{OBS}}_{j\mathrm{2}}-{\mathrm{OBS}}_{j\mathrm{1}}\right)\\ & -\left({\mathrm{Pred}}_{i\left(j\mathrm{2}\right)}-{\mathrm{Pred}}_{i\left(j\mathrm{1}\right)}\right),\end{array}\end{array}$$

Model performance was first evaluated through ten-fold cross-validation, revealing high predictive skill between observed and predicted MDA8 ozone levels during 2015–2023 in the BTH regions (Fig. S2). The index of agreement (IOA) ranged from 0.96 to 0.97, with correlation coefficients (R) between 0.93 and 0.95. Root mean square errors (RMSE) and normalized mean bias (NMB) varied from 16.9 to 21.9 µg m⁻³ and 8 % to 25 %, respectively, indicating high model accuracy. Nationally, the model yielded R values of 0.88–0.91 and IOA of 0.93–0.95, with errors remaining within acceptable ranges (Tables S3–S8). To assess uncertainty stemming from interannual model training variability, we applied a matrix-based resampling approach (see Sect. S2). As shown in Fig. 2, the relative difference in residuals ranged from −9 % to 3 %, and remained within ±12 % for all regions – supporting the robustness of the FEA method. Notably, inclusion of time-related variables could reduce model uncertainty compared to simulations excluding these predictors. The average uncertainty decreased by approximately 2 %–4 % at the regional-mean level (Fig. S3).

Figure 2 Uncertainty assessment of the FEA method. The uncertainty for the FEA method is calculated using the approach described in Sect. S2. The x-axis represents the years used for model training, and the y-axis represents the years predicted by the trained model. The diagonal line in each sub-panel represents the changes in the residuals of the models.

2.3 Weather normalization analysis

To compare the FEA method with other commonly used statistical approaches, we also applied the widely adopted meteorological normalization technique based on the RF algorithm (Grange et al., 2018; Vu et al., 2019). This approach constructs a regression model that relates air pollutant concentrations to meteorological parameters and emission surrogate indicators (i.e., time variables such as unix time, day of year, day of month, and hour of day) (Grange et al., 2018; Vu et al., 2019). Once the model is trained, air pollutant concentrations are predicted by randomly resampling meteorological variables from long-term historical meteorological datasets, thereby generating a new ensemble of predictions (Vu et al., 2019). These predictions are made under consistent meteorological conditions, enabling the isolation of meteorological influences from anthropogenic emission effects on air pollutant trends. The resulting weather-normalized pollutant concentrations (Fig. 1) represent the levels expected under average meteorological conditions, thus reflecting the impact of emission changes alone. In this study, the meteorological normalization follows this established framework, with meteorological variables randomly sampled from the long-term dataset spanning 1970–2023. Each normalization process involves 1000 iterations, and the arithmetic mean of these iterations' simulated values was adopted as the final normalized result. The alignment between FEA-based and weather-normalized trends (Fig. S4) affirms the robustness of the FEA framework.

2.4 Quantification of extreme weather-driven changes in ozone

An unprecedented and persistent heatwave struck central and eastern China during the summer of 2022, with the YRD experiencing the most severe impacts (Wang et al., 2023b; Zhang et al., 2023). This event has been identified as the longest-lasting and most intense heatwave since at least 1961 (Mallapaty, 2022). In contrast, the Yangtze-Huaihe region is climatologically prone to sustained extreme precipitation, where prolonged rainfall episodes frequently occur during the East Asian summer monsoon (Yin et al., 2020). Together, the extreme heatwave (HW) in 2022 and recurrent prolonged rainfall (PR) events provide unique and physically realistic atmospheric conditions to investigate the impacts of typical weather extremes on surface ozone.

To quantify the contributions of extreme meteorological conditions to ozone variability, we applied the SHapley Additive exPlanations (SHAP) method (Lundberg et al., 2020) to interpret predictions from the random forest model. SHAP assigns an importance value to each input feature k, representing its marginal contribution to the model-predicted MDA8 ozone. The PR period was defined as 15 June to 15 July for each year, while the remaining period from June to August was classified as the non-prolonged rainfall (NPR) period. The HW event in 2022 was defined as 16 July to 31 August, with the same calendar period in other years designated as non-heatwave (NHW) conditions.

SHAP values were calculated for all input features during the PR and NPR periods, as well as during the HW and NHW periods, respectively. The relative changes in SHAP values (ΔSHAP) between these conditions were used to assess the responses of MDA8 ozone to the rainy season or the 2022 heatwave weather conditions, as per the following equations:

$$\begin{array}{}\text{(8)}& {\displaystyle }\mathrm{\Delta }{\mathrm{SHAP}}_k={\mathrm{SHAP}}_{\mathrm{PR}}-{\mathrm{SHAP}}_{\mathrm{NPR}}\text{(9)}& {\displaystyle }\mathrm{\Delta }{\mathrm{SHAP}}_k={\mathrm{SHAP}}_{\mathrm{HW}}-{\mathrm{SHAP}}_{\mathrm{NHW}}\end{array}$$

2.5 FEA-based assessment of climate change impacts on ozone

To evaluate the long-term impact of climate change on surface ozone concentrations across China from 1970 to 2023, we extended the framework of our machine learning-based FEA method. The core idea of this analysis is to isolate the influence of long-term meteorological variability on ozone while assuming fixed anthropogenic emissions. Given the availability of relatively complete and continuous hourly ozone observations and meteorological data from 2015 to 2023, this period was used to construct nine emission baseline scenarios. Following the modeling protocol described in Sect. 2.2, nine independent random forest models were trained for each city and scenario, with each year from 2015 to 2023 serving as an emissions reference. Model inputs included hourly ozone observations, key meteorological predictors, and time-related variables (hour of day and month of year). The trained models were then applied to historical meteorological reanalysis data from 1970 to 2023 to simulate ozone trends under fixed emissions (Fig. 1), yielding nine independent ozone trajectories, each reflecting the influence of long-term meteorological variability under a different fixed-emissions assumption. While the choice of emission baseline may affect the absolute magnitude of simulated ozone, it does not alter the primary objective: assessing the sensitivity of surface ozone to meteorological drivers over multidecadal timescales (Lecœur et al., 2014; Leung et al., 2018; Wang et al., 2022b). This approach could capture the climate-induced ozone signal while adopting the common assumption that emissions are not themselves influenced by climate change – a simplification consistent with prior attribution studies (Dang and Liao, 2019; Leung et al., 2018; Shen et al., 2017; Wang et al., 2022b). For comparison, we also estimated the impact of anthropogenic emission changes on ozone concentrations during 2015–2023 using the same FEA methodology and the complete hourly dataset for model training. This dual-track analysis enables a clear distinction between the impacts of climate variability and emission mitigation on observed ozone trends.

Figure 3 Spatial distribution of summertime MDA8 ozone, surface NO₂, and TROPOMI NO₂, HCHO across major city clusters in China. The panels represent the average MDA8 ozone, surface NO₂, and TROPOMI NO₂, HCHO column concentrations for 354 cities in China during the summertime (June–August) from 2018 to 2023. The corresponding five regions includes BTH (37–41° N, 114–118° E); YRD (30–33° N, 118.2–122° E); SCB (28.5–31.5° N, 103.5–107° E); PRD (21.5–24° N, 112–115.5° E) and FWP (33–35° N, 106.25–111.25° E, and 35–37° N, 108.75–113.75° E).

To examine the sensitivity of urban ozone pollution to climate variability under different potential atmospheric conditions (e.g., oxidation capacity) and its possible evolution under global warming, we defined three representative regional scenarios based on typical ozone pollution characteristics in China (Fig. 3a): a high-pollution scenario for BTH (BaseBTH), a moderate-pollution scenario for YRD (BaseYRD), and a low-pollution scenario for Pearl River Delta (PRD) (BasePRD). These scenarios allow assessment of ozone trends and climate impacts under fixed emissions across three distinct atmospheric conditions.

3 Results and Discussion

3.1 Spatiotemporal variation of summertime ozone

Figure 3 presents the spatial distribution of the average summertime (2018–2023) maximum daily 8 h average (MDA8) ozone, surface NO₂, and TROPOMI NO₂, HCHO column concentrations across China, along with the locations of the country's five typical city clusters: BTH, Fenwei Plain (FWP), YRD, Sichuan Basin (SCB), and PRD. Across these five city clusters, the average summer ozone concentrations ranged from 88.9 to 161.3 µg m⁻³ – substantially exceeding the 43.0 µg m⁻³ threshold associated with ecosystem productivity loss (Gong et al., 2021) and the World Health Organization (WHO, 2021) recommended peak seasonal average of 60 µg m⁻³. TROPOMI satellite observations of NO₂ column concentration show notably elevated concentrations over the five city clusters, particularly in the BTH, YRD, and FWP, which align with surface NO₂ distribution patterns and confirm the scale of anthropogenic NO_x emissions in these regions (Zheng et al., 2021). TROPOMI satellite observations of HCHO column concentrations similarly reveal these city clusters as hotspots for VOC emissions (Fig. 3d). These concurrent high levels of NO₂ and HCHO suggest a strong photochemical ozone pollution potential, as the abundant precursors in these urban clusters could drive substantial ozone production during the summer months.

Figure 4 Interannual trends of summertime MDA8 ozone across China (2013–2023). (a) The seasonal variations of MDA8 ozone during the summer months (June, July, and August) over China. (b–f) The average trend across the five city cluster regions in China: BTH, FWP, YRD, SCB, and PRD, respectively. The summer months are defined according to meteorological seasonality, encompassing June, July, and August. In the violin plots, hollow diamond markers denote the mean, while solid diamond markers represent the median. The Mann-Kendall test and Sen's slope estimator were employed to assess the statistical significance and rate of change in the monthly average MDA8 ozone concentrations. The p value represents the significance level from the Mann-Kendall test, which is used to determine the statistical significance of the trend in the data.

Figure 4 presents the interannual variations in MDA8 ozone concentrations during summertime across China, with a focus on the five city cluster regions. During the Phase I (2013–2017), the average nationwide MDA8 ozone increased significantly (p < 0.01), rising from 95.5 to 118.0 µg m⁻³. This growth was especially pronounced in the BTH and FWP regions, where the concentrations increased by approximately 38 % and 41 %, respectively. In contrast, ozone increases were more modest in the YRD (∼ 11 %), SCB (∼ 15 %), and PRD (∼ 16 %) regions, respectively. These results were consistent with the previous studies (Li et al., 2021; Liu and Wang, 2020a, b; Wang et al., 2023a). Corresponding to the implementation of more stringent emission controls on NO_x and VOCs emissions during the Phase II (Geng et al., 2024; Liu et al., 2023), a moderate national decrease in MDA8 ozone was observed, with concentrations declining to 109.0 µg m⁻³ from 2017 to 2020. The declines during this period were most notable in FWP (−16 %) and YRD (−15 %), while BTH (−6 %), SCB (−11 %), and PRD (−4 %) also showed reductions compared to their concentration peaks observed in 2017. However, the MDA8 ozone rebounded, reaching 118.4 µg m⁻³ in 2023 – comparable to its 2017 peak – with a particularly sharp increase during the summer of 2022. From 2021 to 2023, MDA8 ozone concentrations rose by 2.8 µg m⁻³ in BTH, 3.1 µg m⁻³ in FWP, 16.1 µg m⁻³ in YRD, and 18.5 µg m⁻³ in SCB, respectively.

Figure S1 further illustrates the spatiotemporal evolution of summertime MDA8 ozone in China from 2013 to 2023. On average, approximately 68 % of the cities exceeded the WHO air quality guideline of 100 µg m⁻³ for the MDA8 ozone. Elevated ozone levels were primarily observed in densely populated and economically developed regions. Spatially, ozone hotspot regions expanded between 2013 and 2017 (Fig. S1a–e), followed by contraction during 2018–2020 (Fig. S1f–i). However, this progress stalled in 2021. A sharp reversal was observed in 2022, with widespread increases in MDA8 ozone (Fig. S1k). These changes could be closely linked to emission control measures and meteorological conditions, which will be further discussed in Sect. 3.2 and 3.3.

3.2 Anthropogenic emission drivers of ozone trends

To disentangle the relative impacts of anthropogenic emissions and meteorological variability on observed ozone trends, we applied the machine learning-based FEA framework described in Sect. 2.2. As illustrated in Fig. 5, anthropogenic emissions were the dominant driver of ozone increases during 2013–2017, contributing an average rise of approximately 23.2 ± 1.1 µg m⁻³ nationwide. The most pronounced increases occurred in the FWP and BTH (45.0 ± 2.0 and 42.1 ± 2.0 µg m⁻³, respectively), whereas the PRD exhibited a relatively modest enhancement (13.4 ± 1.6 µg m⁻³), reflecting its predominantly NO_x-limited photochemical regime versus VOC-limited regimes in other regions (Ren et al., 2022). As shown in Fig. S5, the precursor gases NO₂ and CO exhibited regionally distinct decreasing trends, partially explaining the spatial heterogeneity of ozone changes. The MDA8 ozone decreased by 10.5 ± 2.0 µg m⁻³ in BTH and 10.4 ± 3.0 µg m⁻³ in FWP, with smaller declines in YRD (−4.8 ± 3.8 µg m⁻³), SCB (−2.8 ± 2.4 µg m⁻³), and PRD (−6.6 ± 1.4 µg m⁻³) between 2017 and 2020 (Fig. 5). These trends were overall consistent with those derived using independent statistical approaches (Wang et al., 2023a). The COVID-19 pandemic (January–April 2020) introduced an unprecedented perturbation to anthropogenic activity, leading to sharp declines in industrial production, energy consumption, and transportation (Shi and Brasseur, 2020; Zheng et al., 2021). National emissions of SO₂, NO_x, PM_2.5, and VOCs were estimated to have decreased by 0.37 Tg (12 %), 0.87 Tg (13 %), 0.25 Tg (10 %), and 1.07 Tg (12 %), respectively, relative to the same period in 2019 (Geng et al., 2024). Despite these reductions, MDA8 ozone concentrations increased by 1.7–2.3 µg m⁻³ across BTH, FWP, YRD, and SCB, while a slight decrease occurred in PRD (Fig. S6). Overall, ∼ 79 % of monitored cities experienced ozone increases, with a national mean enhancement of 2.1 ± 1.3 µg m⁻³ (Fig. S7). In the post-pandemic period (2020–2023), concentrations of NO₂, CO, and PM_2.5 stabilized or declined more gradually (Fig. S5), and the contribution of anthropogenic emissions to ozone variability weakened considerably (Fig. S8). Regionally, emission-driven changes ranged from −1.2 to +2.6 µg m⁻³ in BTH, −1.6 to +4.0 µg m⁻³ in FWP, −4.7 to +7.4 µg m⁻³ in YRD, −3.6 to +3.0 µg m⁻³ in SCB, and −3.8 to +7.7 µg m⁻³ in PRD. These results indicate that while emission controls initially yielded substantial ozone mitigation benefits during the Phase II, their effectiveness has gradually diminished, underscoring the need for more targeted and region-specific emission control strategies under evolving photochemical regimes.

Figure 5 Anthropogenic and meteorological drivers of ozone trends from 2013 to 2023. Changes in summertime MDA8 ozone concentrations were decomposed into contributions from anthropogenic emissions and meteorological variability using the FEA framework. Results reflect ensemble estimates based on multiple baseline years (2015–2023) for emissions. The interquartile range, with values in parentheses denoting the 25th and 75th percentiles across all baseline scenarios.

Figure 6 Trends in the distributions of ozone production sensitivity regimes. Fractions of VOC-limited, NO_x-limited, and transitional ozone sensitivity regimes across five key regions during the summertime (June to August) from 2018 to 2023, based on the FNR analysis. (a–e) The trend across the five city cluster regions in China during the summer months (June, July, and August): BTH, FWP, YRD, SCB, and PRD, respectively. (f) The overall trends for all five regions.

Satellite retrievals of tropospheric NO₂ and HCHO from TROPOMI (Figs. S9–S10) further reveal evolving ozone production chemistry. NO₂ columns exhibited strong east-west gradients, with eastern China maintaining levels five times higher than the west. Between 2018 and 2023, NO₂ columns over the North China Plain (NCP) declined from 4.13 × 10¹⁵ to 3.85 × 10¹⁵ molec. cm⁻², while HCHO remained stable until 2021, followed by a sharp increase in 2022. The spatial pattern of temperature anomalies between heatwave and non-heatwave periods (Fig. S11) reveals strong positive differences in the YRD and SCB, consistent with enhanced biogenic and anthropogenic VOC emissions under extreme heat (Qin et al., 2025; Tao et al., 2024). By 2023, HCHO concentrations had returned to pre-heatwave levels. To diagnose the evolving chemical sensitivity of ozone production, we examined the spatiotemporal evolution of the HCHO $/$ NO₂ ratio (Sect. S1). Figure S12 shows that this ratio exhibited regionally distinct transitions from 2018 to 2023, reflecting shifts in photochemical regimes. Figure 6 summarizes the relative contributions of VOC-limited, NO_x-limited, and transitional regimes across the five key regions. In BTH, NO_x-limited areas accounted for ∼ 82 % of the domain in 2018 and remained above 45 % thereafter, while VOC-limited regions declined from ∼ 14 % to ∼ 2 %. In FWP, summertime ozone formation was largely controlled by NO_x-limited and transitional regimes. The YRD underwent a notable shift from VOC- to NO_x-limited chemistry, with VOC-limited fractions decreasing from ∼ 35 % in 2018 to ∼ 22 % in 2023, particularly during 2022 when extreme heat amplified VOC emissions and photochemical activity (Qin et al., 2025; Tao et al., 2024). The SCB region consistently exhibited strong NO_x limitation (> 75 %), whereas the PRD showed a gradual expansion of the transitional regime alongside a modest contraction of VOC-limited regions. These shifts in photochemical sensitivity correspond well with the ozone decrease observed during the Phase II emission reductions. Spatial distributions of ozone formation sensitivity during the COVID-19 lockdown (Fig. S13) reveal that most of China was in a transitional regime, with major urban clusters remaining VOC-limited and only limited areas in southern China being NO_x-limited. This spatial pattern aligns with the observed widespread ozone increases during the lockdown (Fig. S7). These findings highlight that ozone production chemistry in China was shaped by the complex interplay between emission reduction efforts and the rising frequency of meteorological extremes under a warming climate.

3.3 Meteorological impact on ozone variation

Figure 5 shows the interannual meteorological contributions to summertime MDA8 ozone across different emission-control phases. During the Phase I, meteorology exerted relatively weak influences on ozone variability, with contributions ranging from −4.8 to +3.9 µg m⁻³ – far smaller than those from anthropogenic emission changes. In contrast, meteorological anomalies became a decisive factor from 2017 to 2020, driving substantial ozone reductions. Ozone decreases attributable to meteorology reached −14.4 ± 3.0 µg m⁻³ in the FWP, −15.9 ± 3.8 µg m⁻³ in the YRD, and −11.1 ± 2.4 µg m⁻³ in the SCB, explaining 58 ± 12 %, 77 ± 18 %, and 80 ± 17 % of the total ozone decline, respectively. A notable shift occurred from 2020 to 2023, when the influence of extreme meteorological events increasingly dominated ozone variability. In the summer of 2022, persistent heatwaves across eastern and southern China (Mallapaty, 2022; Wang et al., 2023b) triggered sharp ozone increases of +20.8 ± 3.6 µg m⁻³ in the YRD and +22.1 ± 3.2 µg m⁻³ in the SCB, reflecting the enhanced photochemical activity under high-temperature and intense solar radiation conditions. The following summer (2023) featured anomalously heavy rainfall, resulting in sharp ozone suppression (−17.8 ± 2.3 µg m⁻³ in the YRD and −9.7 ± 3.3 µg m⁻³ in the SCB). This reduction coincided with a remarkable increase in precipitation, i.e., 102 % in YRD and 35 % in SCB (Fig. S14), indicating that rainy meteorological conditions may have suppressed ozone production.

To identify the dominant meteorological drivers, we analyzed Gini importance scores derived from the RF model across 18 predictor variables (Fig. S15). Temperature (T) emerged as the most influential predictor in the BTH and FWP regions, while shortwave solar radiation (SR), relative humidity (RH), and 850 hPa zonal wind (u850) were most important in the YRD. In the PRD, ozone variability was primarily governed by temperature and transport-related indices, including meridional winds at different altitudes. These findings are consistent with the climatological contrast between northern continental and southern coastal regimes: in northern China, stagnant anticyclonic conditions (Gong and Liao, 2019) and strong solar radiation promote photochemical production (Bao et al., 2025), whereas in southern regions, high humidity and convective rainfall could tend to suppress ozone by reducing actinic flux and enhancing removal of precursors(Lu et al., 2019).

Partial dependence analysis (Fig. S16) further illustrates the nonlinear responses of ozone to key meteorological factors (T, RH, and SR) for representative cities in each cluster, revealing clear regional contrasts. In Beijing (BTH), ozone concentrations show the strongest positive response to T (Fig. S16a), consistent with the enhancement of reaction kinetics and biogenic VOC emissions under hot conditions. This behavior reflects the thermodynamic coupling between surface heating, boundary-layer expansion, and photochemical production. In Nanjing (YRD), ozone was more sensitive to solar radiation than to temperature (Fig. S16c), highlighting the dominant role of actinic flux in controlling radical production during warm and dry conditions. Consistent with these findings, Yang et al. (2024) reported that high-temperature and low-RH conditions over the NCP and YRD could enhance photochemical ozone formation, with chemical production dominating during peak pollution periods. In the SCB, both T and RH exhibited strong influences, while ozone variability was shaped primarily by T and large-scale circulation patterns in the PRD associated with subtropical maritime flow and typhoon incursions from the Northwest Pacific (Chen et al., 2024; Wang et al., 2024a, 2022a).

Figure 7 Meteorological impact on predicted ozone concentrations under heatwave and rainy weather conditions. (a) The differences in SHAP values (ΔSHAP) between heatwave and non-heatwave periods in the Yangtze-Huaihe region during summer 2022. (b) The ΔSHAP between prolonged rainfall periods and non-prolonged rainfall periods in the same region from 2015 to 2023. Box plots show the distribution of ΔSHAP across cities; the center line indicates the median, boxes denote the interquartile range (25th–75th percentiles), and the whisker line extends to one standard deviation.

To further quantify these relationships, we applied SHAP (SHapley Additive exPlanations) analysis to evaluate the meteorological influence of the HW and the PR events in the Yangtze-Huaihe region between 2015 and 2023 (Sect. 2.4). As shown in Fig. S17, the HW events were associated with strong positive SHAP values in southeastern coastal areas, especially the YRD and SCB, driven by elevated SR and T. Mean SR during the HW periods was substantially higher than during the NHW periods (Fig. S18), increasing photochemical activity through increased radical generation and faster reaction rates. Conversely, PR events produced consistent negative SHAP contributions across all regions (Fig. S19), reflecting the combined effects of reduced photolysis, increased humidity, and efficient wet scavenging on ozone production (He and Carmichael, 1999). A multi-year comparison (Fig. 7) highlights the opposing effects of key meteorological variables – including RH, T, boundary layer height (BLH), total precipitation (TP), and surface pressure (SP) – on MDA8 ozone. The trend in ΔSHAP values under high-humidity conditions from 2015 to 2023 (Fig. S20) further confirms the model's ability to capture the suppressive effects of wet weather conditions on ozone formation.

3.4 Reshaping distributions of ozone by climate change and emission controls

To assess the long-term influence of climate change on surface ozone concentrations, we applied the FEA framework to simulate summertime ozone trends over the period 1970–2023. In this analysis, anthropogenic emissions were held constant at their 2015–2023 summertime levels, while interannual variations in meteorological variables were introduced using historical reanalysis data. This design isolates the climate-driven component of ozone trends while assuming that emission trajectories are independent of climate change – a simplification aligned with prior attribution frameworks (Wang et al., 2022b). The impact of anthropogenic emission controls was estimated by comparing observed ozone concentrations with FEA-predicted values during 2015–2023, thereby quantifying the residual effect of emissions under fixed meteorology.

Figure 8 Impact of climate change and emission controls on ozone trends. Left panels show ozone trends attributable to long-term climate change from 1970 to 2023, simulated under fixed emission scenarios using the FEA framework. Right panels depict ozone trends from 2015 to 2023, reflecting the impact of anthropogenic emission controls. Each trajectory represents results based on a distinct emissions baseline year. Shaded grey areas indicate the interquartile range (25th–75th percentiles), solid red lines denote trend estimates, and light red shading marks the 5th–95th percentile confidence intervals. Statistical significance and trend slopes were assessed using the Mann-Kendall test. More details of the sensitivity tests are provided in Fig. S18.

As shown in Fig. 8, under the 2015–2023 emission levels, climate change has exerted a statistically significant (p < 0.05) positive influence on urban summertime ozone concentrations across China, resulting in a nationwide increase of approximately 0.06 µg m⁻³ yr⁻¹ since 1970. All five major urban regions displayed upward trends, with the most pronounced increase observed in the BTH and SCB at 0.12 µg m⁻³ yr⁻¹. Three sensitivity simulations (see Sect. 2.5 and Fig. S21) confirm this robustness: trend slopes range from 0.11–0.14 µg m⁻³ yr⁻¹ in BaseBTH (high-pollution scenario), 0.05–0.10 µg m⁻³ yr⁻¹ in the BaseYRD (moderate-pollution scenario), and 0.03–0.10 µg m⁻³ yr⁻¹ in the BasePRD (low-pollution scenario). Despite regional differences in chemical regimes or pollution levels, the consistent upward tendencies underscore the pervasive climatic amplification of ozone formation. These results emphasize that climate change acts as a systematic driver of ozone growth across diverse atmospheric environments, reinforcing the need to embed climate resilience within emission control strategies. Spatial correlations between climate-driven ozone increases and temperature changes (Fig. S22) further confirm that warming is the dominant contributor to long-term ozone enhancement. In particular, the correlation coefficients between ozone trends and temperature anomalies reached 0.90 (BTH), 0.89 (FWP), 0.72 (YRD), and 0.93 (SCB), indicating a strong temperature dependence of climate-induced ozone formation in these regions. The PRD showed a weaker correlation, likely due to its unique subtropical maritime climate and higher humidity and cloud cover, which tend to suppress photochemical ozone production (Yang et al., 2019). The right panel of Fig. 8 depicts summertime ozone trends from 2015 to 2023 under the combined influence of anthropogenic emissions, derived from the FEA method. Ozone concentrations rose across all regions between 2015 and 2018, declined modestly during 2018–2020, and rebounded thereafter in most regions except the PRD.

These findings are consistent with future projections that anticipate more frequent high-ozone episodes under continued warming (Li et al., 2023). Recent analyses (Yang et al., 2024) show that the frequency of high-temperature and low-humidity conditions during 2000–2019 was markedly higher than in 1980–1999, suggesting that ozone pollution in both the NCP and YRD has intensified under historical climate change. Indeed, while national emission controls curbed ozone growth after 2018, a post-2020 rebound has emerged, implying that the climatic penalty on ozone is beginning to offset emission gains. The extreme 2022 heatwave exemplified this effect, substantially enhancing photochemical activity and shifting ozone sensitivity from VOC-limited to transitional or NO_x-limited regimes. Although reductions in anthropogenic precursor emissions have improved ozone control efficiency, the warming-induced enhancement of ozone increasingly interferes with – and in some regions may partially offset – the air-quality gains achieved through emission reduction efforts.

4 Conclusions and implications

In this study, we developed a machine learning-based FEA framework to disentangle and quantify the respective roles of anthropogenic emissions and meteorological drivers in shaping ozone trends during 2013–2023. With a national-level prediction uncertainty of approximately 6 %, the FEA method provides a computationally efficient and scalable tool for diagnosing atmospheric variability across large spatial and temporal domains. However, the current model framework did not explicitly resolve grid-scale spatial heterogeneity, vegetation feedbacks, or land-use dynamics, which may influence the ozone prediction. In addition, the sensitivity of the results to spatial resolution need further investigation through coupled applications of machine learning and chemical transport models.

Our results revealed that increased anthropogenic emissions were the dominant driver of the sharp rise in summertime MDA8 ozone concentrations during the Phase I, contributing an average increase of 23.2 ± 1.1 µg m⁻³. In contrast, the strengthened clean air actions during the Phase II – particularly the synergistic control of NO_x and VOCs – led to measurable reductions in MDA8 ozone, with national-average declines of 4.6 ± 1.5 µg m⁻³ from 2017 to 2020. These improvements were especially evident in regions such as BTH and FWP, where ozone formation was highly sensitive to VOCs. However, the impact of emission reductions diminished considerably during the recent period (2021–2023). This stagnation underscores the urgent need for more targeted, region-specific emission control strategies that address the shifting photochemical sensitivity of ozone formation regimes.

Using the SHAP attribution analysis, we further quantified the influence of meteorological extremes on ozone variability. Record-breaking heatwaves in 2022 enhanced ozone concentrations by up to +5.8 µg m⁻³, while prolonged pluvial episodes, particularly during the East Asian monsoon season, suppressed ozone by as much as −15.2 µg m⁻³. These results highlight the dominance of short-term meteorological extremes in shaping ozone air quality under a warming climate. Complementary satellite-based FNR diagnostics revealed that from 2018 to 2023, summertime ozone formation was predominantly influenced by NOx-limited and transitional regimes, while VOC-limited areas experienced a decline.The 2022 heatwave induced regime shifts in regions such as the YRD, where intensified VOC emissions and elevated temperatures drove transitions toward NO_x-limited conditions. These findings emphasize the need for dynamic, region-specific assessments of ozone sensitivity to guide effective mitigation strategies.

To assess the climate penalty on ozone, we extended the FEA framework to simulate long-term trends from 1970 to 2023 by fixing emissions and allowing meteorological variables to evolve with observed climate trends. Our findings show that climate change has contributed to a significant upward trend in urban summertime ozone, averaging 0.06 µg m⁻³ yr⁻¹, with particularly strong increases in the BTH and SCB. Good correlations between the modelled ozone and surface temperature (r = 0.72–0.93) across major urban clusters indicated that climate warming exerts a persistent control on the long-term evolution of ozone. While reductions in precursor emissions have improved ozone control efficiency, the direct enhancement of ozone by rising temperatures increasingly interferes with, and in some regions may partially offset, the air-quality benefits achieved through emission mitigation. Together, these findings highlight that effective ozone management in a warming world will require integrated strategies that jointly address emission reductions and climate adaptation.