Identifying the spatiotemporal variations of 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 VOCs to NO x . Different ranges of HCHO/NO 2 indicate three formation regimes, including VOCs-limited, transitional and NO x -limited regimes. Due to the unstable interactions between a diversity of precursors, the range of transitional regime, which plays a key role in identifying ozone formation 15 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 HCHO/NO 2 across China based on ground observations and remote sensing datasets. The ranges of transitional regime estimated by polynomial simulation and CCM were [1.0, 1.9] and [1.0, 1.8]. Since 2013, ozone formation regime has changed to the transitional and NO x -limited regime all over China, indicating ozone concentrations across China were mainly controlled by NO x . However, 20 despite the NO 2 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 of VOCs became the major driver for the soaring ozone pollution across China. For an effective management of ozone pollution across China, the emission-reduction of VOCs and NO x should be equally considered.


Introduction 25
With the significant improvement of PM2.5 pollution, surface ozone has become a major airborne pollutant across China since 2017 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 (Liu et al., 2018a;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., 2019cb;Cheng et al., 2019;Chen et al., 2020) and anthropogenic 30 emissions (Chen et al., 2019b;Cheng et al., 2018;Li et al., 2019a;Li et al., 2020) on ozone concentrations have been 2 massively studied. However, due to the highly complicated ozone formation regime, effective ozone control remains challenging. Different from PM2.5, whose main precursor are NOx, VOCs and SO2, the formation and decomposition of ozone are closely related to two types of precursors, VOCs and NOx. There is a diversity of reactions between VOCs and NOx under different 35 meteorological conditions and concentration scenarios (Wang et al., 2017). Since VOCs and NOx can either promote or restrict ozone production, VOCs/NOx is crucial for surface ozone concentrations. However, the thresholds at which VOCs/NOx may promote or restrict ozone production remain unclear (Jin et al., 2017;Schroeder et al., 2017). For instance, under a specific VOCs/NOx scenario, the reduction of NOx may conversely increase surface ozone concentrations (Sillman et al., 1990;Kleinman, 1994). Furthermore, given the large variations of meteorological conditions and ozone level across 40 China, the effects of VOCs/NOx on surface ozone concentrations also demonstrate notable spatiotemporal patterns. In this case, a comprehensive understanding of how the variations of VOCs and NOx could influence ozone concentrations under different VOCs/NOx circumstances is crucial for setting effective emission-reduction policies accordingly in different regions.
To examine the complicated non-linear relationship between ozone concentrations and multiple precursors, a large body of 45 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, NO2 and HCHO proved to be effective proxies for NOx and VOCs (Sillman et al., 1990;Martin et al., 2004). Since NO2 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;Zhang et al., 2020). Cheng et al. (2018) proved that NO2/NO presented a good consistence with long-term ozone concentrations in Beijing. However, NO was not an 50 easily recordable precursor based on satellite observations and not applicable in large-scale monitoring. Cheng et al. (2019) suggested that satellite retrieved HCHO/NO2 was strongly correlated with surface ozone concentrations in Beijing. Different HCHO/NO2 indicates distinct ozone formation regimes, including VOC-limited, transitional and NOx-limited regimes. For VOC-limited (NOx-saturated) regime, the control of VOCs emissions leads to the reduction of organic radicals (RO2), the RO2-NOx reactions and thus ozone concentrations (Milford et al., 1989). On the contrary, the decrease of NOx promotes 55 VOCs-CO reaction, leading to the increase of ozone concentration (Kleinman, 1994). For NOx-limited regime, the reduction of NOx slows down NO2 photolysis, which products free oxygen atoms for ozone formation, and reduces ozone concentrations. The variations of VOCs exert limited influences on ozone concentrations for this regime (Kleinman, 1994).
For transitional (VOCs-NOx mixed) regime, both VOCs and NOx impose positive influences on ozone concentrations. Since the transitional regime divides VOC-limited and NOx-limited regimes, the estimation of the transitional regime range plays a 60 key role to identify 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 3 al. (2020) employed a polynomial model and calculated the transitional regime range over U.S. urban areas as [3.2, 4.1] based on decades of remote sensing and ground observation data. However, given the notable difference of meteorological 65 conditions, ozone levels and the composition of precursors across different countries, whether the transitional regime range extracted in 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 , and may lead to large uncertainties. Consequently, ozone-precursors sensitivity, especially the transitional regime range across China, requires further in-depth analysis. 70 To this end, this research attempts to investigate the spatiotemporal variations of ozone formation regimes across China and identify the transitional regime range of HCHO/NO2 based on the cross-verification of multiple models. Firstly, long-term variations of HCHO and NO2 across China were analyzed. Next, the datasets of HCHO, NO2 and ozone were examined using a polynomial model and a causality model respectively to reveal the crucial thresholds of HCHO/NO2 that separates the NOx-limited, VOCs-limited, and transitional regimes. Specifically, due to the large area of China and potential spatial 75 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 Figure 1), to explore the spatiotemporal variations of ozone formation regimes. Meanwhile, we also compared the ozone formation regimes in urban and rural areas. This research sheds useful lights for better modeling complicated ozone-precursors relationship, understanding the major drivers for enhanced 80 ozone pollution and implementing specific emission-reduction measures to mitigate ozone pollution across China.

Data sources
In this study, OMI HCHO/NO2 datasets were employed for exploring the spatiotemporal variations of HCHO and NO2 in 85 China and calculating HCHO/NO2. We connected surface ozone network data to HCHO, NO2 and HCHO/NO2, which served as the input data for running third-polynomial model and Convergent Cross Mapping (CCM). MODIS land cover product provided the spatial distribution of urban areas, which was employed for identifying urban and rural pixels.
In this study, we employed 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 0.1° × 0.1° 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 95 validation report suggested that the error of this product was effectively controlled within 30% over polluted areas (Gonzá lez Abad et al., 2015), and validated for detecting long-term variations of HCHO columns (Zhu et al., 2017;Shen et al., 2019b).
The daily level-3 gridded OMI NO2 product (OMNO2d), provided by NASA's Goddard Space Flight Center, were 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 100 absorption spectroscopy (DOAS) was employed for retrieving the NO2 slant columns, which were successively transformed into tropospheric and stratospheric vertical columns through AMFs (Bucsela et al., 2013). The OMI NO2 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 NO2 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 105 this study (Kroon et al., 2011;Zhu et al., 2014;Krotkov et al., 2017). The May-to-September OMI HCHO and NO2 products were acquired from NASA's Goddard Earth Sciences Data and Information Services Center (https://disc.gsfc.nasa.gov/).

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/). The unit of surface ozone concentrations in this dataset is μg/m 3 . 110 The network had 1633 monitoring stations, which were distributed in 330 cities across China in 2019. We used the observation data from 13:00 to 14:00 at local time to match the overpass time of OMI. This dataset has been employed in many studies to investigate the variations of surface ozone concentrations in China Shen et al., 2019a;.

MODIS land cover product 115
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). MCD12C1 product was obtained from NASA's Earth System Data and Information System (https://earthdata.nasa.gov/). 120

Data pre-processing
Due to different spatial resolution of OMI HCHO, OMI NO2 and MCD12C1, bilinear interpolation method was used for resampling all above-mentioned products to the same spatial size (0.25° × 0.25°). Meanwhile, we also calculated mean hourly surface ozone concentrations on the 0.25° × 0.25° grid (Figure 1).

Methods 125
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 NOx. However, there were large biases in estimating the range of transitional regime based on chemical transport models (Jin et al., 2017;Jin et al., 2020), due to the uncertainties of emission inventory and the setting of model parameters. Employing observation data alone could effectively overcome these limitations, and the 130 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 PM2.5 concentrations Chen et al., 2019b2019c;Chen et al., 2020), which is a promising tool for investigating the relationships between ozone and its precursors. To increase the reliability of estimated range of transitional regime, both the polynomial 135 model and CCM were employed in this research. We employed the third-order polynomial model for fitting surface ozone concentrations to the indicator of HCHO/NO2. CCM was employed for quantifying the influences of HCHO and NO2 on surface ozone concentrations, and Wilcoxon test (Gehan, 1965) was used for examining whether the differences between the causality of HCHO and NO2 on ozone concentrations at different ranges of HCHO/NO2 was significant. Since the algorithms of the two models are quite different, their cross-verification provides useful reference for their reliability. Meanwhile, 140 Mann-Kendall (M-K) test (Kendall, 1948) was employed for exploring the spatiotemporal variations of HCHO, NO2 and ozone formation regimes in China. Furthermore, we extracted all urban and rural areas in China and compared the differences of ozone formation regimes over these two types of areas. The workflow of the models employed in this study is shown in Figure 2.

Estimating the transitional range of ozone formation regime using polynomial simulation
HCHO and NO2 are considered as proxies for VOCs and NOx, respectively. HCHO/NO2, 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;Cheng et al., 2019;Jin et al., 2020). 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-150 attainment events (surface ozone concentrations exceeding 200 μg/m 3 ) in total events at a given range of HCHO/NO2: where Eventsattainment and Eventsnon-attainment denote the attainment and non-attainment events, respectively (Pusede and Cohen, 2012;Jin et al., 2020).

6
In this study, we used a third-order polynomial model (Jin et al., 2020) to explore the quantitative relationships between 155 HCHO/NO2 and ozone exceedance probability. There were 174868 paired observations of surface ozone concentrations and HCHO/NO2 from 2014 to 2019. The peak of fitting curve highlights the turning point of VOC-limited and NOx-limited regimes (Jin et al., 2020). The range of HCHO/NO2, 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 was of spatial stratified heterogeneity (SSH), as 160 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 q-statistic to quantify SSH and the equation is summarized as follows: Where N and σ 2 denote the number of samples and the variance of population, h is the number of stratifications. The range 165 of q-statistic is [0, 1], and t. 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 above-mentioned stratification, we wouldcould apply the polynomial model in each strata, separately.

Estimating the transitional range of 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 170 the influences of other factors such as meteorological conditions (Chen et al., 2019b2019c;Chen et al., 2020), to extract the causal influences of HCHO and NO2 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 . Therefore, we employed CCM to compare the sensitivity of ozone to HCHO and NO2 at different ranges of HCHO/NO2. CCM utilizes convergent maps to demonstrate the bidirectional coupling between the time 175 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 of two variables X and Y. {X} generates the shadow manifold MX.
Following this, the location of lagged-coordinate vector on MX, x(t) is determined, and then E + 1 nearest neighboring points Where ωi stands for a weight calculated based on the distance between X(t) and its ith nearest neighboring point. Y(ti) stands for contemporaneous value of Y. CCM calculates cross map skill (ρ value) that explains the quantitative relationships.
Number of dimensions for the attractor reconstruction (E), time lag (τ) and number of nearest neighbors to use for prediction 7 (b) are required parameters for CCM. According to previous studies Chen et al., 2019b2019c), E, τ and b 185 was 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, NO2 time series and corresponding surface ozone time series. CCM was implemented using "pyEDM" package in Python. Wilcoxon test (Gehan, 1965) was used to examine whether the differences of ρ values between HCHO and NO2 were significant at the given HCHO/ NO2. No significant difference was regarded as the transitional regime, while significant difference indicated 190 the VOC-limited or NOx-limited regime.

Trend analysis
Mann-Kendall (M-K) (Kendall, 1948) test, which has been used in recent studies on HCHO and NO2 (Cheng et al., 2019;Wang et al., 2019;Zeb et al., 2019), was employed to estimate the significance of trends. M-K test is capable of processing samples with random distributions and mitigating the effects of outliers. Z value is calculated as follow: 195 where S denotes the statistic to be tested, () 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 Z , indicating the trends of samples pass the tests at 90%, 95% and 99%, respectively. 200

Comparison of ozone formation regimes in urban and rural areas in China
To compare the differences of 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 (Figure 3). We set the size of each buffer as 27.75 km, which was close to the size of the 0.25° × 0.25° grid (27.75 km ≈ 0.25°). The first and second buffers 205 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 of meteorological conditions (Yao et al., 2019). 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 ( Figure 6). The largest and smallest increasing trends were shown in NCP and 220 SCB, with a mean value of 0.136 × 10 15 molec/cm 2 year -1 and 0.046 × 10 15 molec/cm 2 year -1 . The increasing trend of YRD and PRD were 0.066 molec/cm 2 year -1 and 0.058 molec/cm 2 year -1 , respectively. Meanwhile, reversed trends were detected for NO2 during the two periods ( Figure 5), which was consistent with previous studies (Jin and Holloway, 2015;Li et al., 2019a). From 2005 to 2012, the averaged NO2 was 2.027 × 10 15 molec/cm 2 and the annually mean increasing trend was 0.098 × 10 15 molec/cm 2 year -1 . Thanks to the implementation of Clean Air Action, the averaged NO2 were reduced to 1.

Transitional range of ozone formation regime 235
According to HCHO/NO2, We divided the paired observations into 200 bins for the whole country and 100 bins for these megacity clusters, and the ozone exceedance probability was calculated for each bin. The third-order polynomial was employed for fitting ozone exceedance probability to HCHO/NO2. As shown in Figure 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 geographical detector to examine the SSH of annual May-to-September mean surface ozone concentrations in China. As 240 shown in Table 1, all the q-statistics from 2014 to 2019 were greater than zero, which indicated the surface ozone concentrations in China were of SSH. As suggested by Chen et al. (2020), meteorological factors including temperature, 9 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, iIn addition to the regime range extracted at 245 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 transitional regime for 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 the range of 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 was 250 not available. Statistical bootstrapping was used for estimating the uncertainty of 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 two standard deviations from the peak of the fitting curve. The uncertainty for the third-polynomial model was 0.4, indicating a significant 255 nonlinear relationship between ozone exceedance probability and HCHO/NO2.
Due to the limited data used for running CCM, we set the bin size of HCHO/NO2 as 0.2 for collecting sufficient ρ values to conduct Wilcoxon test. As shown in Figure 8b, there were no significant difference between ρ of HCHO and NO2 when HCHO/NO2 ranged from 0.9 to 1.9, which indirectly defined the range of transitional regime. For HCHO/NO2 < 0.9, ρ of HCHO was notably higher than that of NO2, and this range was regarded as VOC-limited regime. Similarly, HCHO/NO2 > 260 1.9 suggested the NOx-limited regime. Through the cross verification, it was an important finding that the range of transitional ozone formation regime estimated using the third-order polynomial model and CCM was highly close, indicating the reliability of the extracted range.

Ozone formation regimes in China 265
NO2 demonstrated a significant downward trend since 2013, while HCHO kept the increasing trend during the entire study period. Consequently, HCHO/NO2 increased in a majority of regions across China. Specifically, the annually increasing trend of HCHO/NO2 in NCP, YRD and PRD was 0.035 year -1 , 0.023 year -1 and 0.034 year -1 , respectively. Meanwhile, there were no significant trends in SCB during this period (Figure 9). The variations of HCHO/NO2 indicated the shrinkage of VOC-limited regime and the expansion of transitional and NOx-limited regimes. Since the range of transitional regime 270 estimated by 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 NCP, YRD and PRD. The proportions of areas with the VOC-limited regime in 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 SCB. Areas with the transitional regime occupied 60%, 50%, 14% and 20% 275 10 in NCP, YRD, PRD and SCB. NOx-limited regime dominated other areas (Figure 10a). In 2019, areas with the VOCslimited regime decreased significantly, and was simply found in the fringe areas of NCP and YRD. The proportion of the VOCs-limited regime in NCP and YRD was 2% and 9%, respectively. The transitional regime was widely distributed NCP, YRD and SCB, and occupied the 71 %, 56% and 36% of the total areas. The NOx-limited regime still spread over a majority of China (Figure 10a). We calculated the annually mean ρ of HCHO and NO2 over those megacity clusters from 2014 to 280 2019 (Figure 10b). For all megacity clusters, the ρ of NO2 was higher than HCHO, indicating that NO2 was the dominant factor for surface ozone concentrations. Both models suggested that NO2 played a more important role in affecting surface ozone concentrations than HCHO. In the past several years, NOx-oriented emission-reduction has been conducted across China, leading to the continuous decrease of NOx concentrations. Since both VOCs and NOx imposed positive influences on surface ozone concentrations under the transitional and NOx-limited ozone formation regime, the upward trend of HCHO 285 across China might explain recent soaring ozone concentrations across China Lu et al., 2020). It is noted that the difference between the ρ of NO2 and HCHO decreased notably in NCP and YRD. This may be attributed to the

Variations of ozone formation regimes in urban and rural areas
Previous studies suggested that the differences of ozone formation regimes existed between urban and rural areas (Tong et 295 al., 2017;Liu et al., 2018b;Cheng et al., 2019). We extracted HCHO and NO2 columns in urban and rural pixels in those megacity clusters, and calculated the annually averaged HCHO/NO2 (Figure 11). For NCP, HCHO/NO2 in urban areas was higher than 1.0 since 2015, indicating a transformation from VOC-limited to transitional regime. The increase of HCHO/NO2 was attributed to the reversed variation trends of HCHO and NO2. The rising HCHO resulted from the increase of anthropogenic emissions and biogenic volatile organic compound (BVOC) (Shen et al., 2019b;Wang et al., 2020), while 300 the implementation of Clean Air Action imposed notable influences on the decrease of NO2 (Chen et al., 2019a).
HCHO/NO2 in rural areas was in the range of [1.0, 1.9], indicating rural areas were occupied by the transitional regime from 2005 to 2019. For YRD, which was occupied by the transitional regime, no variation of ozone formation regime was found in urban areas. In rural areas, HCHO/NO2 temporally exceeded the threshold of 1.9 from 2016 to 2018, indicating the ozone formation regime changed from transitional to NOx-limited. This phenomenon was attributed to the slight decline of HCHO, 305 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 of NO2 concentrations, the urban and rural areas in PRD was dominated by transitional regime and NOx-limited regime. For SCB, HCHO/NO2 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. China as [1.3, 2.8], which was notably higher than the range across China. Jin et al. (2020) calculated the range of transitional regime in several major regions in the US using QA4ECV dataset, whose spatial resolution was 0.125°, and the 320 output [3.2, 4.1] was much larger than the averaged range of transitional regime across US. One reason could be the severe ozone pollution in mega cities, leading to different ranges of transitional regime. Meanwhile, the calculated range of transitional regime is closely related to the spatial resolution of employed HCHO and NO2 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;Jin et al., 2020). In addition to the generally consistent outputs, some advances of this research are listed as 325 follows. First, only a few parameters are required for polynomial model and CCM, which effectively reduced the uncertainties of 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. 330 Despite a generally reliable output, some uncertainties exist. First, the accuracy of the estimated range of transitional regime might be influenced by the scaling biases between station-based observations of surface ozone and space-based HCHO and NO2. Since ozone monitoring stations are mainly distributed in urban areas, and a 0.25° × 0.25° grid might cover both the urban and rural areas, the surface ozone concentrations of a grid may be overestimated. In addition to the locations and the 335 spatial resolution of data sources Second, the uncertainties of OMI HCHO and NO2 datasets might impose negative influences on the estimation of transitional regime range (Duncan et al., 2010;Jin et al., 2017;Schroeder et al., 2017;Jin et al., 2020). FirstlyOn one hand, errors exist in the retrieval of HCHO and NO2 vertical columns. SecondlyOn the other hand, vertical mixing was not homogeneous, weakening the capability of using HCHO and NO2 vertical columns to explore the near-surface ozone-precursors sensitivity. Therefore, future improvement of earth observation techniques and the 340 spatiotemporal resolution of HCHO and NO2 products can further enhance the accuracy of the estimated range of 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 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 of transitional regimes in different regions, when conducting small-scale research, the range of [1.0, 1.9] may be adapted accordingly based 345 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 models to cross-verify the range of transitional regimes. Despite a relatively high fitting accuracy in terms of uncertainties, the findings from these studies could 350 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 more precise description of ozone-precursors 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 NO2 and HCHO affect ozone formation 355 regime can be better explained in a theoretical environment. With more reliable experimental reference data, the modelbased analysis on the range of transitional regime at local, regional and national scale can be further improved accordingly.
According to the temporal variations of OMI NO2 concentrations across China, a notable decreasing trend was observed in three major megacity clusters, NCP, YRD and PRD. These regions were heavily polluted by PM2.5 and the notable decrease of NO2 was mainly attributed to the national Clean Air Action since 2013 (Zheng et al., 2018), which aimed to reduce PM2.5 360 concentrations by cutting NOx emissions. The influence of Clean Air Action on the reduction of PM2.5 concentrations and NOx 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 NOx significantly accelerated since 2013.
Moreover, Zhang et al. (2020) employed random forest algorithm to remove the effects of meteorological conditions, and evaluated the impacts of Clean Air Action. The results demonstrated that the deweathered NO2 concentrations in winter 2007 365 and 2017 were 70.3 μg/m 3 and 59.1 μg/m 3 , 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 of NO2 and HCHO led to the increase of HCHO/NO2, and the increase of transitional areas and NOx-limited regime areas. From 2013-2019, all these regions were dominated by the transitional or NOx-limited regimes. Attributed to the long-term variation of formation regimes, a more complicated and 370 fragmented spatial pattern was observed across China. Consequently, for an effective control of ozone pollution, the emission-reduction of both NOx and VOCs is required. Especially for NCP and YRD, where the NOx has been remarkably reduced, effective approaches for controlling VOCs emissions are essential for preventing ozone pollution. This finding was 13 consistent with previous studies (Li et al., 2019b), which recommended the simultaneous reduction of NOx and VOCs for mitigating the composite airborne pollution in China. Admittedly, compared with NOx-reduction, the VOCs-reduction is 375 more complicated and the output of anthropogenic VOCs reduction is more unpredictable. In this case, reducing biogenic VOCs emissions can also be a potential solution. VOCs emitted by vegetation takes 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 in summer . 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 380 waves and successive ozone pollution (e.g. summer, 2017).
The large spatial variations of HCHO/NO2, especially the rapid increase of transitional regime areas across China, indicates that a unified NOx-VOCs reduction strategy is not feasible for the entire country. Instead, to effectively reduce ozone concentrations, the specific proportion of NOx and VOCs reduction should be carefully set according to local HCHO/NO2.
Meanwhile, due to the large differences in vehicle and industrial emissions (Cheng et al., 2019), the concentration of NOx is 385 notably higher in urban areas. Therefore, the further reduction of NOx emissions exerts a stronger influence on ozone reduction in rural areas compared to urban areas.

Conclusions
To better understand the spatiotemporal variations of ozone formation regimes across China, we employed the third-order polynomial model and CCM to estimate the range of transitional regime from 2005 to 2019, the results of which were [1.0, 390 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 transitional regime in four megacity clusters and found the range in 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 NO2 led to the increase of HCHO/NO2, indicating China was dominated by the transitional and NOx-limited regimes in recent years. We also found that the ρ of NO2 was higher than HCHO at all megacities, suggesting 395 that the reduction of NOx emissions would become more effective in controlling surface ozone concentrations. Meanwhile, given the rising VOCs emissions, the simultaneous reduction of NOx and VOCs would be more effective than the sole reduction of NOx in mitigating ozone pollution. Finally, the comparison of ozone regimes in urban and rural areas suggested that the reduction of NOx emissions would impose stronger impacts on the control of ozone pollution in rural areas.            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 p-value indicate statistical significance at α = 0.05, 0.01 and < 0.001 level, respectively.