Global total ozone recovery trends 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 from the trends reported in Weber et al. (2018) using the 2 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 4 (OMPS) satellite instruments (1978–present) as well as the Global Ozone Monitoring Experiment (GOME)-type Total Ozone (GTO-ECV) and GOME-SCIAMACHY-GOME-2 (GSG) merged datasets (both 1995–present), mainly comprising satellite 6 data from GOME, SCIAMACHY, OMI, GOME-2A, -2B, and TROPOMI. The fifth dataset consists of the annual mean zonal mean data from ground-based measurements collected at the World Ozone and UV Radiation Data Center (WOUDC). 8 Trends were determined by applying a multiple linear regression (MLR) to annual mean zonal mean data. The addition of four more years consolidated the fact that total ozone is indeed on slowly recovering in both hemispheres as a result of phasing 10 out ozone depleting substances (ODS) as mandated by the Montreal Protocol. The near global ozone trend of the median of all datasets after 1996 was 0.5±0.2 (2σ) %/decade, which is in absolute numbers roughly a third of the decreasing rate of 12 1.4±0.6 %/decade from 1978 until 1996. 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 1996 which confirms the success of the 14 Montreal Protocol. The observed trends are also in very good agreement with the median of 17 chemistry climate models from CCMI (Chemistry Climate Model Initiative) with current ODS and GHG (greenhouse gas) scenarios. 16 The positive ODS related trends in the NH after 1996 are only obtained with a sufficient number of terms in the MLR accounting properly for dynamical ozone changes (Brewer-Dobson circulation, AO, AAO). A standard MLR (limited to solar, 18 QBO, volcanic, and 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. 20 1 https://doi.org/10.5194/acp-2021-1058 Preprint. Discussion started: 5 January 2022 c © Author(s) 2022. CC BY 4.0 License.

trends (either independent or piecewise linear). The choice of two inflection points were chosen from fits having minimum fit root mean square (rms) errors. As stratospheric halogens are declining steadily since the middle 1990s the interpretation of the 56 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 58 approach by Bloomer et al. (2010), first applied to surface ozone and temperature data at selected stations in the US. Without using any proxy data the regression estimates trends of the seasonality expressed as Fourier series. Attribution of physical and 60 chemical processes to the long-term changes are therefore not possible as also stated by the authors. Latitude and longitude dependent total ozone trends are reported by Coldewey-Egbers et al. (2021) derived from the ESA/DLR 62 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, while in the tropics and NH they are mostly insignificant. Consequently, 64 they only reported significant zonal mean positive trends in the SH.
In Section 2 the updates in the five merged datasets are briefly discussed. In Section 3 the multiple linear regression (MLR) 66 as used in our trend analysis is described. Section 4 presents the total ozone trend results in broad zonal bands: near-global, southern and northern hemispheric extratropics, and tropics. In Section 5 latitude dependent annual mean total ozone trends 68 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 Section 6. In Section 7 a summary and final remarks are given. 70 2 Total ozone datasets Five merged total ozone datasets are used in this study of which one dataset is based upon ground observations. All others are 72 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 74 other two merged datasets are based in large part 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 76 GTO-ECV datasets). These datasets start in 1995.
The ground based dataset is the monthly mean zonal mean data from the network of ground-based Brewers, Dobsons, SAOZ 78 (Système d'Analyse par Observations Zénithales), and filter instruments collected at the World Ozone and UV Data Center (WOUDC) (Fioletov et al., 2002). In addition a brief description of the model data from the CCMI initiative is given. The 80 sources of observational data are listed in Table 1 and brief descriptions of the datasets are given in the following. Annual mean timeseries of all five merged datasets are in very good agreement with each other (see Fig. 2.58 in Weber et al. (2021)).  (Kramarova et al., 2022). Version 8.7 also incorporates an updated a-priori 92 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 94 to adjust each instrument record to an equivalent measurement time of 1:30pm (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 96 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., 98 2013Kramarova et al., 98 , 2022. For merging, data are averaged during periods with multiple operational instruments. The Version 8.7 MOD data contains monthly zonal mean ozone profiles in mixing ratio on pressure levels and in Dobson units on layers. The total ozone 100 is then provided as the sum of the layer data.

102
The NOAA COH (cohesive) dataset is a simple extension in time of the dataset appearing in W18. The data includes v8.6 SBUV on Nimbus 7, v8.6 SBUV/2 from NOAA 9, 11, 16 to 19, and v2r3 OMPS Nadir Profiler (NP) on Suomi-NPP as available from may not always reproduce monthly zonal mean fluctuations well. However, seasonal (and longer) averages can be estimated with a precision comparable with satellite-based data sets (∼1%).

144
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) are used (Eyring et al., 2013). An overview of the models, together 146 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).

148
Here we have used median total column ozone from 17 models taking part in the REF-C2 experiment, an internally consistent seamless simulation from the past into the future between 1960 and 2100.

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

156
Like in our earlier study, the GSG and GTO-ECV timeseries 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  158 period. The NOAA data was here preferred over the NASA data, as the former has shorter data gaps after the major volcanic eruption from Mt Pinatubo in 1991 and subsequent years.
(1) y(t) is the annual mean zonal mean total ozone timeseries and t the year of observations. The coefficients b 1 and b 2 are the linear 166 trends before and after t 0 . In order to make both trends independent of each other (or disjoint), two y-intercepts (a 1 and a 2 ) are added. The multiplication of the independent variable t with X i (t) in the first four terms of Eq. 1 describes mathematically 168 that the first two terms only applies to the period before and the third and fourth terms to the period after the turnaround year.
X 1 (t) and X 2 (t) are given by and 172 respectively. The independent trends before and after t 0 are favored over the use of piecewise linear trends or the use of EESC 174 as a proxy timeseries (see detailed discussions in W18) The maximum of the effective equivalent stratospheric chlorine (EESC) was reached at about the year t 0 = 1996 (Newman et al., 2007) and some years later (t 0 ∼ 2000) in the polar regions (Newman 176 et al., 2006(Newman 176 et al., , 2007. Therefore t 0 was set to 1996 globally, except for the polar regions, where t 0 =2000 was selected. The contributions from the QBO, 11-year solar cycle, and stratospheric aerosols are standard in total ozone MLR analyses (e.g. 178 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 180 a shift by one year) so that no further additional auto-regression term as commonly used for monthly mean ozone timeseries is needed (e.g. Dhomse et al., 2006;Vyushin et al., 2007).

182
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 198 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: 200 In W18 the AAO term was not included. Table 2 summarises the sources of the proxy data used here. The calculation of the 202 BDC proxy from the monthly mean eddy heat fluxes is described in detail in W18. In this study the eddy heat flux data come from the ERA-5 reanalysis (Hersbach et al., 2020).

204
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 206 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 208 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 behavior in the latitude dependent ozone response.

210
The various proxy time series, in particular the atmospheric dynamics related ones, are partially correlated. One way to improve upon this is the possibility to orthogonalize them. Doing so will not change the MLR results, but some contributions 212 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 214 term which makes attribution impossible. The 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.  respectively. The lower value in the SH is due to the influence from the spring Antarctic ozone hole, which exhibits the largest 238 local ozone depletion and leads to mixing of ozone depleted air into the middle latitudes (Atkinson et al., 1989;Millard et al., 2002). Recovery trends are +0.5(0.5) and +0.7(0.6) %/decade in the NH and SH, respectively. Within the trend uncertainty, 240 the 1-to-3 ratio in the linear trends before and after the ODS peak in 1996 are close to the ratio of the rate change in the EESC in both hemispheres.

242
In the tropics the linear trend after 1996 is close to zero and insignificant ( Fig. 2 and Coldewey-Egbers et al., 2021). Table   3 summarises the MLR results in the broad zonal bands from the individual datasets and the median timeseries as well as the 244 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 while in the NH the recovery trends remain unchanged up to the highest latitudes shown. In the tropics recovery trends are close to zero. One notable change from W18 is that the tropical trends during the ODS rising phase are now more negative 264 Table 3. 1979-1996 and 1997-2020 annual mean total ozone trends in various broad zonal bands. Uncertainties are given as 2σ and trends in bold are have an absolute magnitude equal or larger than 2σ. r 2 is the square Pearson correlation between timeseries of observations and MLR and χ the residual defined as χ 2 = i (obsi − modi) 2 /(n − m), where obsi are the observations and modi the MLR, n, the number of data (years) in the timeseries, and m, the number of parameters fitted. All results are obtained using the full MLR. which is most likely caused by larger uncertainties due to the sparsity of ground data at these latitudes. The trend uncertainties are generally larger for the early period before 1996, which in part may be caused by the different lengths of the periods before The dashed pink line shows the expected recovery trends when applying the 1-to-3 ratio (corresponding to the rate change 272 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 recovery trends are 274 slightly positive while the MLR regressions suggest rather near zero trends, but they still agree within their uncertainties.
In order to elucidate further on the interpretation of the independent linear trends after 1996 as recovery trends, we repeated 276 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 recovery 278 trends are nearly unchanged in the SH, the NH recovery trends are reduced to zero in the NH extratropics. On the other hand the tropical trends before 1996 are closer to zero. The expected recovery trends (from the 1-to-3 EESC ratio) have become 280 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 1996 in the NH being close to zero now clearly deviate from the expected 1-to-3 ratio. It appears 282 that the additional atmospheric dynamics terms in the regression balance the positive recovery trends from the full MLR which explains why total ozone in the NH appears more or less stable during the last two decades (panel a of Fig. 2).

284
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 286 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.  The comparisons of trends from the standard and full MLR reveal that the NH 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 290 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.
292 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 timeseries).

294
Well-known factors like solar activity and QBO show the expected behaviour, i.e. more ozone during solar maximum at all latitudes (see e.g. W18) and the opposite sign in the QBO response between inner tropics and extratropics (Bowman, 1989;296 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., In the NH the BDC and and AO mostly contribute to ozone variability. Interestingly, there is an influence from the BDC from 300 one hemisphere to the other in both directions. BDC-N 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 302 region (higher ozone) (e.g. Randel et al., 2002;Weber et al., 2011). The correlation of ozone anomalies in the NH winter/spring to SH total ozone was reported by Fioletov and Shepherd (2003) and is believed to explain the positive response in SH total 304 ozone. Somewhat surprising is the impact of the SH BDC on NH ozone with a negative ozone response, for which we have no explanation.

306
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   (Newman et al., 2007). 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 Earlier signs of ozone recovery has been observed in September above Antarctica (Solomon et al., 2016, W18). Now with four more years of data this recovery of about 11 %/decade remains robust (see panel b of Fig. 7). During September the Antarctic 332 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 334 the polar cap (panel c in Fig. 7). Despite the lack of recovery in the late ozone hole season, 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 336 ozone hole has been shifted to later dates despite the larger than average ozone holes observed in recent years (e.g. 2015 and 2020).

338
Panel a of Fig. 7 shows the March ozone timeseries above the Arctic with ozone recovery trends not statistically different from zero. The trend results in the polar regions as shown in Fig. 7 basically confirm the results from W18 and Langematz et al.

340
(2018). Table 4 summarises the polar trends for the individual datasets and the mean timeseries. Within the trend uncertainties all datasets are in very good agreement.

7 Summary and conclusions
We derived globally total ozone recovery trends from five merged total ozone datasets using a multiple linear regression with 344 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 346 trends may be interpreted as recovery trends related to changes in the stratospheric halogens as a response to the Montreal Protocol and Amendments phasing out ozone depleting substances.

348
For the near-global average we see small recovery 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 350 ratios in stratospheric halogens.
In the tropics recovery trends are not statistically different from zero. In line with earlier observations (Solomon et al., 2016, 352 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 354 of early signs of recovery.
Although we showed that ozone recovery is evident at NH middle latitudes the total ozone levels in the NH extratropics have  Table 4. 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 2 is the square Pearson correlation and χ the residual (see caption of Table 3 to the larger uncertainties of satellite observations in this altitude region. From chemistry climate-models it is expected that the BDC and ozone will increase as a result of greenhouse gases and models so far cannot explain the extratropical decline in Another point which will be important is to show consistency between total ozone trends and both stratospheric and tropo-366 spheric ozone trends. The ozone satellite datasets still show significant differences and opposite signs in trends (Gaudel et al., 2018).