Long-term changes of regional ozone in China: implications for human health and ecosystem impacts

Ozone is a ubiquitous trace gas in the stratosphere and troposphere, with the stratosphere and troposphere containing about 90% and 10% of atmospheric ozone, respectively. Stratospheric ozone protects life on the Earth’s surface from harmful solar ultraviolet radiation. Tropospheric ozone is one of the key greenhouse gases (IPCC, 2013; Monks et al., 2015) and very important in sustaining the atmospheric oxidation capacity (Monks et al., 2005; Lelieveld et al., 2008). Surface ozone is a pollutant detrimental to human health and vegetation (REVIHAAP, 2013; CLRTAP, 2016). Atmospheric ozone can only be formed in photochemical reactions. In the stratosphere, ozone is produced and sustained by the Chapman mechanism and catalytically destroyed by active species, such as OH, NO, and halogenated species (Solomon, 1999). A simulation study by Lelieveld and Dentener (2000) suggests that more than 85% of tropospheric ozone originates from photochemical oxidation of volatile organic compounds (VOCs), CO, and CH₄, catalyzed by NO_x, with the rest from stratosphere-to-troposphere exchange (STE). However, the stratospheric contribution to tropospheric ozone has been debated over the years (e.g., Tarasick et al., 2019a and references therein). It has large regional and seasonal variations. For example, Tarasick et al. (2019a) estimated that relative contributions of stratospheric ozone to the total tropospheric ozone column over eighteen North American sites varied in the range of 14.2%–30.4%, with an average of 19.5%. They also indicated that the STE contribution to tropospheric ozone was much larger in winter and spring than in other seasons and influenced less the ozone level in 0–1 km.

Increased emissions of ozone precursors (NO_x, VOCs, CO, CH₄, etc.) since the industrial revolution have caused a large increase in tropospheric ozone, exerting an additional +0.40 (±0.20) W/m² of radiative forcing on the troposphere (IPCC, 2013; Tarasick et al., 2019b; Yeung et al., 2019). Of particular concern is the impact on surface ozone from increased emissions of anthropogenic ozone precursors. Rapid increases in surface ozone were observed different time windows at many urban, rural, and background sites in North America, Europe, and East Asia (e.g., Vingarzan, 2004; Tang et al., 2009; Wang et al., 2009; Cooper et al., 2012; Cooper et al., 2014; Simon et al., 2015; Strode et al., 2015; Ma et al., 2016; Sun et al., 2016; Xu et al., 2016; Gaudel et al., 2018; Jaffe et al., 2018; Wang et al., 2019). Considering the impacts of ozone on human health, ecosystem, and climate, it is necessary to know the spatiotemporal variations of ozone. Long-term observation is an inevitable way of obtaining such spatiotemporal variations. So far, well designed and relatively dense networks for long-term monitoring of surface ozone have been setup mainly in North America, Europe, and part of East Asia (e.g., Japan). In other regions of the world, however, such networks are either nonexsitent or only recently established (e.g., China). Available long-term measurements of surface ozone at many European and North American sites showed rapid increases until about 2000 and slight decreases or a leveling off after 2000 (Derwent et al., 2010;Cooper et al., 2012; Cooper et al., 2014; Cristofanelli et al., 2015; Chang et al., 2017; Lin et al., 2017; Boleti et al., 2018; Derwent et al., 2018; Yan et al., 2018a, 2018b). With effective reductions in NO_x emissions, the high-end ozone values have been clearly decreasing in North America. For example, data from US EPA show that annual 4th highest daily maximum 8-h (A4MDA8) ozone levels declined in the 1980’s, leveled off in the 1990’s, and showed a notable decline after 2002 (https://www.epa.gov/air-trends/ozone-trends). Some other ozone concentration metrics associated with higher ozone concentrations also showed declining trends (Lefohn et al., 2017; Simon et al., 2015; Strode et al., 2015; Jaffe et al., 2018). The reduction of NO_x emissions also caused increases of lower ozone concentrations. Decreasing trends of high ozone concentrations generally occurred during the summer season in less urbanized areas, while increasing trends of low ozone concentrations occurred in winter in more urbanized areas (Simon et al., 2015). In some regions of East Asia, however, continuous increases in the concentrations of surface ozone have been observed also after 2000 (Tang et al., 2009; Wang et al., 2009; Ma et al., 2016; Sun et al., 2016; Xu et al., 2016; Lu et al., 2018). It is noteworthy that the calculations of surface ozone trends are sensitive to subsets of data from a single site or networks. Daytime ozone better represents the well-mixed boundary layer condition and is less subject to nearby sources/sinks, while nighttime ozone is highly impacted by chemical sinks like NO titration and dry deposition hence less representative (Tarasick et al., 2019b). Ozone trends in different seasons can be substantially different as emissions of ozone precursors, photochemical conditions, transport, dry deposition, etc., all may have large seasonal variations. For example, many sites in Europe and North America show opposite trends in daytime ozone for winter (December–February) and summer (June–August) during 2000–2014, with positive and negative trends being found for winter and summer, respectively (Gaudel et al., 2018). Focusing on daytime ozone in summertime (April–September) 2000–2014, Chang et al. (2017) showed that the regional mean daytime ozone decreased significantly (–0.30 ppb/yr, p < 0.01) in eastern North America, changed little (–0.04 ppb/yr, p = 0.09) in Europe, and increased rapidly (0.40 ppb/yr, p < 0.01) in East Asia (Japan and South Korea). Jaffe et al. (2018) made an assessment of background ozone over the USA and found that summertime ozone decreased at many rural sites in the western USA during 2000–2016. This result differs from those of Cooper et al. (2012), who focused on somewhat different period (1990–2010). Changing emissions of ozone precursors are thought to be responsible for the differences in ozone trends (Zhang et al., 2016) though climate variability also has important impact on the variation of surface ozone (Barnes et al., 2016; Okamoto et al., 2018; Xu et al., 2018; Sun et al., 2019). In addition, variations of aerosol level may exert influences on surface ozone by changing photolysis rate and heterogeneous reactions (Qu et al., 2018). A simulation study by Li et al. (2019) suggests that the reduction of fine particulate matter (PM_2.5) was the most important cause of the ozone increase observed during 2013–2017 in the NCP.

The above-mentioned studies on ozone trends are based on surface ozone data dating back to the 1990s or 1980s. Earlier measurements of surface ozone are far more sporadic and often subject to large uncertainties. Recently, Tarasick et al. (2019b) reports results of careful, synthesized analyses of global surface and free tropospheric ozone. They summarized global ozone observations and techniques used from 1870s through 2016 and developed and applied a set of four criteria to select data for historical reconstruction. Although the reconstructed data are still associated with large uncertainties, they found robust evidence for increases in daytime surface ozone in the temperate and polar regions of the north hemisphere of 30–70% between the period of historic observations (1896–1975) and the modern period (1990–2014) and no evidence for a significant increase in surface ozone in the southern hemisphere (Tarasick et al., 2019b). In addition, these authors found no evidence of a change in rural background ozone during 1896–1975, suggesting that the ozone increases in northern temperate and polar regions occurred mainly after the middle or later 1970s. This is in sharp contrast with some previous studies (e.g., Low et al., 1990; Staehelin et al., 1994;Staehelin et al., 2001). Note that most of the historic observations were made in Europe, the origin of the industrial revolution.

To better answer the scientific questions related with tropospheric ozone the International Global Atmospheric Chemistry Project (IGAC) developed theTropospheric Ozone Assessment Report (TOAR): Global metrics for climate change, human health, and crop/ecosystem research (http://www.igacproject.org/activities/TOAR). This is the first international initiative aiming to provide the research community with an up-to-date scientific assessment of tropospheric ozone’s global distribution and trends from the surface to the tropopause. After a few years of efforts, TOAR has established a database of global surface ozone observations at Forschungszentrum Jülich in Germany (https://join.fz-juelich.de/accounts/login/) and has published to date its findings in the TOAR special feature of Elementa: Science of the Anthropocene, TOAR-Database (Schultz et al., 2017), TOAR-Metrics (Lefohn et al., 2018), TOAR-Health (Fleming et al., 2018),TOAR-Vegetation (Mills et al., 2018), TOAR-Climate (Gaudel et al., 2018), TOAR-Model Performance (Young et al., 2018), andTOAR-Observations (Tarasick et al., 2019b).

In the TOAR publications, measurements of surface ozone in Asia are mainly from Korea and Japan, where dense monitoring networks for surface ozone and other pollutants have been operated for two decades or longer. Long-term observations of surface ozone have been lacking in other parts of Asia. Even in mainland China, a densely populated world region with high intensity of pollutants emissions and agricultural activities, ozone data from a dense air pollution monitoring network were not available until 2013 (Lu et al., 2018). Because of the shorter time-series of ozone data most of the Chinese sites were not included in previous assessments of long-term trends in ozone metrics in the set period 2000–2014, except for several sites from Hong Kong and Taiwan and one site from mainland China (Chang et al., 2017;Fleming et al., 2018; Mills et al., 2018; Gaudel et al., 2018). On the other hand, there have been some sporadic reports on long-term changes of surface ozone at individual sites in China (Shao et al., 2006; Xu et al., 2008; Tang et al, 2009; Wang et al., 2009; Zhang et al., 2014; Ma et al., 2016; Sun et al., 2016; Xu et al., 2016; Gao et al., 2017; Wang et al., 2019) and of tropospheric ozone over the North China Plain (NCP) region (Ding et al., 2008; Wang et al., 2012; Wang et al., 2017). Most of these studies are based on shorter or poor continuity ozone time-series and provide only estimates of changes in ozone concentrations rather than results about common ozone metrics relevant to human health and vegetation. Recent work by Lu et al. (2018) contains trend estimates of some human health and vegetation metrics based mainly on ozone measurements during 2013–2017 from about 1600 urban and suburban sites in the China National Environmental Monitoring Center (CNEMC) Network. This work revealed striking results that China has been experiencing rapid increases in surface ozone concentrations and that ozone air quality nowadays has been as poor as it was in the United States in the early 1980s. However, the evolution of surface ozone before this recent period and regional differences remain unknown.

In this study, we present results from an integrated analysis of long-term measurements of surface ozone at six World Meteorological Organization’s Global Atmosphere Watch (WMO/GAW) sites in mainland China and one rural and one urban sites in the NCP, one of the most polluted regions in China. In addition to average concentrations, we include plentiful results about present-day values and long-term trends of ozone metrics, which can be referenced in the assessment of ozone impacts on human health and vegetation. We determine the metrics and their trends following the TOAR methodologies and data completeness criteria. We have not been able to adopt the fixed period 2000–2014 used by TOAR for trend calculations, and instead accept what is available from the time-series. We also compare the average seasonal and diurnal cycles and present-day metrics values of ozone among regions and sites. Our goal is to present information about trends and present-day results for China that are not available in the current TOAR or other publications.

Ozone measurements analyzed in this study are mainly from six background stations in mainland China. The long-term observations of surface ozone at the stations have been maintained by China Meteorological Administration (CMA) as a part of contribution to the GAW programme of WMO. The names of the six stations are Mt. Waliguan (WLG), Shangdianzi (SDZ), Lin’an (LA), Longfengshan (LFS), Xianggelila (XGLL), and Akedala (AKDL). As comparison, ozone measurements from a rural site (Gucheng, GCH) and an urban site (CMA campus, CMA) are also included in the analysis. The locations of all the sites are shown in Figure 1 and more details about the sites are given in Table 1. Figure 1 also shows trend slope values of tropospheric NO₂ column during 2005–2019. Additionally, averages and trends of tropospheric NO₂ column over grids (0.5° × 0.5°) covering different sites and some other metadata are given in Table S1 in the Supplementary Materials. The WLG station is located on top of Mt. Waliguan in Qinghai Province, western China, and is one of the 31 global baseline GAW stations (http://www.wmo.int/pages/prog/arep/gaw/measurements.html). SDZ, LA, LFS, XGLL, and AKDL are all regional background GAW stations. SDZ is located about 100 km northeast of urban Beijing. LA is located in the Yangtze River Delta (YRD) region, and about 50 km west of Hangzhou and 210 km southwest of Shanghai. LFS is located in Heilongjiang Province, Northeast China, and about 140 km southeast of Harbin. XGLL is located in Yunnan Province, southwest China, and about 450 km northwest of Kunming. AKDL is located the Gobi area of Xinjiang, Northwest China. CMA is a typical urban site in Beijing, located between the 2nd and 3rd ring roads of the city. GCH is a polluted rural site located in Hebei Province and about 110 km southwest of Beijing. More details about the sites can be found in publications (Lin et al., 2008, 2009a,2009b, 2010, 2011; Ma et al., 2014; Xu et al., 1998, 2008, 2016, 2018; Jin et al., 2016). All GAW sites are far away from direct anthropogenic impacts. However, the GAW sites in east part of China (SDZ, LA and LFS) are more subject to impacts of regional air pollution than those in west part of China (WLG, XGLL and AKDL), as indicated by the tropospheric NO₂ data in Table S1.

Surface ozone at all sites has been observed using ultraviolet photometric technique. Details of instruments used at the sites are given in Table 2. Parallel measurements of ozone were made at WLG using dual ozone analyzers. Two ozone analyzers (TE 49, Thermo Environmental Instruments Inc.) were deployed at the beginning of the ozone observation at WLG (Xu et al., 2016). One of the TE 49 analyzers was replaced with a new model (TE 49i, Thermo Fischer Scientific Inc.) on 21 May 2011. The ozone analyzers at WLG have been calibrated every 3 months using an ozone calibrator (TE 49PS, Thermo Environmental Instruments Inc.) maintained at the station. Ozone at other sites has been monitored using different models of analyzers. TE 49C ozone analyzers (Thermo Environmental Instruments Inc.) were used at SDZ (Lin et al., 2008), LA (Xu et al., 2008), GCH (Lin et al., 2009a) and CMA (Lin et al., 2011), while EC 9810B ozone analyzers (Ecotech PTY LTD) at LFS (Lin et al., 2009b), XGLL (Ma et al., 2014) and AKDL (Lin et al., 2010). These analyzers have been calibrated at least every 3–6 months against a travelling standard (ozone calibrator, TE 49CPS, Thermo Environmental Instruments Inc.) maintained by Key Laboratory for Atmospheric Chemistry of China Meteorological Administration (KLAC/CMA). The ozone analyzers and calibrator used at WLG were audited and calibrated in 1994, 1995, 2000, 2004, and 2009 by experts from the WMO World Calibration Centre for Surface Ozone and Carbon Monoxide, EMPA Dübendorf, Switzerland. The audit results indicate excellent or good agreement between the instruments at WLG and transfer standard from EMPA (Zellweger et al., 2000, 2004, 2009). The ozone calibrator at KLAC/CMA was also compared with the EMPA transfer standard and excellent agreement between both instruments was achieved. Therefore, ozone measurements from our sites are traceable to the Standard Reference Photometer (SRP) maintained by EMPA.

Raw data of ozone measurements are stored as 5 min averages, from which hourly averages are calculated in post-processing. However, only those hourly data with at least 45 min measurements are considered valid and used in this paper. This fulfils the data capture criterion of 75% for TOAR (Lefohn et al., 2018). The ozone time-series analyzed in this study cover the periods from observation start times to the end of 2016 for the background sites and to the end of 2018 for GCH and CMA. The time-series of valid ozone data from SDZ has the best data completeness (97.6%), followed by those from LA (94.5%), LFS (90.8%), GCH (89.6%), WLG (89.2%), XGLL (88.8%), AKDL (82.9%), and CMA (81.0%).

Twelve ozone metrics are calculated from the ozone measurements. These metrics allow for characterizing the spatial differences and temporal variations of ozone from different aspects. Detailed descriptions of these metrics are given in Table 3. The annual average (AAVG) and seasonal average (SAVG) represent the arithmetic means of all available hourly ozone concentrations over a year and a season, respectively. Annual daytime (0800–1959 local time) average ozone (AdAVG) and seasonal daytime average ozone (SdAVG) are used to represent annual and seasonal average ozone levels under well-mixed boundary layer conditions, not influenced by localized nighttime deposition and destruction processes (Chang et al., 2017; Gaudel et al., 2018; Lu et al., 2018). These metrics can be used in the assessment of ozone impacts on climate change and in model evaluation. The annual maximum daily 8-h average ozone (AmaxMDA8), and the annual 4th highest daily maximum 8-h ozone (A4MDA8) are metrics related to high ozone values. They are important metrics for assessing human health acute effects from ozone exposure (Lefohn et al., 2017, 2018; Fleming et al., 2018). A4MDA8 corresponds about to the 99th percentile of the MDA8 values in a year and is used in the US for determining compliance with its National Ambient Air Quality Standards for Ozone (US Federal Register, 2015). AmaxMDA8 focuses on the highest MDA8 value in a year and can be used to show the changes of extreme ozone values at different sites. This ozone metric, however, is subject to annual variability of meteorological and emission conditions. Therefore, shorter period variations of AmaxMDA8 should not be over-interpreted. The annual sum of the positive differences between the daily maximum 8-h ozone value and 35 ppb (SOMO35) is a metric for assessing cumulative health impacts from chronic exposure to medium and high ozone levels (Lefohn et al., 2017, 2018; Fleming et al., 2018). In addition to annual SOMO35, the seasonal sum of the positive differences between the daily maximum 8-h ozone value and 35 ppb (SSOMO35) were also calculated to show seasonal differences of cumulative ozone exposure. The annual number of days of MDA8 > 90 ppb (NDGT90) and annual number of days of MDA8 > 70 ppb (NDGT70) are metrics of ozone exceedances (Fleming et al., Lefohn et al., 2018). They are health-related ozone metrics and often used in policy making. The sum of all hourly ozone values over 40 ppb for all daylight hours over 3-months (AOT40) and the sigmoidally weighted sum of hourly ozone values over 3-months and daytime (W126) are two established metrics for assessing the ozone impacts on vegetation (Lefohn et al., 2018; Mills et al., 2018). In the calculations of AOT40 and W126, we did not focus on the growing season of certain type of vegetation in a particular climatic zone, which was done in the systematic global assessment (Mills et al., 2018). Instead, we present AOT40 and W126 values for all seasons (i.e., March–May, June–August, September–November, and December–February) so that researchers can get AOT40 and W126 results for three months, six months, and the whole year, which correspond approximately to the growing periods of interested vegetation types. It is suggested that the TOAR database and tools are used to obtain accurate matches of ozone metrics for vegetations types.

All quality-controlled hourly data available from our sites are included in the metrics calculations. However, ozone data capture is usually more or less incomplete. Consequently, the cumulative ozone metrics (i.e., SOMO35, SSOMO35, AOT40, and W126) calculated from the available ozone measurements may be underestimated. The underestimation of the metrics may vary with the extent of undersampling, leading to inconsistent comparisons between periods and sites. To minimize the effect of this problem, the metrics values are divided by the corresponding data completeness to obtain modified or “corrected” metrics values, i.e.,

where, M_mod and M are metrics values after and before modification, respectively; P_t and P_r are theoretic and real data points within the period (a season or year). After such modification the cumulative metrics values should be more close to the real ones. However, the data completeness in some cases was poor due to technical problems so that the modified cumulative metrics (and other metrics) may be with high uncertainties. Therefore, all metrics values associated with data completeness <50% are not used in further analysis. Note that here our threshold of data capture is relaxed from the TOAR requirement (75%, Lefohn et al., 2018). The main purpose of the relaxation is to include as many reliable data points as possible in the trends analysis. This is a compromise between the accuracy of each data point and the number of data points used in trends calculations. The compromise is necessary since most of our ozone time-series are not as long as 15 years, which is the time scale, required by TOAR in long-term trends calculations (Fleming et al., 2018; Mills et al., 2018). On the other hand, calculated metrics are discarded if large biases are introduced by the relaxation of data completeness. For example, all values of annual ozone metrics of CMA for the year 2016 (with data completeness about 52%) are not used in this study since the measurements were mainly from the low ozone seasons (fall and winter).

To study the long-term changes in the ozone metrics, trend analysis is performed on the annual and seasonal metrics. Following the approaches of TOAR (Lefohn et al., 2018), the Mann-Kendall (M-K) test is used to find the trends in AAVG, AdAVG, SAVG, SdAVG, AmaxMDA8, A4MDA8, SOMO35, SSOMO35, AOT40, and W126. The M-K test needs no assumption of data distribution and is resistant to the influences of outliers (Kendall, 1955). The M-K test is performed using a self-written IGOR macro, which estimates the trends using the Theil-Sen estimator (Sen, 1968) and computes Kendall’s τ, p value, etc. As NDGT90 and NDGT70 are count metrics, for which problematic results may be produced by the M-K test (Lefohn et al., 2018), the linear trend analysis is performed on these two metrics using IGOR’s linear fit function.

As mentioned in section 2.4, the cumulative metrics values are corrected for data capture and those associated with data completeness lower than 50% are abandoned in the analysis. In addition, the annual ozone metrics in 2016 for the CMA site are also abandoned because the majority of the measurements were from the cold seasons though the data completeness in this year was 51.9%.

To see the annual variability of ozone at different sites, 95-, 75-, 25 and 5-percentiles as well as averages and medians were calculated from hourly ozone values in each year and are displayed in Figure 2. Note that some annual statistics are not shown in Figure 2 due to low data completeness (<50%). As shown in Figure 2, the longest ozone record is from WLG, with the annual statistics covering 1995–2016. The annual average (median) of ozone at WLG varied from 47.9 ppb (42.1 ppb) in 1997 to 54.1 ppb (52.1 ppb) in 2016. These represent the highest annual ozone levels among all sites. The WLG site is located on top of the high mountain (3816 m a.s.l.) and subject to influence of the free troposphere (Ma et al., 2002; Xu et al., 2016). And the site is within a region with very low emissions of ozone precursors, as indicated by the very low tropospheric NO₂ column (0.025 DU, Table S1). Therefore, the ozone level at the site is mainly due to high ozone background rather than local or regional anthropogenic activities. The annual average (median) of ozone at XGLL (a high altitude site in southwest China, 3580 m a.s.l.) varied from 31.8 ppb (35.1 ppb) in 2010 to 42.9 ppb (44.3 ppb) in 2008. This site is also significantly impacted by air mass downward transport (Ma et al., 2014). Ozone observation at AKDL started in November 2009. Annual statistics show that the smallest and largest annual average (median) was 29.2 ppb (28.9 ppb) in 2015 and 39.7 ppb (40.3 ppb) in 2010, respectively. Among the three background sites in eastern part of China, SDZ has observed the largest annual ozone averages (33.1 ppb – 40.4 ppb). This reflects the fact that the NCP, in which SDZ is located, has been most impacted by photochemical pollution (Lu et al., 2018). However, much lower annual average ozone was observed at the urban site (CMA, 23.0 ppb – 33.8 ppb) and the rural site (GCH, 18.8 ppb – 31.3 ppb) in the NCP. This can be explained by the large fractions of lower ozone values (Figure2g and 2h) mainly caused by NO titration (Lin et al., 2009; Lin et al., 2011). Tropospheric NO₂ column densities over CMA and GCH were the highest (0.663 DU) and second highest (0.560 DU) among all the sites (Table S1), suggesting strong consumption of surface ozone by NO titration at these sites. Annual ozone averages observed at LA (31.1 ppb – 35.4 ppb) in the YRD and LFS (29.4 ppb – 39.0 ppb) in Northeast China were at intermediate levels.

Compared with the sites in western China (WLG, AKDL and XGLL) and Northeast China (LFS), the sites in the NCP (SDZ, GCH and CMA) and the YRD (LA) are more often impacted by photochemical pollution and NO emission, as indicated by the very high 95th percentiles and very low 5th percentiles of ozone value in Figure2e–2h. The values of tropospheric NO₂ column (Table S1) over these sites were much higher than those over sites in western and northeastern China, indicating stronger influences of local and regional anthropogenic emissions. In addition, the SDZ, GCH, CMA, and LA sites also show larger differences between the average and median values, suggesting that there is a significant deviation of ozone frequency distribution from the normal distribution. Such deviation can be attributed to the strong titration influences on surface ozone at these sites. Previous studies of ozone from SDZ, GCH, and CMA all showed that ozone frequency distributions could be fitted with two Lorentz curves, with one having higher peak values and being more close to normal distribution, and the other having much lower peak values (close to zero) and deviating considerably from a normal distribution (Lin et al., 2008; Lin et al., 2009; Lin et al., 2011). These results suggest that the low ozone data can be important in generating skewed ozone distributions. These low ozone levels were mainly observed during nighttime and in cold seasons, when the NO level was high and the boundary layer was shallow.

All valid monthly ozone averages for each site were averaged for different months to obtain the average seasonal cycle of ozone at the site. Figure 3 shows seasonal cycles for all sites. On average, ozone at WLG peaks in July (61.8 ppb) and has a minimum in January (43.0 ppb). Similar seasonal cycles are observed at LFS, SDZ, GCH, and CMA. These sites are all located at higher latitudes (higher than 36°N). AKDL shows a quite different seasonal pattern of ozone, with a rapid increase from January to March and a slow decrease from March to December. Further inspection shows that the peak of monthly ozone at AKDL varied largely (from February to May) from year to year. Therefore, including more measurements (i.e., those after 2016) from this site in the calculation may produce an average seasonal ozone pattern different from that shown in Figure 3. The seasonal patterns of ozone at the sites in south part of China (LA and XGLL) have common features with those in north part of China but also show their uniqueness. Similar to ozone at the sites in the north, ozone at XGLL and LA shows an increase after January and peaks in May and June, respectively. The level of ozone at XGLL decreases rapidly in summer to its annual minimum (24.0 ppb), stay low from August to October, and increases after October. The level of ozone at LA also shows a decrease in summer but it stays higher in summer and fall than that at XGLL. The lower ozone values observed in summer and fall at XGLL and LA are mainly caused by the Asian summer monsoon, which transports cleaner marine air masses to the sites and is associated with meteorological conditions unfavourable for photochemical reactions (Wang et al., 2001; Xu et al., 2008; Ma et al., 2014).

It is mentioned in section 3.1 that WLG has the highest annual ozone levels among all sites. Data in Figure 3 indicate that among all sites, WLG has the highest average concentrations of ozone in all months, with only one exception (XGLL in April). The monthly average ozone levels at the NCP sites (SDZ, GCH and CMA) are substantially enhanced during summer. Nevertheless, they are lower than those at WLG. In addition to the much lower altitudes of these sites, stronger consumption of ozone by NO titration during dark and/or cold periods, and larger dry deposition of ozone over better vegetated areas, should be the main reasons of lower average ozone levels. These ozone sinks, together with photochemical production, cause larger diurnal (see section 3.3) and seasonal variations of ozone at the low altitude sites. In other word, ozone produced in local boundary layer is partially consumed by ozone sinks. The ozone source (photochemical formation) and sinks (NO titration, dry deposition, etc.) dominate in terms of short-term fluctuations of ozone at level around its background values, and contribute to longer-term changes when they are not balanced with each other. As a high-mountain site, WLG is largely impacted by free tropospheric air and local sources and sinks contribute minor to the ozone level there (Ma et al., 2002). Therefore, the monthly and annual average concentrations of ozone at WLG are mainly determined by background ozone, which is more influenced by the free troposphere. Free tropospheric ozone makes also large contributions to monthly averages of ozone at XGLL, another high-mountain site, particularly in winter and spring (Ma et al., 2014). In summer and early fall Asian Monsoon exerts strong impacts on ozone at XGLL and LA, causing lower ozone levels.

If only the daytime measurements of ozone are considered, the situation looks different. A comparison between the left and right plots in Figure 3 indicates that monthly averages of daytime ozone at SDZ, GCH, CMA, and LA are much larger than those of entire day ozone though the differences for WLG, XGLL, AKDL, and LFS are minor or small. On average, the daytime levels of ozone at GCH, CMA, LN, SDZ, LFS, AKDL, and XGLL are respectively 45.5%, 32.9%, 24.8%, 22.6%, 11.2%, 9.8%, and 4.1% higher than those for entire-day ozone, while the daytime level of ozone at WLG is about 1.6% lower than the entire-day level. Compared with the monthly values for WLG, SDZ, CMA, GCH, LA. and XGLL all have some monthly values larger than those for WLG. This occurs mainly in the warmer months (May–September). As daytime ozone is influenced more by photochemical production and less by NO titration, it reflects better the impacts of photochemical pollution than entire day ozone. Therefore, the results here confirm that ozone at sites in eastern part of China (GCH, CMA, LN, SDZ, and LFS) is much stronger influenced by than photochemical pollution than those in western China (WLG, AKDL, and XGLL). This is consistent with the distribution of tropospheric NO₂ (Table S1).

The patterns of the seasonal cycles of ozone at SDZ and GCH look very similar to each other though the monthly ozone values of GCH are lower than those of SDZ (10.5 ppb on average if entire-day measurements are considered). The average ozone values of CMA are very close to those of GCH in most months but there are larger CMA-GCH differences in May (5.3 ppb) and September (8.3 ppb). The sites GCH, CMA, and SDZ, align along the southwest-to-northeast line (Figure 1), which is a typical route for the transport of air pollutants in the NCP. Since SDZ is located on the northeast edge and downwind of the main NCP region, air masses transported from the polluted upwind region may exert large impacts on air quality at the site. This is particularly true for ozone, which can be significantly formed on the way of transport (Lin et al., 2008; Xu et al., 2011; Ge et al., 2012). A study by Lin et al. (2008) shows that pollution transport contributes about 22 ppb yearly to surface ozone at SDZ and about 29 ppb in summer. On the other hand, the strong titration reaction with NO can remove much of ozone at CMA and GCH, particularly during dark and cold periods (Lin et al., 2009; Lin et al., 2011). Therefore, it is not surprising that the average ozone levels at SDZ are mostly higher than those at CMA and GCH.

It should be noted that the multi-year averaged seasonal cycles shown in Figure3 do not necessarily agree with those for individual years. The large differences in meteorological conditions and/or precursor emissions may cause shifts in the ozone maximum and minimum, as shown in Xu et al. (2008) and Ma et al. (2014).

Multi-year average diurnal cycles of ozone are calculated from the valid hourly measurements from the sites and shown in Figure 4. The most striking phenomenon is the contrast between the diurnal patterns for WLG and XGLL and those for other sites. Ozone at WLG and XGLL displays flat diurnal patterns with very small variations. However, ozone at other sites shows much larger amplitudes, particularly at SDZ, CMA, GCH, and LA. Both WLG and XGLL are high mountain sites that are less impacted by anthropogenic emissions (see Table S1). Transport of air masses from the free troposphere and even the lower stratosphere makes significant contributions to surface ozone at both sites (Ding et al., 2006; Ma et al., 2014; Xu et al., 2016, 2018). At WLG, the mountain-valley breeze plays an important role in the diurnal-seasonal cycles of ozone, with daytime and nighttime ozone being mainly influenced by boundary layer air and free-tropospheric air, respectively (Xu et al., 2016). A chemical budget analysis by Ma et al. (2002) indicates that the net production of ozone at WLG is small and can even be negative. Therefore, the day-to-night variation in ozone concentration is mainly determined by the difference of ozone level between boundary layer and free-tropospheric air. Consequently, slightly lower ozone concentrations are observed at WLG during daytime. Being a regional background site, XGLL is not as remote as WLG and is subjected to more influences from the emissions in the surrounding areas. The station is practically in a forest with pine and broad-leaved trees. The Xianggelila County has a population of about 0.15 millions. Daytime photochemical production of ozone can be sustained by natural and anthropogenic emissions of VOCs, NO_x, etc. As shown in Figure 4, the ozone level at XGLL is slightly enhanced during the daytime, on average by 7.7 ppb from 0800 to1400.

The amplitudes of the diurnal cycle of ozone at the two northernmost sites, AKDL and LFS, are 10.6 ppb and 15.6 ppb, respectively. This indicates that net photochemical formation of ozone at these sites is significant but not very intense. These sites are more impacted by air pollution but are located far from the pollution centres. All other sites are located either directly in (CMA and GCH) or just on the borders of the heavily polluted regions (SDZ and LA). These are well reflected by tropospheric NO₂ column data in Table S1. The average amplitudes of ozone concentrations for GCH, CMA, LA, and SDZ are 41.6 ppb, 38.0 ppb, 34.2 ppb, and 33.0 ppb, respectively. They all are very large, consistent with the significant impacts of photochemical pollution over these sites during daytime and influences of NO titration and dry deposition on nighttime ozone in the shallow nocturnal boundary layer.

Trends in annual metrics, i.e., AAVG, AdAVG, AmaxMDA8, A4MDA8, SOMO35, NDGT90, and NDGT70 (see detailed description in Table 3), were calculated for each site from all valid annual metrics values using the M-K test. Table 4 lists all the trends, with the trend values associated with p < 0.05 or p < 0.10 being in bold and marked with either “**” or “*”. To show the relative change rates, percent trends in different annual ozone metrics at different sites are presented in Table S2. Note that the trends in annual metrics represent average tendencies of the ozone time-series. As can be seen in Table 4, the trend in a metrics can change from positive at one site to negative at the other, and for some sites, the sign of the trends can also change from one metric to the other.

Since the ozone measurements from different sites do not necessarily cover the same time period, it is not expected that the trend values of a metric are comparable among sites. On the other hand, it is surprising that the same ozone metrics increased at some sites and decreased at the other sites. For example, AAVG and AdAVG increased at WLG, SDZ, and CMA but decreased at the other sites even though all the sites share common observation period 2010–2016. This phenomenon implies that the long-term changes in emissions of ozone precursors, chemical formation and removal of ozone, and meteorological conditions were quite different at different sites.

It can be seen in Table 4 that at some sites (i.e., XGLL, LFS, GCH, and CMA), different ozone metrics changed in different directions. For example, AAVG, AdAVG, and NDGT70 at XGLL changed in negative direction, while the other annual metrics changed in positive direction. At CMA, however, AAVG, AdAVG, SOMO35, NDGT90, and NDGT70 changed in positive direction, while the other metrics changed in negative direction. Such phenomenon was studied by Lefohn et al. (2017) and can be attributed to varying responses of different metrics to changes in ozone concentration distribution. Some ozone metrics are more influenced by high ozone values, such as AmaxMDA8, A4MDA8, NDGT90, etc., others more influenced by moderate and high ozone values, such as AAVG, SOMO35, AOT40, W126, etc. Different shifting patterns in the ozone concentration distribution can cause different changes in these metrics and even in different directions, as seen in Table 4. On the other hand, some of the trends have higher uncertainties as suggested by the large p-values. These uncertainties should be considered in the interpretation of the difference between different metrics.

Trends of AdAVG for SDZ and CMA are obviously larger than those of AAVG, suggesting that the trend of ozone at these sites were mainly resulted from increasing photochemical pollution, which becomes severer during daytime. The trend of AdAVG for WLG is only slightly larger than that of AAVG, reflecting the minor role of local photochemistry. Interestingly, AdAVG of AKDL decreased at –1.30 ppb/yr, which is a more rapid decline than AAVG (–1.00 ppb/yr). Trends of AdAVG for other sites are either more negative (GCH and LA) or slightly less negative (LFS and XGLL) than those of AAVG. For these sites the trends values are quite uncertain as indicated by the large p-values.

The WLG global baseline station has the longest ozone time-series among the eight sites. Figure 5 displays the long-term variations of AdAVG and SOMO35 for WLG during 1995–2016. The annual value in 1998 is not shown due to the low data completeness (48.6%) in this year. Both AdAVG and SOMO35 have some multi-year variations (Figure5). In previous studies we found that there were periodicities of ozone at WLG during 1994–2013, with 2–4-year, 7-year, and 11-year cycles (Xu et al., 2018), which were linked to the periodic variations of the quasi-biennial oscillation (QBO), the East Asian summer monsoon (EASM), etc. The multi-year oscillations in the ozone concentration exert significant impacts on ozone tendencies for shorter time-periods. Due to the length of ozone observation at WLG, we obtained relative stable trends in ozone metrics. In Xu et al. (2016) we obtained a Theil-Sen (T-S) trend in AAVG of 0.21 ppb/yr (p = 0.005) from WLG ozone measurements during 1994–2013. Adding new data, we obtained a T-S trend in AAVG of 0.20 ppb/yr (p = 0.00) and a trend in AdAVG of 0.21 ppb/yr (p = 0.00) from WLG ozone measurements during 1995–2016. These trends are virtually the same considering the slightly different periods. They indicate that the ozone concentration at WLG increased at a moderate rate since the observation began. The rates of change in AAVG and AdAVG were 0.39% and 0.42%, respectively (Table S2). As for the causes of the increase, our study confirms that both anthropogenic emissions in South East Asia and climate variability played roles in the long-term increase (Xu et al., 2018). In response to the increase of ozone concentration, the cumulative metric SOMO35 at WLG also increased (Figure 5b). The rate of increase of SOMO35 is 78.8 ppb-day/yr or 1.09%/yr (p = 0.00) (Tables 4 and S2).

The SDZ site has also experienced a significant increase in ozone concentration since the observation began in 2003. All annual ozone metrics studied show substantial increases. The AAVG, AdAVG, and SOMO35 values for SDZ increased during 2004–2015 at 0.45 ppb/yr (p = 0.01), 0.67 ppb/yr (p = 0.01), and 255.5 ppb-day/yr (p = 0.02), respectively (Figure 6, Table 4). These trends are more than twice as large as those for WLG. And the relative rates of increases in AAVG, AdAVG, and SOMO35 at SDZ reach 1.2%/yr, 1.49%/yr, and 3.17%/yr, respectively (Tables 4 and S2). Besides AAVG, AdAVG, and SOMO35, some metrics mainly influenced by high ozone values also show large increases. As can be seen in Tables 4 and S2, AmaxMDA8 and A4MDA8 at SDZ increased significantly, at a rate of 2.52 ppb/yr (or 1.77%/yr, p = 0.08) and 2.15 ppb/yr (or 1.56%/yr, p = 0.06), respectively. In response to the enhancement of high ozone values, the ozone exceedance metrics NDGT90 and NDGT70 increased at 3.40 days/yr (or 7.5%/yr, p = 0.00) and 3.71 days/yr (or 4.04%/yr, p = 0.00), respectively.

Long-term change of MDA8 ozone at Shangdianzi from October 2003 to June 2015 was analyzed by Ma et al. (2016). A Kolmogorov–Zurbenko (KZ) filter method was applied in the analysis to separate short-term, seasonal, and long-term variations in the MDA8 ozone data. After such treatment, they obtained an average increasing rate of 1.13 ppb/yr (R² = 0.92). This trend is less than half of our AmaxMDA8 trend (2.52 ppb/yr). This difference is reasonable because the MDA8 trend is derived from all the valid MDA8 data, while the AmaxMDA8 trend is calculated from the annual largest MDA8 values.

It is shown that the ozone concentration distribution at SDZ shifted from lower values towards middle and higher values during 2004–2015 (Lefohn et al., 2017). Such a change in ozone concentration distribution is consistent with the trends in ozone metrics at SDZ (Table 4). As to large long-term trend in MDA8, increases in emissions, particularly those of VOCs, are believed to be the major cause (Ma et al., 2016). And the same applies to the increases of surface ozone over some other parts of central eastern China, according to the modelling study of Sun et al. (2019).

Increases in AAVG, AdAVG, and SOMO35 have also observed at CMA, an urban site in the megacity Beijing (Table 4). The AAVG and AdAVG values for CMA have increased at a rate of 0.83 ppb/yr (or 2.96%/yr, p = 0.06) and 0.88 ppb/yr (or 2.32%/yr), respectively (Tables4 and S2). Such positive trends are about 84% (AAVG) and 31% (AdAVG) greater than those for SDZ though the two sites are only about 100 km distant from each other. SOMO35 for CMA has increased rapidly as well (276.4 ppb-day/yr, 3.64%/yr, p = 0.10). These increases in AAVG, AdAVG, and SOMO35 at CMA are the largest among all sites, suggesting that ozone air quality in the urban areas has been most rapidly deteriorated. The metrics related with high ozone values at CMA have changed in inconsistent directions, with AmaxMDA8 and A4MDA8 showing decreases, and NDGT90 and NDGT70 showing increases (Tables 4 and S2). This inconsistency suggests that although the high-end values of ozone at this urban site have been declining, the numbers of hourly ozone values over 90 ppb and 70 ppb have increased. However, all these decreasing and increasing trends have lower certainty, as suggested by the higher p-values.

The ozone level at the AKDL site has been low even during the daytime (see Figure 4). There has been no ozone exceedance since the observation commenced so that the trends in NDGT90 and NDGT70 cannot be calculated. However, all other annual ozone metrics show clear and rapid decrease trends (Tables 4 and S2). Such phenomenon is more or less unexpected. Among all the sites, AKDL had the least number of years of monitoring data available. Each of the annual metrics at this site is based on only six annual values. Therefore, the decreasing trends should be interpreted with caution. On the other hand, the negative trends in these metrics are very clear and with low p-values. As shown in Figure 7, the gradual decrease of ozone at AKDL caused substantial reductions in AdAVG and SOMO35 during 2010–2016. The SOMO35 value, for example, declined at a rate of about 426 ppb-day/yr (p = 0.00) or 17.4%/yr (Table S2) and had even a more than 30% drop within the seven years. At present, it is unknown what kinds of changes in the environment have caused the trends in ozone metrics at AKDL. The annual ozone metrics at LFS, LA, XGLL, and GCH show positive or negative trends but none of these trends are with higher certainty as suggested by the larger p-values (Tables 4 and S2).

We calculated trends in four seasonal metrics, i.e., SAVG, SdAVG, SSOMO35, AOT40, and W126 (see detailed description in Table 3) for all sites and seasons using the M-K test. All the trends in seasonal ozone metrics are given in Table 5, with the trend values associated with p < 0.05 or p < 0.10 being in bold and marked with either “**” or “*”. Percent trends in different seasonal ozone metrics at different sites are presented in Table S3. We can see in Table 5 that seasonal ozone trends with p < 0.10 occurred at all sites except XGLL. However, the seasonal metrics of ozone at LFS, LA, and GCH changed not always in the same direction, leading to annual trends with lower certainty (larger p-values) at these sites (see Table 4).

At WLG, all seasonal ozone metrics showed increasing trends. SAVG, SdAVG, SSOMO35, and AOT40 increased all at the largest rate in fall (0.26 ppb/yr or 0.57%/yr, 0.29 ppb/yr or 0.65%/yr, 25.3 ppb-day/yr or 1.95%/yr, and 233 ppb-h/yr or 3.06%/yr), while W126 showed the largest increase trend in summer (288 ppb-h/yr or 1.30%/yr). Figure 8 presents the variations of SdAVG and AOT40 as well as the calculated seasonal trends in four seasons. The data points of SdAVG and AOT40 in fall are most closely aligned near the trend lines, consistent with the lowest p-values (0.000006 and 0.0002) among all seasons. The trends in SdAVG and AOT40 in spring were 0.21 ppb/yr (0.38%/yr) and 205 ppb-h/yr (1.32%/yr), respectively, only next to the counterparts in fall. The trends in SdAVG and AOT40 in summer were 0.18 ppb/yr (0.31%/yr) and 139 ppb-h/yr (0.66%/yr), respectively, with the absolute rates of change ranking third among the seasons. In addition, the summer trends in SdAVG and AOT40 are with largest p-values (0.03 and 0.10) among the seasons. Our results are consistent with those in Xu et al. (2016), who analyzed the 1994–2013 ozone data from WLG and found that the trends in SAVG and SdAVG were largest and most significant in fall (0.32 ppb/yr, p = 0.00005 and 0.28 ppb/yr, p = 0.00006) and smallest and least significant in summer (0.13 ppb/yr, p = 0.20 and 0.04 ppb/yr, p = 0.59). Further study by Xu et al. (2018) indicates that the effect of pollution transport from South East Asia played a key role in the ozone trend in fall and stratosphere-to-troposphere transport contributed mostly to the ozone trend in spring.

As can be seen in Tables 5 and S3, although the SAVG values for SDZ increased in spring, summer, and fall, only is the trend in summer (1.04 ppb/yr or 2.03%/yr) with low p-value (0.07). SdAVG for this site increased clear in spring and fall, with a rate of 0.65 ppb/yr (1.23%/yr, p = 0.09) and 0.36 ppb/yr (0.93%/yr, p = 0.09), respectively. Large and small increases in SSOMO35 occurred in summer and fall, with a rate of 120.2 ppb-day/yr (3.07%/yr, p = 0.00) and 36.7 ppb-day/yr (2.78%/yr, p = 0.09), respectively. Although SAVG in spring and fall increased as well, the rates of increase (0.32 ppb/yr and 0.11 ppb/yr) are much smaller than that in summer, and the trends are with larger p-values. As shown in Figure 9 SdAVG for SDZ increased mainly in summer and spring. The summer trends of SdAVG reached 1.32 ppb/yr (2.09%/yr, p = 0.11). This trend is consistent with the 1–2 ppb/yr summertime trend estimated for the period 2003–2015 by Sun et al. (2016) based on ozone measurements from Mt. Tai (1465m asl), a mountain-top site in the NCP. SAVG, SdAVG, and SSOMO35 in winter showed negligible changes. These results suggest that the trends in AAVG, AdAVG, and SOMO35 for SDZ (Tables 4 and S2) were mainly from the increase of ozone in summer, with smaller contributions from ozone in spring and fall. As shown in Table 5, SAVG and SSOMO35 also had largest increases in summer and spring at the rural site (GCH), and in fall and winter at the urban site (CMA).

The vegetation metric AOT40 at SDZ increased at a large rate (858 ppb-h/yr, 3.92%/yr, p = 0.04) in spring (Tables 5 and S3). The increasing trend in summer (692 ppb-h/yr, 2.01%/yr) is also relatively large though with slightly larger p-value (0.11, Figure 9b). Another vegetation metric W126, which gives higher ozone values larger weight (Lefohn et al., 2018), showed clear trends in spring, summer, and fall, with rates of 1086 ppb-h/yr (5.08%/yr, p = 0.05), 1786 ppb-h/yr (4.53%/yr, p = 0.09), and 426 ppb-h/yr (3.85%/yr, p = 0.04), respectively (Tables 5 and S3). It is noteworthy that the summer trend in W126 at SDZ is the largest among all W126 trends in Table 5, and much larger than the reported global and regional trends in 3-month W126 for 1995–2014 (Mills et al., 2018). At GCH, the spring and summer trends in AOT40 were also very large, with rates of 757 ppb-h/yr (3.57%/yr, p = 0.14) and 502 ppb-h/yr (1.34%/yr, p = 0.34), and the summer trend in W126 (536 ppb-h/yr, 1.35%/yr, p = 0.58) is the fifth largest among all W126 trends in Table 5. The rapid increases in AOT40 and W126 at the rural and regional background sites in the NCP during growing seasons endanger crops and other vegetations in the region.

Nearly all seasonal ozone metrics at AKDL showed decreasing trends, with those in spring and summer being at large declining rates (from –4.59%/yr to –34.3%/yr) and with p-values less than 0.10 (see Tables 5 and S3, Figure 10). The decreases in the metrics in spring and summer are prominent and should be responsible for the observed annual trends at the site (see Table 4). As already mentioned, the trends in ozone metrics at AKDL should be viewed with caution because of the shorter observation period. As shown in Figure 10, the numbers of data points used for calculating the seasonal trends are either 5 or 6, which may not allow for obtaining robust trends. On the other hand, the contrast in the trends between spring/summer and fall/winter implies that the negative trends in spring and summer seemed not to be due to any random variations or technique problems. At present, it is not known why the spring and summer ozone metrics at AKDL declined at such large rates.

There were no significant trends in annual ozone metrics at LFS and LA during 2006–2016 (see Table 4). However, some of the seasonal ozone metrics at these sites did change at rates with higher certainties. At LFS, the spring values of SAVG, SdAVG, SSOMO35, and W126, and the winter values of SdAVG, SSOMO35, AOT40, and W126 all showed clear decreasing trends, while the fall values of SSOMO35, AOT40, and W126 showed clear and relatively large increasing trends. At LA, the spring values of SAVG, SdAVG, SSOMO35, and W126 decreased at considerable rates. No seasonal ozone metric at XGLL had significant trends.

Overall, seasonal ozone metrics increased moderately in all seasons at WLG and largely in summer at the regional and rural sites in the NCP (SDZ and GCH). Clear decreases in ozone metrics were observed in spring and summer at AKDL, and in colder seasons at LFS and LA. These increases in summer ozone metrics at the NCP sites are considerably rapid. For example, the summer trend of SdAVG at SDZ was 1.32 ppb/yr (2.09%/yr, p = 0.11) during 2004–2016, which is much larger than the summertime trends of daytime ozone at global non-urban sites (Gaudel et al., 2018). Although Figure 13 in Gaudel et al. (2018) shows a trend for increased daytime ozone at some sites in Korea and Japan during 2000–2014, none of these rates of increase are over 1 ppb/yr. Estimates of Chang et al. (2017) show that regional average trends of daytime ozone in summer in East Asia were 0.23 ppb/yr (rural) and 0.51 ppb/yr (urban) for the 2000–2014 period, all much smaller than the SDZ trend. These results suggest that the summertime increase trend of daytime ozone at SDZ belongs probably to one of largest trends found for non-urban and low-altitude sites during the same period (if it is not the largest one). The rapid increases in ozone concentrations have caused large increases in vegetation and human health relevant metrics. The summer and spring increasing trends of AOT40 and W126 for SDZ (Table 5) are greater than the 3-months trends of AOT40 and W126 in Mills et al. (2018). The trends of A4MDA8, SOMO35, and NDGT70 for SDZ (Table4) are all larger than those reported by Fleming et al. (2018) for the period 2000–2014. All these indicate a very serious deterioration of ozone air quality in the NCP region, in which SDZ is located.

The long-term changes in ozone metrics at our sites cannot always be directly explained by changes of ozone precursors (e.g., NO₂). While the increases in ozone at WLG are consistent with the positive trend of tropospheric NO₂, those at AKDL are inconsistent with the small increase trend of tropospheric NO₂ (Table S1). Increases and decreases of ozone were found at LA respectively in cold (winter) and warm (spring-fall) periods there. This may probably be related to the clear reduction of tropospheric NO₂ (–0.079 DU/yr, p = 0.04, Figure 1 and Table S1). Ozone values for LFS increased during warmer period (summer and fall) and decreased in colder period (winter and spring). It is unclear whether or not this was a result of the minor increase of tropospheric NO₂ (0.009 DU/10yr). Tropospheric NO₂ has decreased over the sites in the NCP, with a clear declining trend over CMA (–0.142 DU/yr, p = 0.00) (see Figure1 and Table S1). Since the NCP has been one of most polluted regions of the world and with high levels of tropospheric NO₂ (Table S1), even the regional background site (SDZ) could be within the area where ozone production is more sensitive to VOCs (Ge et al., 2012; Ma et al., 2016). If this is true, the increases of ozone at SDZ in warm seasons does not contradict the small declining rate of tropospheric NO₂ (–0.013 DU/yr, Table S1). The rapid decrease in NO₂ over the urban site CMA (0.142 DU/yr, p = 0.00) (see Table S1 and Figure 1) should be consistent with less NO titration with surface ozone hence the clear increases of the ozone levels in winter (Table 5).

In this study, average values for ozone metrics during 2012–2016 are considered as present-day levels and compared to see the regional differences in the metrics. Note that the time frame 2012–2016 for present-day in this study has a two-year shift from that defined in TOAR (2010–2014) (Schultz et al., 2017). The main purpose of this shift is to select a period as present-day that is not too far from now. In addition, some of the 2010–2014 ozone metrics for background sites in mainland China and other sites over the world are presented in TOAR-Health (Fleming et al., 2018),TOAR-Vegetation (Mills et al., 2018), and TOAR-Climate (Gaudel et al., 2018). To keep this paper concise, we do not compare present-day values of all ozone metrics given in Table 3. Instead, we choose three annual and three seasonal ozone metrics for comparison.

Figure 11 shows the present-day values of AdAVG, SOMO35, and AmaxMDA8, which are influenced by the whole range of daytime ozone values, the medium and high range, and the high end of ozone values, respectively. WLG has the largest present-day AdAVG (52.6 ppb), followed by SDZ (48.1 ppb), LA (41.1 ppb), CMA (40.9 ppb), and XGLL (40.5 ppb). AKDL, LFS, and GCH have similar present-day AdAVG values (33.5–35.9 ppb). Although GCH and CMA are sites located in polluted areas in the NCP, the present-day AdAVG values at these sites are about in the middle. As mentioned in section 3.2, the annual ozone concentrations for WLG are mainly determined by background ozone, which is subject to free tropospheric and stratospheric influences (Ma et al., 2002; Xu et al., 2018). XGLL is also significantly influenced by stratosphere-to-troposphere transport (Ma et al., 2014). Therefore, the high present-day AdAVG values at WLG and XGLL are mainly due to higher levels of ozone background. Regional pollution should have played an important role in sustaining relatively high AdAVG at SDZ and LA. The relatively lower present-day AdAVG values at GCH and CMA should be mainly caused by suppressed ozone concentrations in cold seasons (see Figure 3) under high NO conditions (Lin et al., 2009; Lin et al., 2011). The high NO_x levels over the NCP sites (see Table S1) caused not only lower AdAVG at GCH and CMA but also high AmaxMDA8 at SDZ, GCH, and CMA. As shown in Figure 11a, the present-day AmaxMDA8 is largest at CMA (158.9 ppb), second largest at SDZ (149.7 ppb), and third largest at GCH (137.8 ppb). These large AmaxMDA8 values indicate that the NCP has been the region with severest photochemical pollution in China. This is consistent with the results based mainly on ozone measurements from urban sites (Lu et al., 2018). The present-day AmaxMDA8 values at LA and LFS are 125.9 ppb and 107.1 ppb, respectively. They are much lower than those in the NCP but still significantly higher than current national standard for MDA8 ozone (about 75 ppb).

The present-day SOMO35 is the highest at CMA (9054 ppb-day), followed by SDZ (8905 ppb-day) and WLG (8051 ppb-day). CMA is an urban site in megacity Beijing and SDZ is a regional background site downwind of urban Beijing. Both sites are located in the NCP, one of the most polluted regions of world. Therefore, it is not surprising to find the highest and second highest present-day SOMO35 at these sites. Due to the high background, the majority of ozone values from WLG are higher than 35 ppb (see Figure 2a), leading to high SOMO35 at this global baseline station. Although GCH is also located in the NCP, its present-day SOMO35 value (6379 ppb-day) is significantly smaller than those at CMA and SDZ. The main reason for this is the extremely skewed distribution of ozone values. As can been seen in Figure 2g, although the high end ozone values (e.g., the 95-percentiles) at GCH can be quite high, most of the ozone values are below 35 ppb. Therefore, fewer numbers of ozone values from GCH can be included in the calculation of SOMO35 compared with those from CMA and SDZ. LA is another site showing large present-day SOMO35 (7015 ppb-day). Being a regional background site in the YRD, the second largest polluted region in China, LA experienced an elevation of high end ozone values during 1991–2006 (Xu et al., 2008). Photochemical pollution over LA was not worsened during 2006–2016, nor was it significantly mitigated, as suggested by the data in Table 4. Therefore, both AmaxMDA8 and SOMO35 show relatively high present-day values (Figure 11). The present-day SOMO35 values at LFS, XGLL, and AKDL are all small and close to those found at most sites in North America and Europe (Fleming et al., 2018).

To see the regional differences in present-day ozone metrics and their seasonal dependence, we calculated for every site the averages of SdAVG, AOT40, and W126 during 2012–2016. The results are presented in Figure12. The present-day SdAVG at WLG shows only small seasonal variations and its summer value is close to those at the NCP sites (Figure 12a). The present-day SdAVG values at SDZ, CMA, and GCH rank first, second, and third in summer, respectively, but they are very low in colder seasons, particularly in winter. This suggests that surface ozone levels at these polluted sites can be readily enhanced under favourable photochemical conditions though they are suppressed in colder seasons presumably by high NO_x levels (Simon et al., 2015).

The present-day values of AOT40 and W126 show very large differences among sites and seasons (Figure 12b and12c). The present-day AOT40 values at the NCP sites (i.e., CMA, GCH, and SDZ) are prominently high in warmer seasons, particularly in summer, while they are very low in winter. The present-day AOT40 at LA is close to those at the NCP sites in spring and fall, while it is significantly lower than those at the NCP sites in summer. The present-day AOT40 values at CMA, GCH, SDZ, and LA are all over 20000 ppb-h in the warmer seasons, much higher than the threshold (15000 ppb-h) set for the group of most polluted sites (Mills et al., 2018). And the summer AOT40 at CMA (46079 ppb-h) is even more than three times as large as the threshold. As shown in Figure12b, higher present-day AOT40 values occur also at WLG in summer and spring, at LFS in summer, and at XGLL in spring though they are much smaller than those at the NCP sites.

The weighting function for calculating W126 can influence the results from sites and seasons associated with higher and lower ozone levels. The weighting effect can be clearly seen by comparing the W126 data in Figure12c with the AOT40 data in Figure12b. The present-day AOT40 values at CMA, GCH, and SDZ in summer are protruding, and the corresponding present-day W126 values are even more outstanding. The summer W126 values at these sites are all significantly higher than their AOT40 counterparts, and the differences among them are smaller than those among their AOT40 counterparts. In addition, the summer W126 at WLG is also enlarged compared with the summer AOT40 at the site. All other present-day W126 values become smaller than their AOT40 counterparts. Compared with the W126 threshold (25000 ppb-h) set for the group of most polluted sites (Mills et al., 2018), the summer values of present-day W126 at CMA, GCH, and SDZ are 73%–94% higher. The W126 values at WLG in summer and at SDZ and CMA in spring are slightly higher than the threshold. The high warm season values of present-day AOT40 and W126 at some sites, particularly those in the NCP and YRD suggest that the ozone levels at these sites are very detrimental to local vegetation. Potential impacts on crops and forests need to be assessed. Recently, Feng et al. (2019) published an assessment of economic losses caused by ozone impacts on human health, forest productivity, and crop yield in China. They used one year (2015) ozone measurements from more than 1497 sites in mainland China and found significant ozone-related losses. Including more data, particularly those from long-term observation at background sites may gain more systematic results.

In this work, we present for the first time an integrated analysis of long-term measurements of surface ozone at six background sites (WLG, SDZ, LA, LFS, XGLL, and AKDL), a rural site (GCH), and an urban site (CMA) in mainland China. We characterize the long-term variations in annual statistics and the seasonal and diurnal cycles, discuss the long-term trends of some annual and seasonal ozone metrics relevant to human health, vegetation, and climate change, and compare present-day levels of some ozone metrics among the sites. The ozone time-series analyzed range from 7 to 22 years, with most of them being over 11 years.

High annual average values of entire-day ozone are observed at WLG (47.9–54.1 ppb) and XGLL (31.8–42.9 ppb). Both sites are far from anthropogenic emissions but significantly impacted by ozone-rich air from the free troposphere and the lower stratosphere due to high altitudes (>3.5 km). Among the sites in eastern China, the highest annual average ozone (40.4 ppb) occurs at SDZ, a regional background site in the NCP. At the rural site GCH and urban site CMA in the NCP, however, the lowest (18.8–31.3 ppb) and second lowest (23.0–33.8 ppb) ranges of annual average ozone are observed due to very frequent low ozone values presumably caused by strong NO titration. After the titration effect is minimized by taking only daytime ozone into account, the average levels of ozone at GCH and CMA become respectively 45.5% and 32.9% higher than the averages of entire-day ozone.

We find clear seasonal cycles at all sites studied, with average ozone peaking around July (at WLG, LFS, SDZ, GCH, and CMA), June (at LA), May (at XGLL), and March (at AKDL). Located in the southern part of China, LA and XGLL are strongly influenced by Asian summer monsoon hence show an ozone decrease in summer. Ozone at the sites in the polluted regions (NCP and YRD) show very large diurnal variations, with peaks around 15:00 or 16:00, indicating strong photochemical ozone production under the daylight condition and large reduction at night due to NO titration and dry deposition. The ozone levels at XGLL, ADKL, and LFS have minor or small rises during daytime, while there is a small valley in ozone level around noon at WLG, which is attributable to ozone-poor air brought to the site by valley breeze. The impact of local photochemistry on ozone at these sites is small as they are remote or only slightly polluted sites.

Following the definitions and methodologies of TOAR, we have calculated the long-term trends and present-day values of some annual and seasonal ozone metrics relevant to human health, vegetation, and climate change. Values of seven annual ozone metrics (AAVG, AdAVG, SOMO35, AmaxMDA8, A4MDA8, NDGT90, and NDGT70) and five seasonal metrics (SAVG, SdAVG, SSOMO35, AOT40, and W126) are obtained for each site based on the ozone measurements fulfilling the data capture criterion for TOAR. We find the trends in the ozone metrics using the M-K test and estimate the trend values using the Theil-Sen estimator. This allows for comparability of our results with those reported in TOAR.

We find clear and very large increase trends in the annual ozone metrics (AAVG, AdAVG, SOMO35, AmaxMDA8, A4MDA8, NDGT90, and NDGT70) and in the seasonal metrics (SAVG, SdAVG, SSOMO35, AOT40, and W126) in summer at the background site SDZ in the NCP. In addition, some spring and fall values of the seasonal ozone metrics also show large increases. These indicate that photochemical pollution has intensified in warm seasons and caused large positive trends in many annual and seasonal ozone metrics at this site. The circumstance looks different at the rural and urban sites in the NCP. We only find clear increases in fall and winter SAVG as well as winter SdAVG at CMA. The summer and spring AOT40 and W126 as well as summer SAVG at GCH show relatively large increases but with lower certainty. For the sites used in our analyses, this implies that trends in annual ozone metrics can be masked by stronger disturbances at the sites close to anthropogenic impacts and better resolved using observations at regional background sites. The global baseline site WLG shows clear moderate increasing trends in AAVG, AdAVG, SOMO35, and NDGT70, and in all seasonal metrics values calculated. Changes in emissions of ozone precursors in South East Asia and climate variability are responsible for the increase trends in ozone at WLG (Xu et al., 2018). Annual ozone metrics for AKDL show large negative trends because of the large downward trends of ozone in spring and summer. Although robust ozone trends from AKDL cannot be obtained at present from the shorter ozone time-series, the clear decreases of many ozone metrics values for this site suggest that something interesting might have occurred in the local/regional environment, in particular in spring and summer months. We find no clear trends in the annual ozone metrics at other sites. However, some values of seasonal ozone metrics at LFS and LA show large positive or negative trends.

For obtaining representative present-day values of ozone, we consider the period 2012–2016 as present-day and have made averages of three annual metrics (AdAVG, AmaxMDA8 and SOMO35) and three seasonal metrics (SdAVG, AOT40 and W126). The sites located in the polluted NCP and YRD (e.g., SDZ, GCH, CMA, and LA) have quite high present-day values of annual AmaxMDA8 and SOMO35, while the global baseline site WLG has the highest present-day values of AdAVG and high SOMO35. The phenomenon can be explained by the large differences in ozone concentration distributions of these sites. The ozone concentrations at WLG are distributed mainly around medium and higher levels, while those at the polluted sites in the range from very low to very high due to stronger photochemical pollution in warm seasons and more NO titration in the dark and cold periods. The same reasons can explain the large seasonal differences in SdAVG, AOT40, and W126 at the polluted sites. As indicated by the tropospheric NO₂ values in Table S1, the sites in the NCP (i.e., SDZ, GCH, and CMA) are most polluted. The present-day for ozone metrics at these sites, particularly those in summer, are much larger than the highest thresholds set in TOAR for metrics relevant to climate change, human health, and vegetation (Gaudel et al, 2018; Fleming et al., 2018; Mills et al., 2018). The summer AOT40 for CMA is even larger than triple of the threshold. The AOT40 values for LA in spring-fall are also much larger than the threshold. Consequently, ambient air over the NCP and YRD is very unhealthy for human and toxic for vegetation in warm seasons, particularly in summer.

Results from this work, together with others, extend our knowledge about the spatiotemporal variations of surface ozone in China. Background ozone shows a very rapid long-term increase in the NCP but no significant trends in the YRD, Northeast China, and Southwest China. Ozone has increased strongly at SDZ, especially with respect to ozone metrics relevant to human health. Comparisons with the results from TOAR (Gaudel et al, 2018; Fleming et al., 2018; Mills et al., 2018) indicate that the summer concentrations of surface ozone in the NCP are some of highest seen world-wide, and the summer values of several metrics for SDZ have been increasing at rates rarely observed at other background sites over similar periods. Reducing the observed levels of these metrics would require reductions of regionally emitted ozone precursor gases. Although ozone metrics at WLG show significant increases in all seasons, those at AKDL, a site in the Gobi area of Xinjiang (about 1700 km northwest of WLG), show large decreases, particularly in spring and summer. The reasons for the different trends of ozone at these two sites are currently unknown.

We also find large negative trends in some seasonal ozone metrics but no clear trends in the annual metrics at LFS and LA. It is worthy to mention that the high-end values of ozone were found to increase during 1991–2006 at LA, the background site in East China (Xu et al., 2008). For the similar period (1994–2007), Wang et al. (2009) found an 0.58 ppb/yr trend of ozone at Hok Tsui (a coastal site in Hong Kong), which is largely impacted by air masses from East China. Our results in this paper indicate that the increase of ozone at LA leveled off during 2006–2016. Recently, Wang et al. (2019) revealed that the increase of ozone at Hok Tsui occurred mainly in 1994–2006 and leveled off in 2007–2018. They found a large ozone increase in the maritime air from Southeast Asia but stabilization in the impact from East Asia in the recent decade. The consistency between our results and those from Wang et al. (2019) suggests that the reductions of ozone precursors in East China (see Figure 1) have stopped increases of ozone at some background or rural sites in this region and in the air masses transported out of this region. This can also be confirmed by the study of Gao et al. (2017), who analyzed long-term ozone measurements from Shanghai, East China. For the 2006–2015 period, they found no clear trend of ozone at the background site (Dongtan) but a large increasing trend of ozone at the urban site (Xujiahui) in the city. The increase of ozone at the urban site was mainly caused by the reduction of NO_x emission (Gao et al., 2017). Although tropospheric NO₂ has been declining over the major air pollution regions (Figure 1), strong increases of ozone at many of urban sites in China were observed during 2013–2017 (Lu et al., 2018). These imply that China’s policy of emissions reduction in past years was unfavorable for controlling ozone levels at urban sites though it has stopped the increases of ozone at some background sites and effectively lowered PM_2.5 concentrations. As the relations of ozone to NO_x, VOCs, PM_2.5, etc., are very complex, more in-depth studies are needed to better understand the different long-term changes of ozone in different regions and assess the impacts of changing ozone levels on human, ecosystem, and climate.

The ozone data analyzed in this work are available from the TOAR database (http://doi.org/10.34730/19cdb0ca3ad046bd88449f38b426542b).

We thank many individuals that are not listed as coauthors for their contributions to the ozone observations at all stations. These contributions include but not limited to routine observations and station management. We also thank Drs. Christoph Zellweger, Peter Hofer, and Jörg Klausen from Swiss/WCC-Empa for their calibrations of our instruments. We acknowledge Martin Schultz and Sabine Schröder from Research Center Jülich, Germany for their support in metadata preparation and DOIs creation for our datasets. Finally, we thank David Tarasick, Allen S. Lefohn, and two anonymous reviewers for their constructive comments.

The observation systems and the operation costs at the GAW stations (WLG, SDZ, LA, LFS, XGLL, and AKDL) are provided by CMA. WLG received support from WMO/GEF in the instrumentations. The long-term observations at the GCH and CMA sites are supported by the Natural Science Foundation of China (No. 41330422, 40775074, and 41505107), China Special Fund for Meteorological Research in the Public Interest (No. GYHY201206015), and the Basic Research Fund of CAMS (No. 2013Z005, 2016Z001, and 2017Z011).

The authors have no competing interests to declare.

• Contributed to conception and design: XX, WL

• Contributed to acquisition of data: XX, WL, JJ, XZ, YW, WX, GZ, ZM, QM, DY, YD, ZL, DW, HZ

• Contributed to analysis and interpretation of data: XX, WL, WX

• Drafted and/or revised the article: XX, WL, WX, JJ, XZ, YW, GZ, ZM, QM, DY, YD, ZL, DW, HZ

• Approved the submitted version for publication: XX, WL, JJ, XZ, YW, WX, GZ, ZM, QM, DY, YD, ZL, DW, HZ