Analysis and attribution of total column ozone changes over the Tibetan Plateau during 1979–2017

Various observation-based datasets have confirmed positive zonal mean column ozone trends at midlatitudes as a result of the successful implementation of the Montreal Protocol. However, there is still uncertainty about the longitudinal variation of these trends and the direction and magnitude of ozone changes at low latitudes. Here, we use the extended Copernicus Climate Change Service (C3S) dataset (1979–2017) to investigate the long-term variations in total column ozone (TCO) over the Tibetan Plateau (TP) for different seasons. We use piecewise linear trend (PWLT) and equivalent effective stratospheric chlorine loading (EESC)-based multivariate regression models with various proxies to attribute the influence of dynamical and chemical processes on the TCO variability. We also compare the seasonal behaviour of the relative total ozone low (TOL) over the TP with the zonal mean at the same latitude.

Both regression models show that the TP column ozone trends change from negative trends from 1979 to 1996 to small positive trends from 1997 to 2017, although the later positive trend based on PWLT is not statistically significant. The wintertime positive trend starting from 1997 is larger than that in summer, but both seasonal TP recovery rates are smaller than the zonal means over the same latitude band. For TP column ozone, both regression models suggest that the geopotential height at 150 hPa (GH150) is a more suitable and realistic dynamical proxy compared to a surface temperature proxy used in some previous studies. Our analysis also shows that the wintertime GH150 plays an important role in determining summertime TCO over the TP through persistence of the ozone signal. For the zonal mean column ozone at this latitude, the quasi-biennial oscillation (QBO) is nonetheless the dominant dynamical proxy.

We also use a 3-D chemical transport model to diagnose the contributions of different proxies for the TP region. The role of GH150 variability is illustrated by using two sensitivity experiments with repeating dynamics of 2004 and 2008. The simulated ozone profiles clearly show that wintertime TP ozone concentrations are largely controlled by tropics to midlatitude pathways, whereas in summer variations associated with tropical processes play an important role. These model results confirm that the long-term trends of TCO over the TP are dominated by different processes in winter and summer. The different TP recovery rates relative to the zonal means at the same latitude band are largely determined by wintertime dynamical processes.

1 Introduction

The Tibetan Plateau (TP), also known as the third pole, is an area very sensitive to global climate change. It exerts important thermal and dynamical effects on the general circulation and climate (Yanai et al., 1992; Ye and Wu, 1998). Furthermore, climate changes over the TP have a significant impact on the distribution of stratospheric ozone.

There is well-established observational evidence of a persistent total column ozone low (TOL) centred over the TP (e.g. Zhou et al., 1995; Zheng et al., 2004; Bian et al., 2006; Tobo et al., 2008). Zou (1996) found that the largest ozone deficit over TP occurs in May, while the smallest deficit occurs in wintertime. Ye and Xu (2003) proposed that the high topography and the elevated heating source associated with thermally forced circulations are the two main reasons for its occurrence. Other studies have also suggested that the thermal–dynamical forcing of the TP, for example, by air expansion, uplifting of the tropopause, thermal convection and monsoon circulation, makes a dominant contribution to the TOL, especially in summer (e.g. Tian et al., 2008; Bian et al., 2011; Guo et al., 2012, 2015). However, the exact coupling pathways between the thermal–dynamical forcing and long-term total column ozone (TCO) changes during different seasons are still not well established.

TP column ozone trends can be significantly affected by internal variability. Zou (1996) reported strong negative TCO trends over Tibet for the 1979–1991 time period and, in subsequent studies, analysed the effects of the quasi-biennial oscillation (QBO) and the El Ninõ–Southern Oscillation (ENSO) (e.g. Zou et al., 2000, 2001). Zhou and Zhang (2005) presented decadal ozone trends over the TP using the merged Total Ozone Mapping Spectrometer/solar backscatter ultraviolet (TOMS/SBUV) ozone data for the period 1979–2002 and found that the downward trends are closely related to the long-term changes of temperature and geopotential height. Zhou et al. (2013) found substantial downward ozone trends over the TP in the longer 1979–2010 TOMS/SBUV record during the winter–spring seasons. They also showed that long-term ozone variations are largely correlated with thermal–dynamical proxies such as lower stratospheric temperature, with its contribution reaching around 10 % of the total ozone change. Zhang et al. (2014) indicated that the TOL over the TP in winter deepened during the period 1979–2009 and that thermal–dynamical processes associated with the TP warming may account for more than 50 % of the TCO decline in this region.

While previous studies have demonstrated the contributions of dynamical processes to the long-term ozone variation over the TP (e.g. Zhou et al., 2013; Zhang et al., 2014), a better proxy is needed to explain the dynamical influence for this region. The geopotential height in the free atmosphere is an important thermal–dynamical proxy that not only conveys information about the thermal structure of the atmosphere but also serves as an indicator of synoptic circulation changes (Christidis and Stott, 2015). The natural and anthropogenic contributions to the changes in geopotential height (GH) establish the coherent thermal–dynamical nature of externally forced changes in the regional climate system, which provides the basis for the validation of climate models. In this study, the GH at 150 hPa over the TP is used as a new thermal–dynamical proxy which incorporates coupling between the local TP circulation and various tropospheric teleconnection patterns and represents the tropospheric dynamical influence more realistically.

With the extended Copernicus Climate Change Service (C3S) TCO time series available from 1979 to early 2018, the aim of this paper is to study the long-term trend and variability in ozone over the Tibetan region. Based on statistical regression analysis of C3S ozone data and TOMCAT/SLIMCAT three-dimensional (3-D) chemical transport model (CTM) simulations, the contributions of different influencing variables, including the local thermal–dynamical proxy (GH), are diagnosed to help understand the long-term ozone variability in different seasons and over different areas.

The layout of the paper is as follows. Section 2 introduces the C3S ozone dataset and TOMCAT/SLIMCAT model used for the analysis of the total ozone variability. The long-term TCO time series and TOL over the TP region are presented in Sect. 3. Regression models as well as analysis of the contribution of different proxies to the total ozone variations in different seasons are given in Sect. 4. Section 5 discusses our 3-D model sensitivity experiments and is followed by our summary and conclusions in Sect. 6.

2 Data

2.1 Ozone dataset from C3S

High-quality observational-based datasets are necessary for better quantification of decadal TCO trends. This is because interannual variability can cause variations of up to 20 %, whereas ozone trends are generally less than half a percent. As the lifetime of most satellite instruments is less than two decades, merged satellite datasets are widely used to determine long-term ozone trends. These datasets are created by combining total ozone measurements from different individual instruments to provide global coverage over several decades (e.g. Frith et al., 2014). SBUV provides nearly continuous satellite-based measurements of total ozone to analyse trends. The variations from all the instruments are within 2 % relative to the ground-based data at all latitudes (Labow et al., 2013). SBUV-merged data are obtained from https://acd-ext.gsfc.nasa.gov/Data_services/merged/instruments.html (last access: 16 April 2020). However, this merged satellite dataset is available as zonal mean values at 5^∘ latitude resolution and therefore is not well suited to study relatively small geographical areas such as the TP.

Hence, here we use the total column ozone from the C3S which is produced by the European Centre for Medium-Range Weather Forecasts (ECMWF). For a detailed description and data availability, see https://cds.climate.copernicus.eu/cdsapp#!/dataset/satellite-ozone?tab=overview (last access: 10 July 2020). In brief, these are monthly mean gridded data that span from 1970 to present. They are created by combining total ozone data from 15 satellite sensors including GOME (1995–2011), SCIAMACHY (2002–2012), OMI (2004–present), GOME-2A/B (2007–present), BUV-Nimbus4 (1970–1980), TOMS-EP (1996–2006), SBUV-9, -11, -14, -16, -17, -18, -19 (1985–present) and OMPS (2012–present). The horizontal resolution of the assimilated product after January 1979 is ${\mathrm{0.5}}^{\circ }\times {\mathrm{0.5}}^{\circ }$. The document describing the methodology adopted for the quality assurance in the C3S-Ozone procurement service, with detailed information about the ground-based measurements used to verify satellite observations, the specific technical project implemented to compare the gridded (level-3) and assimilated (level-4) data, and the metrics developed to associate validation results with user requirements, can be downloaded from https://cds.climate.copernicus.eu/cdsapp#!/dataset/ozone-monthly-gridded-data-from-1970-to-present?tab=doc (last access: 13 November 2019). The strength of this dataset is the long-term stability of the total column monthly gridded average product that is below the 1 % per decade level. Systematic and random errors in these data are below 2 % and 3 %–4 %, respectively, hence making it suitable for long-term trend analysis.

We use four different area-weighted total ozone time series during 1979–2017: TP (27.5–37.5^∘ N, 75.5–105.5^∘ E), zonal TP (full zonal mean for 27.5–37.5^∘ N) as well as zonal mean for latitude bands to the south (10–20^∘ N) and north (40–50^∘ N) of the TP region. These regions represent the tropics and midlatitudes with the TP and zonal TP in the critical zone. We choose them to compare the contribution of different dynamical proxies to their ozone variations especially over the TP region. In this paper, we also use the direct ozone observations from the SBUV series of satellites to validate the results based on C3S.

2.2 TOMCAT/SLIMCAT model

Chemistry–transport models are important tools to investigate how past and present-day ozone-depleting substance (ODS) and greenhouse gas (GHG) concentrations have influenced the ozone layer (e.g. Shepherd et al., 2014; Zvyagintsev et al., 2015). In combination with observed ozone time series, simulations allow the attribution of ozone changes, thus encapsulating our understanding of the fundamental physics and chemistry that controls ozone and its variations (e.g. Chipperfield et al., 2017). TOMCAT/SLIMCAT (hereafter SLIMCAT) is a 3-D offline chemical transport model (Chipperfield et al., 2006), which uses winds and temperatures from meteorological analyses (usually ECMWF) to specify the atmospheric transport and temperatures and calculates the abundances of chemical species in the troposphere and stratosphere. The model has the option of detailed chemical schemes for various scenarios with different assumptions of factors affecting ozone (e.g. Feng et al., 2011; Grooß et al., 2018), including the concentrations of major ozone-depleting substances, aerosol effects from volcanic eruptions (e.g. Dhomse et al., 2015), and variations in solar forcing (e.g. Dhomse et al., 2013, 2016) and surface conditions. For this study, the model has been forced by ECMWF ERA-Interim reanalysis (Dee et al., 2010) and run from 1979 to 2017 at a horizontal resolution of $\mathrm{2.8}{}^{\circ }\times \mathrm{2.8}{}^{\circ }$ with 32 levels (up to around 60 km).

We perform control and sensitivity simulations based on the SLIMCAT CTM to elucidate the impact of dynamical changes on the total ozone variations over the TP region. The control experiment R1 uses standard chemical and dynamical parameters for the time period of 1979–2017, which is identical to the control run of Chipperfield et al. (2017). To understand the special dynamical influences (e.g. GH) on ozone variations over the TP, two sensitivity experiments (R2 and R3) were performed with all configurations the same as R1 except the simulations used annually repeating meteorology for the years 2004 and 2008, respectively. These years were chosen because the 150 hPa GH in wintertime is substantially different, while other dynamical proxies are almost the same for the two years. We also take a 5-year average from model dates in 2004–2008 for each sensitivity experiment to exclude the influence from other time-dependent changes (e.g. chemical processes).

3 TCO and TOL over the TP

Figure 1 shows the TCO time series averaged for December–January–February (DJF) and June–July–August (JJA) seasons during 1979–2017 over the north-TP (40–50^∘ N), south-TP (10–20^∘ N), zonal-TP region (27.5–37.5^∘ N) and the TP region (27.5–37.5^∘ N, 75.5–105.5^∘ E). Zonal mean values from SBUV-merged total ozone data for the same latitude band are also shown in Fig. 1a and b to compare with the C3S datasets. Differences between C3S and SBUV are less than 2 %–3 % throughout the data record and are shown in the Supplement Fig. S1, confirming that there is no long-term drift in the C3S data. As shown in Fig. 1, the magnitudes of interannual variations increase with latitude, with amplitude of DJF ozone variations being much larger than JJA. Besides, Fig. 1c and d show the TP and zonal-TP ozone time series, highlighting much smaller difference in DJF (< 5 DU) compared to about 20 DU difference in JJA. This is consistent with previous studies (e.g. Ye and Xu, 2003; Zhang et al., 2014).

Figure 1 C3S-based total column ozone (TCO) time series averaged for December–January–February (DJF) and June–July–August (JJA) seasons during 1979–2017 over (a, b) the north-TP (40–50^∘ N) and south-TP (10–20^∘ N) regions, and (c, d) the zonal-TP (27.5–37.5^∘ N) and the TP regions (27.5–37.5^∘ N, 75.5–105.5^∘ E). Panels (a) and (b) also show the satellite-based observations from SBUV (solid black lines).

To illustrate TOL characteristics, we calculate the zonal deviations by subtracting the zonal mean total ozone for each latitude band from the TCO at each grid point (Fig. 2). The negative zonal deviations suggest that the TOL centred over the TP exists for all the seasons. As shown in Fig. 1c and d, TOL over the TP is most discernible in summer (JJA), followed by spring (March–April–May, MAM) and autumn (September–October–November, SON), while it is weakest in winter (DJF). The TOL centre also moves from the northwest in spring (MAM) to the south in winter (DJF). The mechanisms for these seasonal differences over the TP are very different in winter and summer. In wintertime, the plateau geographic effect is somewhat less effective in modifying the lower stratospheric circulation as the subtropical jet moves southwards (e.g. Luo et al., 2019). During summertime, the TP is an elevated heating source causing thermally forced anticyclonic circulation. The upper-level Asian summer monsoon anticyclone coupling with deep convection over the TP can potentially transport ozone-poor air from the boundary layer upward into the upper troposphere and lower stratosphere (Liu et al., 2003; Gettelman et al., 2004; Randel and Park, 2006; Bian et al., 2011). Seasonal variations in TCO over the TP and zonal-TP region are shown in Fig. S2. The wintertime ozone buildup and steady summertime ozone decline are evident over both regions. However, the high topography of the TP causes an earlier phase (about 1 month) and smaller amplitudes in TCO variability over the TP. The different TOL magnitudes in different seasons could be associated with the fact that wintertime ozone concentrations are largely controlled by large-scale dynamical processes, while photochemical loss is the only dominant process in summer. Thus, it is necessary to analyse the influences of the chemical and dynamical processes (e.g. EESC, solar, QBO and the local thermal–dynamical proxy) on the total ozone variability under different atmospheric conditions.

Figure 2 Latitude–longitude cross section of the zonal ozone deviations for (a) March–April–May (MAM), (b) June–July–August (JJA), (c) September–October–November (SON) and (d) December–January–February (DJF) seasons based on C3S total ozone dataset for the 1979–2017 time period. The solid and dashed contours represent the positive and negative zonal deviations. The contour interval is 5 DU. The TP region (27.5–37.5^∘ N, 75.5–105.5^∘ E) is marked by the white rectangle.

4 Multivariate linear regression based on C3S

4.1 Regression models

Multivariate linear regression models are widely used to assess the long-term total ozone trends. In these models, proxies are included to separate the influences of important short- and long-term processes on trend determination. Typically, multivariate linear regression models use equivalent effective stratospheric chlorine (EESC) or piecewise linear trend (PWLT) terms for long-term ozone trends (e.g. Reinsel et al., 2002; Nair et al., 2013; Chehade et al., 2014). EESC is a measure of the total inorganic chlorine and bromine amounts in the stratosphere, which drive chemical ozone depletion. Previous studies have indicated that EESC is a main contributor to the long-term global ozone decline and the trend changes after the end of 1990s (e.g. Newman et al., 2004; Fioletov and Shepherd, 2005; Randel and Wu, 2007). We use this method to study the effect of EESC on the long-term ozone variations over the TP and the other zonal regions. A PWLT-based regression method is used to statistically analyse robustness of decreasing and recovery trends in the total ozone before and after the EESC peak in 1997. Our aim is to clarify statistical significance of the key processes responsible for the total column ozone variations over the TP in different seasons using two different regression models.

Traditional explanatory proxies to account for influence of chemical and dynamical processes, include the F10.7 solar flux for the 11-year solar cycle, quasi-biennial oscillation (QBO) at 30 and 10 hPa (QBO30 and QBO10), and ENSO (e.g. Baldwin et al., 2001; Camp and Tung, 2007). Some studies also include aerosol optical depth at 550 nm, to account for ozone loss due to volcanically enhanced stratospheric aerosol loading after El Chichón (1982) and Mt. Pinatubo (1991) eruptions. To account for dynamical variability typical indices are wind near vortex, Arctic oscillation (AO) index, Eliassen–Palm flux or eddy heat flux (e.g. Chehade et al., 2014 and references therein). Due to unique nature of TP orography, the local thermal–dynamical forcing, e.g. the geopotential height at 150 hPa (GH150) and the surface temperature (ST), are also considered as dynamical proxies. We calculate the GH150 and ST over the TP and zonal latitude bands from the ECMWF ERA-Interim reanalysis dataset obtained via https://apps.ecmwf.int/datasets/data/interim-full-daily/ (last access: 10 January 2020). Radiosonde-based GH150 data from a nearby Lhasa station (http://weather.uwyo.edu/upperair/seasia.html, last access: 30 December 2019) are also used for comparison with ECMWF data (Fig. S3). The statistically significant correlation (0.96) validates our use of the ECMWF GH150 data for the TP region.

Table 1 Correlation coefficients for the DJF mean TCO and explanatory variables over the TP during 1979–2017.

${}^{***}$ 99 % confidence level. ${}^{**}$ 95 % confidence level.

Due to the large differences in scales and units of the explanatory variables, we have standardized all the time series to ensure each factor contributes approximately proportionately to the final ozone variations. The transformation does not change the correlation and fitting results. Another important criterion for multivariate regression model is that explanatory variables should not be highly correlated with each other. Table 1 shows the correlation values for the DJF mean TCO (over the TP region) and explanatory variables during 1979–2017 (a similar analysis for JJA is presented in Table S1 in the Supplement). The local thermal–dynamical proxy (GH150 or ST over the TP) is de-trended before being used in the regression models. As shown from the correlation analysis, the DJF mean TCO has significant negative correlations with EESC, QBO and GH150. The solar variability proxy (F10.7 index) is strongly correlated with the AO (0.398) time series. Also, the GH150 time series shows relatively stronger correlation with the ENSO (−0.468), AO (0.374) and ST (0.618) time series. We also find that aerosol and ENSO are correlated (0.372). Hence, to avoid any aliasing effects, we omit the data after the El Chichón (1982, 1983) and Mt. Pinatubo (1991, 1992) volcanic eruptions. AO is also removed, as it shows strong correlation with the solar and GH150 proxies. As for the other partially correlated proxies (ENSO, ST and GH150), we make three groups of independent variables to analyse the TCO variations and compare the corresponding regression results under different situations:

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}\mathrm{TCO}\left(t\right)& =C_{\mathrm{0}}+C_{\mathrm{1}}\cdot \mathrm{EESC}\left(t\right)+C_{\mathrm{2}}\cdot \mathrm{solar}\left(t\right)+C_{\mathrm{3}}\\ & \cdot \mathrm{QBO}\left(t\right)+C_{\mathrm{4}}\cdot \mathrm{ENSO}\left(t\right)+\mathit{\epsilon }\left(t\right)\end{array}\text{(2)}& {\displaystyle }\begin{array}{rl}\mathrm{TCO}\left(t\right)& =C_{\mathrm{0}}+C_{\mathrm{1}}\cdot \mathrm{EESC}\left(t\right)+C_{\mathrm{2}}\cdot \mathrm{solar}\left(t\right)+C_{\mathrm{3}}\\ & \cdot \mathrm{QBO}\left(t\right)+C_{\mathrm{4}}\cdot \mathrm{ENSO}\left(t\right)+C_{\mathrm{5}}\cdot \mathrm{ST}\left(t\right)+\mathit{\epsilon }\left(t\right)\end{array}\text{(3)}& {\displaystyle }\begin{array}{rl}\mathrm{TCO}\left(t\right)& =C_{\mathrm{0}}+C_{\mathrm{1}}\cdot \mathrm{EESC}\left(t\right)+C_{\mathrm{2}}\cdot \mathrm{solar}\left(t\right)+C_{\mathrm{3}}\\ & \cdot \mathrm{QBO}\left(t\right)+C_{\mathrm{4}}\cdot \mathrm{GH}\mathrm{150}\left(t\right)+\mathit{\epsilon }\left(t\right),\end{array}\end{array}$$

where t is a running index corresponding to the years during the period 1979–2017, excluding the 4 years due to the volcanic aerosol loading. QBO herein is equivalent to ($a\times \mathrm{QBO}\mathrm{30}+b\times \mathrm{QBO}\mathrm{10}$). C₀ is a constant for the long-term average. C₁–C₅ represent the time-dependent regression coefficients of each proxy and ε is the residual. In the PWLT regression model, the C₁×EESC(t) term is replaced by ($c\mathrm{1}\times t\mathrm{1}+c\mathrm{2}\times t\mathrm{2}$) in Eqs. (1)–(3) with linear trends (Trend1 and Trend2) in the periods 1979–1996 and 1997–2017, respectively.

4.2 Regression analysis

We apply the multivariate linear regression models to the seasonal mean TCO time series to determine long-term ozone changes over the TP, zonal-TP, south-TP and north-TP zonal bands, respectively. Table 2 lists the adjusted determination coefficients (adj. R²) based on the PWLT regression model for DJF mean TCO time series with three groups of independent explanatory variables over four different regions. Compared to the regression results based on Eq. (1), the additional consideration of the ST proxy in Eq. (2) improves the adjusted R² over all these regions, especially over the TP. By replacing ENSO in Eq. (1), GH150 in Eq. (3) improves the regression fit more significantly for the TP and zonal-TP time series compared to the ST. However, similar improvements are not visible for the non-TP zonal time series. This seems more feasible, as the changes in GH150 represent locally relevant dynamical variability that is modulated by the orography and local circulations over the TP. EESC-based regression results with adjusted determination coefficients are also shown in the Table S2 and are consistent with PWLT-based regression results in Table 2.

Table 2 Adjusted determination coefficients of the PWLT-based regression model for DJF mean TCO over different regions with different proxies.

Using the PWLT-based regression model, we analyse the TCO trends for 1979–1996 and 1997–2017. The fitted signals of the TCO anomalies and explanatory terms in Eq. (3) for DJF and JJA are shown in Fig. 3, and corresponding regression coefficients along with 2σ standard deviations are listed in Table 3. As the summer/autumn ozone variability is much weaker compared to seasonal ozone buildup during winter and spring, the long-term ozone anomalies as well as the contributions from different explanatory variables show much weaker contribution in JJA. Hence, the adjusted coefficient of determination in JJA mean TCO regression (0.61) is also much smaller. As expected, linear trends in both DJF and JJA show a decline over the TP during 1979–1996 (Trend1) and a recovery starting from 1997 (Trend2). Furthermore, the upward trend starting from 1997 in JJA is relatively weaker than that in DJF. EESC-related ozone trends over the TP and zonal-TP region in both seasons are given in Table S3. The TCO trends over the TP, compared to those over the zonal-TP region, show relatively smaller decline and recovery rates before and after 1997. These differences indicate the zonal asymmetry in ozone trends due to longitudinal variations. The comparison between the EESC and PWLT trends shows good agreement, except that EESC trends are statistically significant within 2σ due to the full data record; however, the positive trends (Trend2 term) in PWLT are always non-significant, highlighting complexities in determination of ozone trends at low latitudes.

Except for the linear trends, all the other explanatory proxies (solar cycle, QBO and GH150) contribute significantly to the ozone variations in DJF (above the 99 % confidence level), especially combined contribution from three dynamical proxies (QBO30, QBO10 and GH150) which adds up to 40 DU. As shown later JJA ozone concentration are largely controlled by photochemical ozone loss, contribution from GH150 drops sharply (∼0.27 DU). Hence, the main contributors to the JJA mean TCO variations are linear trends (7.86 DU), solar cycle (4.61 DU) and QBO at 10 hPa (6.56 DU). Results obtained from the EESC-based regression model (not shown) are very similar to those shown in Table 3, confirming the robustness of the results.

Figure 3 (a) PWLT regression results with contributions from linear trends for the 1979–1996 and 1997–2017 time periods, solar cycle, QBO at 30 and 10 hPa, and the GH at 150 hPa in DJF based on C3S during 1979–2017 over the TP region. (b) Similar to panel (a) but with all factors averaged in JJA.

Table 3 PWLT-based regression coefficients and standard deviations for the DJF and JJA mean TCO over the TP during 1979–2017.

To describe quantitatively the contributions of different explanatory proxies to the DJF and JJA mean total ozone variability over different regions, we calculate the percentage ozone change for comparison, as shown in Fig. 4. These contributions using percentage ozone change are represented by Eq. (4):

$$\begin{array}{}\text{(4)}& \mathrm{\Delta }\mathrm{TCO}\left[\mathit{\%}\right]={\displaystyle \frac{max\left(X\right[\mathrm{DU}\left]\right)-min\left(X\right[\mathrm{DU}\left]\right)}{\mathrm{mean}\left(\mathrm{TCO}\right[\mathrm{DU}\left]\right)}}\times \mathrm{100}\mathit{\%},\end{array}$$

where X means the contribution of one proxy (in DU) to the long-term ozone variability. In Fig. 4, the percentage contribution with an error bar indicates the statistical significance within 2σ. During the wintertime (DJF), dynamical proxies (QBO and GH150) exert a significant effect on the total ozone variability over the TP (about 8 % each), while QBO dominates over the zonal-TP region (up to 7 %). However, in summertime (JJA), contributions from dynamical proxies are much smaller, although the contribution from the QBO10 remains above 2 %, the contribution from GH150 almost disappears.

Figure 4 Peak contributions of various explanatory variables to variability in the total ozone column (in %) in (a) DJF and (b) JJA over the TP and the zonal-TP region based on C3S data during 1979–2017. The hatched bars without error bars indicate the contribution is not significant within the 2σ level.

Previous studies have found that changes in GH150 associated with an enhanced South Asian high (SAH) results in significant TCO deviations at 150–50 hPa over the TP (Tian et al., 2008; Bian et al., 2011; Guo et al. 2012). From April onwards, as the SAH advances over the TP, summertime GH150 starts increasing (Fig. S4). Between the TP and zonal-TP region, the GH150 contribution shows a maximum difference in May when the negative TOL is also strongest (Fig. S2), with a correlation coefficient of −0.86 within the 0.001 significance level. Thus, the amplitude of SAH imposes an important constraint on the formation of the summertime TOL over the TP. However, here we find that GH150 makes a major contribution to the TCO variability in wintertime but not in summertime. The sharp contrast between the contributions of the 150 hPa GH in DJF and JJA is an interesting feature and a possible explanation for those differences is discussed below.

The seasonal variability in TCO over the TP (Fig. S2) indicates a marked seasonal cycle with a buildup of total ozone through the winter and a decline through the summer. The correlation of the DJF mean TCO with the subsequent JJA means over the TP during 1979–2017 is 0.44, which is statistically significant above the 95 % confidence level. This significant positive correlation indicates that negative or positive wintertime TCO anomalies over the TP appear to persist through the summer period (as shown in Fig. 5). Table 4 lists the correlation coefficients of TCO variations in a given season of the year with those in subsequent seasons. The correlation decreases from the buildup in winter to the end of summer, and there exists a sharp drop between the summer (JJA) and autumn (SON) which may reflect that dynamical variability is nearly absent during summer months and ozone simply drops off photochemically in a predictable way (Fioletov and Shepherd, 2003). Detailed analysis of the correlation between subsequent months of the year is provided in Table S4.

Figure 5 (a) Correlation map of the DJF and JJA mean TCO based on C3S during 1979–2017. Correlation values in the stippled area are statistically significant above the 95 % confidence level. The white rectangle represents the TP region. (b) Correlation fit between the DJF and JJA mean ozone anomalies (DU) during 1979–2017 over the TP region.

Table 4 Correlation coefficient between ozone values in a given season and the subsequent season.

One lag is 3 months; bolded numbers are statistically significant within 2σ.

Fioletov and Shepherd (2003) highlighted the seasonal persistence of midlatitude total ozone anomalies and indicated that seasonal predictability is applicable for latitudinal belts or large regions only. The seasonal persistence of ozone anomalies over the TP also implies a causal link between the wintertime ozone buildup due to planetary-wave-induced transport and the subsequent chemical loss. The ozone buildup in wintertime when transport dominates is largely modulated by QBO (Holtan and Tan, 1980). However, GH150 represents large part of wintertime variability in the ozone transport. In summertime, as expected, photochemical processes become more important, while dynamical impact from QBO decreases and almost disappears for GH150. Seasonal persistence in TCO anomalies shows that if there is more transport in DJF as represented by GH150 changes, higher ozone values will persist for at least 6 months, even though there is little correlation between summertime ozone anomalies and GH150. This analysis clearly highlights dynamical influence of the wintertime GH150 on the summertime (JJA) ozone concentrations.

5 Model sensitivity simulations

To investigate the role of wintertime GH150 on ozone transport, we use the SLIMCAT 3-D chemical transport model to understand its role under different conditions. The simulated TCO time series obtained from the control experiment R1 are shown in Fig. S5. Overall, modelled TCO is consistent with the C3S-based TCO data although they are biased low. By applying the PWLT regression model in Eq. (3) to the simulated TCO time series, the percentage ozone change from each explanatory proxy is shown in Fig. S6. The simulation results are similar to the C3S regression results, although contributions from most explanatory proxies are larger except for the GH150. This difference is probably due to the coarse model resolution and the inhomogeneities in ERA-Interim data, especially before 2000 (e.g. Dhomse et al., 2011, 2013; McLandress et al., 2014). The contribution of the GH150 proxy to the simulated TCO variations over the TP is statistically significant in DJF but not in JJA. To further elucidate the role GH150 plays in the total ozone variability over the TP, we performed two sensitivity experiments (R2 and R3) with repeating dynamics from the years 2004 and 2008, respectively. For the two years, the wintertime difference in QBO is modest but GH150 is significantly different. We then take a 5-year average based on a time-slice simulation during 2004–2008 for each sensitivity experiment to ensure that other chemical factors (EESC, solar cycle, etc.) are the same between the simulations. Thus, the model settles down with the GH150 as the main proxy that influences the ozone variations over the TP.

Figure 6 (a) Pressure–latitude and (b) pressure–longitude cross sections of DJF mean ozone differences (colours in DU) between the 5-year averaged SLIMCAT sensitivity experiments (R2 and R3) and the GH differences (contours in gpm) between the years 2004 and 2008. (c, d) Similar to panels (a, b) but averaged for JJA. Positive ozone and GH differences are shown with red colours and solid contours, whereas blue colours and dashed contours indicate negative differences. The shaded area shows the TP region averaged over the 75.5–105.5^∘ E longitude band (a, c) and the 27.5–37.5^∘ N latitude band (b, d). The arrows in panels (a) and (c) indicate the TP GH differences influenced by those from the high and low latitudes; the dashed blue and red boxes in panel (b) indicate the negative and positive ozone anomalies over the TP and the Pacific Ocean.

A caveat is that none of the dynamical processes are independent. The GH150 proxy represents the overall tropospheric dynamical influence somewhat realistically as it incorporates coupling between various tropospheric teleconnection patterns and the local TP circulation. To better understand the zonal and meridional pathways, the vertical DJF mean GH differences between the years 2004 and 2008 as well as the 5-year averaged ozone differences (2004–2008) based on the SLIMCAT sensitivity simulations are represented by the contours and colours in Fig. 6. The shaded area shows the TP region. In DJF (Fig. 6a and b), a positive anomaly centre of the GH difference occurs near the 150 hPa pressure level, co-located with a negative ozone anomaly. In JJA (Fig. 6c and d), there are no such clear anomaly centres for mean GH and ozone differences over the TP.

By comparing the GH variation with latitude in Fig. 6a and c, we find that the DJF mean GH differences over the TP are mainly influenced by those over the high latitudes, and in JJA they are mainly influenced by those from low latitudes (as shown by the arrows therein). This may be because the TP lies near the boundary between the tropics and midlatitudes in the troposphere. Due to the movement of the Intertropical Convergence Zone (ITCZ), the TP in wintertime is located in midlatitude band where ozone variability is determined by the tropopause height or folds in the lower stratosphere, while in summer, the TP lies in the tropical band where ozone variability is largely determined by QBO (and QBO-induced circulation) in the mid-stratosphere (Baldwin et al., 2001).

The GH variation with longitude in Fig. 6b and d suggests a tropospheric coupling between the local TP circulation and some tropospheric teleconnection patterns (e.g. ENSO or Walker circulation). As the TP is an elevated heat source, the differences in heat distribution between the plateau and ocean will cause air motions in the zonal and vertical direction. In the normal condition, the pressure gradient force that results from a high-pressure system over the eastern Pacific Ocean and a low-pressure system over the TP will cause the global general circulation (such as the Walker circulation) and therefore affect the ozone distribution. A correlation analysis shows that the GH150 proxy over the TP is in a strong, negative relation to ENSO in DJF, which means during an El Niño event GH150 near the TP also increases, thereby increasing tropopause height, leading to a decrease in TCO over the TP. The positive–negative vertical band-like features in DJF mean ozone differences shown in Fig. 6b seem to closely resemble Walker-circulation-type anomalies (Hu et al., 2016). They also explain why the ozone differences over the TP and the Pacific Ocean are opposite in sign, as indicated by the dashed blue and red boxes therein. Thus, we suggest that wintertime GH fluctuations associated with ENSO events or Walker circulation may play an important role in controlling the TCO variability over the TP. In JJA, however, there are no distinctive features of GH and ozone differences near the TP. As the summertime ozone is less controlled by the dynamical processes (especially GH150), there would not exist such a clear correlation as that in wintertime. Overall, the model results support the hypothesis that wintertime TP ozone variations are largely controlled by tropical to high-latitude transport processes, whereas summertime concentrations result from the combined effect of photochemical decay and tropical processes.

6 Summary and conclusions

In this study, we have analysed the variations and trends of the total column ozone and the relative total ozone low over the Tibetan Plateau in different seasons during the period of 1979–2017. The most recent C3S datasets based on model assimilation of meteorological and ozone observations are used and compared with merged SBUV satellite observations. We use the PWLT- and EESC-based multivariate regression models to analyse the contributions and trends associated with the dynamical and chemical processes that modify the total ozone changes over the TP and zonal areas. In addition to conventional regression proxies (EESC, solar cycle, QBO, ENSO, etc.), we also use the local thermal–dynamical proxy (ST or GH150) to account for the dynamical influence on the wintertime and summertime ozone changes over the TP. Based on the SLIMCAT 3-D model, we have performed sensitivity experiments to explore the role 150 hPa GH plays in the DJF mean ozone variations over the TP.

Our main conclusions are as follows:

• The comparison of the C3S ozone dataset with the merged SBUV satellite-based observations has verified the feasibility of using assimilated C3S data to study long-term variations over the relatively small TP region.

• With the C3S data extended up to early 2018, the long-term variations of TCO and TOL averaged in different seasons are compared over 1979–2017. The TOL over the TP compared to the zonal mean at the same latitude band exists throughout the year, though the magnitude and the centre location change with season. Both PWLT and EESC-based multivariate regression models show a change in TCO trends from the pre-1997 decline to the post-1997 recovery, although the positive trend based on PWLT is not statistically significant. Compared to the zonal mean trend over the same latitude band, the TP ozone trend shows a relatively smaller rate of increase after 1997, which highlights the zonal asymmetry in ozone recovery.

• Overall, regression results based on three groups of independent explanatory variables show that the GH150 proxy improves the regression especially for the TP region and is more significant than the ST proxy. By comparison of the contributions of different proxies in DJF and JJA, dynamical proxies (QBO and GH150) dominate the wintertime TCO variations over the TP, with statistical significance at 99 % confidence level, but in summertime photochemical processes dominate and the dynamical process decays (QBO at 10 hPa persists but GH150 disappears). The positive correlation between the DJF and JJA TCO over the TP indicates the seasonal persistence of total ozone variations from the ozone buildup in winter to the decreasing period in summer. Our analysis clearly highlights the influence of wintertime GH150 variations on summertime TCO trends.

• Results from the SLIMCAT control experiment (R1) reproduce the TCO time series and regression results for the TP region, and are consistent with the C3S-based results. Sensitivity experiments (R2 and R3) are performed to explore the significant contribution of the GH150 proxy to the wintertime Tibetan ozone variations. The composite analysis shows that GH150 fluctuations play a key role in controlling the DJF mean TCO variability over the TP, which may be associated with ITCZ, ENSO events or the Walker circulation.

Overall, our results show that stratospheric ozone recovery due to the impact of the Montreal Protocol is not expected to behave similarly at all longitudes within a certain latitude region. In the specific case of the Tibetan Plateau, other local factors, which vary with season, will affect column ozone variations. Given the impact of dynamical proxies described above, column ozone over the TP will be subject to long-term changes beyond halogenated ozone-depleting substances and needs careful monitoring.