Long-term trends of surface ozone and its influencing factors at the Mt Waliguan GAW station , China – Part 1 : Overall trends and characteristics

Tropospheric ozone is an important atmospheric oxidant, greenhouse gas and atmospheric pollutant at the same time. The oxidation capacity of the atmosphere, climate, human and vegetation health can be impacted by the increase of the ozone level. Therefore, long-term determination of trends of baseline ozone is highly needed information for environmental and climate change assessment. So far, studies on the long-term trends of ozone at representative sites are mainly available for European and North American sites. Similar studies are lacking for China and many other developing countries. Measurements of surface ozone were carried out at a baseline Global Atmospheric Watch (GAW) station in the north-eastern Tibetan Plateau region (Mt Waliguan, 3617 N, 10054 E, 3816 m a.s.l.) for the period of 1994 to 2013. To uncover the variation characteristics, long-term trends and influencing factors of surface ozone at this remote site in western China, a two-part study has been carried out, with this part focusing on the overall characteristics of diurnal, seasonal and long-term variations and the trends of surface ozone. To obtain reliable ozone trends, we performed the Mann–Kendall trend test and the Hilbert–Huang transform (HHT) analysis on the ozone data. Our results confirm that the mountain-valley breeze plays an important role in the diurnal cycle of surface ozone at Waliguan, resulting in higher ozone values during the night and lower ones during the day, as was previously reported. Systematic diurnal and seasonal variations were found in mountain-valley breezes at the site, which were used in defining season-dependent daytime and nighttime periods for trend calculations. Significant positive trends in surface ozone were detected for both daytime (0.24± 0.16 ppbv year) and nighttime (0.28± 0.17 ppbv year). The largest nighttime increasing rate occurred in autumn (0.29± 0.11 ppbv year), followed by spring (0.24± 0.12 ppbv year), summer (0.22± 0.20 ppbv year) and winter (0.13± 0.10 ppbv year), respectively. The HHT spectral analysis identified four different stages with different positive trends, with the largest increase occurring around May 2000 and October 2010. The HHT results suggest that there were 2–4a, 7a and 11a periodicities in the time series of surface ozone at Waliguan. The results of this study can be used for assessments of climate and environment change and in the validation of chemistry–climate models.


Introduction
Ozone (O 3 ) is one of the key atmospheric species and is closely related to climate change and environmental issues (IPCC, 2013).The stratospheric ozone layer protects living organisms at the Earth's surface against the harmful solar UV radiation, while tropospheric ozone is an important greenhouse gas and governs oxidation processes in the Earth's atmosphere through the formation of OH radical (Monks et al., 2015;The Royal Society, 2008;Young et al., 2013;Lelieveld and Dentener, 2000;Levy, 1971).In the surface layer, ozone W. Xu et al.: Long-term trends of surface ozone and its influencing factors -Part 1 is also one of the toxic gases for human beings and vegetation.
Data on the spatiotemporal variations of ozone are greatly needed for assessing the impacts of ozone on human health, ecosystems, and climate.Since ozone is a secondary pollutant, with a lifetime of 22 days (Young et al., 2013), its mixing ratios are influenced both by local photochemistry and by transport of ozone and its precursors (Wang et al., 2006a;Lal et al., 2014).Deep convection and stratosphereto-troposphere exchange (STE) events can also bring down ozone-rich air from above and influence surface ozone mixing ratios at high-elevation sites (Stohl et al., 2003;Lefohn et al., 2012;Jia et al., 2015;Ma et al., 2014;Langford et al., 2015;Lin et al., 2015a).In the troposphere, particularly in the surface layer, ozone is highly variable in space and time due to the large variabilities of its dominant sources and sinks, which are impacted by anthropogenic activities and meteorological conditions.So far, there has been no better way than networked monitoring to obtain the spatial distribution and long-term temporal variation of surface ozone.
The Global Atmosphere Watch (GAW) program of the World Meteorological Organization (WMO) has been one of the key international initiatives in long-term monitoring of the chemical and physical properties of the atmosphere.Many GAW stations have been set up to monitor atmospheric compositions including surface ozone due to its importance and due to the urgent need to evaluate the trends of baseline ozone.Based on data from some GAW sites and other sources, past trends in surface baseline ozone have been reported for Europe and North America (Cui et al., 2011;Gilge et al., 2010;Oltmans et al., 2013;Vingarzan, 2004;Parrish et al., 2012;Logan et al., 2012), which mostly reveal strong increases in ozone before 2000 and slow or even no growth afterwards.Data from some important regions, e.g.East Asia and South America, are very scarce, which makes them even more valuable.China, as one of the rapidly developing countries, is contributing increasing ozone precursor emissions to the atmosphere and transport from the South and East Asian sector was thought to be most responsible for the increase in ozone in the western United States (Cooper et al., 2010), though other studies suggest that STE events had an equivalent important role in causing high-ozone events at western US alpine sites during spring (e.g.Langford et al., 2009;Ambrose et al., 2011;Lin et al., 2012Lin et al., , 2015a)).A recent study by Lin et al. (2015b) found that although rising Asian emissions contributed to increasing springtime baseline ozone over the western US from the 1980s to the 2000s, the observed western US ozone trend over the short period of 1995-2008 previously reported by Cooper et al. (2010) was strongly biased by meteorological variability and measurement sampling artefacts, resulting in an overestimate of the trend.Extending the analysis to 1995-2014, a weaker ozone trend of 0.31 ± 0.21 ppbv year −1 was found from observations and a similar one from model results.Nevertheless, the impact of Asian pollution outflow events on western US surface ozone is evident (Lin et al., 2012(Lin et al., , 2015a)).
Besides the impact of Asian pollution outflow on the surface ozone in other regions, it is at least equally important to know how the level of surface ozone in Asia, particular in China, has been changing.Long-term changes in ozone in China, however, have only been reported in a few publications.Ding et al. (2008) studied the tropospheric ozone climatology over Beijing based on data from the MOZAIC (Measurement of Ozone and Water Vapor by Airbus In-Service Aircraft) program and found a 2 % year −1 increase of boundary layer ozone from the period of 1995-1999 to 2000-2005 over Beijing in the North China Plain (NCP) region and a weaker increasing trend of free-tropospheric ozone.Wang et al. (2012) reported a similar increasing trend of lower tropospheric ozone and larger ozone increases in the middle and upper troposphere for the period of 2002-2010 based on ozonesonde measurements over Beijing.Ma et al. (2016) found an increase of 1.1 ppbv year −1 for the period 2003-2015 in the maximum daily average 8 h concentration of ozone at Shangdianzi, a background site in the NCP.Xu et al. (2008) reported positive trends of extreme values and increased variability in six periods of ozone measurements from 1991 to 2006 at Lin'an, a background site in the Yangtze River Delta (YRD) region.Wang et al. (2009) found a significant increasing trend of 0.58 ppbv year −1 during 1994-2007 at a coastal site of Hong Kong in the Pearl River Delta (PRD) region, which was caused by rapid increases in ozone precursor emissions in the upwind source regions.The above studies all focus on the most polluted regions in the eastern part of China, i.e. the NCP, YRD and PRD, aiming to study the impact of growing precursor emissions on ozone trends.The trend of ozone over remote regions in China still remains to be studied based on long-term observations.
Continuous long-term observations of surface ozone have been made only at a few representative sites in China, among which is the Mt Waliguan (WLG) GAW station, which has the longest ozone measurement record in China.It is a pristine high elevation site located downwind of the European, Central Asian and Indian outflow, representative of the baseline of Eurasia.A few studies have already been performed on short-term measurements of ozone at WLG. Surface ozone at the site has been proved to be highly representative of free-tropospheric ozone (Ma et al., 2002b) and hence is subjected to the influences of STE events (Ding and Wang, 2006;Zhu et al., 2004).Air masses from the west are dominant at WLG and were found to be associated with the highest ozone mixing ratios (Wang et al., 2006b).Only in summer a substantial portion of the airflows comes from the eastern sector and exposes the surface ozone mixing ratio to some regional anthropogenic influences (Wang et al., 2006b;Xue et al., 2011).Previous studies of ozone at WLG, based on short-term measurements and modelling results, clarified the causes for certain episodes or for the diurnal and seasonal cycles of ozone (Ma et al., 2002a(Ma et al., , 2005;;Zhu et al., 2004).The overall variation characteristics and long-term trend of ozone at WLG have not yet been studied.Considering the geographical representativeness of the WLG site, results on the long-term variations of ozone at WLG may improve understanding of ozone changes in the northern mid-latitudes, particularly in the rural and remote regions of Eurasia.
The most common method used in the detection of ozone trends is the linear least squares method (Tarasova et al., 2009;Cui et al., 2011;Wang et al., 2009;Xu et al., 2008).Other studies directly compared mean ozone levels of different periods to detect possible trends (Ding et al., 2008;Lin et al., 2014).Oltmans et al. (2013) determined trends and their significance in the W126 and W_Low metrics using the Theil-Sen estimate together with the Mann-Kendall's tau test (Kendall, 1955;Sen, 1968).The non-linear variation of ozone mixing ratio with season and many other meteorological factors can introduce uncertainties into the linear trend analysis.Wang et al. (2012) de-seasonalised the monthly data by subtracting the average of all monthly data for a given month from the original data of the same month before performing a linear regression analysis.Oltmans et al. (2006Oltmans et al. ( , 2013) ) first performed an autoregressive model fitting incorporating explanatory variables (that are known sources of ozone variability) and a cubic polynomial fit to better represent the long-term variations of ozone, then used a bootstrap method to determine the trends of ozone.However, surface ozone typically is influenced by many factors, which makes it hard to determine which to incorporate.The seasonal Mann-Kendall test, which is a modified version of the non-parametric Mann-Kendall trend test, can account for the seasonal variation within the data (Hamed and Ramachandra Rao, 1998).It has been widely applied in hydrology and seldom in atmospheric chemistry.The Hilbert-Huang transform (HHT) analysis, which has been widely applied on the analysis of meteorological data sets and not yet on that of atmospheric composition data, is a precise and adaptive spectral analysis method that can divide the signal into various oscillation modes and study the anomaly and periodicity within the data (Rao and Hsu, 2008).Applications of HHT on temperature, wind, rainfall and solar radiation data have proved that the HHT method is capable of capturing synoptic and climatic features, revealing known diurnal, seasonal, annual and inter-annual cycles (Huang, 2014).
In this paper, we present the first part of an analysis on 20year surface ozone mixing ratio at WLG, focusing mainly on the overall, diurnal, seasonal and long-term variations characteristics and the trends of surface ozone.We apply a linear regression as well as a seasonal Mann-Kendall test together with the Theil-Sen estimate to calculate the overall trend of ozone.The HHT spectral analysis method is used for the first time to investigate the ozone trends during different periods and the underlying anomalies and periodicities within the ozone data.A detailed discussion on the influenc- ing factors contributing to the ozone variation at WLG will be presented in a companion paper.
2 Data and methodology

Site and measurements
The Mt Waliguan site (WLG, 36 WLG is situated at the northeast edge of the Qinghai-Tibetan plateau and surrounded by highland steppes, tundra, deserts, salt lakes, etc. (Fig. 1).With very low population (about 6 persons km −2 ) and hardly any industry within 30 km, the WLG site is far from major anthropogenic sources of ozone precursors.However, some impact of long-range transport of anthropogenic pollutants from the NE-SE sector cannot be excluded, particularly from the major cities Xining (about 90 km northeast of WLG, population ∼ 2.13 millions) and Lanzhou (about 260 km east of WLG, population ∼ 3.1 millions).Such impact, if any, may be significant only during the warmer period (May-September), as suggested by previous air-mass trajectory studies (Zhang et al., 2011).The vehicle emission on the highway running from the northwest to the northeast of WLG may also be another source for anthropogenic ozone precursors (Wang et al., 2015).The WLG baseline station was established in 1994.The long-term monitoring program for surface ozone began in August 1994.The mixing ratios of surface ozone have been measured using two ozone analysers (Model 49, Thermal Environmental Instruments; one of the analysers was replaced with a Model 49i ozone analyser in 2011) at a sampling height of 7 m.The analysers have been automatically zeroed alternatively every second day by introducing ozone-free air for 45 min.Seasonal multipoint calibrations (8 points) have been conducted using an ozone calibrator (Model 49PS, Thermal Environmental Instruments).The analysers have been checked weekly for changes in instrument parameters.The inlet filters have been replaced weekly.Maintenance on the observation system has been performed yearly and whenever it was necessary.The yearly maintenance includes cleaning of absorption tubes, pumps, inlet tubing and other connecting parts, and checking of the inlet loss.In the years 1994, 1995, 2000, 2004, and 2009, the ozone calibrator and analysers at WLG were compared with the transfer standard from the WMO World Calibration Centre for Surface Ozone and Carbon Monoxide, EMPA Dübendorf, Switzerland.Intercomparison results show excellent or good agreement between the WLG instruments and the transfer standard (Zellweger et al., 2000(Zellweger et al., , 2004(Zellweger et al., , 2009)).The two instruments performed parallel measurements, recording surface ozone data as 5 min averages, which were corrected annually based on the zero checks and multipoint calibrations.If the observed ozone values from the two analysers agreed within 3 ppb, average values were calculated and included in the final data set.Otherwise, causes for the differences were searched by the principal investigator and only data from the well-performing analyser were included in the data set.95 % of the data pairs show discrepancies within ±1.0 ppb and the difference between two instruments shows a nearly random distribution around zero. 5 min averaged ozone mixing ratios from August 1994 to December 2013 were then averaged into hourly data and used in this study.In the trend analysis, monthly average ozone mixing ratios were acquired by first calculating the daily average ozone values and then performing a monthly averaging.A data completeness of 75 % was required for each averaging step.
Meteorological observations have been made at the site using automatic weather stations (AWS) installed on the ground level and on an 80 m tower at 2, 10, 20, 40 and 80 m height.These observations provide meteorological parameters such as temperature, pressure, precipitation, and wind speed/direction in 5 min resolution.Additionally, the vertical velocity is measured at the 80 m platform.The 10 m horizontal wind and 80 m vertical wind data from August 1994 to December 2013 are used in this study and have been accordingly averaged into hourly data, which meet a data completeness requirement of 75 %.

Determination of daytime and nighttime
Past research has already revealed that the surface ozone at WLG is governed by different air masses during daytime and nighttime (Ma et al., 2002b).The WLG station experiences upslope winds during the day and is influenced by boundary layer (BL) air, while during the night, winds are downslope and the site is controlled by free tropospheric (FT) air.The boundary layer air is largely influenced by local photochem-istry and may contain pollutants transported from nearby areas, while the free-tropospheric air may sometimes contain signals of long-range transport or STE events.Hence, it is of necessity to differentiate between daytime and nighttime ozone mixing ratios in order to study the trends associated with different air masses.
In the previous study (Xu et al., 2011), daytime and nighttime were defined as fixed time ranges (e.g.11:00-16:00 LT for daytime and 23:00-04:00 LT for nighttime).However, the actual well-developed daytime and nighttime ranges vary with season, and so does the local wind.Figure 2a-c respectively show the season-diurnal variation characteristics of 10 m zonal (u) and meridional (v) wind velocity and the 80 m vertical (ω) wind velocity.Due to the local topography, the WLG station is under the influence of mountainvalley breezes and all three wind vectors exhibit distinct diurnal variation characteristics.The height difference to the west of Mt WLG is much larger than that to the east; hence valley breezes during daytime come from the west accompanied by upward drafts, resulting in a diurnal maximum u and w vector between noontime and middle afternoon depending on the season.The v vector changes from southern to northern winds around noontime.Mountain breezes during the night come from the south-east sector accompanied by subsiding air flows, resulting in low u and w and high v during the night.The dominant air flow at WLG is westerly during the cold seasons, which enhances the westerly valley breeze during the day and cancels out the easterly mountain breeze during the night.During the warm seasons, easterly winds gain in frequency, which sometimes cancels out the daytime valley breeze and enhances the nighttime mountain breeze.The distinct diurnal variation of the wind can be used to define a daytime and nighttime range that varies with season.The white dots in Fig. 2 represent the monthly average occurrence hour of the diurnal maximum u.In this study, a 6 h time range that is centred around the white dots is used as the daytime range (white dashed lines in Fig. 2).The nighttime window also covers 6 h and is considered to be offset by 12 h to the daytime window.

Trend analysis
The trend analysis is performed using both the Spearman's linear trend analysis and the modified Mann-Kendall trend test.The Mann-Kendall test is performed using a Fortran program developed by Helsel et al. (2006).Here, a brief description on the modified Mann-Kendall test is given.The Mann-Kendall test is a non-parametric test commonly used to detect trends.Hamed and Ramachandra Rao (1998) modified the test, so that it can be used on data with seasonality.
For two sets of observations X = x 1 , x 2 , . . ., x n and Y = y 1 , y 2 , . . ., y n , the rank correlation test as proposed by Kendall (1955) is performed as the following: where If Y is replaced with the time order T = 1, 2, . . ., n, the test becomes a trend test and S = i<j a ij .The significance of the trend is tested by comparing the standardised test statistic Z = S/ √ var(S) to the standard normal variate at a given significance level (Z α ).Here, a modified var(S) is given by the following: where n n * S represents a correction for the autocorrelation that exists in the data and can be obtained by an approximation to the theoretical values. (4) Here ρ s (i) is the autocorrelation function of the ranks of the observations.If |Z| > Z 1−α/2 , then the data are non-stationary, a positive Z indicates a positive trend and a negative Z suggests a declining trend.If |Z| ≤ Z 1−α/2 , then the data are stationary.Here we use α = 0.05, hence the corresponding critical Z 1−α/2 = 1.96.A non-parametric method (Theil-Sen estimate) is then used to estimate the magnitude of the trend; details can be found in Sen (1968).

The Hilbert-Huang transform analysis
The Hilbert-Huang transform (HHT) analysis is a combination of the empirical mode decomposition (EMD) and the Hilbert spectral analysis proposed by Huang et al. (1998).It is often used to analyse the time-frequency variation of non-linear and non-stationary processes.The EMD acts as a time-frequency filter; it decomposes the data into several oscillation modes with different characteristic timescales.The HHT method has been proved to be an efficient and precise method in investigating the periodicity, long-term oscillations and trends that are embedded within the data (Huang and Wu, 2008).So far, it has been widely applied in the studies of meteorology and climate, including wind field, temperature, radiation and rainfall analysis (Rao and Hsu, 2008;Lundquist, 2003;El-Askary et al., 2004), but it has not been used on atmospheric composition data.Here we give a brief description of the HHT method.
First, the EMD is performed on the data, to decompose the data into n intrinsic mode functions (IMFs), c 1 , c 2 , . . ., c n , and one residual r n , which are ordered from the smallest to the largest variation timescale (Huang et al., 2003).
Then the Hilbert transform is applied to each IMF using Eq. ( 6), where P is the Cauchy principal value.An analytical signal is then obtained with Eq. ( 7), where The instantaneous frequency ω can be calculated as the following: Thus, Eq. ( 5) can be transformed into the following expression: where R is the real part of the complex number.
To obtain the Hilbert amplitude spectrum H (ω, t), we assign for each time t, the calculated amplitude a j (t) to the associated ω j (t).An integration of H (ω, t) over the frequency span yields the instantaneous energy (IE), which represents the time variation of the energy.An integration along the time span yields the marginal Hilbert spectrum h(ω), which provides information on how the frequency is distributed over the entire span.
The degree of stationarity DS(ω) is often used to investigate the stationarity and periodicity of the data; it is defined as follows: where T is the entire time span.The volatility V (t, T ) is defined as the ratio of the sum of certain IMF components S h (t) to the original signal S(t).
Here we use the summation of the residual and all the IMFs except for the first one as S h (t): where n is the number of IMFs.

The gap-filling of the monthly average ozone data
To perform the HHT analysis, a complete, even-spaced data series is required.Hence we need to fill the gaps in the monthly average surface ozone mixing ratio data.In our monthly ozone time series, gaps of 1-6 months can be found in 1997, 1998, 1999 and 2002.If the gap is small and occurs in between the ozone seasonal low and peak value, then a spline interpolation suffices.However, this is not the case for some gaps.In 1997 and 1998, the gaps occurred during summertime, when the seasonal peak of ozone mixing ratio was expected.In 2002, the gap continued on to winter, when the lowest ozone mixing ratio was expected.A simple spline interpolation would have underestimated the seasonal peak value and overestimated the seasonal low.Hence, we applied the following method to fill the gaps.First, the monthly mean ozone time series from 1994 to 2013 is shaped into an array O 3 (i, j ) of the size [20 years × 12 months], where i = 1994, . .., 2013 and j = 1, . .., 12.
The gaps in O 3 (i, j ) are filled by applying a spline interpolation on each row of the array: O 3,spline (1994, . . ., 2013, j ) = spline(O 3 (1994, . .., 2013, j )), j = 1, . .., 12. (13) In this way, both the average value of ozone mixing ratio at a certain month and the overall ozone trend will be considered.A complete data set of average monthly ozone mixing ratio can then be recreated by using interpolated data only on months of missing observation data: Our method could yield a reasonable interpolated time series with both seasonal low and peak values occurring at the right time of year.

Season-diurnal variation characteristics of ozone
The average season-diurnal variation of surface ozone during 1994-2013 is displayed in Fig. 3 with the monthly average local times associated with the diurnal minimum ozone and maximum zonal wind.The seasonal maximum ozone occurs during summer, with an average peak in June-July, while the minimum is found in winter (Fig. 3a), which will be discussed in detail in Sect.3.2.Daily maximum ozone usually occurs during nighttime, while the daily minimum ozone is found around noontime, on average at 12 a.m., Beijing local time (Fig. 3c).Ma et al. (2002b) suggest that the WLG station is mostly influenced by boundary layer (BL) air that is brought up through an upslope flow during the day, while a downslope flow brings down free tropospheric (FT) air during the night.The BL air masses are typically characterised by lower ozone mixing ratios in comparison with FT air masses; hence the occurrence of the daily ozone minimum value indicates the time when the BL is fully developed and the air within is well mixed.
From Fig. 3b it can be denoted that the occurrence time of the daily minimum ozone mixing ratio (red dots) shows a significant annual variation similar to that of the maximum zonal wind velocity (white dots), with the former occurring 1-2 h earlier than the later.Due to the annual variation of the BL development, the daily minimum ozone should occur earlier in the day during warm seasons and later in the day during cold seasons.This phenomenon can indeed be confirmed by Fig. 3b; however, the ozone minimum of June-August seems to occur later than expected.This phenomenon is not seen in the season-diurnal variation of horizontal or vertical wind speeds, indicating that it is not caused by the BL development.A possible explanation might be that the photochemical production of ozone is enhanced at early noon during summertime, leading to a delayed noontime minimum.The in situ ozone production/destruction in different seasons is not well quantified at the moment.Previous studies focused on modelling the photochemical net production in winter and summer, reporting net ozone production in winter and destruction in summer (Ma et al., 2002b).Observation results, however, suggest that there should be net photochemical production during summertime (Wang et al., 2006b).Hence there is a need for more investigation into the cause for such a phenomenon.

Season-annual variation characteristics of ozone
Figure 4 displays the season-annual variation of surface ozone during 1994-2013.Again, the ozone mixing ratios peak in summer and are lowest during winter (Fig. 4b), with an average seasonal peak occurring in June during 1994-2013 (Fig. 4c).Previous studies reported the same seasonal ozone pattern, but attributed the summertime peak to differ- ent causes, e.g. more frequent STE events (Ding and Wang, 2006;Tang et al., 2011), enhanced vertical convection (Ma et al., 2005), long-range transport from eastern-central China, central-southern Asia or even Europe during summer (Zhu et al., 2004) and stronger cross boundary transport and vertical convection during the East Asian summer monsoon season (Yang et al., 2014).From Fig. 2c it can be noted that nighttime subsiding wind is indeed strongest in summer, which supports the hypothesis of downward transport of ozone.Zheng et al. (2011) argued that STE reaches maximum strength in spring and shows a decline in late spring based on 10 Be / 7 Be measurements, indicating that the continuous ozone increase in summer is caused by the photochemical production.The long-term variation of the annual average ozone exhibits a clear increasing trend (Fig. 4a).A 2-4-year cycle seems to exist within the long-term variation of surface ozone.Previous study has shown that there is a quasibiannual oscillation (QBO) within the total ozone column density over the Tibetan Plateau, which is in antiphase with the QBO of the tropical stratospheric winds, exhibiting a 29month cycle (Ji et al., 2001).The influence of the QBO could extend to WLG station at the 3.8 km altitude via STE.Thus, surface ozone at WLG might also have a QBO with a similar periodicity, which is related to that of the total ozone column.The periodicity within the surface ozone data will be further discussed in Sect.3.4.

Long-term trends of ozone
The trends of monthly average all-day, daytime and nighttime ozone during 1994-2013 are displayed in Fig. 5a1 3) summer (JJA), ( 4) autumn (SON) and ( 5) winter time average all day (a), daytime (b) and nighttime (c) surface ozone mixing ratio during 1994-2013 (black solid curves or black circles) and its trend (red lines: dotted line stands for the linear variation and solid line stands for the Theil-Sen trend estimates).
Table 1.The linear slope, 95 % confidence interval (in ppbv year 1 ) and the p values (in parentheses) of all-year and seasonal average surface ozone mixing ratio for the all-day data and for the daytime and nighttime data subsets during 1994-2013.respectively.Ozone data in Fig. 5b1 and c1 are the subsets of data from the daytime and nighttime ranges determined in Sect.2.2 based on the zonal wind information.The increase in surface ozone in the past two decades is evident in all three data subsets, with a slightly stronger increase in the nighttime data.The linear trends for all-day, daytime and nighttime ozone mixing ratios reached 0.25 ± 0.17, 0.24 ± 0.16 and 0.28 ± 0.17 ppbv year −1 , respectively, while the Theil-Sen trend estimates reached 0.18, 0.17, and 0.19 ppbv year −1 , respectively.The Theil-Sen trend estimate is smaller than the linear regression slope, mainly because the linear regression method is influenced by the seasonality within the data.However, both methods yielded statistically significant increasing trends.
To further investigate the trend of ozone in different seasons, the trend of seasonal average ozone during 1994-2013 was calculated and is shown in Fig. 5a-c (2-5).After elim-inating the seasonality in the data, the linear least squares fitting slopes and Theil-Sen trend estimates yielded very similar results; thus only the linear slopes and p values are listed in Table 1.The strongest increase in surface ozone is found in autumn (SON), followed by spring (MAM), respectively, reaching 0.28 ± 0.11 and 0.24 ± 0.11 ppbv year −1 in the seasonal average of all-day ozone mixing ratios.In comparison, summer (JJA) and winter (DJF) both show much weaker increasing trends, with rates of 0.15 ± 0.19 and 0.14 ± 0.09 ppbv year −1 , respectively, and the summertime trend cannot even reach a confidence level of 95 %.In summer the daytime increasing rate is significantly lower than the nighttime one, respectively reaching 0.07 ± 0.18 and 0.22 ± 0.20 ppbv year −1 .The nighttime slope reaches the confidence level of 95 %, while the daytime slope is statistically insignificant.
Previous investigations on the air-mass origin of WLG have shown that WLG is mostly governed by western and northwestern air masses, air masses coming from the eastern sector takes up only 2, 5 and 8 % in winter, spring and autumn, respectively (Zhang et al., 2011).However, a significant percentage (30 %) of air masses come from the eastern direction during summertime.Since the two major cities in the vicinity of WLG are both in the east, summertime is believed to be the season in which WLG is most influenced by nearby anthropogenic activities.From the diurnal variation of the horizontal wind speeds (Fig. 2a-b) it can be discerned that daytime winds are weak northerly winds, while nighttime winds are rather strong north-easterly winds, which are more in favour of transporting anthropogenic pollution to WLG.
As already mentioned before in Sect.3.2, some researchers believe that STE is also most frequent in summer at WLG (Ding and Wang, 2006).During the night the WLG site is governed by downwards winds, which may bring down air with high ozone mixing ratios from above.Hence, an increase in the frequency of STE events would also result in increasing nighttime ozone mixing ratios in summer.Whether it is anthropogenic activities or rather meteorological factors that have led to the distinct daytime and nighttime ozone variation slopes in summer still needs further investigations.
The seasonal peak of ozone in the northern midlatitudes typically occurs in summer over populated continental areas, due to local and regional photochemical production and in late spring for remote continental areas, due to both enhanced stratospheric input and photochemical production in that season (Monks, 2000;Parrish et al., 2013).Recently, a shift in the seasonal peak towards an earlier time in year has been observed at several high elevation sites (Cooper et al., 2014;Parrish et al., 2013).Unlike other remote sites in the northern mid-latitudes, the seasonal ozone peak at WLG occurs in summer.However, the largest increase in ozone mixing ratio is found in autumn rather than in summer.Lin et al. (2014) also reported significant increasing ozone trends in autumn rather than spring at the Mauna Loa Observatory in Hawaii in the past 4 decades and attributed this phenomenon to strengthened ozone-rich air mass transport from Eurasia.The fact that we observed the largest ozone increase in autumn is possibly linked to changes in atmospheric circulation.Details will be discussed in the companion paper.
Here we present a comparison of our ozone trends in different seasons with those at other high altitude (> 1200 m a.s.l.) sites in the Northern Hemisphere, which have been reported in literature and are based on time ranges similar to that of WLG (Table 2).The stations have been sorted by latitude.The low-latitude sites, Mauna Loa and Izaña, both show increasing trends (0.31 ± 0.07 and 0.14 ± 0.05 ppbv year −1 ) during 1991-2010 (Oltmans et al., 2013).Lin et al. (2014) compared the ozone levels at the Mauna Loa site in Hawaii during the period of 1995 to 2011 to that of 1980 to 1995, and discovered a strong increase dur-ing summer and autumn.The mid-latitude stations exhibit inconsistent trends.Significantly positive trends were detected at Mt Happo, Japan (0.65 ± 0.32 ppbv year −1 , Cooper et al., 2014), at the Rocky Mountains National Park site, USA (0.33 ± 0.05 ppbv year −1 , Oltmans et al., 2013), at Lassen Volcanic National Park, USA (0.27 ± 0.13 ppbv year −1 , Cooper et al., 2014.) and at Jungfraujoch, Switzerland (0.32 ± 0.18 ppbv year −1 , Cui et al., 2011).Tarasova et al. (2009) found evidence for increased stratospheric contribution to surface ozone at Jungfraujoch.The strongest increase at Jungfraujoch was detected in winter and the weakest in summer.Gilge et al. (2010) also reported increased wintertime ozone at two other alpine sites in central Europe during 1995-2007. Cooper et al. (2014) reported significant daytime increasing trends at Lassen Volcanic National Park during spring (0.39 ± 0.15 ppbv year −1 ) and winter (0.21 ± 0.14 ppbv year −1 ).Lin et al. (2015b) reported an increasing trend of 0.31 ± 0.21 ppbv year −1 in springtime free-tropospheric ozone over western North America during 1995-2014; however, by shutting off North American emissions in the model and focusing on the subset of ozone associated with Asian influence (also possibly mixed with stratospheric intrusions), the baseline ozone revealed a more significant increasing rate of 0.55 ± 0.14 ppbv year −1 during 1992-2012.No significant trends were found at Gothic, USA, Pinadale, USA and Zugspitze, Germany.Negative trends were revealed at Kislovodsk, Russia (−0.37 ± 0.14 ppbv year −1 , Tarasova et al., 2009) and the Whiteface Mountain Summit site, USA (−0.22 ± 0.06 ppbv year −1 , Oltmans et al., 2013).Tarasova et al. (2009) attributed the strong decrease in ozone at Kislovodsk to control measures of Europe and the breakdown of the former USSR.Both the strong increasing and decreasing trends at Jungfraujoch and Kislovodsk were mostly caused by the variation in ozone mixing ratios in the 1990s.The positive trend at Jungfraujoch during the 1990s was strongest in spring and weakest in summer and autumn, while the reduction at Kislovodsk was strongest in summer and weaker in autumn and winter (Tarasova et al., 2009).After 2000, the eastern US showed a significant decrease due to the implementation of NO x emission control measures, while ozone mixing ratios at the other sites in the northern mid-latitudes have entered a steady stage with either slow or no growth (Tarasova et al., 2009;Oltmans et al., 2013).
In comparison, WLG shows a continuous rise of ozone mixing ratio throughout the past two decades and the most significant positive trends appear in autumn and spring, unlike the other mid-latitude stations.

Hilbert-Huang spectral analysis of surface ozone at WLG
The long-term variation of surface ozone may be the result of changes in emissions of ozone precursors, but may also be caused by year-to-year fluctuations or multiyear oscillations To unravel the potential oscillations on different timescales in the ozone time series, we performed an HHT analysis on the ozone data from WLG using the method described in Sect.2.4.Our effort is the first time that the HHT method has been applied in the analysis of atmospheric composition data.The first step of this analysis was the EMD filtering of the time series of monthly average ozone mixing ratio.The results of the EMD are shown in Fig. 6.The monthly average ozone signal could be decomposed into five IMFs with different characteristic timescales.The lowest-order IMF (c1) shows an oscillation with the highest frequency.The second IMF (c2) shows the seasonal variation in the ozone signal.Panel c3 reveals 3-4-year oscillations, c4 shows 7-year oscillations and the highest-order IMF (c5 in Fig. 6f) shows the longest oscillations pattern, with a quasi-11-year periodicity.
A segmentation analysis was performed by finding the local extrema of c5.The total time span can be separated into 4 segments, as indicated by the dotted lines in Fig. 6a.The slopes of the segments of c5 indicate whether the value is increasing or declining.To determine the significance of the trend, the modified Mann-Kendall trend test was performed on each segment and the results are given in Table 3.The first segment lasted 3 years (from August 1994 to June 1997) and revealed no significant trend (z = 1.42), with an increas-  ing slope of 0.27 ppbv year −1 .The second segment lasted for 5 years (from July 1997 to May 2002) and displayed a significant upward trend (z = 3.66), with an increasing slope of 0.42 ppbv year −1 .Afterwards the increase of the ozone mixing ratio at WLG slowed down in segment 3, lasting 6 years (from June 2002 to April 2008), with an increasing slope of 0.30 ppbv year −1 ; however, the increasing trend remained significant (z = 3.57).In the last segment (from May 2008 to the end of July 2013), the significant upward trend continued (z = 3.65) with a larger increasing slope (0.36 ppbv year −1 ) than that in segment 3.
Overall, surface ozone mixing ratio at WLG increased continuously from 1997 to 2013. Figure 7a shows the anomaly of the interpolated monthly average ozone during 1994-2013, its overall trend (represented by c5 + r in Fig. 6) and its variation on a scale of 7-year or longer (represented by c4 + c5 + r in Fig. 6).The corresponding variation slopes of the overall trend and the 7-year or longer variation are depicted in Fig. 7b.The overall trend confirms the continuous increase since January 1997.The two largest slopes are respectively detected in May 2000 and October 2010.The 7-year or longer trend line displays a rise in ozone after August 1996, which reaches a maximum increasing speed in September 2003.Afterwards, the increase slows down and turns into a decreasing trend in September 2005.After January 2009, ozone mixing ratios went up again, reaching a maximum increasing speed in December 2010.
The Hilbert energy spectrum is depicted in Fig. 8d, along with the volatility, instantaneous energy (IE) and the degree of stationarity (DS) (Fig. 8b, c and e).Both the volatility and the IE reflect the variation of energy with time.Compared to the mean IE, which represents the temporal variation of the frequency averaged energy, volatility rather focuses on the ratio of the variation of certain signals to the total signal.Peaks in the mean IE can be found in 1994-1995, 2000-2001, 2003, 2008and 2013 (Fig. 8c (Fig. 8c), which correspond to the high ozone mixing ratio values in the data.High values of volatility are found around 2003, 2008 and 2012 (Fig. 8b), which mostly agree with those of the IE.The cause for these high anomalies still needs to be investigated.
The DS corresponding to each frequency, as displayed in Fig. 8e, can provide information on the underlying periodicity within the original signal.The smaller the DS is, the more stationary the data are at this frequency.Lower DS values are observed in the low frequency part.A local minimum at the frequencies between 0.08 and 0.12 can be found, which corresponds to the annual cycle of ozone.Other local minima are found at even lower frequencies, corresponding to 2.5a, 3.5a, 7a and 11a cycles.Among all the known atmospheric factors that have an impact on the ozone mixing ratio at WLG, QBO, ENSO, etc., could be responsible for these periodicities.
Overall, the HHT analysis is able to detect variations in surface ozone trends during different periods, and is successful in finding the anomalies and periodicities within the data.Results of this analysis can further facilitate the attribution of the variations of surface ozone at WLG to the influencing factors, which will be discussed in the companion paper.

Summary
In this paper we present the characteristics, trends and periodicity of surface ozone mixing ratio at a global baseline GAW station in the eastern Tibetan Plateau region (Mt Waliguan) during the past two decades.The trends and periodicity of ozone are investigated using a modified Mann-Kendall test and an adaptive method (Hilbert Huang transform) that is suited for analysing non-stationary and nonlinear natural processes.
While confirming the reported diurnal and seasonal characteristics of surface ozone at WLG, our study reveals a relationship between the seasonality in mountain-valley breeze and the seasonal shift in the occurrence time of daily maximum and minimum ozone at the site.Based on this relationship, season-dependent daytime and nighttime periods are defined for separately analysing the daytime and nighttime trends of surface ozone.Both daytime and nighttime surface ozone have increased significantly at WLG. Autumn and spring revealed the largest increase rates, while summer and winter showed relatively weaker increases.A significant daytime and nighttime difference in trend could only be found in summer, where nighttime ozone was significantly increasing and daytime ozone had no significant trend.Results of the HHT spectral analysis confirm the increasing trends in surface ozone mixing ratio and further identify four different stages with different increasing rates.The overall trend indicates that the largest increase occurred around May 2000 and October 2010.The ozone signal can be decomposed into five intrinsic mode functions with different timescales.2-4year, 7-and 11-year periodicities are found within the data, the cause of which still needs further investigation.
The results obtained in this work are valuable for related climate and environment change assessments of western China and surrounding areas, and can be used in the validation of chemistry-climate models.As WLG is a high altitude mountain-top site in a remote region, measurements of surface ozone and other species can well represent a largescale situation.Previous air mass origin and modelling studies (Zhang et al., 2011;Li et al., 2014) suggest that WLG is mostly under the influence of transport from the north-west direction; hence the upward trend in ozone might be an indication of impact of transport from that direction.Since eastern China is downwind of WLG, our results imply that under rising baseline ozone conditions, even greater effort needs to be applied to reducing ozone precursors in eastern China in order to improve ozone air quality.In the second part of our study, influencing factors or potential causes of the observed long-term trends of surface ozone at WLG will be addressed and discussed.

Figure 1 .
Figure 1.The location of the Mt Waliguan GAW site and the two major cities in its vicinity.The shading stands for the topographic height.

Figure 2 .
Figure 2. The average season-diurnal variation of surface zonal (a), meridional (b) and vertical (c) wind velocity on top of Mt Waliguan during 1995-2013.The monthly average hour associated with the diurnal maximum zonal wind speed is given by the white dots, the daytime range is provided by the white dashed lines, which covers 6 h centred around the white dots.

Figure 3 .
Figure 3.The average seasonal variation (a), season-diurnal variation (b) and diurnal variation (c) of ozone during 1995-2013.The red and white dots indicate the monthly average local times associated with the diurnal minimum ozone and the diurnal maximum zonal wind (U max ), respectively.

Figure 6 .
Figure6.The interpolated monthly average ozone mixing ratio at WLG from 1994 to 2013 (the interpolated data given in dashed lines, a) and its intrinsic mode functions c1-c5 (b-f, from the lowest to the highest-order IMF) and its residue, r (g).The time segments in (a) were determined by the slope of the c5.The red slashed lines are the Theil-Sen trends and the numbers are the Theil-Sen trend estimates (in ppbv year −1 ).
Figure 7. (a)The anomaly of the interpolated monthly average ozone (black line, with the dashed line segments representing values interpolated using the method in Sect.2.5), the sum of last IMF and the residual (c5 + r, red line), and the sum of the last two IMFs and the residual (c4 + c5 + r, blue line); (b) the slope of the sum of last IMF and the residual (c5 + r, red line) and the sum of the last two IMFs and the residual (c4 + c5 + r, blue line).

Figure 8 .
Figure 8.The interpolated monthly average ozone mixing ratio signal at Mt WLG during 1994-2013 (a), the volatility (b), the normalised mean value of the instantaneous energy (red lines: ±2σ ) (c), Hilbert energy spectrum (d) and the degree of stationarity (e).

Table 2 .
The linear slopes (in ppbv year −1 ) and the 95 % confidence intervals of all-year and seasonal average surface ozone mixing ratio at WLG and other North Hemispheric high-altitude sites.

Table 3 .
Modified Mann-Kendall trend test on segments based on the last IMF.