Multi-Temporal Scale Response Analysis and Exploration of Ozone Pollution over East Asia

1. Introduction

Since the implementation of the Air Pollution Prevention and Control Action Plan in China in 2013, national PM_2.5 concentrations have achieved substantial reductions. However, in the early stages, ground-level O₃ pollution received relatively limited attention. To address this issue, the Chinese government has progressively strengthened the coordinated reduction of VOCs and NOx since 2017. The State Council successively issued the Three-Year Action Plan for Winning the Blue Sky Defense War (2018) and the Action Plan for Continuous Improvement of Air Quality (2023), explicitly incorporating ozone control into the national air quality strategy and emphasizing regional joint prevention and control of PM_2.5 and O₃[1,2,3,4,5].

The formation mechanism of O₃ is highly complex, governed by the nonlinear coupling of precursor concentrations, meteorological conditions, and photochemical reactions. This uncertainty renders research spanning from macroscopic ozone variations to microscopic chemical mechanisms all of significant value[6,7,8,9,10,11,12,13,14,15,16]. Ozone concentration variations manifest across multiple scales, including short-term meteorological perturbations, seasonal cycles, and long-term trends. Current research has established an analytical framework covering multiple temporal scales of ozone, encompassing diverse approaches such as time series decomposition, interpretable machine learning, spatial modeling, multifractal cross-correlation analysis, and climate driving factor attribution. He et al. (2025) utilized STL decomposition to quantify the differentiated responses of long-term declining trends and seasonal/short-term perturbations of ozone in the Yangtze River Delta. Zhong et al. (2025) revealed the differentiated climate drivers of interannual variability in heat-ozone compound extreme events between North China and South China through EOF analysis and climate attribution. Ma et al. (2025) and Wu et al. (2026) deepened the mechanistic understanding of short-term extreme ozone events from the perspectives of OVOC photolysis-driven differences in winter and summer formation pathways and nocturnal low-level jet residual layer entrainment processes, respectively. Sowmya et al. (2025) revealed the regulatory role of boundary layer dynamic processes on diurnal and seasonal ozone variations from the perspectives of ventilation coefficient and multifractal cross-correlation[17,18,19,20,21,22,23,24,25,26]. There are many additional studies focusing on ozone at various scales, and we can greatly advance our understanding and knowledge of ozone from previous research.

However, we consider that certain gaps remain in the existing body of work. We suggest that, based on the close chemical linkages between ozone and its precursors, it is worthwhile to attempt treating them as parallel objects for analysis. Furthermore, ozone pollution is essentially a process characterized by time-varying coupling and multi-scale interactions. Nevertheless, existing studies have mostly focused on individual and relatively static relationships, with limited discussion on the dynamic evolution of response relationships across different temporal scales.

This study aims to reveal the evolutionary characteristics of O₃ pollution in East Asia from a multi-temporal scale perspective and its responses to emission reduction policies in China. First, O₃ and its major precursors are incorporated into a unified covariance matrix for EOF decomposition to extract principal component time series reflecting the co-evolution of pollution, thereby overcoming the limitations of traditional single-pollutant analysis. The principal components are then decomposed into ST, SN, and LT temporal scale components through Gaussian low-pass filtering, and statistical tests are applied to each component to identify potential anthropogenic emission reduction signals. Second, two complementary pathways, OLS and CNN-BiLSTM combined models, are employed to explore the cross-scale response relationship β from larger scales to smaller scales and its time-varying characteristics. A PINN model is further incorporated to conduct exploratory analysis of the potential physical features underlying each mode, thereby enhancing the interpretability of cross-scale response relationships. This study provides a new analytical perspective for deepening the understanding of multi-scale coupling mechanisms of O₃ pollution, identifying key response pathways, and evaluating the effectiveness of emission reduction policies across different temporal scales.

2. Materials and Methods

2.1. Data

This study focuses on East Asia (105° to 125°E, 20° to 42.5°N) as the study area. The primary datasets include the following: (1) Three-dimensional field data (0.75°) of O₃, NO₂, HCHO, and isoprene. Previous studies have supported the reliability of these species as representative components of the ozone pollution system[27,28,29,30,31,32,33,34]; (2) Emission data (0.1°) of NOx and VOCs, categorized into four sectors: agriculture, industry, residential, and transportation, along with natural source emission data (0.1°). Both datasets (1) and (2) are derived from the ECMWF EAC4 reanalysis data and global emission inventory data, with the vertical level set at the near-surface 1000 hPa. (3) Six meteorological variables closely associated with ozone: 2-meter temperature, relative humidity, wind speed (WS, synthesized from u and v components), boundary layer height (BLH), daily mean precipitation, and daily surface net upward thermal radiation[35,36,37,38,39]. These data are sourced from the ECMWF ERA5 reanalysis dataset (0.1°). To ensure consistency with the four chemical species for subsequent analysis, area-weighted averaging was applied to regrid meteorological and natural source emission inventory data to 0.75° resolution. The temporal coverage of all datasets spans from January 1, 2017 to December 31, 2025 (daily). Grid-point-wise Z-score standardization was applied to each pollutant, where ${\overline{\mathrm{X}}}_{\mathrm{i},\mathrm{j}}$ and σ_i,j denote the temporal mean and standard deviation at grid point (i,j), respectively. Standardization eliminates the effects of dimensional differences among variables and spatial concentration gradients. Subsequent results primarily reflect anomalous variations of each variable relative to its own historical background.

$$Z_{i,j}\left(t\right)={\displaystyle \frac{{\mathrm{X}}_{\mathrm{i},\mathrm{j}}\left(\mathrm{t}\right)-{\overline{\mathrm{X}}}_{\mathrm{i},\mathrm{j}}}{{\mathsf{\sigma }}_{\mathrm{i},\mathrm{j}}}}$$

(1)

2.2. Gaussian Filter Decomposition

This study employs a dual-layer Gaussian low-pass filter to decompose the standardized pollutant fields along the temporal dimension. This approach possesses smooth frequency response characteristics and favorable time-domain localization performance. The specific procedure is as follows. (1) A one-dimensional Gaussian kernel with standard deviation σ_SN=10 days is applied along the time axis as a low-pass filter, yielding the smoothed field G_SN(t); (2) A one-dimensional Gaussian kernel with standard deviation σ_LT=183 days is applied as a second low-pass filter, yielding the long-term smoothed field G_LT(t).The short-term component(ST):X_ST(t) = X − G_SN(t), captures weather-scale fluctuations. The seasonal component(SN):X_SN(t) = G_SN(t) − G_LT(t), captures seasonal cycles. The long-term component(LT):X_LT(t) = G_LT(t), captures multi-year trends. ST, SN, and LT constitute a non-orthogonal decomposition, and the variance proportions of the three components do not necessarily sum to 100%. The sample sizes across the three datasets remain consistent. The Gaussian filter performs frequency decomposition rather than partitioning the dataset into separate blocks.

X(t) = X_ST(t) + X_SN(t)+X_LT(t)

(2)

2.3. Joint EOF Decomposition

The formation of O₃ involves complex underlying causes. By incorporating O₃ and its precursor components at each temporal scale after Gaussian filtering into a unified covariance matrix for decomposition, the information is efficiently compressed into a small number of comprehensive modes. While achieving dimensionality reduction, the PC series serves as a composite evolution signal of the ensemble, whose variational characteristics reflect the overall evolution, consistent with the actual behavior of the atmospheric chemical system. The three-dimensional data of each pollutant are reshaped into a “time × space” matrix and then concatenated along the feature dimension to form the joint matrix[40,41,42]:

$${\mathrm{X}}_{\mathrm{s}}=\left[{\mathrm{X}}_{{\mathrm{O}}_{3,}\mathrm{s},}{\mathrm{X}}_{\mathrm{Iso},\mathrm{s},}{\mathrm{X}}_{\mathrm{N}{\mathrm{O}}_2,\mathrm{S},}{\mathrm{X}}_{\mathrm{HCHO},\mathrm{S}}\right]$$

(3)

where s denotes the three temporal scales. The temporal mean is removed from each spatial feature of the joint matrix, and truncated singular value decomposition is applied:

$${\mathrm{X}}_{\mathrm{S}}={\mathrm{PC}}_{\mathrm{s}}{\bullet \mathrm{M}}_{\mathrm{s}}^{\mathrm{T}}$$

(4)

Where PC_s is the principal component time series and${\mathrm{M}}_{\mathrm{s}}^{\mathrm{T}}$ is the corresponding spatial mode. The spatial modes are fixed and can be understood as stationary concentration distribution patterns manifested by the collective activity of ozone and precursors within the study area. Mode 1 (M1) corresponds to the most frequently occurring distribution pattern, and so forth. Some modes possess clear physical interpretations, while others may represent data residuals. The PC series reflects the temporal activity of each mode. For example, $\mathrm{P}{\mathrm{C}}_{{\mathrm{ST}}_{\mathrm{M}1}}$ represents the PC series of M1 at the ST scale, and so forth. The explained variance of the m-th mode relative to the total sum of squares of the joint matrix is:

$$\mathrm{E}{\mathrm{V}}_{\mathrm{m}}={\displaystyle \frac{{\mathrm{S}}_{\mathrm{m}}^2}{{\sum }_{\mathrm{t},\mathrm{k}}{\mathrm{X}}_{\mathrm{t},\mathrm{k}}^2}}$$

(5)

The explained variance reflects the proportion of variance accounted for by the corresponding mode within the entire concentration dataset. The spatial modes from EOF decomposition exhibit alternating positive and negative loading patterns. Positive loading regions (red) should not be simply interpreted as high pollution areas, nor negative loading regions (blue) as low pollution areas. The scientific meaning is derived from the synchronous variation between spatial signs and the corresponding PC values. A grid point is associated with a high ozone concentration component only when the spatial loading and the corresponding PC value share the same sign at a given time. The first two modes at each of the three temporal scales are extracted for subsequent analysis, based on the following rationale. (1) The cumulative explained variances of the first two modes for SN and LT reach 85.40% and 83.48%, respectively, which sufficiently characterize the major spatial variability structures. (2) ST fluctuations are complex and mutually independent, with individual mode variance contributions below 6%. Selecting the first two representative modes that also yield satisfactory deep learning performance reduces uncertainty (SI Section1.2). On this basis, the following analytical indicators are established for the three scales:

(1) For the PC series of each ST mode, the 90th percentile is used as a threshold to identify short-term joint pollution mode anomaly events (adjacent events separated by ≤1 day are merged). Three indicators are computed annually: event frequency (events per year), mean peak intensity (exceedance above threshold), and mean duration (days).

(2) For each SN mode PC series, a harmonic model containing the fundamental frequency and the second harmonic is fitted annually:

PC_SN(t) = a₀ + a₁cos(ωt) + b₁cos(ωt) + a₂cos(2ωt) + b₂cos(2ωt)

(6)

Where ω = 2π/365. The following interannual indicators of the seasonal joint pollution mode are extracted: amplitude $\mathrm{A}=\sqrt{{\mathrm{a}}^2+{\mathrm{b}}^2}$, peak phase $\mathsf{\Phi }=\arctan \left({\displaystyle \frac{{\mathrm{b}}_1}{{\mathrm{a}}_1}}\right)/\mathsf{\omega }$, and mean PC value during the May to September high ozone pollution season.

(3) For the annual mean series of each LT mode PC, the following trend tests are applied. (i) Nonparametric Mann-Kendall test: tests the significance of monotonic trends, reporting the two-sided p-value. (ii) Theil-Sen slope estimate: the median of all pairwise slopes are taken as the robust trend estimator. (iii) Pettitt change point detection: identifies abrupt mean shifts based on the cumulative deviation sequence ($\mathrm{CUSUM}={\displaystyle {\sum }_{\mathrm{i}=1}^{\mathrm{t}}}({\mathrm{PC}}_{\mathrm{i}}-\overline{\mathrm{PC}}))$, with the test statistic K = max|S_t|，and significance computed as P = 2exp(−6K²/(n³ + n²)).

2.4. Multivariate Linear Regression and β Series Calculation (OLS)

To further explore the cross-scale response characteristics of the ozone pollution system, two pathways are adopted. The first pathway employs conventional multivariate linear regression (OLS), establishing dynamic linear responses of target scales to driving scales within a 365-day sliding window. A cross-scale response framework from larger to smaller scales is adopted, whereby long-term trends and seasonal cycles serve as relatively stable background states that influence the occurrence frequency and intensity of more rapidly varying short-term events. This assumption is based on the objective of simplifying model structure and avoiding circular dependencies between fast and slow scales that would lead to research complexity and difficulties in interpretation. Two tasks are formulated:

TaskA: PC_ST(t) = β₀ + β₁PC_SN,M1 + β₂PC_SN,M2 + β₃PC_LT,M1 + β₄PC_LT,M2 + ε

(7)

$$\mathrm{TaskB}:{\mathrm{PC}}_{\mathrm{SN}}\left(\mathrm{t}\right)={\mathsf{\beta }}_0+{\mathsf{\beta }}_1{\mathrm{PC}}_{\mathrm{LT},\mathrm{M}1}+{{\mathsf{\beta }}_2\mathrm{PC}}_{\mathrm{LT},\mathrm{M}2}+\mathsf{\epsilon }$$

(8)

Where PC corresponds to the PC series of each mode. Tasks A and B correspond to the response of ST to LT and SN, and the response of SN to LT, respectively, yielding a total of 12 relationship pairs. Each sliding window requires a minimum of 300 valid sample days. The resulting β series describes the cross-scale dynamic linear response intensity in the sense of conditional correlation. A larger |β| indicates stronger explanatory power of the driving mode over the variation of the target mode. A positive β indicates that the response mode strengthens synchronously when the driving mode intensifies, while a negative β indicates that the response mode weakens when the driving mode intensifies. A positive trend in β indicates that the driving effect strengthens over the years, while a negative trend indicates that the driving effect weakens or reverses direction over the years. All indicators characterize the evolution of relationships rather than the ozone concentration changes resulting from those relationships.

2.5. Cross-Scale Response Modeling

The second pathway employs a CNN+BiLSTM[43,44,45,46,47,48] neural network combination to learn the variational characteristics of all PC series. The model architecture consists of multi-branch one-dimensional CNN and BiLSTM. The CNN contains five parallel convolution branches with kernel widths of 1, 3, 5, 9, and 15 days, respectively, designed to extract local features. The outputs of each branch are processed through GroupNorm and GELU activation before concatenation. The BiLSTM further models the long-range temporal dependencies of multi-scale features. Its output is processed through LayerNorm, residual connections, and a two-layer fully connected network to generate per-timestep predictions. Simultaneously, the CNN produces local response predictions through an independent 1×1 convolution prediction head, which are summed with the BiLSTM output to form a dual-pathway residual output, improving gradient propagation and training stability. For Task A, inputs include PC_SN，PC_ST series along with six meteorological factors. For Task B, inputs include PC_LT series along with six meteorological factors. Additionally, for Task A, first−order differences, third−order differences, 7−day rolling standard deviations, and 7−day rolling deviation measures are computed separately for meteorological variables, PC_SN，PC_ST series as derived features to capture the high-frequency variational characteristics at the ST scale, constituting a specialized scheme (SI Section1.3). The model uses PC data from January 1, 2017 to December 31, 2024 for training and leverages the advantages of neural networks to fit the β series, possessing the capability to explore nonlinear cross-scale relationships (in contrast to the analytical solution of OLS). The results are ultimately compared with the β series from the OLS pathway. The year 2025 serves as an independent validation period for predicting that year’s PC and β series, providing a comprehensive assessment of model reproduction and learning capabilities (SI Section 2.1 and Section 2.2).

2.6. PINN Physical Constraint Design

The multi-scale response relationships involve a total of 12 groups of research objects. The secondary variables from EOF decomposition, spatial modes ${\mathrm{M}}_{\mathrm{s}}^{\mathrm{T}}$ and the corresponding PCs describing each mode, are abstract concepts. It is necessary to make informed judgments and prior assessments regarding the actual physical features that each mode may represent. This study incorporates PINN into the existing neural network architecture, constructing a hybrid attribution framework that integrates data-driven learning with physical knowledge constraints. The PINN embeds six representative fundamental laws of O₃-related atmospheric chemistry and physics directly into the loss function as soft penalty terms. During training, these physics-based regularization terms guide the model to search for optimal solutions only within the physically plausible solution space, effectively constraining model behavior without requiring additional labeled data. The PINN first learns physically consistent residual corrections and generates physically regularized feature representations. These physics-enhanced features[49,50,51,52,53,54] are then fed into the neural network for ensemble learning and attribution analysis (SI Section 3.1 and Section 3.2). The specific constraint mechanisms are as follows:

Feature 1: Arrhenius Temperature Factor. Ozone photochemical production involves a series of thermochemical reactions whose rate constants k(T) follow the Arrhenius law:

$$\mathrm{k}\left(\mathrm{T}\right)=\mathrm{A}.\exp \left(-{\displaystyle \frac{{\mathrm{E}}_{\mathrm{a}}}{\mathrm{R}.\mathrm{T}}}\right)$$

(9)

Where A is the pre−exponential factor, E_a is the activation energy, R is the ideal gas constant, and T is the absolute temperature. This study uses the relative rate factor with respect to a reference temperature T_ref(298.15K):

$${\mathrm{f}}_{\mathrm{Arrhenius}}\left(\mathrm{T}\right)=\exp \left[-{\displaystyle \frac{{\mathrm{E}}_{\mathrm{a}}}{\mathrm{R}}}\left({\displaystyle \frac{1}{\mathrm{T}}}-{\displaystyle \frac{1}{{\mathrm{T}}_{\mathrm{ref}}}}\right)\right]$$

(10)

Feature 2: Formaldehyde-to-NO₂ Ratio (FNR). The sensitivity of ozone production to precursors depends on the chemical regime. Here, ε is a prescribed constant of 10⁻¹² to prevent division by 0, and 50 serves as the upper clipping limit for the FNR ratio to prevent physically unrealistic extreme values.

$$\mathrm{FNR}=\mathrm{clip}\left({\displaystyle \frac{\left|\mathrm{HCHO}\right|}{\left|\mathrm{N}{\mathrm{O}}_2\right|+\mathrm{ℇ}}},0,50\right)$$

(11)

Feature 3: O₃ Photochemical Production Rate. Here, max(radiation,0) represents the photolysis driving force, corresponding to the radiation value at the relevant grid point. f_arrhenius(T) is the thermochemical reaction acceleration factor corresponding to Feature 1. VOCs_proxyrepresents the total precursor emissions at the relevant location, comprising both anthropogenic and biogenic emission values. The net photochemical production rate of tropospheric ozone P(O₃) can be expressed in simplified form as:

$${\mathrm{P}}_{{\mathrm{O}}_3}^{\mathrm{proxy}}=\max \left(\mathrm{radiation},0\right)\times {\mathrm{f}}_{\mathrm{arrhenius}}\left(\mathrm{T}\right)\times \left(\mathrm{VOC}{\mathrm{s}}_{\mathrm{proxy}}+\mathsf{\epsilon }\right)$$

(12)

Feature 4: O₃-NOx Titration Effect. In high-NOx environments, freshly emitted NO directly consumes ozone through titration reactions. The intensity of this process is proportional to emission intensity and inversely proportional to mixing volume and boundary layer height. Here, ε=10⁻¹² also serves as numerical protection to prevent division by zero. The constraint formula is designed as:

$${\mathrm{L}}_{{\mathrm{O}}_3}^{\mathrm{titration}}={\displaystyle \frac{{\mathrm{E}}_{\mathrm{industry}}+{\mathrm{E}}_{\mathrm{residential}}+{\mathrm{E}}_{\mathrm{transport}}+{\mathrm{E}}_{\mathrm{agriculture}}}{{\mathrm{BLH}}_{\mathrm{clipped}}+\mathsf{\epsilon }}}$$

(13)

Feature 5: Boundary Layer Exchange. Here, E_industry + E_residential + E_transport + E_agriculture represents the total emissions from all sources. This constraint characterizes the effective concentration of all emissions within the boundary layer, explaining the general dilution effect of the boundary layer, which influences all precursor concentrations. The pollutant concentration released by emission sources into the boundary layer is approximated as:

$${\mathrm{R}}_{{\displaystyle \frac{\mathrm{E}}{\mathrm{BLH}}}}={\displaystyle \frac{{\mathrm{E}}_{\mathrm{total}}}{\max \left(\mathrm{BLH},50\right)+\mathsf{\epsilon }}}$$

(14)

Feature 6: Ventilation Coefficient. The horizontal transport capacity of atmospheric pollutants depends on wind speed and boundary layer height. The former determines horizontal advective diffusion, and the latter determines the thickness of the air column participating in transport. Here, $\sqrt{{\mathrm{u}}^2+{\mathrm{v}}^2}$ represents the horizontal wind speed scalar (m/s), and max(BLH,50) represents the clipped boundary layer height (m). Their product represents the ventilation coefficient.

$${\mathrm{V}}_{\mathrm{wind}}=\sqrt{{\mathrm{u}}^2+{\mathrm{v}}^2}\times \max \left(\mathrm{BLH},50\right)$$

(15)

3. Results

3.1. Potential Linkages Between Anthropogenic Emission Reductions and Multi-Scale O₃ Pollution System

The Gaussian filtering results for O₃ and its precursors are provided in SI Section1.1. Figure 1(a, c) shows that the two ST modes cumulatively explain only 9.94% of the total variance, indicating that the spatial structure at the ST scale is highly fragmented. Different weather events and regional characteristics activate different spatial patterns, with no single dominant mode. A total of 144 and 187 anomaly events were detected for ST_M1, ST_M2, respectively. The interannual variations in frequency, intensity, and duration are all statistically insignificant (|P|>0.2), suggesting that the intensification of ozone occurrence frequency and intensity at the ST scale has not been sustained. However, it is worth noting that the duration trends of ST_M1, ST_M2(−0.016 days/year and −0.024 days/year), although insignificant (P = 0.208), suggest that the duration of major ST modes is gradually shortening, albeit slowly. Figure 1(b) shows that the ST modes exhibit large interannual fluctuations overall, precluding the detection of monotonic trends. The low explained variance of ST scale modes indicates limited representativeness. The two SN modes cumulatively explain 85.4% of the variance (Figure 1(d, f)). The MK test results for SN_M1, SN_M2 indicate no significant trends in any of the three indicators ((|P|>0.3). The ozone component at the SN scale across eastern China has remained largely stable over the nine−year period. This is broadly consistent with the understanding that O₃ background concentrations are rising while seasonal peaks have remained relatively stable[55,56,57,58,59,60]. The phase shift of SN_M2 is a signal worthy of attention. Although statistically insignificant, if the peak of SN_M2 advances by an average of 1.4 days per year, this may reflect an earlier onset of the ozone season in certain regions. Overall, the activity characteristics of the ozone pollution system at the SN scale also exhibit fluctuating patterns without significant trends (Figure 1(e)). The two LT modes cumulatively explain 83.48% of the variance, indicating that long-term changes are characterized by two distinct spatial mode features. LT_M1(TS=+10.23/yr，MK P=0.048) represents the only marginally significant trend in the entire analysis, while LT_M2(TS=+5.40/yr，MK P=0.602) shows no significant trend. The gradually rising positive-phase trend of the LT_M1 PC series, combined with the uniformly positive loadings of LT_M1 across the entire domain shown in Figure 1, indicates that both share the same sign. This suggests that the ozone component corresponding to LT_M1 is undergoing overall enhancement during the study period. The gradually rising trend of the LT_M2 PC series, combined with the negative spatial loadings of LT_M2 across most regions shown in Figure 1(i), suggests that the ozone component corresponding to LT_M2 is undergoing overall weakening during the study period. Pettitt change point detection identifies change points at 2019 and 2022. Notably, these coincide closely with both the two years of extreme ozone pollution and the two years when ozone control policies were issued during the study period. The aforementioned ozone activity at both the ST and SN scales exhibits no trends but rather fluctuating patterns. Combined with the fact that LT_M2 reflects the only partial upward trend in ozone pollution and LT_M2 reflects the only declining trend in ozone components, as well as the strong correspondence between change point information and policy promulgation periods, we are inclined to consider that LT_M2 has a certain association with anthropogenic emission reduction activities.

3.2. Evolution of Multi-Scale Response Relationships in Ozone Pollution

The multi-scale response relationships involve a total of 12 groups of research objects. Therefore, establishing a prior understanding of the modes at each scale before analysis is necessary. The PINN conducted exploratory analysis of the physical scenarios that each mode may represent, and the model demonstrates satisfactory training performance (SI Section3.3). The two modes at the LT scale point to different primary driving sources. Behind LT_M1,temperature, radiation, and natural sources occupy extremely high SHAP explanatory power, potentially reflecting ozone variations associated with long-term temperature change trends and vegetation emission responses. This is broadly consistent with the statistical results in Section 3.1. LT_M2reflects extremely high SHAP explanatory power from anthropogenic emissions (emission contribution of 76.4%), while the contribution of meteorological elements in this mode is minimal, indicating that LT_M2 has effectively isolated the emission trend. This is also broadly consistent with the statistical results in Section 3.1. The mode characteristics at the SN and ST scales are more complex. From the SHAP perspective, emission and meteorological factors alternately appear among the top-ranked contributors in SN_M1.Agricultural and residential emissions exhibit pronounced seasonal cycles. The winter heating season and the agricultural planting season appear to correspond well with these patterns. Additionally, radiation and temperature rank prominently, reflecting the seasonal variation of photochemical production rates. In SN_M2, the contribution gaps among factors are small, indicating that SN_M2 reflects a seasonal mode governed by multi-factor synergistic effects, with no single factor occupying a dominant position. In ST modes, the top four contributors are all dynamic and diffusion-related meteorological elements. The combined contribution rates of wind(U) and wind(V) reach 32% and 52.2%, respectively, sufficiently demonstrating the importance of wind field transport and boundary layer diffusion as typical meteorological processes for ozone pollution activity at this scale. Emission factors contribute limitedly to ST scale fluctuations. In summary, the ST scale reflects local wind field transport and boundary layer diffusion activity, with emission and natural source contributions being relatively minor. The seasonal scale exhibits strong seasonality, with pronounced characteristics of the monsoon and high-temperature seasons, accompanied by partial contributions from human activities. The long-term scale reflects two opposing types of background concentration changes (Figure 2). The SHAP results and the EOF decomposition reasoning in Section 3.1 provide mutual verification. Although attribution conclusions cannot be drawn solely from SHAP results, based on the above understanding, these scenario classifications possess a certain degree of credibility. The following discussion on β coefficients can also establish a degree of prior knowledge bridging modal language and physical contexts.

The model learning R² values are presented in Table 2. The performance of Task A is generally inferior to that of Task B. As previously discussed, the ST scale fails to extract effective signals despite EOF dimensionality reduction. In Task A, the influence of SN on ST is primarily manifested as a deepening modulation of the secondary seasonal mode on short-term dominant weather events (Slope = 0.323), while the remaining three relationship pairs show no significant trends. In comparison, the response relationships from LT to ST exhibit stronger variation amplitudes (Mean β values are generally larger) and more significant trends. Whether the global warming background signal driven primarily by temperature (LT_M1, temperature contribution 24.5%) or the anthropogenic emission trend signal(LT_M2, combined contribution 76.4%), their negative modulation on short-term weather processes dominated by wind field transport and boundary layer diffusion(ST_M1, combined contribution of wind field and BLH approximately 59%) is deepening (trends of −0.05 and −0.41 yr⁻¹, respectively, P<0.001). As long-term warming and emission patterns continue to evolve, conventional short-term scale activity dominated by wind field and diffusion processes is being attenuated. Concurrently, the positive modulation of these two long-term signals on the secondary short-term mode characterized by zonal wind and vertical motion (ST_M2， combined contribution of wind(U) and vertical wind 52.2%) is strengthening (trends of +0.081 and +0.25 yr⁻¹, respectively, P< 0.001). This suggests that long-term changes may be activating the redistribution of ozone through weather processes such as frontal passages and convective transport (Figure 2, Table 2). Task B reveals a deep restructuring of the seasonal scale driving structure by long-term signals. The driving efficiency of the global warming background signal (LT_M1) on the dominant seasonal feature is declining (LT_M1 to SN_M1,Slope=−2.82), while simultaneously shifting toward the secondary seasonal feature(LT_M1 to SN_M2，Slope=+2.98). Considering that SN_M1 is primarily driven by the seasonal cycles of agricultural and residential emissions (combined 27.0%), with radiation and temperature as secondary drivers, while SN_M2 exhibits a multi-factor synergistic pattern with relatively uniform contributions from 13 factors and industrial emissions ranking slightly first, this opposing evolution implies that sustained background warming is attenuating the classical seasonal ozone pattern dominated by winter heating emission peaks or summer photochemical production peaks, and instead enhancing more complex seasonal response pathways co-modulated by industrial emissions, meteorological conditions, and biogenic sources. Meanwhile, the modulation of the anthropogenic emission trend signal (LT_M2) on the entire seasonal scale is strengthening comprehensively in the positive direction (Slope = +0.5/+4.28, P< 0.001). This indicates that the long-term evolution of emission patterns, particularly the structural changes in agricultural and residential sources, is systematically amplifying its influence on the ozone seasonal cycle (Table 2). Synthesizing the results of both tasks, the most prominent feature during the historical period is that responses generally weaken or remain stable when major modes (M1) serve as targets, while responses generally strengthen when secondary modes (M2) serve as targets. Combined with the physical driving characteristics behind each mode, this pattern carries clear physical implications. O₃ variation patterns traditionally dominated by single primary factors (such as the agricultural emission cycle at the seasonal scale or wind field transport at the short-term scale) are yielding to secondary patterns driven by multi-factor synergy. Emission activities are also involved in this restructuring process. This suggests that the spatiotemporal structure and chemical composition of the future ozone system may deviate from historical characteristics and past experience, and that anthropogenic emission reduction measures may play an increasingly significant role. For ozone forecasting and control strategies, this implies that empirical judgments relying solely on historical dominant mode characteristics may face the risk of systematic bias.

4. Discussion

Among the 14 trend indicators examined in Section 3.1, only LT_M1 reaches the significance level, and the only declining trend in ozone components is found in LT_M2. In contrast, the ST and SN scales exhibit predominantly fluctuating characteristics, with emission reduction effects not manifested at these scales. Change points are detected in LT_M1, LT_M2 (2019 and 2022). M1 corresponds to an upward trend in the ozone component, consistent with the prevailing understanding that ozone background concentrations are rising under the current climate warming context, and the change point years coincide with regional extreme O₃ pollution events dominated by high temperatures[61,62,63,64]. M2 exhibits a declining trend, with change points coinciding with the Three-Year Action Plan for Winning the Blue Sky Defense War (2018) and the Action Plan for Continuous Improvement of Air Quality (2023). The 2019 inflection point also echoes the findings of Wang et al[65]. SHAP provides additional evidence that LT_M2 can be reasonably associated with a portion of ozone component decline attributable to anthropogenic emission reductions. Within the analytical logic of this study, the value of β curves lies in their characterization of the strength variations in coupling relationships between features of the ozone pollution system across different temporal scales. The strengthening or weakening of any single relationship pair cannot be directly translated into changes in ozone concentration. Results indicate that responses generally weaken or remain stable when major modes (M1) serve as targets, while responses generally strengthen when secondary modes (M2) serve as targets. This is consistent with recently observed phenomena including the increasing complexity of short-term ozone pollution events, seasonal extension, and the growing frequency of high-value events during non-traditional pollution months. The occurrence of these phenomena involves the participation of emission activities. The above findings yield two implications. First, attribution analysis of ozone pollution needs to be refined to the “mode” level. Broadly discussing “the influence of long-term trends on short-term events” or “the modulation of extreme events by seasonal signals” may obscure fundamentally different or even opposing modulation directions between modes. Second, the evolution of response relationships itself may constitute a source of system variability equally important as emission trends. For example, even if the long-term concentration rising trend is curbed, short-term extreme events may still increase rather than decrease due to the effects of coupling relationships (LT_M1 to ST_M1, Slope = -0.053，Mean β=-0.03).

This study provides a new perspective for understanding cross-temporal scale coupling relationships in O₃ pollution. The history of O₃ governance in China has progressed from an early focus on coal-combustion type pollutants, to a mid-stage focus on VOCs reduction, and subsequently to a renewed focus on NOx control. Throughout this series of emission reduction practices, the chemical and physical background of O₃ formation has also been continuously changing. Mode A, which frequently appeared and dominated O₃ pollution activity in the study area during the early period, may be replaced by Mode B in the new period due to phased governance measures and uncertain variations in meteorological conditions. This indicates that emission reduction strategies designed based on historically dominant pollution characteristics may not effectively address emerging characteristic modes in the future. Although this study has identified potential background factors for relevant modes based on SHAP, the evidence remains incomplete. Future work should investigate these as potential cross-scale risk amplification pathways through in-depth mechanistic studies.

Finally, several limitations should be noted. First, the historical period sample spans only nine years, and trends at the long-term scale may be insufficient to support definitive conclusions regarding climate-scale long-term changes. Second, this study does not discuss in depth the specific physical mechanisms underlying cross-scale responses between modes. In other words, we have only discussed the possible existence of cross-scale responses without investigating what specific physical process mechanisms the 12 relationship pairs correspond to. Future work advancing in this direction could explore this aspect in greater depth.

5. Conclusions

The only declining trend in O₃ components is identified exclusively at the long-term scale, indicating that the effects of anthropogenic emission reductions are predominantly manifested in long-term trends. O₃ concentration components at the short-term and seasonal scales exhibit fluctuating patterns without significant trends.

This study reveals that future short-term pollution events and seasonal characteristics may no longer conform to the regionally dominant or compositionally characteristic patterns observed during the historical period, but instead gradually shift toward secondary features or different event types. Control priorities need to be dynamically adjusted based on these evolving relationship characteristics.

The evolution of cross-scale response relationships itself may constitute a source of system variability equally important as the trends of individual scales. In policy evaluation, focusing solely on total emission quantities or annual mean concentration trends while neglecting the evolving relationships across different temporal scales may lead to systematic misestimation of the transmission efficiency of event risks.