Identifying the spatiotemporal variations in ozone formation regimes across China from 2005 to 2019 based on polynomial simulation and causality analysis

Ozone formation regimes are closely related to the ratio of volatile organic compounds (VOCs) to NO_x. Different ranges of $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ indicate three formation regimes, including VOC-limited, transitional, and NO_x-limited regimes. Due to the unstable interactions between a diversity of precursors, the range of the transitional regime, which plays a key role in identifying ozone formation regimes, remains unclear. To overcome the uncertainties from single models and the lack of reference data, we employed two models, polynomial simulation and convergent cross-mapping (CCM), to identify the ranges of $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ across China based on ground observations and remote sensing datasets. The ranges of the transitional regime estimated by polynomial simulation and CCM were [1.0, 1.9] and [1.0, 1.8]. Since 2013, the ozone formation regime has changed to the transitional and NO_x-limited regime all over China, indicating that ozone concentrations across China were mainly controlled by NO_x. However, despite the NO₂ concentrations, HCHO concentrations continuously exert a positive influence on ozone concentrations under transitional and NO_x-limited regimes. Under the circumstance of national NO_x reduction policies, the increase in VOCs became the major driver for the soaring ozone pollution across China. For an effective management of ozone pollution across China, the emission reduction in VOCs and NO_x should be equally considered.

1 Introduction

With the significant improvement of PM_2.5 pollution, surface ozone has become a major airborne pollutant across China since 2017 (Li et al., 2019a; Lu et al., 2020). Due to its severe threat to public health even during a short-period exposure, ozone pollution has received growing emphasis from governments and scholars (H. Liu et al., 2018; Xie et al., 2019). In the past several years, spatiotemporal distribution of ozone concentrations (Wu and Xie, 2017; Shen et al., 2019a) and the influence of meteorological conditions (Chen et al., 2019c; Cheng et al., 2019, 2020) and anthropogenic emissions (Chen et al., 2019b; Cheng et al., 2018; Li et al., 2019a, 2020) on ozone concentrations have been massively studied. However, due to the highly complicated ozone formation regime, effective ozone control remains challenging.

Different from PM_2.5, whose main precursors are NO_x, volatile organic compounds (VOCs), and SO₂, the formation and decomposition of ozone are closely related to two types of precursors, VOCs and NO_x. There is a diversity of reactions between VOCs and NO_x under different meteorological conditions and concentration scenarios (Wang et al., 2017). Since VOCs and NO_x can either promote or restrict ozone production, the $\mathrm{VOCs}/{\mathrm{NO}}_x$ ratio is crucial for surface ozone concentrations. However, the thresholds at which $\mathrm{VOCs}/{\mathrm{NO}}_x$ may promote or restrict ozone production remain unclear (Jin et al., 2017; Schroeder et al., 2017). For instance, under a specific $\mathrm{VOCs}/{\mathrm{NO}}_x$ scenario, the reduction in NO_x may conversely increase surface ozone concentrations (Sillman et al., 1990; Kleinman, 1994). Furthermore, given the large variations in meteorological conditions and the ozone level across China, the effects of $\mathrm{VOCs}/{\mathrm{NO}}_x$ on surface ozone concentrations also demonstrate notable spatiotemporal patterns. In this case, a comprehensive understanding of how the variations in VOCs and NO_x could influence ozone concentrations under different $\mathrm{VOCs}/{\mathrm{NO}}_x$ circumstances is crucial for setting effective emission reduction policies accordingly in different regions.

To examine the complicated nonlinear relationship between ozone concentrations and multiple precursors, a large body of studies has been conducted (Duncan et al., 2010; Choi et al., 2012; Pusede and Cohen, 2012; Chang et al., 2016; Jin et al., 2020). Through small-scale experiments, NO₂ and HCHO proved to be effective proxies for NO_x and VOCs (Sillman et al., 1990; Martin et al., 2004). Since NO₂ and HCHO can be monitored using remote sensing data, the two precursors have been increasingly considered in ozone–precursor sensitivity research (Jin et al., 2020; X. Zhang et al., 2020). Cheng et al. (2018) proved that ${\mathrm{NO}}_{\mathrm{2}}/\mathrm{NO}$ presented a good consistence with long-term ozone concentrations in Beijing. However, NO was not an easily recordable precursor based on satellite observations and not applicable in large-scale monitoring. Cheng et al. (2019) suggested that satellite-retrieved $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ was strongly correlated with surface ozone concentrations in Beijing. Different $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ indicates distinct ozone formation regimes, including VOC-limited, transitional, and NO_x-limited regimes. For the VOC-limited (NO_x-saturated) regime, the control of VOC emissions leads to the reduction in organic radicals (RO₂), the RO₂–NO_x reactions and thus ozone concentrations (Milford et al., 1989). In contrast, the decrease in NO_x promotes VOC–CO reaction, leading to the increase in ozone concentration (Kleinman, 1994). For the NO_x-limited regime, the reduction in NO_x slows down NO₂ photolysis, which produces free oxygen atoms for ozone formation and reduces ozone concentrations. The variations in VOCs exert limited influences on ozone concentrations for this regime (Kleinman, 1994). For the transitional (VOC–NO_x mixed) regime, both VOCs and NO_x impose positive influences on ozone concentrations. Since the transitional regime divides VOC-limited and NO_x-limited regimes, the estimation of the transitional regime range plays a key role in identifying different ozone formation regimes.

Duncan et al. (2010) calculated the transitional regime range as [1.0, 2.0] using the Community Multiscale Air Quality Modeling System (CMAQ) model, whose uncertainties may influence the estimation accuracy (Schroeder et al., 2017). Jin et al. (2020) employed a polynomial model and calculated the transitional regime range over US urban areas as [3.2, 4.1] based on decades of remote sensing and ground observation data. However, given the notable difference in meteorological conditions, ozone levels, and the composition of precursors across different countries, whether the transitional regime range extracted in the US is applicable to other countries remains unclear. Furthermore, the polynomial model may ignore the complicated inner interactions between multiple precursors, meteorological factors, and ozone concentrations in the atmospheric environment (Chen et al., 2020) and may lead to large uncertainties. Consequently, ozone–precursor sensitivity, especially the transitional regime range across China, requires further in-depth analysis.

To this end, this research attempts to investigate the spatiotemporal variations in ozone formation regimes across China and identify the transitional regime range of $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ based on the cross-verification of multiple models. Firstly, long-term variations in HCHO and NO₂ across China were analyzed. Next, the datasets of HCHO, NO₂, and ozone were examined using a polynomial model and a causality model, respectively, to reveal the crucial thresholds of $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ that separate the NO_x-limited, VOC-limited, and transitional regimes. Specifically, due to the large area of China and potential spatial variations in ozone formation regimes, we respectively investigated ozone formation regimes in several major regions, including the North China Plain (NCP), Yangtze River Delta (YRD), Pearl River Delta (PRD), and Sichuan Basin (SCB) (the geographical locations of four megacity clusters were shown in Fig. 1), to explore the spatiotemporal variations in ozone formation regimes. Meanwhile, we also compared the ozone formation regimes in urban and rural areas. This research sheds useful light for better modeling the complicated ozone–precursor relationship, understanding the major drivers for enhanced ozone pollution, and implementing specific emission reduction measures to mitigate ozone pollution across China.

Figure 1 The May-to-September mean hourly surface ozone network data from 2014 to 2019. Mean hourly surface ozone concentrations are calculated on the $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ grid. Purple, blue, green, and red outlines indicate the boundaries of North China Plain (NCP), Yangtze River Delta (YRD), Pearl River Delta (PRD), and Sichuan Basin (SCB), respectively.

2 Materials and methods

2.1 Data sources

In this study, Ozone Monitoring Instrument (OMI) $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ datasets were employed for exploring the spatiotemporal variations in HCHO and NO₂ in China and calculating $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$. We connected surface ozone network data to HCHO, NO₂, and $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$, which served as the input data for running third-polynomial model and convergent cross-mapping (CCM). The MODIS land cover product provided the spatial distribution of urban areas, which was employed for identifying urban and rural pixels.

2.1.1 OMI $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$

The Ozone Monitoring Instrument (OMI), on board the Aura satellite, monitors global solar backscatter in the UV–vis domain (270–500 nm). The OMI provides daily global observations and crosses the Equator at 13:38 LT (Levelt et al., 2006). In this study, we employed the daily level-3 gridded OMI HCHO product (OMHCHOd) from the Smith Astrophysical Observatory (SAO) (González Abad et al., 2015). The HCHO vertical columns are the weighted mean values for the $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$ grid. Backscattered solar radiation, ranging from 328.5–356.5 nm, was used for fitting HCHO slant columns. Air mass factors (AMFs) were employed for converting HCHO slant columns to vertical columns (González Abad et al., 2015). The validation report suggested that the error in this product was effectively controlled within 30 % over polluted areas (González Abad et al., 2015) and validated for detecting long-term variations in HCHO columns (Zhu et al., 2017; Shen et al., 2019b). The daily level-3 gridded OMI NO₂ product (OMNO2d), provided by NASA's Goddard Space Flight Center, was utilized in this study (Bucsela et al., 2013; Lamsal et al., 2014). The spatial resolution of OMNO2d is 0.25^∘, and each grid is generated as the weighted average of the corresponding level-2 data pixels (Krotkov et al., 2017). Differential optical absorption spectroscopy (DOAS) was employed for retrieving the NO₂ slant columns, which were successively transformed into tropospheric and stratospheric vertical columns through AMFs (Bucsela et al., 2013). The OMI NO₂ column product agrees well with other satellite products, and its overall uncertainties range from 30 %–60 % (Bucsela et al., 2013; Lamsal et al., 2014). To reduce uncertainties, we only selected those OMI HCHO and NO₂ data that (1) passed quality checks, (2) had a cloud coverage less than 30 %, (3) had a solar zenith angle less than 60^∘, and (4) were not affected by row anomalies for this study (Kroon et al., 2011; Zhu et al., 2014; Krotkov et al., 2017). The May-to-September OMI HCHO and NO₂ products were acquired from NASA's Goddard Earth Sciences Data and Information Services Center (https://disc.gsfc.nasa.gov/, last access: 1 September 2021).

2.1.2 Surface ozone network data

The May-to-September hourly surface ozone concentrations from 2014 to 2019 were obtained from the China Ministry of Ecology and Environment (MEE) (https://quotsoft.net/air/, last access: 1 September 2021). The unit of surface ozone concentrations in this dataset is µg m⁻³. The network had 1633 monitoring stations, which were distributed among 330 cities across China in 2019. We used the observation data from 13:00 to 14:00 LT to match the overpass time of the OMI. This dataset has been employed in many studies to investigate the variations in surface ozone concentrations in China (Li et al., 2019a; Shen et al., 2019a; Lu et al., 2020).

2.1.3 MODIS land cover product

The annual MODIS land cover product (MCD12C1) with a spatial resolution of 0.05^∘ from 2005 to 2019 was employed for extracting urban and rural areas. The urban and water pixels from the International Geosphere–Biosphere Program (IGBP) classification layer were employed for the following processing. The land cover product was generated based on a decision tree algorithm with boosting techniques, and its overall accuracy was about 75 % (Palmer et al., 2015; Bajocco et al., 2018). The MCD12C1 product was obtained from NASA's Earth System Data and Information System (https://earthdata.nasa.gov/, last access: 1 September 2021).

2.1.4 Data pre-processing

Due to the different spatial resolution of OMI HCHO, OMI NO₂, and MCD12C1, a bilinear interpolation method was used for resampling all abovementioned products to the same spatial size ($\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$). Meanwhile, we also calculated mean hourly surface ozone concentrations on the $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ grid (Fig. 1).

2.2 Methods

Chemical transport models, such as the global chemical transport model (GEOS-Chem) (Jin et al., 2017; Li et al., 2019a) and the Community Multiscale Air Quality Modeling System (CMAQ) (Duncan et al., 2010), have been frequently employed for exploring the ozone sensitivity to VOCs and NO_x. However, there were large biases in estimating the range of the transitional regime based on chemical transport models (Jin et al., 2017, 2020) due to the uncertainties in the emission inventory and the setting of model parameters. Employing observation data alone could effectively overcome these limitations, and the relationships between ozone and its precursors were fitted using linear and polynomial models (Sun et al., 2018; Jin et al., 2020). Meanwhile, convergent cross-mapping (CCM) (Sugihara et al., 2012), as a robust causality analysis model, has been widely employed for quantifying the influences of meteorological factors on surface ozone and PM_2.5 concentrations (Chen et al., 2018, 2019c, 2020). It is a promising tool for investigating the relationships between ozone and its precursors. To increase the reliability of the estimated range of the transitional regime, both the polynomial model and CCM were employed in this research. We employed the third-order polynomial model for fitting surface ozone concentrations to the indicator of $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$. CCM was employed for quantifying the influences of HCHO and NO₂ on surface ozone concentrations, and the Wilcoxon test (Gehan, 1965) was used for examining whether the differences between the causality of HCHO on ozone and NO₂ on ozone at different ranges of $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ was significant. Since the algorithms of the two models are quite different, their cross-verification provides useful reference for their reliability. Meanwhile, the Mann–Kendall (M–K) test (Kendall, 1970) was employed for exploring the spatiotemporal variations in HCHO, NO₂, and ozone formation regimes in China. Furthermore, we extracted all urban and rural areas in China and compared the differences in ozone formation regimes over these two types of areas. The workflow of the models employed in this study is shown in Fig. 2.

Figure 2 The workflow of the polynomial simulation and the causality analysis.

2.2.1 Estimating the transitional range of the ozone formation regime using polynomial simulation

HCHO and NO₂ are considered to be proxies for VOCs and NO_x, respectively. $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$, as an effective indicator, has been widely employed for determining ozone formation regimes (Duncan et al., 2010; Jin and Holloway, 2015; Jin et al., 2017, 2020; Cheng et al., 2019). Pusede and Cohen (2012) suggested that ozone exceedance probability (OEP) was an effective indicator to interpret the ozone sensitivity to its precursors. The indicator is defined as the proportion of non-attainment events (surface ozone concentrations exceeding 200 µg m⁻³) in total events at a given range of $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$:

$$\begin{array}{}\text{(1)}& \mathrm{OEP}={\displaystyle \frac{{\text{Events}}_{\text{non-attainment}}}{{\text{Events}}_{\text{attainment}}+{\text{Events}}_{\text{non-attainment}}}},\end{array}$$

where Events_attainment and Events_{non-attainment} denote the attainment and non-attainment events, respectively (Pusede and Cohen, 2012; Jin et al., 2020).

In this study, we used a third-order polynomial model (Jin et al., 2020) to explore the quantitative relationships between $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ and ozone exceedance probability. There were 174 868 paired observations of surface ozone concentrations and $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ from 2014 to 2019. The peak of fitting curve highlights the turning point of VOC-limited and NO_x-limited regimes (Jin et al., 2020). The range of $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$, which corresponded to the top 10 % ozone exceedance probability, was defined as the transitional regime. Since we aimed to apply a global model to determine the transitional range, it was necessary to examine whether the surface ozone concentrations in China were of spatially stratified heterogeneity (SSH), as suggested by Wang et al. (2016). We employed the geographical detector (Wang et al., 2010) to measure the SSH of surface ozone concentrations. The geographical detector calculates the q statistic to quantify SSH, and the equation is summarized as follows:

$$\begin{array}{}\text{(2)}& q=\mathrm{1}-{\displaystyle \frac{{\sum }_{h=\mathrm{1}}^L{N}_h{\mathit{\sigma }}_h^{\mathrm{2}}}{N{\mathit{\sigma }}^{\mathrm{2}}}},\end{array}$$

where N and σ² denote the number of samples and the variance of population, and h is the number of stratifications. The range of the q statistic is [0, 1]. The larger the q statistic is, the stronger the SSH is. In this study, the boundaries of four megacity clusters served as strata. If the SSH is detected based on the abovementioned stratification, we could apply the polynomial model in each strata separately.

2.2.2 Estimating the transitional range of the ozone formation regime using convergent cross-mapping

We also employed a causality model named convergent cross-mapping (CCM) (Sugihara et al., 2012), which could reduce the influences of other factors such as meteorological conditions (Chen et al., 2019c, 2020), to extract the causal influences of HCHO and NO₂ on surface ozone concentrations. Thanks to its capability of detecting weak coupling, CCM is advantageous for reliably comparing the influences of different meteorological factors on surface ozone concentrations (Chen et al., 2020). Therefore, we employed CCM to compare the sensitivity of ozone to HCHO and NO₂ at different ranges of $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$. CCM utilizes convergent maps to demonstrate the bidirectional coupling between the time series of two variables. A convergent curve indicates that one variable imposes influences on the other variable, whilst a non-convergent curve denotes no causality between two variables. The main idea of CCM is summarized as follows. Firstly, CCM defines {X} and {Y} as the temporal variations in two variables X and Y. {X} generates the shadow manifold M_X. Following this, the location of the lagged-coordinate vector on M_X, x(t) is determined, and then E+1 nearest neighboring points of x(t) are extracted. Finally, the cross-mapped estimate of Y(t), Y(t)|M_X is calculated as follows:

$$\begin{array}{}\text{(3)}& Y\left(t\right)|M_X=\sum _{i=\mathrm{1}}^{E+\mathrm{1}}{\mathit{\omega }}_{i}Y(t_i),\end{array}$$

where ω_i stands for a weight calculated based on the distance between X(t) and its ith nearest neighboring point. Y(t_i) stands for the contemporaneous value of Y. CCM calculates cross-map skill (ρ value), which explains the quantitative relationships. Number of dimensions for the attractor reconstruction (E), time lag (τ), and number of nearest neighbors to use for prediction (b) are required parameters for CCM. According to previous studies (Chen et al., 2019c, 2020), E, τ, and b were set as 3, 2, and 4, respectively. Since the existence of missing values imposes negative impacts on CCM results, only the consecutive time series were retained for this research. There were 1660 observation records of HCHO time series, NO₂ time series, and corresponding surface ozone time series. CCM was implemented using the “pyEDM” package in Python. The Wilcoxon test (Gehan, 1965) was used to examine whether the differences in ρ values between HCHO and NO₂ were significant at the given $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$. No significant difference was regarded as the transitional regime, while significant difference indicated the VOC-limited or NO_x-limited regime.

2.2.3 Trend analysis

The Mann–Kendall (M–K) (Kendall, 1970) test, which has been used in recent studies on HCHO and NO₂ (Cheng et al., 2019; Wang et al., 2019; Zeb et al., 2019), was employed to estimate the significance of trends. The M–K test is capable of processing samples with random distributions and mitigating the effects of outliers. The Z value is calculated as follows:

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }Z=\left[\begin{array}{c}\frac{S-\mathrm{1}}{\sqrt{\text{Var}\left(S\right)}}(S>\mathrm{0})\\ \frac{S+\mathrm{1}}{\sqrt{\text{Var}\left(S\right)}}(S<\mathrm{0})\end{array}\right],\end{array}$$

where S denotes the statistic to be tested, and Var(S) stands for the variance of S. The sign and absolute value of Z indicate the direction and significance of trends, respectively. Specifically, the positive and negative values of Z indicate the upward and downward trend; 1.28, 1.64, and 2.32 are the threshold values of $\left|Z\right|$, indicating that the trends of samples pass the tests at 90 %, 95 %, and 99 %, respectively.

2.2.4 Comparison of ozone formation regimes in urban and rural areas in China

To compare the differences in ozone formation regimes in urban and rural areas in China, the key step is to extract urban and rural pixels, respectively. Urban pixels were used for buffer analysis (Imhoff et al., 2010) to identify rural pixels. Following Peng et al. (2018), two buffers were set for urban pixels to extract candidate rural pixels (Fig. 3). We set the size of each buffer as 27.75 km, which was close to the size of the $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ grid ($\mathrm{27.75}\mathrm{km}\approx \mathrm{0.25}{}^{\circ }$). The first and second buffers were determined as the urban fringes and candidate rural areas, respectively. Water pixels were firstly removed from candidate rural areas to avoid following uncertainties. Consequently, rural areas were regarded as buffers of 27.75–55.50 km surrounding urban areas. The use of two buffers not only assisted a complete separation of the urban and rural areas but also minimized the uncertainties in meteorological conditions (Yao et al., 2019).

Figure 3 The geographical locations of urban area, urban fringe, and rural area.

3 Results

3.1 Spatial and temporal variations in HCHO and NO₂

Given the national Clean Air Action implemented in 2013, we set this year as a break point to explore the spatial and temporal variations in HCHO and NO₂ in 2005–2012 and 2013–2019, respectively. Figure 4 shows the spatial distribution of HCHO in the two periods. The mean HCHO values during the period of 2005–2012 and 2013–2019 were 4.335×10¹⁵ and 4.845×10¹⁵ $\mathrm{molec}.{\mathrm{cm}}^{-\mathrm{2}}$, characterized by a 12 % increase. Both periods presented an increasing trend of HCHO, and the averaged values during the two periods were 0.164×10¹⁵ and 0.213×10¹⁵ $\mathrm{molec}.{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (Fig. 5). A faster increasing trend was detected during the period of 2013–2019. The variation trend of HCHO agreed well with previous studies (Jin and Holloway, 2015; Shen et al., 2019b). We also calculated the overall linear trends of HCHO in four megacity clusters from 2005 to 2019 (Fig. 6). The largest and smallest increasing trends were shown in the NCP and SCB, with a mean value of 0.136×10¹⁵ and 0.046×10¹⁵ $\mathrm{molec}.{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. The increasing trends of the YRD and PRD were 0.066 and 0.058 $\mathrm{molec}.{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, respectively. Meanwhile, reversed trends were detected for NO₂ during the two periods (Fig. 5), which was consistent with previous studies (Jin and Holloway, 2015; Li et al., 2019a). From 2005 to 2012, the averaged NO₂ was 2.027×10¹⁵ $\mathrm{molec}.{\mathrm{cm}}^{-\mathrm{2}}$, and the annual mean increasing trend was 0.098×10¹⁵ $\mathrm{molec}.{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. Thanks to the implementation of the Clean Air Action, the averaged NO₂ was reduced to 1.900×10¹⁵ $\mathrm{molec}.{\mathrm{cm}}^{-\mathrm{2}}$, with a decreasing trend of $-\mathrm{0.029}\times {\mathrm{10}}^{\mathrm{15}}$ $\mathrm{molec}.{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ from 2013 to 2019. Except for the SCB, all other megacity clusters presented significant downward trends of NO₂ from 2005 to 2019. Amongst these megacity clusters, NO₂ in the YRD demonstrated the largest decreasing trend of 0.104×10¹⁵ $\mathrm{molec}.{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. NO₂ in the NCP and PRD decreased by 0.010×10¹⁵ and 0.092×10¹⁵ $\mathrm{molec}.{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, respectively. A slightly increasing trend of 0.012×10¹⁵ $\mathrm{molec}.{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ was detected in the SCB (Fig. 7).

Figure 4 May-to-September averaged HCHO and NO₂ across China during the period of 2005–2012 and 2013–2019.

Figure 5 The linear trends of May-to-September HCHO and NO₂ across China during the period of 2005–2012 and 2013–2019.

Figure 6 The time series of HCHO columns in the four megacity clusters from 2005 to 2019. Black lines indicate the linear trend of HCHO columns. Red, white, and blue areas stand for VOC-limited, transitional, and NO_x-limited regimes, respectively.

Figure 7 The time series of NO₂ columns in the four megacity clusters from 2005 to 2019. Black lines indicate the linear trend of NO₂ columns. Red, white, and blue areas stand for VOC-limited, transitional, and NO_x-limited regimes, respectively.

3.2 Transitional range of the ozone formation regime

According to $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$, we divided the paired observations into 200 bins for the whole country, and the ozone exceedance probability was calculated for each bin. The third-order polynomial was employed for fitting ozone exceedance probability to $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$. As shown in Fig. 8a, the peak of the fitting curve was 1.4, and the vertical shaded area indicated that the transitional regime over China ranged from 1.0 to 1.9. We employed a geographical detector to examine the SSH of annual May-to-September mean surface ozone concentrations in China. As shown in Table 1, all the q statistics from 2014 to 2019 were greater than zero, which indicated that the surface ozone concentrations in China were of SSH. As suggested by Chen et al. (2020), meteorological factors including temperature, humidity, and sunshine duration imposed great impacts on surface ozone concentration. Moreover, the composition of ozone precursors was closely related to ozone levels (Cheng et al., 2019). Both the meteorological conditions and ozone precursors contributed to the SSH of surface ozone concentrations across China. Therefore, in addition to the regime range extracted at the national scale, we also examined the range of ozone formation regimes in four major megacity clusters. The paired observations of these megacity clusters were divided into 100 bins. The range of the transitional regime for the NCP, YRD, PRD, and SCB was [1.2, 2.1], [1.0, 1.9], [0.9, 1.8], and [1.1, 2.0], respectively, which was generally consistent with the range at the national scale. The small differences between four megacity clusters across China suggested that the range of the transitional regime at the national scale [1.0, 1.9] can be employed to regional- or local-scale research if small-scale data and investigation were not available.

Figure 8 (a) Fitting ozone exceedance probability to $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ through the third-order polynomial model. The curve indicates the fitting result of the third-order polynomial. The vertical line denotes the maximum of the curve, and the shaded area represents the top 10 % ozone exceedance probability. (b) The cross-map skill of HCHO and NO₂ on surface ozone (the skill of using HCHO and NO₂ for predicting surface ozone concentrations) at different ranges of $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$. The symbols and texts above the bars are the results of the Wilcoxon test. ${}^{***}$ and ${}^{**}$ indicate that the difference was significant at the p=0.01 and 0.05 confidence level, respectively. NS suggests non-significant differences.

Table 1 The q statistic and p value calculated by the geographical detector, which indicate the SSH of annual May-to-September mean surface ozone concentrations in China. ^*, ${}^{**}$, and ${}^{***}$ of the p value indicate statistical significance at the α=0.05^*, 0.01${}^{**}$, and <0.001${}^{***}$ level, respectively.

Statistical bootstrapping was used for estimating the uncertainty in the fitting model. Specifically, we iteratively extracted 50 randomly selected subsets from the paired observations to run the model, and the uncertainty was defined as 2 standard deviations from the peak of the fitting curve. The uncertainty for the third-polynomial model was 0.4, indicating a significant nonlinear relationship between ozone exceedance probability and $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$.

Due to the limited data used for running CCM, we set the bin size of $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ as 0.2 for collecting sufficient ρ values to conduct the Wilcoxon test. As shown in Fig. 8b, there was no significant difference between ρ of HCHO and NO₂ when $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ ranged from 0.9 to 1.9, which indirectly defined the range of the transitional regime. For $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}<\mathrm{0.9}$, ρ of HCHO was notably higher than that of NO₂, and this range was regarded as the VOC-limited regime. Similarly, $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}>\mathrm{1.9}$ suggested the NO_x-limited regime. Through the cross-verification, it was an important finding that the range of the transitional ozone formation regime estimated using the third-order polynomial model and CCM was highly close, indicating the reliability of the extracted range.

3.3 Ozone formation regimes in China

NO₂ demonstrated a significant downward trend since 2013, while HCHO kept the increasing trend during the entire study period. Consequently, $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ increased in a majority of regions across China. Specifically, the annually increasing trend of $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ in the NCP, YRD, and PRD was 0.035, 0.023, and 0.034 yr⁻¹, respectively. Meanwhile, there were no significant trends in the SCB during this period (Fig. 9). The variations in $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ indicated the shrinkage of the VOC-limited regime and the expansion of the transitional and NO_x-limited regimes. Since the range of the transitional regime estimated by the third-order polynomial model and CCM was very close, and the former included more reliable observation data, [1.0, 1.9] was employed for identifying different ozone formation regimes. In 2005, areas with the VOC-limited regime were concentrated in the NCP, YRD, and PRD. The proportions of areas with the VOC-limited regime in the NCP, YRD, and PRD were 26 %, 16 %, and 6 %, respectively. Areas with the transitional regime were mainly distributed in the marginal regions of those megacity clusters and scatteredly distributed in the SCB. Areas with the transitional regime occupied 60 %, 50 %, 14 %, and 20 % in the NCP, YRD, PRD, and SCB. The NO_x-limited regime dominated other areas (Fig. 10a). In 2019, areas with the VOC-limited regime decreased significantly; this regime was simply found in the fringe areas of the NCP and YRD. The proportion of the VOC-limited regime in the NCP and YRD was 2 % and 9 %, respectively. The transitional regime was widely distributed throughout the NCP, YRD, and SCB and occupied 71 %, 56 %, and 36 % of the total areas. The NO_x-limited regime still spread over a majority of China (Fig. 10a). We calculated the annual mean ρ of HCHO and NO₂ over those megacity clusters from 2014 to 2019 (Fig. 10b). For all megacity clusters, the ρ of NO₂ was higher than HCHO, indicating that NO₂ was the dominant factor for surface ozone concentrations. Both models suggested that NO₂ played a more important role in affecting surface ozone concentrations than HCHO. In the past several years, NO_x-oriented emission reduction has been conducted across China, leading to the continuous decrease in NO_x concentrations. Since both VOCs and NO_x imposed positive influences on surface ozone concentrations under the transitional and NO_x-limited ozone formation regime, the upward trend of HCHO across China might explain recent soaring ozone concentrations across China (Shen et al., 2019a; Lu et al., 2020). It is noted that the difference between the ρ of NO₂ and HCHO decreased notably in the NCP and YRD. This may be attributed to the following reason. The NCP and YRD are the regions that received severe PM_2.5 pollution, and strict NO_x reduction policies have been conducted since 2013. With the remarkably reduced NO₂ concentrations, the variations in HCHO concentrations plays an increasingly important role in affecting ozone concentrations in the NCP and YRD. The reduction in VOC emissions is key for an effective management of surface ozone pollution in the NCP and YRD.

Figure 9 The time series of $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ in the four megacity clusters from 2005 to 2019. Black lines indicate the linear trend of $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$. Red, green, and blue dots stand for VOC-limited, transitional, and NO_x-limited regimes, respectively.

Figure 10 (a) The spatial distribution of $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ across China in 2005 and 2019. The boundaries of the NCP, YRD, PRD, and SCB are denoted with the purple, blue, yellow, and red bold lines. Red, green, and blue stand for VOC-limited, transitional, and NO_x-limited regimes. (b) The annual mean cross-map skill (ρ value) of four megacity clusters. The red and blue shadow areas indicate the standard deviations.

3.4 Variations in ozone formation regimes in urban and rural areas

Previous studies suggested that the differences in ozone formation regimes existed between urban and rural areas (Tong et al., 2017; Y. Liu et al., 2018; Cheng et al., 2019). We extracted HCHO and NO₂ columns in urban and rural pixels in those megacity clusters and calculated the annually averaged $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ (Fig. 11). For the NCP, $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ in urban areas was higher than 1.0 since 2015, indicating a transformation from the VOC-limited to the transitional regime. The increase in $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ was attributed to the reversed variation trends of HCHO and NO₂. The rising HCHO resulted from the increase in anthropogenic emissions and biogenic volatile organic compounds (BVOCs) (Shen et al., 2019b; Wang et al., 2020), while the implementation of the Clean Air Action imposed notable influences on the decrease in NO₂ (Chen et al., 2019a). $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ in rural areas was in the range of [1.0, 1.9], indicating that rural areas were occupied by the transitional regime from 2005 to 2019. For the YRD, which was occupied by the transitional regime, no variation in ozone formation regime was found in urban areas. In rural areas, $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ temporally exceeded the threshold of 1.9 from 2016 to 2018, indicating that the ozone formation regime changed from transitional to NO_x-limited. This phenomenon was attributed to the slight decline in HCHO, which might be attributed to the restrictions on crop residue burning in this area (Zhuang et al., 2018; Shen et al., 2019b). Due to the large differences in NO₂ concentrations, the urban and rural areas in the PRD were dominated by the transitional regime and NO_x-limited regime. For the SCB, $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ in both urban and rural areas fluctuated around the threshold value of 1.9, and no significant difference between urban and rural areas was found.

Figure 11 The temporal variations in $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ from 2005 to 2019 in the NCP, YRD, PRD, and SCB. The two dashed red lines indicate the threshold values of 1.0 and 1.9, which separate the NO_x-limited, transitional, and VOC-limited regime.

4 Discussion

This research employed CCM and a third-order polynomial model to estimate the transitional regime of ozone formation across China, and the calculated range of $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ was [0.9, 1.9] and [1.0, 1.9], respectively. Our findings were generally consistent with previous studies. For the US, Duncan et al. (2010) and Choi et al. (2012) employed the OMI and GOME-2 data, whose 0.25^∘ resolution was close to this research, and calculated the range of the transitional regime as [1.0, 2.0]. The similar range of the transitional regime in the US and China further proved the reliability of the calculated range [1.0, 1.9] at a national scale. On the other hand, the range of the transitional regime can vary significantly across regions (Schroeder et al., 2017; Jin et al., 2020). Sun et al. (2018) employed station-based data and calculated the range of the transitional regime in Anhui Province, China, as [1.3, 2.8], which was notably higher than the range across China. Jin et al. (2020) calculated the range of the transitional regime in several major regions in the US using the QA4ECV dataset, whose spatial resolution was 0.125^∘, and the output [3.2, 4.1] was much larger than the averaged range of the transitional regime across the US. One reason could be the severe ozone pollution in megacities, leading to different ranges of the transitional regime. Meanwhile, the calculated range of the transitional regime is closely related to the spatial resolution of employed HCHO and NO₂ data, and high-resolution data are more advantageous in extracting the sensitivity of ozone concentrations to precursors at the local scale (Martin et al., 2004; Jin et al., 2017, 2020). In addition to the generally consistent outputs, some advances of this research are listed as follows. First, only a few parameters are required for the polynomial model and CCM, which effectively reduced the uncertainties in model setting. Second, considering the differences between model and satellite-retrieved datasets (Jin et al., 2020), only observation data were employed in this research, which reduced potential data inconsistences and uncertainties. Most importantly, given the lack of actual reference data, this research employed two different models to examine ozone formation regimes, and the close outputs further proved the reliability of this research.

Despite a generally reliable output, some uncertainties exist. First, the accuracy of the estimated range of the transitional regime might be influenced by the scaling biases between station-based observations of surface ozone and space-based HCHO and NO₂. Since ozone monitoring stations are mainly distributed in urban areas, and a $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ grid might cover both the urban and rural areas, the surface ozone concentrations of a grid may be overestimated. Second, the uncertainties in OMI HCHO and NO₂ datasets might impose negative influences on the estimation of the transitional regime range (Duncan et al., 2010; Jin et al., 2017, 2020; Schroeder et al., 2017). On one hand, errors exist in the retrieval of HCHO and NO₂ vertical columns. On the other hand, vertical mixing was not homogeneous, weakening the capability of using HCHO and NO₂ vertical columns to explore the near-surface ozone–precursor sensitivity. Therefore, future improvement of earth observation techniques and the spatiotemporal resolution of HCHO and NO₂ products can further enhance the accuracy of the estimated range of the transitional regime. In general, according to the cross-verification and comparison with previous studies, [1.0, 1.9] from this research is a reliable range for the transitional ozone formation regime across China and can be used as an approximate criterion to follow when implementing national emission reduction policies. On the other hand, given the potential variations in transitional regimes in different regions, when conducting small-scale research, the range of [1.0, 1.9] may be adapted accordingly based on local data.

Previous studies on the range of ozone formation regimes were mainly conducted using statistical models or chemical transport models. For this research, we employed both a statistical and a causality model to cross-verify the range of the transitional regimes. Despite a relatively high fitting accuracy in terms of uncertainties, the findings from these studies could not be effectively compared or interpreted due to the lack of reliable reference data. To this end, as well as numerical models, lab experiments should also be considered to extract a more precise description of the ozone–precursor relationship. With the rapid development of atmospheric science, smog chambers have been increasingly employed to investigate complicated interactions between multiple precursors. By setting specific meteorological conditions (e.g., temperature and humidity) and gradually adjusting the proportion of different precursors, how the proportion of NO₂ and HCHO affects the ozone formation regime can be better explained in a theoretical environment. With more reliable experimental reference data, the model-based analysis on the range of the transitional regime at the local, regional, and national scale can be further improved accordingly.

According to the temporal variations in OMI NO₂ concentrations across China, a notable decreasing trend was observed in three major megacity clusters: NCP, YRD, and PRD. These regions were heavily polluted by PM_2.5, and the notable decrease in NO₂ was mainly attributed to the national Clean Air Action since 2013 (Zheng et al., 2018), which aimed to reduce PM_2.5 concentrations by cutting NO_x emissions. The influence of the Clean Air Action on the reduction in PM_2.5 concentrations and NO_x has been investigated by previous studies. Zheng et al. (2018) employed index decomposition analysis to quantify the contribution of the Clean Air Action and suggested that the decreasing rate of NO_x significantly accelerated since 2013. Moreover, Y. Zhang et al. (2020) employed a random forest algorithm to remove the effects of meteorological conditions and evaluated the impacts of the Clean Air Action. The results demonstrated that the deweathered NO₂ concentrations in winter 2007 and 2017 were 70.3 and 59.1 µg m⁻³, with a decreasing rate of 16 %. Conversely, HCHO concentrations during this period increased remarkably across China due to the combined effects of anthropogenic and biogenetic emissions (Shen et al., 2019b; Wang et al., 2020). The distinct temporal variations in NO₂ and HCHO led to the increase in $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$ and the increase in transitional areas and NO_x-limited regime areas. From 2013–2019, all these regions were dominated by the transitional or NO_x-limited regimes. Attributed to the long-term variation in formation regimes, a more complicated and fragmented spatial pattern was observed across China. Consequently, for an effective control of ozone pollution, the emission reduction in both NO_x and VOCs is required. Especially for the NCP and YRD, where the NO_x has been remarkably reduced, effective approaches for controlling VOC emissions are essential for preventing ozone pollution. This finding was consistent with previous studies (Li et al., 2019b), which recommended the simultaneous reduction in NO_x and VOCs for mitigating the composite airborne pollution in China. Admittedly, compared with NO_x reduction, the VOC reduction is more complicated, and the output of anthropogenic VOC reduction is more unpredictable. In this case, reducing biogenic VOC emissions can also be a potential solution. VOCs emitted by vegetation take up to 50 % of total VOCs in the atmospheric environment, especially in summer. The key factor that may cause enhanced biogenic emissions is summertime high temperature (Chen et al., 2020). Therefore, such projects as wind corridors or contingent artificial precipitation, which can effectively reduce urban heat effects, should be implemented properly to avoid summertime heat waves and successive ozone pollution (e.g., summer, 2017).

The large spatial variations in $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$, especially the rapid increase in transitional regime areas across China, indicate that a unified NO_x–VOC reduction strategy is not feasible for the entire country. Instead, to effectively reduce ozone concentrations, the specific proportion of NO_x and VOC reduction should be carefully set according to local $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$. Meanwhile, due to the large differences in vehicle and industrial emissions (Cheng et al., 2019), the concentration of NO_x is notably higher in urban areas. Therefore, the further reduction in NO_x emissions exerts a stronger influence on ozone reduction in rural areas compared to urban areas.

5 Conclusions

To better understand the spatiotemporal variations in ozone formation regimes across China, we employed the third-order polynomial model and CCM to estimate the range of the transitional regime from 2005 to 2019, the results of which were [1.0, 1.9] and [0.9, 1.9], respectively. The close outputs from two distinct models proved the reliability of the extracted range. At the regional scale, we also investigated the range of the transitional regime in four megacity clusters and found that the range in the NCP, YRD, PRD, and SCB demonstrated limited differences and was generally consistent with the range at the national scale. The reverse trends of HCHO and NO₂ led to the increase in $\mathrm{HCHO}/{\mathrm{NO}}_{\mathrm{2}}$, indicating that China was dominated by the transitional and NO_x-limited regimes in recent years. We also found that the ρ of NO₂ was higher than HCHO in all megacities, suggesting that the reduction in NO_x emissions would become more effective in controlling surface ozone concentrations. Meanwhile, given the rising VOC emissions, the simultaneous reduction in NO_x and VOCs would be more effective than the sole reduction in NO_x in mitigating ozone pollution. Finally, the comparison of ozone regimes in urban and rural areas suggested that the reduction in NO_x emissions would impose stronger impacts on the control of ozone pollution in rural areas.