Global total ozone recovery trends attributed to ozone-depleting substance (ODS) changes derived from five merged ozone datasets

We report on updated trends using different merged zonal mean total ozone datasets from satellite and ground-based observations for the period from 1979 to 2020. This work is an update of the trends reported in Weber et al. (2018) using the same datasets up to 2016. Merged datasets used in this study include NASA MOD v8.7 and NOAA Cohesive Data (COH) v8.6, both based on data from the series of Solar Backscatter Ultraviolet (SBUV), SBUV-2, and Ozone Mapping and Profiler Suite (OMPS) satellite instruments (1978–present), as well as the Global Ozone Monitoring Experiment (GOME)-type Total Ozone – Essential Climate Variable (GTO-ECV) and GOME-SCIAMACHY-GOME-2 (GSG) merged datasets (both 1995–present), mainly comprising satellite data from GOME, SCIAMACHY, OMI, GOME-2A, GOME-2B, and TROPOMI. The fifth dataset consists of the annual mean zonal mean data from ground-based measurements collected at the World Ozone and Ultraviolet Radiation Data Centre (WOUDC).

Trends were determined by applying a multiple linear regression (MLR) to annual mean zonal mean data. The addition of 4 more years consolidated the fact that total ozone is indeed slowly recovering in both hemispheres as a result of phasing out ozone-depleting substances (ODSs) as mandated by the Montreal Protocol. The near-global (60^∘ S–60^∘ N) ODS-related ozone trend of the median of all datasets after 1995 was 0.4 ± 0.2 (2σ) %/decade, which is roughly a third of the decreasing rate of 1.5 ± 0.6 %/decade from 1978 until 1995. The ratio of decline and increase is nearly identical to that of the EESC (equivalent effective stratospheric chlorine or stratospheric halogen) change rates before and after 1995, confirming the success of the Montreal Protocol. The observed total ozone time series are also in very good agreement with the median of 17 chemistry climate models from CCMI-1 (Chemistry-Climate Model Initiative Phase 1) with current ODS and GHG (greenhouse gas) scenarios (REF-C2 scenario).

The positive ODS-related trends in the Northern Hemisphere (NH) after 1995 are only obtained with a sufficient number of terms in the MLR accounting properly for dynamical ozone changes (Brewer–Dobson circulation, Arctic Oscillation (AO), and Antarctic Oscillation (AAO)). A standard MLR (limited to solar, Quasi-Biennial Oscillation (QBO), volcanic, and El Niño–Southern Oscillation (ENSO)) leads to zero trends, showing that the small positive ODS-related trends have been balanced by negative trend contributions from atmospheric dynamics, resulting in nearly constant total ozone levels since 2000.

1 Introduction

The stratospheric ozone layer protects the biosphere from harmful UV radiation. How much UV reaches the surface depends, among other factors like clouds, on the overhead total ozone column. The discovery of the Antarctic ozone hole (Chubachi, 1984; Farman et al., 1985; Solomon et al., 1986) raised awareness of the need to protect the ozone layer that culminated in the 1985 Vienna Convention and a commitment to take actions. One of the actions was the signing of the Montreal Protocol in 1987 that started the phaseout of ozone-depleting substances (ODSs), which are sufficiently long-lived to reach the stratosphere and release active halogens that destroy ozone (e.g. Solomon, 1999). As a consequence of the Montreal Protocol and its later amendments, stratospheric halogens started to decline in the middle 1990s (e.g. Anderson et al., 2000; Solomon et al., 2006). A corresponding ozone increase has been detected from satellite and ground-based observations, particularly in the upper stratosphere (Braesicke et al., 2018, and references therein).

Changes in total ozone column are representative of lower stratospheric ozone changes as the majority of ozone resides in the lower stratosphere (“ozone layer”). Lower stratospheric ozone is sufficiently long-lived to be influenced by transport and circulation changes. The rapid increase in Northern Hemisphere total ozone in the late 1990s (Harris et al., 2008) revealed the important role of ozone transport via the Brewer–Dobson (BD) circulation. These circulation changes also cause large variability on inter- and intra-annual timescales in lower stratospheric ozone and the total column (e.g. Fusco and Salby, 1999; Randel et al., 2002; Dhomse et al., 2006; Harris et al., 2008; Weber et al., 2011) and make detection of ozone recovery challenging. Apart from the observed variability, zonal mean total ozone levels in both hemispheres remained stable since about the year 2000 (e.g. Weber et al., 2018). The success of the Montreal Protocol is nevertheless undisputed as the earlier decline in total ozone was successfully stopped (Mäder et al., 2010; Braesicke et al., 2018).

Global and continuous ozone observations from satellites through 2020 now span a total time period of 42 years, of which 25 years cover the period after the stratospheric halogen peak (around 1995). The added years should help in improving the statistical significance of ozone recovery after the middle 1990s (Weatherhead et al., 2000). This paper reports on updated zonal mean total ozone trends from Weber et al. (2018) (abbreviated to W18 in the following) by adding 4 more years of data (2017–2020) to five merged total ozone datasets. In our earlier study, ozone recovery trends in the extratropics were on the order of +0.5 %/decade. The derived trends depend on the proper treatment of dynamical processes in the multilinear regression. Changes in circulation and ozone transport, in part due to increasing greenhouse gas (GHG) levels, have variability on decadal and longer timescales and can therefore mask ODS-related recovery trends. Longer data records are helpful to further disentangle the various processes responsible for long-term changes in ozone. In this work we focus specifically on trend estimates directly related to ODS changes in order to evaluate the direct impact from the Montreal Protocol.

The main results from our earlier paper (W18) were latitude-dependent annual mean total ozone trends from the middle 1990s to 2016, which were reported to be on average +0.5 %/decade in the extratropics and only significant in the Southern Hemisphere (SH; W18). Since W18 was published, there have been three recent studies on global and regional ozone column trends (Bozhkova et al., 2019; Krzyścin and Baranowski, 2019; Coldewey-Egbers et al., 2022). Krzyścin and Baranowski (2019) derived total ozone column trends from a multivariate linear regression (MLR) applied to the Multi-Sensor Reanalysis-2 (MSR-2) total ozone dataset up to 2017 (van der A et al., 2015). In their MLR, they split the entire period from 1978 to 2017 into three periods with separate trends (either independent or piecewise linear). The choice of two inflection points was chosen from fits having minimum fit root mean square (rms) errors. As stratospheric halogens have been declining steadily since the middle 1990s, the interpretation of the segmented trends is difficult. Trends of the first period (before middle 1990s) are in agreement with W18 and this study.

Bozhkova et al. (2019) applied a regression to TOMS and OMI total ozone at northern hemispheric mid-latitudes using the approach by Bloomer et al. (2010), first applied to surface ozone and temperature data at selected stations in the United States. Without using any proxy data, the regression estimates trends of the seasonality expressed as Fourier series. Attribution of physical and chemical processes to the long-term changes is therefore not possible, as also stated by the authors. Latitude- and longitude-dependent total ozone trends are reported by Coldewey-Egbers et al. (2022) derived from the ESA/DLR GTO-ECV dataset, which is one of the five observational datasets used in this study. They report significant positive linear trends after 1995 over large regions in the extratropical Southern Hemisphere (SH), while in the tropics and Northern Hemisphere (NH) they are mostly insignificant. Consequently, they only reported significant zonal mean positive trends in the SH.

In Sect. 2 the updates in the five merged datasets are briefly discussed. In Sect. 3 the multiple linear regression (MLR) as used in our trend analysis is described and discussed. Section 4 presents the total ozone trend results in broad zonal bands: near-global, southern and northern hemispheric extratropics, and tropics. In Sect. 5, latitude-dependent annual mean total ozone trends are presented and discussed. Polar ozone trends for the months where polar ozone losses are largest (e.g. during ozone hole season) are presented in Sect. 6. In Sect. 7, a summary and final remarks are given.

2 Total ozone datasets

Five merged total ozone datasets are used in this study of which one dataset is based upon ground-based observations. All others are based on satellite observations. Two different merged datasets are derived from the series of SBUV and SBUV-2 satellite instruments (SBUV MOD V8.7 from NASA and SBUV COH V8.6 from NOAA), operating continuously since the late 1970s. The other two merged datasets are based largely upon the series of European satellite spectrometers GOME, SCIAMACHY, GOME-2A, and GOME-2B, with different retrieval and merging algorithms applied (University of Bremen GSG and ESA/DLR GTO-ECV datasets). These datasets start in 1995.

The ground-based dataset is the monthly mean zonal mean data from the network of Brewers, Dobsons, SAOZ (Système d'Analyse par Observations Zénithales), and filter instruments collected at the World Ozone and Ultraviolet Radiation Data Centre (WOUDC) (Fioletov et al., 2002). In addition, a brief description of the model data from the CCMI Phase 1 is given. The sources of observational data are listed in Table 1, and brief descriptions of the datasets are given in the following. Annual mean time series of all five merged datasets are in very good agreement with each other to within a few Dobson units (DU; see also Fig. 2.58 in Weber et al., 2021). All datasets cover the entire earth, except for months and latitudes under polar night conditions.

Table 1 Source of merged total ozone datasets.

2.1 NASA SBUV MOD V8.7

The NASA Merged Ozone Data (MOD) time series is constructed using data from the Nimbus 4 BUV, Nimbus 7 SBUV, and six NOAA SBUV-2 instruments numbered 11, 14, and 16–19, and the Ozone Mapping and Profiler Suite Nadir Profiler (OMPS-NP) instrument aboard the Suomi-NPP satellite (Frith et al., 2014, 2022).The instruments are of similar design, and measurements from each are processed using the same V8.7 algorithm. To maintain consistency over the entire time series, the individual instrument records are analysed with respect to each other, and absolute calibration adjustments are applied as needed based on comparison of radiance measurements during periods of instrument overlap (DeLand et al., 2012).

Version 8.7 uses the same core algorithm as Version 8.6 (Bhartia et al., 2013) but includes new inter-instrument calibration adjustments for instrument records since 2000 (NOAA-16 SBUV/2 through OMPS NP) based on a new approach to radiance intercomparisons across overlapping instruments (Kramarova et al., 2022). Version 8.7 also incorporates an updated a priori ozone profile climatology with improved tropospheric representation based on GMI model output and diurnal adjustments to ensure the a priori profile correctly reflects the local solar time of each measurement (Ziemke et al., 2021). A post-retrieval diurnal correction is applied to adjust each instrument record to an equivalent measurement time of 13:30 (Frith et al., 2020). Remaining offsets between instruments exist (mostly below 5 % for layers, below 1 % for total ozone), but their cause is not understood. We therefore do not make adjustments to the data. Rather we set limitations on the data included in the merged product based on data quality analysis by the instrument team and on comparisons with independent measurements (DeLand et al., 2012; Kramarova et al., 2013, 2022). For merging, data are averaged during periods with multiple operational instruments. The Version 8.7 MOD data contain monthly zonal mean ozone profiles in mixing ratio on pressure levels and in Dobson units on layers. The total ozone is then provided as the sum of the layer data.

2.2 NOAA SBUV COH V8.6

The NOAA COH (cohesive) dataset is a simple extension in time of the dataset appearing in W18. The data include v8.6 SBUV on Nimbus 7, v8.6 SBUV/2 from NOAA 9, 11, 16 to 19, and v3r2 OMPS Nadir Profiler (NP) on Suomi-NPP as available from NESDIS STAR. The merging approach differs from NASA MOD in two important ways. NASA MOD averages data from all relevant satellites in any time period for which the data meet certain quality criteria. NOAA COH uses data from a single “best” satellite in any time period. Which satellite is used depends on known data quality issues, on minimizing the solar zenith angle of the measurement, and on maximizing global coverage. NOAA COH does not shift to an equivalent measurement time (13:30) but performs an adjustment between data from differing satellites. For post-2000 data, where drift of the measurement time is minimized, the data are all adjusted to NOAA 18. For data from 1999 and prior to this year, the inter-satellite overlap is often short, and the satellite drift often significant; we choose only to adjust NOAA 9 to the two branches NOAA 11 prior to and after the NOAA 9 time period. The total ozone is calculated from the sum of the adjusted profile layer data. By vertical integration, many of the layer adjustments to a large extent cancel out such that the final total ozone product is altered by less than 1 %, and in most cases by less than 0.5 %, from the original satellite datasets.

2.3 University of Bremen GSG

The merged GOME, SCIAMACHY, GOME-2A, and GOME-2B (GSG) total ozone time series (Kiesewetter et al., 2010; Weber et al., 2011, 2018) consists of total ozone data that were retrieved using the University of Bremen weighting function DOAS (WFDOAS) algorithm (Coldewey-Egbers et al., 2005; Weber et al., 2005; Orfanoz-Cheuquelaf et al., 2021). The merging of the data has been described in W18. The most recent modification was to replace GOME-2A data after January 2015 with data from GOME-2B (2012–present), which has a better global coverage after changes in the GOME-2A scanning pattern. Latitude-dependent bias corrections for GOME-2B were applied from the overlapping period 2014–2020 with GOME-2A.

2.4 DLR/ESA GTO-ECV

The latest version of the GOME-type Total Ozone Essential Climate Variable (GTO-ECV) data record (Coldewey-Egbers et al., 2015, 2022; Garane et al., 2018) has been generated as part of the European Space Agency’s Climate Change Initiative+ ozone (ESA_CCI+ ozone) project. Total columns from six sensors (GOME, SCIAMACHY, OMI, GOME-2A, GOME-2B, and TROPOMI), retrieved with the GOME Direct Fitting (GODFIT) version 4 algorithm (Lerot et al., 2014; Garane et al., 2018), were combined into a coherent record that covers the period 1995–2020. OMI was used as a reference instrument, and the other sensors were adjusted by means of latitude- and time-dependent correction factors determined from overlap periods.

2.5 WOUDC data

The WOUDC zonal mean dataset (Fioletov et al., 2002) was formed from ground-based measurements by Dobson, Brewer, SAOZ instruments, and filter ozonometers available from the WOUDC. The overall performance of the ground-based network was discussed by Fioletov et al. (2008), and the present state of the network is described by Garane et al. (2019). This dataset is provided by the WOUDC and updated regularly.

First, ground-based measurements were compared with an ozone “climatology” (monthly means for each point of the globe) estimated from satellite data for 1978–1989. Then, for each station and for each month, the deviations from the climatology were calculated, and a zonal mean value for a particular month was estimated as a mean of these deviations. The calculations were done for 5^∘ wide zonal bands. In order to take into account various densities of the network across regions, the deviations of the stations were first averaged over 5 by 30^∘ cells, and then the zonal mean was calculated by averaging these first set of averages over the 5^∘ wide zonal band. Then the zonal averages were smoothed by approximating them using Legendre polynomials.

The WOUDC dataset was compared with merged satellite time series and demonstrated a good agreement (Chiou et al., 2014). Estimates based on relatively sparse ground-based measurements, particularly in the tropics and Southern Hemisphere, may not always reproduce monthly zonal mean fluctuations well. However, seasonal (and longer) averages can be estimated with a precision comparable with satellite-based datasets (∼1 %) (Chiou et al., 2014).

2.6 Chemistry–climate model data

In this study, output from the chemistry–climate models (CCMs) and chemistry-transport models (CTMs) participating in phase 1 of CCMI (Chemistry-Climate Model Initiative) is used (Eyring et al., 2013). An overview of the models, together with details particular to each model and an overview of the available simulations, is given in Morgenstern et al. (2017), along with a detailed description of the full forcings used in the reference simulations (Eyring et al., 2013; Hegglin et al., 2016). Here we have used median total column ozone from 17 models taking part in the REF-C2 experiment, an internally consistent seamless simulation between 1960 and 2100.

2.7 Data preparation

From the zonal mean monthly mean data in 5^∘ latitude steps (all datasets), annual means were calculated. Wider zonal bands (like 35–60^∘ N) were averaged from the 5^∘ data using area weights (see W18). All annual mean zonal mean time series were bias-corrected by subtracting the difference to the mean of all datasets during the 1998–2008 period. The multi-dataset mean was then added back to each dataset, such that all bias-corrected time series are provided in units of the total column amounts (W18). However, the trend results derived from them are identical to those derived using anomaly time series.

Like in our earlier study, the GSG and GTO-ECV time series were extended from 1995 back to 1979 using the bias-corrected NOAA data. This way one ensures that all terms other than the trend terms are determined from the full time (1979–2020) period. The NOAA data were preferred here over the NASA data, as the former have shorter data gaps after the major volcanic eruption from Mt. Pinatubo in 1991 and subsequent years.

3 Multiple linear regression

The standard MLR model is identical to the one used in W18 and includes two independent linear trend terms (before and after the ODS-related turnaround year t₀=1995), two aerosol terms (Mt. Pinatubo 1992 and El Chichón 1983), a solar cycle term, two Quasi-Biennial Oscillation (QBO) terms (50 and 10 hPa), and ENSO (El Niño–Southern Oscillation):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& y\left(t\right)=[a_{\mathrm{1}}+b_{\mathrm{1}}\cdot (t_{\mathrm{0}}-t\left)\right]X_{\mathrm{1}}\left(t\right)+[a_{\mathrm{2}}+b_{\mathrm{2}}\cdot (t-t_{\mathrm{0}}\left)\right]X_{\mathrm{2}}\left(t\right)\\ & +{\mathit{\alpha }}_{\text{sun}}\cdot S\left(t\right)+{\mathit{\alpha }}_{\text{qbo50}}\cdot Q_{\mathrm{50}}\left(t\right)+{\mathit{\alpha }}_{\text{qbo10}}\cdot Q_{\mathrm{10}}\left(t\right)+{\mathit{\alpha }}_{\text{ENSO}}\cdot E\left(t\right)\\ & +{\mathit{\alpha }}_{\text{El Chichón}}\cdot A_{\mathrm{1}}\left(t\right)+{\mathit{\alpha }}_{\text{Pinatubo}}\cdot A_{\mathrm{2}}\left(t\right)+P\left(t\right)+\mathit{ϵ}\left(t\right).\end{array}\end{array}$$

y(t) is the annual mean zonal mean total ozone time series and t the year of observations. The coefficients b₁ and b₂ are the linear trends before and after t₀. In order to make both trends independent of each other (or disjoint), two y intercepts (a₁ and a₂) are added. The multiplication of the independent variable t with X_i(t) in the first four terms of Eq. (1) describes mathematically that the first two terms only apply to the period before and the third and fourth terms to the period after the turnaround year. X₁(t) and X₂(t) are given by

$$\begin{array}{}\text{(2)}& X_{\mathrm{1}}\left(t\right)=\left\{\begin{array}{ll}\mathrm{1}& \text{if }t\le t_{\mathrm{0}}\\ \mathrm{0}& \text{if }t>t_{\mathrm{0}}\end{array}\right.\end{array}$$

and

$$\begin{array}{}\text{(3)}& X_{\mathrm{2}}\left(t\right)=\left\{\begin{array}{ll}\mathrm{0}& \text{if}t\le t_{\mathrm{0}}\\ \mathrm{1}& \text{if}t>t_{\mathrm{0}}\end{array}\right.,\end{array}$$

respectively. The independent trends before and after t₀ are favoured over the use of piecewise linear trends or the use of EESC as a proxy time series (see detailed discussions in W18). The maximum of the effective equivalent stratospheric chlorine (EESC) was reached at about the year t₀=1995 (Newman et al., 2007), except for the polar region ($>\mathrm{60}{}^{\circ }$) where t₀=2000 (Newman et al., 2006, 2007).The contributions from the QBO, 11-year solar cycle, and stratospheric aerosols are standard in total ozone MLR analyses (e.g. Staehelin et al., 2001; Reinsel et al., 2005). ϵ(t) is the residual from fitting the coefficients to match the regression model (right side) to the observations. By using annual mean total ozone, auto-correlation is very low here (below 0.1 in absolute value for a shift by 1 year) so that no further additional auto-regression term as commonly used for monthly mean ozone time series is needed (e.g. Dhomse et al., 2006; Vyushin et al., 2007).

The stratospheric aerosols are dominated by the major volcanic eruptions from El Chichón (1982) and Mt. Pinatubo (1991). Enhanced aerosols in the lower stratosphere lasting for a few years impact both ozone chemistry and transport (Schnadt Poberaj et al., 2011; Dhomse et al., 2015). The stratospheric aerosol optical depth (SAOD) at 550  nm from Sato et al. (1993) is used as the explanatory variable before 1990 (includes the El Chichón event), while newer data from the WACCM model (Mills et al., 2016) are used for the period after 1990 (include Mt. Pinatubo major volcanic eruption and the series of more minor volcanic eruptions from the last decade). Missing years after 2015 were filled with background values from the late 1990s.

As mentioned in W18, there is not a sufficient number of months and/or 5^∘ latitude bands available in the SBUV data records for some years, and for these years annual means were treated as missing data. Annual means were only used in the regression if at least 80 % of the 5^∘ bands of the data were contained in the broad zonal bands and 80 % of months available in that year. If annual means of the years 1982 and 1983 are missing, the “El Chichon” term is not used in the MLR; similarly if missing all years from 1991 to 1994, the “Pinatubo” term is excluded in the MLR.

The MLR equation, Eq. (1), without the P(t) term has been commonly applied for determining trends from ozone profile data (e.g. Bourassa et al., 2014, 2018; Harris et al., 2015; Tummon et al., 2015; Sofieva et al., 2017; Steinbrecht et al., 2017). The extra term P(t) in Eq. (1) accounts for additional factors of dynamical variability that have been used in different combinations and definitions (e.g. accumulated, time-lagged) in the past. It includes contributions from the Arctic (AO) and Antarctic Oscillation (AAO), and the Brewer–Dobson circulation (BDC) (e.g. Reinsel et al., 2005; Mäder et al., 2007; Chehade et al., 2014; Weber et al., 2018). The BDC terms are usually described by the eddy heat flux at 100 hPa that is considered a main driver of the BDC (Fusco and Salby, 1999; Randel et al., 2002; Weber et al., 2011). The term P(t) is given as follows:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}P\left(t\right)=& {\mathit{\alpha }}_{\text{AO}}\cdot \text{AO}\left(t\right)+{\mathit{\alpha }}_{\text{AAO }}\cdot \text{AAO}\left(t\right)\\ +& {\mathit{\alpha }}_{\text{BDCn}}\cdot \text{BDCn}\left(t\right)+{\mathit{\alpha }}_{\text{BDCs}}\cdot \text{BDCs}\left(t\right).\end{array}\end{array}$$

In W18 the AAO term was not included. Table 2 summarizes the sources of the proxy data used here. The BDCn and BDCs are 100 hPa eddy fluxes in the Northern (n) and Southern Hemisphere (s). The calculation of the BDC proxy from the monthly mean eddy heat fluxes is described in detail in W18. In this study, the eddy heat flux data were derived from the ERA-5 reanalysis (Hersbach et al., 2020).

(Snow et al., 2014)(update from Naujokat, 1986)(Wolter and Timlin, 2011)(update from Sato et al., 1993)(Mills et al., 2016)

Table 2 Sources of explanatory variables/proxy time series used in the MLR.

One may argue that the addition of P(t) will lead to some overfitting by the MLR. We justify this addition as it enables us to obtain MLR fits matching the extreme events like very high annual mean ozone in the NH in 2010 and the very large warming events above Antarctica in 2002 and 2019 with unusually high ozone. The better the dynamical variations are represented in the MLR, the more likely we can separate out dynamical trend contributions and the linear trend terms best approximate EESC-related trends. In our previous study, only selected terms from P(t) were used dependent on their significance in specific zonal bands. Retaining all terms in all MLRs leads to smoother behaviour in the latitude-dependent ozone response.

The various proxy time series, in particular the atmospheric-dynamics-related ones, are partially correlated. One way to improve upon this is to orthogonalize them. Doing so will not change the MLR fit results, but some contributions from the original proxy terms will be redistributed among the proxies that were orthogonalized. It is also common to detrend the proxy time series. In that case, all linear changes of the various processes or proxies will be added up in the linear trend term which makes attribution impossible. For these reasons, we do not detrend nor orthogonalize the proxy time series in this study. Our goal here is that linear changes of all the processes as expressed by the various proxy terms shall be excluded from the linear trend terms such that the linear trends can be attributed as close as possible to ODS changes.

Figure 1 Near-global (60^∘ S–60^∘ N) total ozone time series of five bias-corrected merged datasets. The thick orange line is the result from applying the full MLR (Eqs. 1 and 4) to the median time series. The square of the correlation between observations and MLR is given by r². χ² is the sum square of the differences between median observational and MLR time series divided by the degrees of freedom (difference between the number of years, n, and number of parameters used in the MLR, m). The solid lines indicate the linear trends before and after the ODS peak, respectively. The dotted lines indicate the 2σ uncertainty of the MLR trend estimates. Trend numbers are indicated for the pre- and post-ODS peak period in the top part of the plot. Numbers in parentheses are the 2σ trend uncertainty. The dashed orange line shows the mean ozone level from 1964 until 1980 from the WOUDC data. The thick grey line is the median of 17 chemistry–climate models from the CCMI.

4 Total ozone trends in broad zonal bands

Figure 2 Same as Fig. 1 but for broad zonal bands: (a) 35–60^∘ N (Northern Hemisphere), (b) 20^∘ S–20^∘ N (tropics), and (c) 35–60^∘ S (Southern Hemisphere).

Figure 1 shows the near-global mean time series (60^∘ S–60^∘ N) of the five bias-corrected merged datasets. The thick orange line is the MLR time series from applying the full regression model (Eqs. 1 and 4) to the median of the five time series. A total of 94 % of the variability in total ozone is well captured by the full MLR. A positive trend of $+\mathrm{0.4}\pm \mathrm{0.2}$(2σ) %/decade after 1995 is derived. This trend is about one-third of the absolute trend during the phase of increasing ODSs before 1995, which is $-\mathrm{1.5}\pm \mathrm{0.6}$ %/decade. The ratio of trends before and after 1995 is very close to the ratio of rate changes in the effective equivalent stratospheric chlorine (EESC) before and after the middle 1990s (Dhomse et al., 2006; Newman et al., 2007). Therefore, the observed linear trend of roughly half a percent per decade up to 2020 can be attributed to reductions in ODSs following the Montreal Protocol. This ODS-related recovery appears statistically robust (to within 2σ), even though the ozone levels have stayed more or less constant apart from the year-to-year variability since the year 2000. The magnitude of the post-ODS peak trend remained unchanged from W18. The trend results vary only slightly if the turnaround year (1995) of the ODS change is shifted by 1 year back and forward. Even if the MLR fit of the post-ODS peak period is limited to years after 2000, the ODS-related trend remains robust at +0.5(0.3) %/decade.

The current near-global ozone level (2017–2020) is about 2.3 % below the average from the 1964–1980 time period, the latter derived from the WOUDC data (see Fig. 1). Recovery of total ozone to the 1980 level is generally not expected before about the middle of this century (Braesicke et al., 2018). The near-global total ozone time series from the median of the 17 CCMI chemistry–climate models is in very good agreement with the observations from which we conclude that the chemical and dynamical changes in total ozone under current ODS and greenhouse gas (GHG) scenarios are well understood and consistent with observations.

Figure 2 shows the ozone time series in the Northern (NH) and Southern Hemisphere (SH) as well as in the tropics. Again, the current ozone levels are well below the 1964–1980 mean, specifically −3.6 % and −4.7 % in the NH and SH (35–60^∘ latitudes), respectively. The lower value in the SH is due to the influence from the spring Antarctic ozone hole, which exhibits the largest local ozone depletion and leads to mixing of ozone-depleted air in the middle latitudes (Atkinson et al., 1989; Millard et al., 2002). ODS-related trends are +0.5(0.5) and +0.7(0.6) %/decade in the NH and SH, respectively. Within the trend uncertainty, the 1 : 3 ratio in the linear trends before and after the ODS peak in 1995 is close to the ratio of the rate change in the EESC in both hemispheres.

In the tropics the linear trend after 1995 is close to zero and insignificant (Fig. 2 and Coldewey-Egbers et al., 2022). Table 3 summarizes the MLR results in the broad zonal bands from the individual datasets and the median time series as well as the mean and median of the individual trends.

In most cases the results from the individual datasets are highly consistent, in particular for the near-global time series. All datasets indicate significant near-global ODS-related trends of around half a percent per decade. The trend derived from the NASA data is a bit lower at +0.2 %/decade. The median and mean trends of all datasets agree here with the trends of the median time series as shown in Figs. 1 and 2. For the narrower zonal bands, not all datasets show significant trends after 1995. The NASA and GSG datasets show lower trends in the NH (+0.3 and +0.1 %/decade, respectively), while all others are between +0.5 and +0.7 %/decade and significant. In the SH, all trends agree to within 0.001 %/decade (+0.7 %/decade), except for the NOAA dataset showing a somewhat higher trend of +1 %/decade.

Table 3 1979–1995 and 1996–2020 annual mean total ozone trends in various broad zonal bands. Uncertainties are given as 2σ; trends in bold have an absolute magnitude equal to or larger than 2σ. r² is the square Pearson correlation between time series of observations and MLR and χ the residual defined as ${\mathit{\chi }}^{\mathrm{2}}={\sum }_i({\text{obs}}_i-{\text{mod}}_i{)}^{\mathrm{2}}/(n-m)$, where obs_i denotes the observations, mod_i the MLR model, n the number of data (years) in the time series, and m the number of parameters fitted. All results are obtained using the full MLR. The second column shows the percent difference of the current total ozone level (2017–2020) with respect to the period before 1980 derived from the median time series.

Bold numbers: statistical significance at 2σ.

In the tropics, trends are close to zero, with the exception of the GSG and GTO datasets that have very small and barely significant positive trends of +0.3 ± 0.3 %/decade. The variations in the trend results from the different datasets are most likely due to some residual drifts in the datasets that are not accounted for in the data merging. With the use of the full MLR with all terms and with 4 years added in the time series, the ozone trends in the various zonal bands after 1995 remain quite similar to the results reported in W18, but uncertainties are slightly reduced.

Table 4 shows different MLR settings applied to the median total ozone time series in broad zonal bands (as defined in Table 3). Here the results from the standard and full MLR are listed. In addition, we applied an iterative MLR approach where statistically insignificant terms (2σ criterion) from Eq. (4) and the El Niño term are successively excluded before the final MLR run. The inclusion of the dynamical proxies generally improved the MLR fit (r² and χ values). Except for the NH zonal band (35–60^∘ N), the various MLR settings yield nearly the same post-ODS peak trends for all broad zonal bands (Table 4). There are, however, larger changes in the trends before the middle 1990s. In the extratropics the early-period trends are lower in the standard retrieval (−4.0 % versus −1.9 %/decade in the NH and −3.1 % versus −1.9 %/decade in the SH). This means that atmospheric dynamics and transport changes contributed to lower early-period extratropical total ozone trends in the standard regression (due to the lack of these dynamical terms in the MLR). The opposite is the case in the tropics, where the early-period trends in the standard MLR are slightly higher than in the full MLR. This opposite behaviour is consistent with ozone transport patterns due to the Brewer–Dobson circulation.

It appears that the post-ODS trends are in most cases unchanged, regardless of the number of extra terms used in the MLR. The linear trend term is the only low-frequency term in the MLR equations, while the dynamical proxies have some high-frequency contributions. This makes the trend estimates rather robust and less sensitive to the various other terms used in the MLR. The only significant changes in the post-ODS peak trends are seen in the NH extratropics. In the standard MLR this trend is zero, while the full and iterative MLR shows trends of half a percent per decade. The sum of the ODS-related trend (full MLR) and atmospheric dynamics contribution (difference in the trends between full and standard MLR) cancel out to result in a zero trend in the standard MLR. The negative dynamical trend contribution in the NH is further discussed later in the paper. The correlation between regression and observations is substantially lower in the standard retrieval (r²=0.74 versus 0.88), which indicates that the standard MLR seems not to capture all variability and changes of total ozone.

In order to document the changes from the current MLR fits (Table 4) to results from the period up to and including 2016 (as in W18), the different MLR settings were applied to the current data for the shorter period as summarized in Table S1 (see Supplement). Note that the results in Table S1 may differ from W18 as the merged datasets have been updated, and data before 2017 may have changed as well. Results from the shorter time period are nearly identical to those shown in Table 3. There is one notable change. The uncertainties of the NH trends from the full MLR up to 2020 are reduced such that these trends have become barely significant (2σ). The NH post-ODS peak trend of the standard MLR is slightly positive up to 2016 but statistically insignificant and within the uncertainties not different from the current results.

Table 4 Different MLR settings applied to the broad zonal mean median ozone time series. For explanations of terms, see Table 3. Standard MLR is based upon Eq. (1) with P(t)=0. Iterative MLR means that only terms of P(t) and El Niño term are fitted when they are statistically significant (2σ), while full MLR includes all terms in Eq. (1) and P(t).

bold numbers: statistical significance at 2σ.

5 Latitude-dependent total ozone trends

Latitude-dependent trends in steps of 5^∘ are shown from 60^∘ S to 60^∘ N for all five merged datasets (thin lines) in Fig. 3. The two thick blue and red lines are the results before and after 1995 from applying the full MLR to the median time series, including 2σ uncertainties shown as error bars. In the extratropics, the ODS-related trends are on the order of +0.5 %/decade, with 2σ uncertainties of about the same magnitude. In the SH, the trends continuously increase to nearly +1.3 %/decade in the 55–60^∘ S band, while in the NH, the trends remain unchanged up to the highest latitudes shown. In the tropics, trends are close to zero. One notable change from W18 is that the tropical trends during the ODS rising phase are now more negative (down to −1 %/decade), while before they were mainly close to zero. This may be caused by the additional proxy terms used in this study.

After 1995, all trends of all datasets are in good agreement to within ±0.3 %/decade. There are some notable differences in the northern subtropical and northern tropical trends for the WOUDC data (up to +1 %/decade) compared to the other datasets, which is most likely caused by larger uncertainties due to the sparsity of ground-based data at these latitudes. The trend uncertainties are generally larger for the early period before 1996, which is caused by the different lengths of the periods before (17 years) and after 1995 (25 years).

The dashed pink line shows the expected ODS-related trends when applying the 1 : 3 ratio (corresponding to the rate change of the EESC) to the trends before 1996. It agrees quite well in the extratropics with the independent linear trend estimates and therefore give us confidence that ozone is responding to the long-term ODS decline. The expected tropical ODS-related trends are slightly positive, while the MLR regressions rather suggest near-zero trends, but they still agree within their uncertainties. In the NH extratropics, the expected ODS-related trend is slightly higher than the observed trends but also agrees within the uncertainties of the observed trends.

In order to elucidate the interpretation of the independent linear trends after 1995, we repeated the analyses using the standard MLR, which excludes several terms responsible for changes in atmospheric dynamics and transport (Eq. 1 with P(t)=0). The latitude-dependent trends from the standard MLR are shown in Fig. 4. While the observed trends for both MLRs are nearly unchanged in the SH, the NH trends are reduced to zero in the NH extratropics. On the other hand the tropical trends before 1996 are closer to zero. The expected ODS-related trends (from the 1 : 3 EESC ratio) have become larger, with increases to +1.5 %/decade at the higher latitudes now in both hemispheres. The most obvious result is that the independent linear trends after 1995 in the NH being close to zero now clearly deviate from the expected 1 : 3 ratio. It appears that the additional atmospheric dynamics terms in the regression balance the positive trends from the full MLR, which explains why total ozone in the NH has appeared more or less stable during the last 2 decades (Fig. 2a).

The declining trends in the NH before 1996 (Fig. 4a) are stronger in the standard MLR and are comparable to the SH (about −4 %/decade near 60^∘ latitude). On the other hand, ODS-related trends are expected to be somewhat stronger in the SH as the influence from polar ozone losses on mid-latitude ozone is thought to be larger in the SH, since Arctic ozone losses are more sporadic and generally smaller. In that regard, the trends from the full MLR seem to support this notion.

Figure 3 Latitude-dependent ozone trends in steps of 5^∘ from applying the full MLR (Eqs. 1 and 4) to the median time series of the five merged total ozone datasets. Trends and 2σ standard deviations are shown in blue for the time period before 1995 and in red after 1995. The thin lines show the trends of the individual total ozone datasets. The dashed pink lines are the post-ODS peak trends as expected from the 1 : 3 ratio (corresponding to changes in the stratospheric halogen), applied to the median time series' trends before 1996.

Figure 4 Same as Fig. 3 but from applying the standard MLR (Eq. 1 and P(t)=0)).

The comparisons of trends from the standard and full MLR reveal that the NH ODS-related ozone recovery is balanced by long-term changes in atmospheric dynamics (circulation and transport changes), or in other words the near-zero linear post-ODS peak trends are caused by the combination of ODS-related recovery and dynamical changes. These two signals are more clearly separated in the full MLR. Before discussing this further, we will take a look at the contributing factors or terms in the MLR. Figure 5 shows the maximum response of the various terms in Eqs. (1) and (4) as a function of latitude (from the fit to the median time series).

Figure 5 Signed maximum time series contribution from various terms in the full MLR equation applied to the median of the five merged total ozone datasets. Solid line indicates values of fit coefficients that are larger than their 2σ uncertainty. The sign of the BDCs proxy time series is reversed, so that in both hemispheres positive BDC term values correspond to enhanced Brewer–Dobson circulation. Negative values correspond to an anti-correlation of the ozone response to the proxy. For instance, positive solar contributions mean high solar activity, which leads to more ozone.

Well-known factors like solar activity and QBO show the expected behaviour, i.e. more ozone during solar maximum at all latitudes (see, for example, W18) and the opposite sign in the QBO response between inner tropics and extratropics (Bowman, 1989; Baldwin et al., 2001). The solar response is of similar magnitude at all latitudes, which means that the solar effect in the lower stratosphere is mostly indirect via changes in temperature and associated atmospheric circulation changes (e.g. Dhomse et al., 2022).

Figure 6 (a) March NH total ozone (40–90^∘ N) from the five bias-corrected merged datasets (coloured) and the smoothed median time series (thick black line). (b) DJFM Arctic Oscillation (AO) index. Black lines shows the three-point triangular smoothed time series. Note the inverted y scale. (c) 100 hPa winter eddy heat flux September to March average (BDCn proxy), with the black line showing the three-point triangular smoothed time series. (d) Inverted stratospheric halogen time series in parts per billion (ppb), representative of middle latitudes (Newman et al., 2007).

Figure 7 Polar ozone trends derived from the full MLR applied to the median time series. (a) March 60–90^∘ N, (b) September, and (c) October 60–90^∘ S. See Fig. 1 for more details.

In the NH, the BDC and AO mostly contribute to ozone variability. Interestingly, there is an influence from the BDC from one hemisphere to the other in both directions. BDCn results in opposite responses in the tropics and NH extratropics. This is expected from the planetary waves driving the BDC, leading to ascent in the tropics (lower ozone) and descent in the polar region (higher ozone) (e.g. Randel et al., 2002; Weber et al., 2011). The correlation of ozone anomalies in the NH winter and spring to SH total ozone was reported by Fioletov and Shepherd (2003) and is believed to explain the positive response in SH total ozone. Somewhat surprising is the impact of the SH BDC on NH ozone with a negative ozone response, for which we have no explanation.

Table 5 Polar total ozone trends in March (NH), September (SH), and October (SH) before and after 2000. Uncertainties are provided for 2σ, and trends in bold indicate statistical significance. r² is the square Pearson correlation and χ the residual (see caption of Table 3). The results were obtained from the full MLR.

Bold numbers: statistical significance at 2σ.

The major volcanic eruption of Mt. Pinatubo in 1991 had a stronger impact on the NH, reducing ozone for several years after the event, while ozone advection apparently balanced the surface acid particle (aerosol)-related ozone losses in the SH (Schnadt Poberaj et al., 2011; Aquila et al., 2013; Dhomse et al., 2015). The second large major volcanic eruption from spring 1982 led to aerosol-related ozone loss in the tropics and NH, while surprisingly a positive ozone response in the SH is seen, possibly related to some atmospheric circulation changes compensating for chemical effects from the El Chichon eruption. In contrast to Mt. Pinatubo, which spread sulfuric acid particles into both hemispheres, enhanced aerosols from El Chichon were confined to lower latitudes in the NH (McCormick and Swissler, 1983), consistent with the region of negative ozone response shown in Fig. 5.

The main reason for stable ozone levels observed at NH mid-latitudes since 2000 was identified to stem from the balancing of the positive observed ODS-related trend by negative trends due to circulation changes and ozone transport (see Figs. 2 and 3). The change in the BDCn proxy and AO over the last 55 years is shown in Fig. 6, along with March total ozone northward of 40^∘ N. The variability in the extratropical annual mean is usually dominated by the variability in winter and spring, where BDC maximizes in the seasonal cycle. Apart from the strong drop in ozone in the 1990s related to the major volcanic eruption and associated circulation changes, NH total ozone has been steadily declining over the last 55 years (about 25 DU). This decline is coherent with an overall positive shift of the AO index. A weakening of the BDC is also seen but appears less clear than for the AO.

A positive shift in the AO and a weakening of the BDC result in a strengthening of the polar vortex, which is associated with larger polar ozone losses (Lawrence et al., 2020). Hu et al. (2018) linked a recent strengthening of the stratospheric Arctic vortex in part to a warming of sea surface temperatures in the central northern Pacific. A downward trend in extratropical lower stratospheric ozone has been reported by Ball et al. (2018) that could be consistent with the total ozone observations. Other studies with many different ozone profile datasets did not show significant trends in the lower stratosphere due to very large variability and lower accuracy of the satellite data in this altitude region (Sofieva et al., 2017; Steinbrecht et al., 2017; Arosio et al., 2019).

6 Trends in polar spring

Earlier signs of ozone recovery were observed in September above Antarctica (Solomon et al., 2016, W18). Now, with 4 more years of data, this recovery of about +12 %/decade remains robust (see Fig. 7b). During September, the Antarctic ozone hole size usually increases and reaches its maximum in late September and early October. In a typical Antarctic winter, ozone is completely destroyed in the lower stratosphere, which may explain why no recovery is yet observed in October over the polar cap (Fig. 7c). Several diagnostics clearly indicate a healing of Antarctic ozone as a consequence of the Montreal Protocol. Stone et al. (2021) show that the onset of the Antarctic ozone hole has been shifted to later dates despite the larger than average ozone holes observed in recent years (e.g. 2015 and 2020).

Figure 7a shows the March ozone time series above the Arctic, with ODS-related ozone trends not statistically different from zero. The trend results in the polar regions basically confirm the results from W18 and Langematz et al. (2018). Table 5 summarizes the polar trends for the individual datasets and the mean time series. Within the trend uncertainties, all datasets are in very good agreement.

7 Summary and conclusions

We derived globally total ozone trends from five merged total ozone datasets using a multiple linear regression with independent linear trend (ILT) terms before and after the turnaround in stratospheric halogens in the middle 1990s (∼2000 in the polar regions). When properly accounting for dynamical changes via atmospheric circulation and transport, these retrieved trends may be directly related to changes in the stratospheric halogens (and ODSs) as a response to the Montreal Protocol and its amendments phasing out ozone-depleting substances.

For the near-global average we see small ODS-related trends of about +0.5 %/decade, with main contributions from the extratropics in both hemispheres. The ratio of ozone trends after and before the turnaround year is in very good agreement with the trend ratios in stratospheric halogens or ODSs.

In the tropics, trends are not statistically different from zero. In line with earlier observations (Solomon et al., 2016, W18), polar ozone recovery has been only identified in September above Antarctica, which is connected to the observed delay in the onset of the Antarctic ozone hole (Stone et al., 2021). In the Arctic, large interannual variability still prevents the detection of early signs of recovery.

Although we showed that ODS-related ozone recovery is evident at NH middle latitudes, the total ozone levels in the NH extratropics have been more or less stable since about 2000. Our regression results show that the positive ODS-related trend here is balanced by changes in ozone transport. A long-term positive drift in the AO index over the last 55 years is indicative of a strengthening of the Arctic vortex (Hu et al., 2018; Lawrence et al., 2020; von der Gathen et al., 2021) and reduced winter and spring transport of ozone into middle and high latitudes. This result may be consistent with the observed decline in lower stratospheric ozone in the extratropics as reported by Ball et al. (2018) and Wargan et al. (2018). They mainly attribute this decline to enhanced horizontal mixing with the tropical region, where lowermost stratospheric ozone decreases (Thompson et al., 2021). Other studies and datasets, however, have so far not confirmed the long-term extratropical decline in the lower stratosphere (Arosio et al., 2019; Steinbrecht et al., 2017; Sofieva et al., 2017), which may be in part due to the larger uncertainties of satellite observations in this altitude region. From chemistry–climate models, it is expected that with a strengthening of the BDC due to increasing GHGs, tropical ozone declines and extratropical ozone increases in the lower stratosphere (Dietmüller et al., 2021).