Increase in Surface Ozone over Beijing-Tianjin-Hebei and the Surrounding Areas of China Inferred from Satellite Retrievals, 2005-2018

Received: November 21, 2019
Revised: April 10, 2020
Accepted: April 22, 2020

 Copyright The Author(s). This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are cited.

Cite this article:

Li, Y., Cheng, M., Guo, Z., He, Y., Zhang, X., Cui, X. and Chen, S. (2020). Increase in Surface Ozone over Beijing-Tianjin-Hebei and the Surrounding Areas of China Inferred from Satellite Retrievals, 2005-2018. Aerosol Air Qual. Res. 20: 2170–2184. https://doi.org/10.4209/aaqr.2019.11.0603

HIGHLIGHTS

• The long term SOC was retrieved based on satellite data and GWR model.
• The highest and lowest SOCs occurred in the southeast and northeast, respectively.
• The growth of the SOC in cold season had a key contribution for annual growth.
• The SOC increased faster since 2012 because of the rapid decrease of NO₂.

Keywords: OMI ozone profile; Geographically weighted regression model; Surface ozone concentration; Spatio-temporal distribution; Beijing-Tianjin-Hebei and the surrounding areas.

INTRODUCTION

Surface ozone, which is a strong oxidant, with respect to radiative effects, is one of the three most significant greenhouse gas with negative impact on environmental change, public health, vegetation, crop production and environmental conditions (Lee et al., 2014; Monks et al., 2015; Plocoste et al., 2017; Wang et al., 2019a); the compound has recommendations in air quality guidelines in all countries (Kim et al., 2017; Maji et al., 2019). Surface ozone is resulting from number of complicated chemical reactions among the volatile organic compounds (VOCs) and nitrogen oxides (NO_x = NO + NO₂) in sunlight within the wavelength range of 200–300 nm (Lin et al., 2009; Choi et al., 2012; Jin and Holloway, 2015; Wang et al., 2017a; Zhong et al., 2019). It is very difficult to control the surface ozone as the relation between ozone and its precursors is complicated and non-linear (Zoogman et al., 2014; Wang et al., 2019b). Literature reports have indicated surface ozone to be leading summertime air pollutant overtaking fine particles with an aerodynamic diameter of 2.5 µm or smaller (PM_2.5) in 74 Chinese cities (Zhao et al., 2018b; Wang et al., 2019b) and has become an increasingly severe type of air pollution in China (Zhang et al., 2014; Feng et al., 2015; Feng et al., 2018; Qiao et al., 2019).

Ozone was identified as an air pollutant with the 1-h-averaged concentration exceeding 600 ppbv on occasion in Los Angeles in the ’40s and ’50s of the 20^th century (Yu, 2019). To study the spatio-temporal characteristics of surface ozone, various sites for monitoring for long period of time have been established worldwide (Wang et al., 2017b; Zhao et al., 2018b), such as the Photochemical Assessment Monitoring Stations (PAMS) established in the United States (Tsai et al., 2008), the European air quality database network (Airbase) established in Europe (Yan et al., 2019), the Network for the Detection of Atmospheric Composition Change (NDACC) and other ozonesondes and Lidar stations contributing to the World Meteorological Organization’s Global Atmosphere Watch (WMO GAW) (Keppens et al., 2018). Based on the advantages of global coverage and high spatio-temporal resolutions (Cai et al., 2018), atmospheric observation via satellite remote sensing is an effective complement to surface observations. Most satellite ozone instruments provide global coverage and total ozone columns, such as the Total Ozone Mapping Spectrometer (TOMS) (Fishman et al., 2010; Vyas and Saraswat, 2013), Global Ozone Monitoring Experiment (GOME) (Choi et al., 2008), Total Ozone Unit (TOU) (Zhang et al., 2018), Scanning Imaging Absorption Spectrometer for Atmospheric Chartography (SCIAMACHY) (Anton et al., 2013) and Ozone Monitoring Instrument (OMI) (Shukla et al., 2017; Bozhkova et al., 2019). Chen et al. (2014) analyzed the spatio-temporal variations in the total ozone column (TOC) over the Yangtze River Delta from 1978 to 2013 using remote sensing TOC data derived from the TOMS (1978–2005) and OMI (2004–2013) satellite instruments. Some satellite ozone instruments also provide a distribution of ozone profiles, such as the Solar Backscatter Ultraviolet Radiometer (SBUV) (Noreen et al., 2018), Tropospheric Emission Spectrometer (TES) (Kulawik et al., 2006), Ozone Mapping Profiler Suite (OMPS) (Keppens et al., 2018), Microwave Limb Sounder (MLS) (Ziemke et al., 2011), GOME-2 (van Peet et al., 2018) and OMI. Shen and Wang (2012) analyzed the tropospheric ozone observations over China (western, northern and southern) using TES in the time frame of 2005–2010 for different pressure levels (484, 681 and 825 hPa). Since 2005, the ozone data derived from the OMI have been widely used because of its high spatio-temporal resolution and the rich data products supporting ozone research at different atmospheric layers. Especially at the sea level, ozone concentrations based on OMI are similar to the average concentrations observed from monitors in the ground state (Di et al., 2017). Researchers often use OMI alone or in combination with other instruments that can observe ozone profiles (Huang et al., 2017). Grigoras et al. (2016) presented a method for estimating the tropospheric ozone concentrations over Bucharest, Romania, in 2009 using the measured data from the OMI. Algorithm derived by van Peet et al. (2018) for the combinational assimilation of the GOME-2 and OMI ozone profiles, resulted in improvement in the measurement of derived ozone concentrations by measuring at various time frames within a day.

China from 1980s has started building observation stations for continuous surface monitoring for WMO GAW, like at Shangdianzi, Linan, Longfengshan and Waliguan (Shen and Wang, 2012; Zhao et al., 2018a). China started releasing air quality monitoring data in the ground level to the public though a national air quality monitoring network (AQMN) set up in 2012 (Wang et al., 2017b; Maji et al., 2019). Moreover, China’s ambient air quality standards (GB3095-2012) set the mean daily maximum 8-h (DMA8) ozone concentration limit at 100 µg m^–3 for primary and 160 µg m^–3 for secondary levels in 2012. However, compared to those in developed countries, the data of ozone observations in China are still very spatio-temporally limited. Although numerous studies have been conducted on ozone episodes, scarce reports are available on the spatio-temporally variability over longer duration in the surface ozone concentration (SOC) in the Beijing-Tianjin-Hebei (BTH) region because of limited data availability (Wang et al., 2009; Cai et al., 2018). Herein, we used the data of the surface monthly mean DMA8 SOC and ozone profile from OMI to understand the SOC over a long period along with the spatio-temporal variations and the causes of the increasing surface ozone were analyzed over the BTH and its surrounding areas from 2005 to 2018.

DATA AND METHOD

Study Area

The study area is the BTH and the surrounding areas (34°12ʹ–41°03ʹN, 111°30ʹ–119°19ʹE), which is shown in Fig. 1. The abbreviations for the 28 cities of study area are shown in Table S1. The Yanshan Mountains in the north, North China Plain in the south and Bohai Sea in the east show the topographic characteristics of high in the northwest and low in the southeast which make the pollutant difficultly diffuse. Relevant reports have indicated that the study area emits more than 10% of the SO₂ and VOCs and more than 15% of the NO_x and primary particulate matter in less than 3% of the Chinese land area (China’s Ministry of Environment Ecology, 2018). The intensity of atmospheric pollutant emissions is high, and heavy pollution weather occurs frequently in the study area due to economic development and industrial agglomeration. 

Fig. 1. Schematic map of the study area and ground monitoring points.

Ground Monitoring Data

In 2013, the China National Environmental Monitoring Center (CNEMC) conducted 24-h continuous online monitoring of the conventional atmospheric pollution factors and meteorological parameters at 496 national monitoring sites in major 74 cities (Lu et al., 2018). And by 2018, the nationwide observational network included 1436 national monitoring sites covering 338 cities. These measurements provided precise and reliable data supporting the evaluation of air quality levels and for controlling and preventing air pollution and assists environmental protection decision-making. The ultraviolet absorption spectrometry method and differential optical absorption spectroscopy were used in monitoring sites to measure ozone concentration and the CNEMC reported air pollutants at the standard atmospheric condition for temperature (273.15 K) and pressure (1 atm) (Lu et al., 2018; Maji et al., 2019). Obtained data from monitoring sites were refined to remove discrepancy; generally, if less than 6 h of ozone data were available for a certain 8-h period of measurement in each third of the day (i.e., 00:00–08:00, 08:00–16:00, and 16:00–24:00 Beijing Standard Time [UTC+8]), then this DMA8 was assigned as the “missing” value, and no more than two consecutive hourly measurements missing within one day (Zhan et al., 2018). In 2013, only 13 cities in the study area were included in 74 cities providing monitoring data, and the distribution was uneven, which was not enough to represent the whole study area. In addition, the National Urban Air Quality Real-time Publishing Platform (NUAQRPP) of the CNEMC only provided the real-time monitoring data nationwide since May 2014 (http://beijingair.sinaapp.com/). Therefore, the ground monitoring data of this study are the DMA8 SOC from May 2014 to December 2018 from the NUAQRPP of the CNEMC. The distribution of the ground monitoring sites is shown in Fig. 1.

Although the ground monitoring data can best reflect the SOC near the monitoring stations, the ground monitoring points are mostly distributed in urban areas, and their coverage is limited. Therefore, the ground monitoring data were spatially interpolated to obtain the regional ground monitoring data set first. Herein, 80% of the ground monitoring stations in the study area were selected as the training set, while the remaining 20% were selected as the verification set. Moreover, considering that there almost were urban sites in the verification set which was more representative of urban regions, we selected 14 suburban environmental evaluation sites of the study area as the non-urban verification set to verify the applicability of the derived ozone in the non-urban areas. The locations of the urban and suburban sites are shown in Fig. S1. The study used the inverse-distance-weighted (IDW) method in the training set to obtain the monthly mean DMA8 SOC raster data set over BTH and its surrounding areas from May 2014 to December 2018. Based on the station information of the 20% urban sites and 14 suburban environmental evaluation sites of the study area, the interpolated training set was sampled and extracted through geostatistical analysis, and then the results of the geostatistics were tested for correlation with the urban and suburban verification sets. In the verification of interpolated training set and urban sites monitoring data, the coefficient of determination (R²) is 0.97 and the mean percentage of the absolute error (MPAE) (∑[|ozone_I –ozone_G|/ozone_G × 100%]/n) is 6% where ozone_I and ozone_G are the interpolated ozone concentration and ground monitoring ozone concentration, respectively, and n is the number of ground monitoring points in the verification set. In the correlation test of interpolated training set and suburban sites monitoring data, the R² is 0.90 and the MPAE is 10%. The trend line of urban sites was almost coincident with the 1:1 indication line and the suburban sites almost fallen within the 1:2 and 2:1 lines, as shown in Fig. S2. And the linear regression results for interpolated training set and urban/suburban sites monitoring data were shown in Table S2. The verification results showed that the overall correlation was good with small error, and the interpolated results of the study area were reliable for both urban and non-urban regions.

Remote Sensing Data

Compared to ground-based ozone monitors, satellite instruments have the advantages of high spatio-temporal resolutions and global coverage. The OMI is on board the Earth Observing System (EOS) Aura satellite, which was put to orbit in July 2004. With the 2600 km OMI swath, global coverage is achieved in one day under normal operating conditions. The OMI Level 2 product OMO3PR contains the retrieved ozone profile and the corresponding a priori ozone profile, error covariance matrix and averaging kernel. Furthermore, it contains ancillary information produced by the optimal estimation algorithm applied to OMI global-mode measurements. The retrieval algorithm of OMO3PR is based on optimal estimation with climatological a priori information (Rodgers, 2000). Basically the amount of ozone in each atmospheric layer is adjusted such that the difference between the modeled and measured sun-normalized radiance is minimal. As the information content in the measured spectrum is not large enough to determine the ozone amount for all of the layers, a side constraint is used by demanding that the retrieved profile does not differ too much from the climatological average. The measurements are taken from the UV1 channel (270–308.5 nm) and the first part of the UV2 channel (311.5–330 nm). Here two UV2 pixels are combined to obtain the spectrum of a pixel that corresponds to a pixel in the UV1 channel. According to the atmospheric pressure level, the ozone vertical profile is divided vertically into 18 layers from the surface pressure (1000 hpa) to 0.3 hpa, which is about 0–60 km height. The ozone column concentration in each layer is in Dobson units (DU), which can be converted into a volume mixing ratio (ppm) as follows (Liu et al., 2016):

V_mri = 1.2672N_i/∆P_i                                                                  (1)

where layer-wise concentration of ozone column (DU) is designated as N_i, ΔP_i is the pressure difference (hPa) between the top and bottom of the layer, and V_mri is the average volume mixing ratio (ppm).

The OMI data selected in the study is the lowest layer of OMO3PR, which is approximately 0–3 km in height from 700 to 1000 hpa. The factors considered in screening OMI data are OMI row anomaly, averaging kernels, total root mean square error (RMSE), cloud contamination, the altitude of the study area. In the OMI measurements, three row anomalies (first starting June 25, 2007, second starting May 11, 2008, and third starting January 24, 2009) have been observed. And the anomalies have been investigated. The ozone profiles where the field refectanceCostFunction in the OMO3PR product file has a value smaller than 30 can still be used. Hence, this field has been filtered. Fu et al. (2018) reported that the MUlti-SpEctra, MUlti-SpEcies, MUlti-SEnsors algorithm uses the measured radiances from the quality-assured OMI off-nadir pixels and the corresponding collocated AIRS measurements for retrieval since 2009 to solve row anomaly issues. The averaging kernels provide information on the vertical resolution of the retrieved profile and they are affected by ozone concentration in the boundary layer (Worden et al., 2007). The vertical resolution of OMO3PR is 3 km. The OMI Science Requirements Document (SRD) requires the total RMSE of retrieved ozone profile is less than 30%. According to the OMI Algorithm Theoretical Basis Document (ATBD; Version 2), the total RMSE of the 700–1000 hpa layer is 26.5%, which indicates that the quality of the lowest layer of OMO3PR is reliable. Wang et al. (2011) indicated that the error range of OMI-retrieved data is less than 10% by comparing with the Community Multi-scale Air Quality (CMAQ) model operation data and OMI observations at the lowest layer representing the mean values by comparing with Environmental Protection Agency (EPA)’s surface monitoring data. Referring to Liu et al. (2016) study, we removed the OMO3PR data with cloud fraction greater than 0.5 in the process of OMI data retrieval. In addition, according to the statistical analysis of the DEM 500-m image in the study area, the highest altitude of the study area is 2687 m within the height range of the 700–1000 hpa layer. After screening, a single OMO3PR file contains approximately 330 OMI measurements of 30 ground pixels each. The pixel size of OMO3PR is 13 × 48 km² at nadir. The retrieval algorithm and screening method of OMO3PR have not changed in the research period, therefore the sizes and quality of screened OMI data are consistent.

Finally, the study analyzed the spatial correlation between screened OMI data and ground monitoring data. As shown in Fig. S3, the result of the spatial correlation test between the lowest layer of OMO3PR and ground monitoring data set showed that the correlation coefficient (R) was greater than 0.8 over most study area. The total RMSE of the 700–1000 hpa layer, relevant researches and the spatial correlation test of the study all showed that the lowest layer of OMO3PR can be used to characterize surface ozone.

Geographically Weighted Regression Model

The geographically weighted regression (GWR) model was used to retrieve the long-term surface ozone based on space-based and ground-based data.

The GWR model is an extension of the ordinary linear regression (OLR). Based on the OLR single multivariate linear model, the geographic location of the sample data is embedded into the regression parameters to generate an independent regression equation for each grid point in study area (Gilbert and Chakraborty, 2011). Each GWR equation can be expressed as follows:

y_i = β₀(u_i, v_i) + _β1(u_i, v_i) x_i1 + β₂(u_i, v_i)x_i2 + ⋯ β_k(u_i, v_i)x_ik + ε_i (2)

where y_i is the dependent variable, x_ik is the k^th independent variable of grid point i, (u_i, v_i) is the coordinate position of grid point i, β_k is the k^th regression parameter of grid point i, and ε_i is the random error of grid point i.

The GWR can model the relationships with the spatial variation by generating a continuous surface of the regression parameters through local linear regression (Ma et al., 2014). The choice of spatial weight function and bandwidth is the core of GWR. Any grid point in space will have some influence on its neighboring grid points, and this influence will decrease as the distance increases. In this study, the distance threshold method was used to establish the spatial weight function, and an adaptive kernel method for selecting the number of neighboring grid points which are optimal was used based on the Akaike Information Criterion (AIC). This method is considered to be the closest for finding the most suitable neighborhood number for the regression locations (Chen et al., 2018). In the study, the ground monitoring data set and the satellite data set in 2014–2017 were used as the training set for the retrieved-data model. The training set data are used to establish the GWR equation for each grid of the study area as follows:

Y_i = β₀(r_i, n_i) + β₁(r_i, n_i)X_i + ε_i                                                    (3)

where Y_i is the ground monitoring mean monthly matrix of grid point i from May 2014 to December 2017, (r_i, n_i) is the number of rows and columns of grid point i, X_i is the satellite mean monthly matrix of grid point i from May 2014 to December 2017, β₀(r_i, n_i) and β₁(r_i, n_i) are the regression parameters of grid point i, and ε_i is the random error of grid point i.

According to the established GWR parameters, the regional spatial linear relationship was constructed. The mean monthly satellite data in 2005–2018 were retrieved to obtain the monthly average SOC of the BTH and its surrounding areas in 2005–2018.

RESULTS AND DISCUSSIONS

Validation of Retrieved SOC

To verify the retrieved-data accuracy of the GWR model, the study first compared the average annual retrieved SOC spatially with the 2018 ground monitoring data. The retrieved SOC was generally consistent with the concentrations from the ground monitoring stations, and most of the retrieval errors were within 10 µg m^–3, as shown in Fig. 2(a). In addition, the study spatially sampled the mean monthly retrieved SOC for 2018 through geostatistical analysis to test the correlation with the verification set, which was the 2018 ground monitoring data set. As shown in Fig. 2(b), when the ground monitoring data were lower or higher than 89 µg m^–3, the retrieved results were overestimated by 14% or underestimated by 7%, respectively. And the result of linear regression for validation of retrieved SOC was shown in Fig. 2(b). The R² is 0.82 and MPAE is 18%. In Fig. 2(b), the points in the black rectangle that were located in BJ and collected in September were relatively far from the trend line, which may be related to the Beijing branch of the China-Africa Cooperation Forum in 2018 and other key activities. In response to the mismatch between the predicted and observed surface ozone over Beijing in September 2018, the study showed the comparison between ground-based observations and retrieved SOC over Beijing, 2018, as shown in Fig. 3. By comparing the ground monitoring data and satellite-retrieved SOC over Beijing in September 2018 and other months in 2018, it can be seen that the observed SOC in September 2018 plummeted compared with last months. This showed that satellite data is not sensitive to the short-term abrupt changes of near-ground ozone concentration. Special months like September 2018 in Beijing were rare, so the overall correlation between ground monitoring data and satellite-retrieved data was good and the retrieved SOC was reliable.

Fig. 2. Validation of the retrieved SOC: (a) spatial comparison of the retrieved surface ozone and the ground monitoring data in space and (b) scatter plots for verifying the retrieved results of the satellite data. The red solid line is the result of linear regression fitting the scatter plots. The red dash-dot line is the result of 1:1 fitting through 0. The black dotted lines indicate ± factor 1:2 and 2:1.

The study also compared the surface ozone trend recorded at Miyun site in Beijing in 2005–2007 which showed in Cooper et al. (2014) study with the OMI-retrieved SOC trend of the study. The Miyun site is a suburban environmental evaluation site which is located in the northeast about 100 km away from central Beijing. The region surrounding the site is characterized by a mix of agriculture and small villages. The trend of OMI-retrieved SOC at Miyun site in Beijing in 2005–2007 as shown in Fig. S4 was roughly consistent with the trend of Cooper et al. (2014) study. Cooper et al. (2014) indicated that the 2005–2007 monthly 95^th midday ozone percentiles, summertime (June–August) median ozone mixing ratio at Miyun site was 75 ppbv, corresponding to 146 µg m^–3. The retrieved SOC in the study was 137 µg m^–3 at Miyun site in the same time range. The retrieval error of 9 µg m^–3 indicated that the GWR model established using 2014–2017 ground monitoring data applied to the 2005–2013 period. 

Fig. 3. Comparison between ground-based observations and retrieved SOC over Beijing, 2018.

Increasing of Surface Ozone in Spatio-temporal Distribution from 2005 to 2018

The spatial distributions of SOC and growth rate over the study area during the whole year of 2005–2018 are shown in Figs. 4(a) and 4(b), respectively. The mean SOC over the study area ranged from 65 to 104 µg m^–3. The mean concentration in the area of study was 86 µg m^–3, and the standard deviation was 7 µg m^–3. The range and standard deviation of the SOC reflected the difference in the spatial distribution of the long-term sequence of the study area in 2005–2018. The highest concentrations were in the southeast and southwest, including in PY, HZ, JN, CZHI, JC and JZ, and the lowest concentration was in TJ. The SOC of the study area during 2005–2018 was generally on the rise, with an average growth rate of 3.4 µg m^–3 y^–1. The growth rate of SOC over the study area ranged from 0.7 to 5.9 µg m^–3 y^–1. The highest growth rate distributed in the northeast, including LF, TJ, CZ, HS, DZ, BZ, and the lowest growth rate distributed in the southwest, including TY, YQ, CZHI, JC. 

Fig. 4. Spatial distributions of (a) the SOC and (b) the growth rate in the whole year, (c) the SOC and (d) the growth rate in the warm season and (e) the SOC and (f) the growth rate in the cold season over the study area during 2005–2018.

The spatial distributions of SOC and growth rate over the study area in the warm season (April–September) and cold season (October–March of the following year) of 2005–2018 are shown in Figs. 4(c), 4(d), 4(e) and 4(f). In the warm season, the investigated area showed a mean concentration in the range of 102–146 µg m^–3. The average value in the area of investigation was 126 µg m^–3, with 7 µg m^–3 being the standard deviation. The SOC showed “east-west” zoning, with BJ, BD, HS, DZ, JNAN, HZ and JN as the areas in the “high-value belt.” The highest concentrations were in the “high-value belt,” and the lowest concentrations were on the eastern and western sides of the “high-value belt.” The growth rate of SOC over the study area ranged from 0 to 4.2 µg m^–3 y^–1. The highest growth rate distributed in CZ, BZ, and the lowest growth rate distributed in the southwest. In the cold season, the investigated area had concentrations between 25 and 69 µg m^–3, with 45 µg m^–3 and 10 µg m^–3 being the average concentration and standard deviation respectively. The highest concentrations were in the southeast and southwest, including PY, HZ, JN, KF, YQ, CZHI, JC and JZ, and the lowest concentrations were in LF and TJ. The growth rate of SOC over the study area ranged from 1.5 to 8.1 µg m^–3 y^–1. The highest growth rate distributed in the northeast, and the lowest growth rate distributed in the southwest. The spatial distribution in the cold season was consistent with that throughout year. The average growth rates were 2.3 µg m^–3 y^–1 and 4.6 µg m^–3 y^–1, respectively, in the warm season and cold season. The average concentration in the warm season was approximately 2.9 times that in the cold season, while the average growth rate in the cold season was twice that in the warm season. The analysis indicated that the cold season had a great influence on the spatial distribution and annual variation of the whole year, while the warm season contributed greatly to the SOC.

The monthly mean satellite-retrieved SOC in BTH and its locality from 2005 to 2018 and the monthly mean ground-observed SOC from 2014 to 2018 are shown in Fig. 5. And the linear regression results for the trends of the monthly mean retrieved SOC and ground-observed SOC shown in Fig. 5 was shown in Table 1. The R between satellite-retrieved ozone and ground observations was 0.94 indicating that the OMI-retrieved ozone trends were consistent to ground-observed ozone. Averaged over the study area, the mean SOC was 86 ± 7 µg m^–3 and increased at rates of 3.4 µg m^–3 y^–1 from 2005–2018. The growth rate increased in each the monthly maximum, monthly mean and monthly minimum. The month of June represented maximum SOC except in 2009 when the maximum concentration was in July and as for minimum SOC, it was in the months between November to January. Since 2012, the growth of the SOC has been more significant especially in winter. Comparing with previous studies, Lu et al. (2018) discussed the severe surface ozone pollution in China by using the 2013–2017 surface ozone measurements from the Chinese monitoring network and found that all ozone metrics showed continuous increases over China in 2013–2017. Li et al. (2019) also showed increasing ozone trends of 1–3 ppbv y^–1 in megacity clusters of eastern China including BTH that they attributed to changes in anthropogenic emissions. The temporal trends of SOC in the study was consistent with the previous studies. 

Fig. 5. The monthly mean of satellite-retrieved and ground-observed SOC in the study area from 2005 to 2018. The black and blue solid lines are the satellite-retrieved and ground-observed SOC, respectively. The blue dotted line is the trend line of ground-observed SOC. The black, red and green dotted lines are the trend lines of mean, maximum and minimum of satellite-retrieved SOC, respectively. 

Causes for Increasing in Surface Ozone over the Study Area

The study focused on the impact of the VOCs and NO_x which are the primary precursors to the surface ozone in China on the increasing in surface ozone over the study area. Surface ozone and emitted NO react directly as per the reaction known as the “NO_x titration effect” (Wang et al., 2017b). NO₂ forms secondary aerosols and also is a participant in the tropospheric formation of ozone via catalytic pathway, whereas formaldehyde (HCHO) is a proxy for total VOCs and itself one of the most important VOCs (Fu et al., 2007; Wang et al., 2017c). The value of the OMI tropospheric column ratio of HCHO to NO₂ (FNR) has been used to identify the surface ozone formation regime (Castell et al., 2009; Duncan et al., 2010; Witte et al., 2011; Mahajan et al., 2015). From earlier literature reports, FNR < 1.0 was interpreted as indication of a VOCs-limited regime, NO_x-limited regime was indicated by FNR > 2.0, and the value of FNR between 1.0 and 2.0 was indicative of a transitional regime (when both NO_x and/or VOCs reduction is expected to reduce ozone) (Jin and Holloway, 2015). In this study, the mean monthly tropospheric column concentrations of NO₂ in 2005–2018 and HCHO in 2005–2018 from OMI were used. The spatial distributions and annual variation in FNR over the study area during 2005–2018 were shown in Figs. 6(a) and 6(b), respectively. And the linear regression results for the annual variation of FNR in the whole year, summer and winter shown in Fig. 6(b) was shown in Table 2. During 2005–2018, the FNR of most study area was less than 1.0 showing VOCs-limited regime, and the FNR of small area was greater than 1.0 and less than 2.0 showing transitional regime. The FNR was all less than 1.0 whether in the whole year or winter from 2005 to 2018, indicating that the study area had a VOCs-limited regime. The FNR was greater than 1.0 in summer from 2005 to 2018, indicating that the study area had a transitional regime or NO_x-limited regime. The study area had a transition from VOCs-limited to the transitional regime, especially after 2012, which was consistent with Wang et al. (2019b). 

Fig. 6. (a) The spatial distribution of FNR and (b) the annual variation in FNR (means for the whole year, summer and winter) over the study area during 2005–2018. The dotted lines are the corresponding trend lines.

The fluctuations in the SOC, HCHO column and NO₂ column over the study area in the whole year, summer and winter from 2005–2018 are described in Figs. 7(a), 7(b) and 7(c), respectively. And the linear regression results for the annual variations of surface ozone, HCHO column and NO₂ column in the whole year, summer and winter shown in Fig. 7 were shown in Table 3. In the whole year, the SOC and HCHO column were increasing, while the NO₂ column showed an increasing first and then decreasing tendency. The average SOC annually was 73 µg m^–3 in 2005–2011 which became 98 µg m^–3 in 2012–2018, indicating that the growth of the SOC has been more significant since 2012. The growth rate of the HCHO column in 2012–2018 was higher than that in 2005–2011. And the rate of decrease of the NO₂ column in 2012–2018 was much higher than the growth rate in 2005–2011, which may have been related to the 12^th Five-Year Plan for energy conservation and emission reduction issued by China in August 2012. Li et al. (2019) also indicated that anthropogenic NO_x emissions in China were estimated to have decreased by 21% during 2013–2017, whereas VOCs emissions changed little. Moreover, the reduction rate of the NO₂ column was greater than the growth rate of the HCHO column in 2012–2018. It was concluded that effect of the HCHO column on the SOC was stable during the period of the study, while the influence of the NO₂ column on the SOC in 2012–2018 was much greater than that in 2005–2011. In 2005–2018, underlying the continuous increase in the VOCs concentration, the accelerated decrease in the NO2 has led to an increment in the SOC especially since the winter of 2012. The conclusion was consistent with Li et al. (2019) study which indicated that decrease in anthropogenic NOx emissions could increase ozone in urban areas where ozone production is expected to be VOCs-limited and Lu et al. (2018) study which indicated that either decreasing NOx levels or increasing VOCs levels could potentially enhance ozone pollution for those VOCs-sensitive regions. 

Fig. 7. Annual variations in the (a) SOC, (b) HCHO column and (c) NO₂ column over the study area in the whole year, summer and winter during 2005–2018. The break lines are the corresponding mean SOC. The dotted lines are the corresponding trend lines. 

In summer, surface ozone is driven by photochemical production. The average SOC annually was 137 µg m^–3 in 2005–2011 which became 148 µg m^–3 in 2012–2018. The growth rate of SOC in 2012–2018 was slightly higher than that in 2005–2011. HCHO column in 2012–2018 showed an increasing trend compared with 2005–2011 while NO₂ column in 2012–2018 showed a decreasing trend compared with 2005–2011. With the continuous growth of HCHO, more NO₂ was consumed for stronger photochemical reaction to produce higher ozone in 2012–2018. In winter, titration by NO become more important driver of surface ozone when photochemistry production is not too strong. The average SOC annually was 17 µg m^–3 in 2005–2011 which became 52 µg m^–3 in 2012–2018, indicating that the growth of the SOC has been more significant since 2012. NO₂ column in 2005–2011 showed an increasing trend which indicated the NO titration effect was strong to consume more surface ozone in 2005–2011. In 2012–2018, the sharp decreasing of NO₂ weakened the NO titration effect and led to the accumulation of surface ozone.

In addition, the interannual variability of surface ozone could be influenced by meteorological conditions, for example, the hotter and drier weather conditions which favored ozone production and led to higher ozone levels (Lu et al., 2018). Li et al. (2019) indicated that the temperature and ultraviolet (UV) radiation were increased while surface air relative humidity was decreased over BTH during the summer of 2013–2017, which created favorable meteorological conditions for surface ozone growth after 2012. Reduced PM_2.5 and its precursor emissions levels also caused an increase in the level of ozone due to aerosol chemistry particularly in the North China Plain where PM_2.5 concentrations were highest (Lu et al., 2018; Li et al., 2019). Li et al. (2019) indicated that observed average PM_2.5 levels in summer during 2013–2017 decreased by 41% for BTH. The study inferred that the increasing in surface ozone since 2012 may be related to the significant decrease of PM_2.5. Land cover was also an important factor especially vegetation coverage which could increase the total VOCs. Wang et al. (2018) indicated that forest coverage rose from 20.6 to 35.8% from 1998 to 2013 in Beijing which may lead to the increase in the VOCs concentration fitting our study.

Although the reduction of NO₂ and PM_2.5 led to an increase in the SOC to varying degrees, we still need to continue to strengthen the reduction of the NO₂ and VOCs co-emissions to overcome the effect of decreasing PM_2.5 whether in the whole year, summer or winter. Literature reports have indicated that only when the ratio of VOCs/NO_x is reduced by more than 2:1, does the effective reduction of surface ozone pollution is possible in the study area (Wang et al., 2019b). That is, by setting the SOC as the air quality monitor index, control of NO_x emission need to be coherent with faster reduction of VOCs for mitigating surface ozone pollution.

CONCLUSIONS

A GWR model was developed to estimate the SOC between 2005 and 2018 in the BTH region by establishing a regional spatial linear relationship between ground-based and satellite-retrieved spatio-temporal data sets for the DMA8 ozone from May 2014 till December 2017. Based on the modeling results, the increase and the spatio-temporal variations in the SOC were analyzed. The main conclusions are as follows:

1. The spatio-temporal verification of the satellite-retrieved results confirmed that they were accurate.
   
   
2. The SOC during the study period averaged 86 µg m^–3, 126 µg m^–3 and 45 µg m^–3 over the entire year, during the warm season and during the cold season, respectively, with the highest and the lowest concentrations occurring in the southeast and the northeast, respectively. The average SOC growth rates for the study area were 4 µg m^–3 y^–1, 2.3 µg m^–3 y^–1 and 4.6 µg m^–3 y^–1 for the entire year, the warm season and the cold season, respectively, with the highest and the lowest rates occurring in the northeast and the southeast, respectively.
   
   
3. Underlying the continual rise in the VOCs levels, the accelerated decrease in the NO₂ led to an increase in the SOC from 2005 till 2018, especially beginning the winter of 2012. Compared to 2005–2011, the HCHO column for 2012–2018 showed an increasing trend, whereas the NO₂ column for the same period of time showed a decreasing trend. During summer, the production of surface ozone was driven by photochemical production. With the continual increase in the HCHO, a larger amount of NO₂ was consumed in stronger photochemical reactions, which resulted in higher ozone concentrations between 2012 and 2018. During winter, when the photochemical production was weaker, titration by NO become a greater factor in the surface ozone concentrations. The NO₂ column between 2005 and 2011 displayed an increasing trend, indicating a stronger NO titration effect, which reduced surface ozone production during these years. Between 2012 and 2018, the sharp decrease in the NO₂ weakened the NO titration effect and led to the accumulation of surface ozone.
   
   
4. Strengthening reduction measures for VOCs and NO₂ co-emissions can effectively mitigate the surface ozone pollution.

DECLARATION OF INTEREST

There are no conflicts of interest to declare.

FUNDING

This work was supported by the National Natural Science Foundation of China (41601457), National Key R&D Program of China (2016YFC0206001), China Postdoctoral Science Foundation (2018M641531), National research program for key issues in air pollution control (DQGG0101-02).

ACKNOWLEDGEMTNS

This work was supported by the National Natural Science Foundation of China (41601457), National Key R&D Program of China (2016YFC0206001), China Postdoctoral Science Foundation (2018M641531), National Research Program for Key Issues in Air Pollution Control (DQGG0101-02). We thank NASA (https://mirador.gsfc.nasa.gov) for providing the ozone profile data set (OMO3PR 003), and the CNEMC for providing the air quality data.