Ozone Variability and Trend Estimates from 20-Years of Ground-Based and Satellite Observations at Irene Station, South Africa

Abstract

While the stratospheric ozone protects the biosphere against ultraviolet (UV) radiation, tropospheric ozone acts like a greenhouse gas and an indicator of anthropogenic pollution. In this paper, we combined ground-based and satellite ozone observations over Irene site (25.90° S, 28.22° E), one of the most ancient ozone-observing stations in the southern tropics. The dataset is made of daily total columns and weekly profiles of ozone collected over 20 years, from 1998 to 2017. In order to fill in some missing data and split the total column of ozone into a tropospheric and a stratospheric column, we used satellite observations from TOMS (Total Ozone Mapping Spectrometer), OMI (Ozone Monitoring Instrument), and MLS (Microwave Limb Sounder) experiments. The tropospheric column is derived by integrating ozone profiles from an ozonesonde experiment, while the stratospheric column is obtained by subtracting the tropospheric column from the total column (recorded by the Dobson spectrometer), and by assuming that the mesospheric contribution is negligible. Each of the obtained ozone time series was then analyzed by applying the method of wavelet transform, which permitted the determination of the main forcings that contribute to each ozone time series. We then applied the multivariate Trend-Run model and the Mann–Kendall test for trend analysis. Despite the different analytical approaches, the obtained results are broadly similar and consistent. They showed a decrease in the stratospheric column (−0.56% and −1.7% per decade, respectively, for Trend-Run and Mann–Kendall) and an increase in the tropospheric column (+2.37% and +3.6%, per decade, respectively, for Trend-Run and Mann–Kendall). Moreover, the results presented here indicated that the slowing down of the total ozone decline is somewhat due to the contribution of the tropospheric ozone concentration.

1. Introduction

Atmospheric ozone plays a key role by protecting the biosphere against harmful ultraviolet (UV) radiation. About 90% of ozone mass is found in the stratosphere, forming the ozone layer where damaging UV radiation from the sun is filtering out. While stratospheric ozone is important for climate system protection, tropospheric ozone is considered as greenhouse gas contributing directly to global warming [1] and presents a serious environmental concern due to its negative impact on human health and terrestrial ecosystems. Mainly formed in the tropical regions, stratospheric ozone is transported to the subtropics and higher latitude regions following the Brewer–Dobson circulation [2]; its variability and change depend on several dynamical proxies, as well as chemical processes including tropospheric pollutants and the emission of stratospheric ODS (Ozone Depleting Substances).

Irene (25.90° S, 28.22° E, 1524 m above sea level) weather station is located in the southern subtropics (see Figure 1), between two highly industrialised cities (Pretoria and Johannesburg), which are significantly marked by chemical activities. Ozone concentrations in the boundary layer are higher during austral winter due to a systematic increase of ozone precursors from the household combustion for space heating and the concentration of low-level pollutants near the surface [3]; while free tropospheric ozone variability shows a seasonal maximum during austral spring time [4], which is attributed to biomass burning activity in Southern Africa [5]. As a subtropical site, Irene is affected by low and high latitudinal dynamical processes defined by isentropic and latitudinal transport in the UT-LS (Upper Troposphere-Lower Stratosphere) and the upper stratosphere, respectively [6]. The behaviour of ozone variability over Irene is different in comparison with other subtropical stations [5,7,8]. The processes modulating ozone distribution, variability, and trend over Irene are complex due to the mentioned above chemical and dynamical processes. However, an upward trend of ozone over the site was reported in previous studies [8,9]. Thompson et al. [10] indicated a significant increase (~1 ppbv/yr) of ozone in the free troposphere over Irene between 1990 and 2007, especially during winter. A recent report released by SPARC-LOTUS (Long-term Ozone Trends and Uncertainties in Stratosphere) [11] highlights an increase of stratospheric ozone in the southern subtropics above 25 km. However, results differ from region to region and depend on methods and instruments used. By applying a multivariate model called Trend-Run on total columns of ozone (TCO) data obtained from ground-based and satellite instruments, an upward trend of about 1.7% per decade was found over Irene from January 1998 to December 2012 [8]. However, many questions remain as to the sources of this trend. Additional details were necessary to define if the observed trend is due to an increase of stratospheric ozone as a result of a success implementation of the Montreal protocol or to a systematic increase of the tropospheric ozone as consequence of anthropogenic emission of ozone precursors (NOx, VOC, CO). The Irene site in South Africa is one of the sites in the southern tropics equipped to undergo continuous ozone observations by ozonesondes (ozone profiles) and by Dobson spectrometer (total columns of ozone) in routine mode. In this paper, we combine ground-based and satellite ozone measurements over Irene collected for 20 years of continuous observations, with the aim to assess variability and trends in tropospheric and stratospheric ozone columns during the 1998–2017 periods.

Satellite data used are from Earth-Probe (EP) Total Ozone Mapping Spectrometer (TOMS), Ozone Monitoring Instrument (OMI), and Microwave Limb Sounder (MLS) instruments. This choice is due to the high quality of ozone data [12,13,14]. Good agreement between TCO obtained from a combination of TOMS and OMI with Dobson data measured over Irene from 1998 to 2012 were found [9]. OMI observations match well with Dobson data over Irene [15] and have chosen among others to construct long-term TCO data for ozone climatological analysis [16]. TOMS ozone data also have a good agreement with data from Dobson instruments installed worldwide including at Irene station [17]. In this paper, monthly average values of TCO recorded by OMI and TOMS are used to complete Dobson measurements recorded from January 1998 to December 2017.

Regarding the construction of tropospheric partial column time series, a recent product of tropospheric partial column retrieved from OMI/MLS observations recorded between August 2004 and December 2017 is used in order to complete balloon-sondes observations of SHADOZ (Southern Hemisphere Additional Ozonesonde) network. The latter provides ozonesonde profiles over the southern tropics and subtropics at different locations. Their measurement precision is approximately 5% with a vertical resolution between 50 and 100 m [18]. The ozone-soundings have been operating at Irene since 1998 under the SHADOZ programme and were used in several studies to validate satellite data and to determine Irene ozone climatology, seasonal variability, and trend [5,7,19,20,21].

In this paper, the inter-annual variability and the trends analysis were conducted separately on tropospheric and stratospheric columns of ozone using a multivariate regression model called Trend-Run [22,23]. Multivariate regression models are considered as a powerful numerical tool to study ozone variability and trends [19,24,25,26,27,28]. In these models, the indexes of the chosen parameter should be based on atmospheric forcing that has been historically accepted as having an influence on ozone variability [29]. More details about the Trend-Run model can be found in Begue et al. [23]. In the present study, we determined the various forcings contributing to ozone variability using the wavelet analysis technique. The Mann–Kendall test method was also performed on Irene ozone time series for comparison and for detecting trend change points, if any.

2. Instruments and Data

2.1. Ozone from Satellite Instruments

Ozone datasets from satellite used in this study were recorded by EP-TOMS and OMI-Aura instruments. EP-TOMS data considered is the final L2 TCO product available between 1996 and 2005. Note that TOMS instruments were successfully flown aboard different satellites: the Nimbus-7 satellite (1978–1993), Meteor-3 (1991–1994), Advanced Earth Observing Satellite (ADEOS) (1996–1997), and the latest one was aboard the Earth Probe (EP) satellite operational from July 1996 to December 2005. It was the only instrument aboard the EP satellite launched the 2 July 1996 into a polar orbit, initially at 500 km with an inclination angle of 0.98°. Due to the failure of ADEOS spacecraft in June 1997, EP was raised to 739 km of altitude during December 1997, in order to provide complete global coverage of ozone data and other atmospheric pollutant measurements including aerosol index and sulphur dioxide (SO₂). TOMS is a spectrometer with nadir viewing and had vertical and horizontal resolutions of 77 km and 50 km × 50 km, respectively. The instrument uses a single monochromator and a scanning mirror to sample the backscattered solar ultraviolet radiation at 35 sample points at three intervals along a line perpendicular to the orbital plane [12,17]. It measures the incoming solar energy and the backscatter UV radiance in six discrete wavelengths (379.95, 359.88, 339.66, 331.06, 317.35, and 312.34 nm) selected by a chopper wheel in the light path behind the monochromator grating. The reader may refer to Petropavlovskikh et al. [12] for more details about TOMS instrument and its products. TOMS Ozone data employed in this work is the V8 L3 product overpass over Irene station and it is available to download on the following website https://acdisc.gesdisc.eosdis.nasa.gov/data/EarthProbe_TOMS_Level3/TOMSEPOVP.008/. In terms of the data uncertainties, the absolute error is ±3%, the random error is ±2% (though somewhat higher at high latitudes), and the drift after 1.5 years of operation is less than ±0.6% [12]. As EP-TOMS data have a temporal cover of 10 years (from July 1996 to December 2005), TCO from OMI overpass for Irene, from August 2004 to December 2017, was combined with that from EP-TOMS in order to construct ozone dataset having temporal cover of 20 years of observation (from January 1998 to December 2017). As for TCO values from OMI instrument, they were downloaded from the Aura validation data center: http://avdc.gsfc.nasa.gov/pub/data/satellite/Aura/OMI/V03/L2OVP. This is the OMI-TOMS V4 L2 product recorded above the Irene site. Note that the V3 L2 product has good agreement with TCO recorded by the EP-TOMS instrument overpass Irene [9]. OMI-TOMS ozone data are retrieved using two wavelengths: the 317.5 and 331.2 nm are employed under most conditions, while 331.2 and 360 nm are specially employed for conditions of high ozone concentration and high solar zenith angle. The precision of this L2 product is about 3%. Previous studies had shown good agreement between this L2 product and ground-based (Dobson and SAOZ (Système D’Analyse par Observations Zénithales)) measurements in many stations located around the southern tropics and subtropics [8,15,30]. OMI was among the instruments aboard Aura satellite launched on 15 July 2004 into a near-polar helio-synchronous orbit at an altitude of approximately 705 km. It is a compact nadir-viewing instrument operating at a spectral resolution of the order of 0.5 nm and measures trace gas concentration and aerosol index in three broad spectral regions (UV-1, UV-2 and VIS). In terms of spatial coverage, OMI has a viewing angle of 57° and covers a swath width of about 2600 km. The ground pixel size of each scan is 13 × 24 km² in the UV-2 (310–365 nm) and visible (350–500 nm) channels, and 13 km × 48 km for the UV-1 (270–310 nm) channel. For further details about OMI instrument and its operation mode, the reader may refer to the OMI Algorithm Theoretical Basis Document Volume II [13].

MLS is among the four instruments aboard Aura satellite. The MLS instrument was first on board the UARS satellite, which was operational from September 1991 to August 2001, then for the second time aboard the EOS-Aura satellite launched on 15 July 2004 into a near-polar helio-synchronous orbit at an altitude of approximately 705 km. The MLS instrument observes microwave emissions from the Earth’s limb in the “forward” direction of the Aura orbit [14]. The MLS provides around 3500 vertical profiles of ozone, temperature, and other chemical trace gases per day from the troposphere to the mesosphere. The measurements are achieved every 1.5° along the satellite orbit track, which corresponds to a horizontal spacing of ~165 km. The observations are achieved from 82° S to 82° N at two fixed local solar times and cross the tropical region at around 01:30 a.m./p.m. (UT).

A partial tropospheric column of ozone from satellite used here was obtained by a combination of TCO from OMI and MLS profile recorded between August 2004 and December 2017. These are global data and are available at the following link https://acd-ext.gsfc.nasa.gov/Data_services/cloud_slice/new_data.html. Daily OMI/MLS tropospheric ozone was first determined by subtracting co-located MLS stratospheric column ozone from OMI total column ozone. Stratospheric column ozone from MLS was spatially interpolated (2D Gaussian/linear latitude-longitude interpolation) each day to fill in between the actual along-track measurements. Monthly mean fields were then determined by averaging all available daily data within each month. OMI total column ozone was filtered for near clear-sky conditions by including only measurements when coincident OMI reflectivity was less than 0.3. Further details for the OMI/MLS method are given by Smit et al. [31].

2.2. Ozone Data from Ground-Based Instrument

The Irene site is operated by the South African Weather Service (SAWS). It is a part of the SHADOZ network since 1998. Ozonesonde were launched bimonthly from October 1998 to December 2003, after then, frequency has changed to weekly. Indeed, the number of launches can reach 4 per month. A total of 341 profiles was recorded during the study period (1998–2017), with an interruption period between 2007 and 2012. The station uses an electrochemical cell (ECC) ozonesonde device, equipped with radiosonde for temperature, humidity, pressure, and ozone partial pressure measurements from ground to the altitude where the balloon burst occurs (~26–32 km). During the measurements, data are recorded at 2-s intervals and are subsequently averaged over 100-m-height intervals [5]. The release time of the balloon is usually fixed at 8h00 UTC, corresponding at 10h00 (local time) in South Africa. For information regarding ECC ozonesonde instrument validation, operating mode, and algorithm used for data retrieval, the reader may refer to Smit et al. [32] and published articles of Thompson et al. [18,33]. Ozone profiles employed in this work were downloaded from the SHADOZ website http://croc.gsfc.nasa.gov/shadoz/. These profiles were used to calculate the Tropospheric Partial Column (TPC) of ozone. The obtained times-series were completed by using satellite data.

Moreover, TCO observations are made at Irene site using a Dobson Spectrophotometer. The principle of measurement is based on comparing the relative intensities of selected pairs of UV wavelengths, emanating from the sun, moon, or zenith sky. Thus, by measuring the relative intensities of suitably selected pair wavelengths with the Dobson instrument, it is possible to determine how much ozone is present in a vertical column of air extending from ground level to the top of the atmosphere in the vicinity of the instrument. The result is expressed in terms of a thickness of a layer of pure ozone at standard temperature and pressure [34] and is given in Dobson Units (DU). There are currently about 60–70 operational Dobson instruments around the world, which form part of the primary ground-based ozone observational network.

The first South African total column ozone measurements were made during the 1964–1972 period of time with Dobson-089 operating at Pretoria. This was as a result of the International Geophysical Year during the mid-1950s when Earth observations were greatly enhanced over the globe. Due to skills capacity and instrument challenges, Dobson-089 was decommissioned in 1972. After that, reinstating South Africa’s commitment to the Vienna Convention, the SAWS now operates two Dobson ozone spectrophotometers, Dobson-089 at Irene near Pretoria (25.9° S, 28.2° E) since 1989, and Dobson-132 at Springbok (29.7° S, 17.9° E) since 1995. Both instruments have been regularly calibrated with reference to the world standard Dobson-083 instrument at Mauna Loa, Hawaii. The most recent calibration was undertaken during October 2019, at Irene. The Dobson’s calibration remains verified and stable within 1% of the World Standard Reference instruments housed by NOAA, USA, and the European Standard Instrument operated by the Deutscher Wetterdienst (DWD), Germany. Dobson calibration traceability can be found at this link http://www.o3soft.eu/dobsonweb/welcome.html. The Dobson observations (Direct Sun Observations only) are processed for daily average values by software originally developed by Stanek [35].

3. Ozone Time Series Construction

3.1. Construction of TCO Time Series

Ozone product retrieved over Irene station from the above-mentioned instruments were used to construct the TCO time series. Figure 2a shows TCO monthly mean values obtained from Dobson (black line), OMI (red dotted line), and TOMS (blue dotted line). It is clear from Figure 2a that Dobson measurements agree with OMI and TOMS. Ground-based data fit well with satellite data; the correlation coefficients are higher than 0.90 and the mean difference is less than 2 DU. As the Dobson spectrometer has a good spectral resolution in comparison to satellites, we used ground-based measurements for the ozone stratospheric partial column calculation. However, due to the consistency between satellite and ground-based monthly mean values, satellite measurements are specifically used to fill the missing ground-based data in order to construct a complete TCO time series. The resulting TCO time series is presented in Figure 2b over the study period from January 1998 to December 2017.

3.2. Construction of Tropospheric Column Ozone Time Series

Similar to the construction of the TCO time series presented in the previous subsection, tropospheric ozone column time series was built by using ground and satellite observations. Satellite data used are the OMI/MLS tropospheric partial column of ozone (TPC) as mentioned in the earlier Section 2.1, while ground-based data are ozone profiles recorded over Irene from ozonesondes in the framework of SHADOZ network. The example depicted in Figure 3 shows that the ozone profile recorded from ozonesonde was in good agreement with the MLS profile recorded on the same day for the study location.

The TPC was retrieved by integrating the balloon-sonde ozone profile from ground up to the tropopause height, defined here as the lapse-rate tropopause (LRT). LRT is defined as the lowest level at which the lapse rate is found to be less than 2 Kelvin.km⁻¹ and it remains within this level for the next 2 km (World Meteorological Organisation, 1957). Other tropopause definitions, such as the cold-point tropopause (CPT), the ozone tropopause (OT), or dynamic tropopause, may also be used to identify the tropopause height [21]. However, as explained in Ziemke et al. [31], OMI/MLS tropospheric ozone is retrieved using the LRT. In order to harmonise tropospheric ozone dataset, it is found justifiable to determine the LRT height from balloon-sonde temperature profiles as well. Table 1 presents the LRT climatology over Irene. Figure 4 shows the obtained monthly mean climatological temperature and ozone-concentration time-height cross-sections, with LRT and CPT seasonal variations superimposed.

As expected, the LRT heights are lower than CPT and show more pronounced seasonal variations. On average, the LRT height was at 15.27 ± 0.60 km, while the CPT was located about 2 km higher (17.3 ± 0.8 km). This is in agreement with results published by Sivakumar et al. [21], based on 11 years of ozonesonde observations (1998–2008) from SHADOZ stations. They reported the LRT and CPT heights over Irene at 15.6 km and 17.2 km, respectively. Furthermore, it can be seen on Figure 4 that the LRT height was characterized by a well-defined seasonal cycle, high during the austral summer (16.19 km) and low by winter (14.64 km).

Plot (a) of Figure 5 depicts the TPC values obtained from balloon-sonde profiles (blue) together with the OMI/MLS values (in black). The correlation between OMI/MLS and ozonesonde is evaluated at 0.79. However, some differences between both observations are noticed from the Figure. To assess the observed differences, the relative difference (RD) between satellites and ozonesonde data with respect to the ozonesonde data was calculated, following the equation:

(1)${RD}_m=100\times \frac{{\mathrm{SAT}}_m-Sond{e}_m}{Sond{e}_m}$

where “m” index refers to the month where Aura satellite and ozonesonde observations were performed. Figure 5b depicts the RD, which is within ±20% except in July and October 2015 where an unusual increase of ozone is observed by ozonesonde in the upper troposphere. The case of these events will be investigated in a future study.

The mean bias (in percentage) of both observations is taken as the average of RD values as described by [15]. It is found to be negative (−5.07%) indicating a slight underestimation of TPC recorded by satellites when comparing to ground observation. The root mean square (RMS) associated to these two datasets (ozonesonde and satellite) is assessed at 3.11 DU (8.45%). This latter value was used to adjust satellite data with respect to ozonesonde. The RMS represents the mean systematic error between satellite and ground-based values. As the mean bias error is negative, satellite data are adjusted to ground-based data following the equation:

OMI/MLS _adjusted (m) = OMI/MLS (m) + RMS(2)

The OMI/MLS time series was adjusted to radiosonde measurements and together with ozonesonde measurements is superimposed on Figure 6a. It is evident from the Figure that adjusted satellite data is in good agreement with ozonesonde observations. Thus, satellite monthly values adjusted to ozonesonde measurements were used to fill the gaps in the ozonesonde time series. Obtained time series from this process is plotted in Figure 6b.

3.3. Construction of Stratospheric Column Ozone Time Series

The time series of stratospheric partial column of ozone (SPC) is obtained by subtracting the tropospheric column of ozone from the TCO and assuming that ozone concentrations are negligible above the stratopause. The obtained SPC time-evolution data are shown in Figure 7. As expected, the SPC shows similar seasonal variations as TCO (see Figure 2). It is mainly driven by the annual component, with minimum and maximum values around the equinoxes in April and September, respectively.

4. Methods Used for Data analysis

4.1. Wavelet Method

Prior to quantifying the inter-annual variability of tropospheric and stratospheric ozone over Irene, forcing modes were investigated by using wavelet decomposition. Wavelet analysis is a common tool for analyzing localized variations of discrete temporal signal. As explained by [36], once the temporal signal is decomposed into time-frequency space, it is possible to determine the periods of dominant modes and how those modes vary throughout time. Wavelet transform analysis has been used in previous works to study the variability of total column of ozone [37,38] and ozone profiles [39]. Wavelet transform methods were used to determine the different forcing involved in the long-term variability of ozone in southern Brazil [38]. They highlighted a variability in ozone dominated by an annual cycle, a 600-days mode associated to the Quasi-Biennial Oscillation (QBO) and modes associated to the first and second harmonics of the 11-year solar cycle.

The wavelet transform of a discrete data series is defined as the convolution between the series and a scaled translated version of the chosen wavelet function [36]. Note that in order to be “admissible” as a wavelet, this function must have a mean of zero and be localized in both time and frequency space. In the present work, we have used Morlet wavelet function. It consists of a plane wave modulated by a Gaussian function and expressed as followed:

(3)${\phi }_0\left(t\right)={\pi }^{-1/4}·e^{i{w}_{0}t}·e^{-t^2/2}$

where w₀ is a non-dimensional frequency while t is a term of non-dimensional time. Thus, for a giving X_n time series composed of “n” number of observation varying from 0 to N − 1 and spaced by a time interval of δt, the continuous wavelet transform of the discrete time series X_n is defined as followed:

(4)$W_n\left(f\right)={\displaystyle \sum }_{t=0}^{N-1}{X}_t{\phi }^{*}\left[\frac{\left(t-n\right)·\delta t}{f}\right]$

where the (*) indicates the complex conjugate of the wavelet function. The continuous wavelet transform gives the temporal signal behavior along the localized time “$\left(t-n\right)$” for a wavelet scale (frequency) “f”. However, the global wavelet spectrum provides objective and consistent estimation of temporal signal power spectrum and highlights the average intensity of the dominant and significant frequency bands. Note that the power spectrum is the square of amplitude of the wavelet transform ${\left|W_n\left(f\right)\right|}^2$. Thus, the global wavelet spectrum is found by averaging the wavelet transform power spectrum in the temporal scale. It is expressed as follows:

(5)$\overline{W^2}\left(f\right)=\frac{1}{N}{\displaystyle \sum }_{n=0}^{N-1}{\left|W_n\left(f\right)\right|}^2$

4.2. Trend-Run Model

Trend-Run model was applied to the constructed time series of tropospheric and stratospheric ozone partial columns in order to quantify the contribution of different modes to the ozone variability. As explained above, The Trend-Run is a multivariate model adapted at Reunion University to study ozone and temperature inter-annual variability and to estimate the associated trends. In this study, the input parameters are based on monthly mean values covering the same period of time: ozone data (TPC or SPC), annual and semi-annual oscillations, QBO at 30 hPa pressure level, El Niño-Southern Oscillation (ENSO), and the 11-year solar cycle. Overall, the Trend-Run model is designed to simulate the geophysical signal using the main forcings. It thus enables to determine the contributions of the forcings by minimizing the residual term. The latter was used to estimate the linear trend based on the least-square method. For more details on Trend-Run model, the reader may refer to previous works [8,22,23].

In the present study, the model was applied to estimate TPC and SPC trends over Irene by taking into account the 5 above forcings found to contribute to ozone variability, following the time equation:

(6)$\left[O_3\right]\left(t\right)=a_o+ a_1·SAO\left(t\right)+ a_2·AO\left(t\right)+ a_3·QBO\left(t\right)+ a_4·ENSO\left(t\right)+ a_5·Solar\left(t\right)+r\left(t\right)$

where coefficients $a_k$ are computed by the model for each forcing, and r is the residual term, assumed to be consisting of Trend and noise:

(7)$\mathrm{r}\left(\mathrm{t}\right) = \mathrm{Trend}\left(\mathrm{t}\right)+ \mathsf{\epsilon }$

4.3. The Mann–Kendall Method

4.3.1. The Mann–Kendall Test

The MK (Mann–Kendall) non-parametric test is a statistical tool usually used to quantify the significance of time series trends [40]. The strength of the trend depends on the magnitude, sample size, and variations of data series [41]. For X_i (i = 1 to ‘n’) temporal data, the MK test statistics “S” is given as follows:

(8)$S={\displaystyle \sum }_{i=1}^{n-1}{\displaystyle \sum }_{j=i+1}^{n}sign\left(X_j-X_i\right)$

where, X_i and X_j are sequential data in the times series, ‘n’ the length of time series, and

(9)$sgn\left(X_j-X_i\right)=\{\begin{array}{c}+1, \mathrm{if} X_j-X_i>1\\ 0, \mathrm{if} X_j-X_i=0\\ -1, \mathrm{if} X_j-X_i<1\end{array}$

A downward or upward trend of data is defined by positive or negative of the S value. When the length of time series is more or equal to ten (10), S behaves as a normal distribution and MK test is characterized by normal distribution with a mean µ (S) = 0 and variance defined as followed:

(10)$Var\left(S\right)=\frac{n\left(n-1\right)\left(2n+5\right)-{\sum }_{k=1}^m{t}_k\left(t_k-1\right)\left(2t_k+5\right)}{18}$

where t_k defines the ties of the kth value, and m represents the number of the tied value [42,43,44]. The standardized test statistics “Z” is computed as followed:

(11)$$\begin{gathered} Z=\frac{\mathrm{S}-1}{\sqrt{Var\left(S\right)}} \mathrm{if} \mathrm{S} > 0 \\ Z=0 \mathrm{if} \mathrm{S} = 0 \\ Z=\frac{\mathrm{S}+1}{\sqrt{Var\left(S\right)}} \mathrm{if} \mathrm{S} < 0 \end{gathered}$$

The significance of trend is assessed from the Z value. Positive or negative value of Z indicates an upward or downward trend of time series. In this paper, the MK test with significance level p = 5% was applied in order to detect if the inter-annual trend of ozone is statistically significant. It is worth noting that at the 5% significance level, the null hypothesis of no trend is rejected if the absolute value of Z is higher than 1.64 [42].

4.3.2. The Sequential Mann–Kendall Test

The sequential Mann–Kendall (SQ-MK) test is usually used to determine the approximate time when a significant trend began. This test sets up two series (u(t) and u’(t)): one is progressive and the other one is backward with time. If they cross each other and diverge beyond the specific threshold value, then there is a statistically significant trend [45]. The point where the two time series are crossing indicates the approximate time where the trend started to be significant. In this paper, the threshold values used are ±1.96 (p = 5%), with the crossing time estimating the month at which the trend started to be statistically significant. For details about the sequential Mann–Kendall test assessment, the reader may refer to Zarenistanak et al. [45].

5. Results and Discussion

5.1. Dominant Forcing Modes

The aim here was to identify the dominant forcings in the variability of TCO, SPC, and TPC ozone time series. In that regard, the wavelet transform was applied to the above-mentioned ozone time series. The obtained wavelet power spectra are showed in Figure 8. It is clear from the power spectra of all the ozone time series that the annual forcing (12-month period) has a high intensity (then it is the prevailing mode) expressed by a strong peak within the 95% significance level contours (black contours). While all the ozone time series present a signature of the 11-year solar cycle, the SPC time series (Figure 8b) seems to indicate a strong 11-year solar mode periodicity. In decreasing order of intensity, it can be noticed from the TCO and SPC wavelet spectrum the occurrence of a mode in the band of 20–32 months, and another one with a 6-month period (see Figure 8a). It is also important to mention that all the time series seem to experience some periodic occurrence periodicities that are less than 4 months.

In general, stratospheric ozone variability over Irene appears to be modulated by 4 main modes, with periods of 6 months, 12 months, 20–32 months, and 11 years corresponding to the annual and semi-annual seasonal forcings, the QBO, and the 11-year solar cycle. As for tropospheric ozone variability, in addition to the two modes of 12 months and 11 years, the wavelet spectrum shows an 80-month (7-year period) mode and a weak 6-months mode (see Figure 8b). The 7-year mode could be associated to the influence of the ENSO (El-Niño Southern Oscillation), an oceanic dynamical proxy having a periodicity varying from 2 to 8 years [36,46].

In summary, taking into account results from wavelet analysis applied to time series of stratospheric and tropospheric columns of ozone, it can be noted that the variability of ozone is expected to be meanly driven by five forcings, i.e., semi-annual (SAO) and annual oscillations (AO), QBO, ENSO, and 11-year solar cycle (Solar).

5.2. Variability Analysis

The variability analysis was conducted using the multivariate Trend-Run model described above. This model consists of a multi-regression least-square method to compute the coefficient vector $\overrightarrow{a} \left(a_1,a_2,a_3,a_4,a_5\right)$ by minimizing the sum of the residual squares following the equation:

(12)${\left[r\right]}^2={\left[O_3\left(t\right)-\overrightarrow{F}.\overrightarrow{a}\right]}^2$

where $\overrightarrow{F}$ is the forcing vector $\overrightarrow{F}\left(SAO, AO, QBO, ENSO, Solar\right)$.

This was applied to TPC and SPC time series. The latter were smoothed by a three-month moving average filter to remove possible inter-seasonal variability. The obtained coefficients are reported in Table 2.

Forcing contributions to ozone variability were computed based on methods adopted by Toihir et al. [8]. The obtained results are presented in percentage on Table 2. Overall, we found similar values as with the wavelet analysis. The annual oscillation appears to be the dominant forcing for both the TPC and SPC time series. It is also found that the QBO and the 11-year solar cycle are more prevalent in the stratosphere. This was in agreement with previous findings over subtropical regions [8,47]. Furthermore, it should be noted that the regression coefficients obtained for the QBO are negative for both time series of TPC and SPC. This indicated that the QBO was in an opposite phase compared with ozone time series over the study region. Similar results were reported by Vaz Peres et al. [37] over a subtropical site located in the south of Brazil. Moreover, this was in agreement with the study by Chehade et al. [48]. By analysing the ozone time series, they obtained negative values for the regression coefficients for the QBO forcing between 15° in latitude and polar regions. Regarding the 11-year solar cycle response, negative and positive values were found in tropospheric and stratospheric ozone time series, respectively. This means that unlike the stratosphere, the ozone variability in the troposphere is in phase opposition with the solar cycle, indicating that the increasing phase of solar index is associated to a decrease of tropospheric ozone. Similar results were reported by [49]. They found that the annual mean solar cycle response for troposphere and stratosphere columns of ozone were −2.29 ± 1.01 DU and 6.64 ± 1.53 DU, respectively, for an increase of 100 units of 10.7 cm solar flux.

In terms of percentage contributions, it can be seen from Table 2 that ~79.11% and ~69.37% of total variability of ozone is explained by SAO, AO, QBO, ENSO, and solar forcings, respectively, in tropospheric and stratospheric columns of ozone, while the corresponding residual terms show contributions in the range of 20.89% and 30.63%, respectively. Furthermore, it is evident from Table 2 that the percentage contribution of the QBO to stratospheric ozone (~8.23%) is much higher than that of tropospheric ozone (0.27%). This result was in agreement with previous published works. Based on ozone profiles recorded at Irene site for the 1998–2013 period, Toihir et al. [8] showed that the maximum contribution of QBO on ozone variability was in the 20–28 km altitude range. Previous studies [19,27,29,50] have also reported significant QBO contribution on ozone variability in the stratosphere and especially in the 20–27 km altitude range, which corresponds to the stratospheric ozone layer where the maximum of ozone is observed in tropics and subtropics. Note that QBO modulation on ozone variability is linked to the downward propagation of zonal wind throughout the time. The associated QBO downward disturbances start in the upper stratosphere and occur as far down as the lower stratosphere. They gradually attenuate and are not observable at lower altitudes when atmospheric pressure is higher than 100-hPa [47]. As reported in Table 2, this is in agreement with the low-percentage contribution of QBO in the TPC time series (0.27%). The SAO have significant influence on ozone variability for both tropospheric and stratospheric columns of ozone (see Figure 8), with a higher contribution in the stratosphere (~3.94%) in comparison with the troposphere (2.79%). According to previous works by earlier researchers [51,52], the SAO is generated by the zonal wind with its maximum influence at the stratopause, from where it propagates downward into the stratosphere and to the troposphere by decreasing intensity. This agrees with percentage of contribution found here in the present study. Regarding the ENSO forcing, it has a higher contribution to the tropospheric ozone in comparison to stratospheric ozone. In fact, unlike QBO and SAO, ENSO is a forcing generated at the sea surface due to ocean-atmosphere interactions. It is worth recalling that the wavelet analysis presented above (Figure 8) also revealed that the ENSO forcing has greater contribution in the tropospheric than in the stratospheric columns of ozone. The solar cycle contribution appeared to be greater in the stratosphere ozone (18.43%) than in the troposphere ozone (1.01%). It was reported that the solar cycle response on ozone is positive and statistically significant in the lower and upper stratosphere [53]. Moreover, the percentage contribution of the solar cycle to tropospheric ozone obtained here for the study site ($1.01\%$) is slightly less than the averaged one ($2.5\%$) reported at global scale [54], with the maximum nearby the equatorial zone and the minimum over the mid-latitude regions.

5.3. Trend Analysis

5.3.1. Obtained Results from Trend-Run

The trend analysis was performed using the Trend-Run model. A long-term linear function is generated by the model to characterize the trend index. This fuction was used among the model parameters and the trend value was computed based on its rate (regression coeficient) evolution. Model and observation time series are shown in Figure 9. Here, the regression model fits with the observations and the obtained correlation between model and observations was significant (82% and 91% for tropospheric and stratospheric ozone, respectively).

The linear function linked to ozone time-evolution obtained by the model is presented in red line. The lines reveal an increasing trend in the tropospheric column of ozone (+0.76 ± 0.05 DU per decade) and an increasing trend in the stratospheric column of ozone (−1.27 ± 0.21 DU per decade). The calculated trends are reported in Table 3.

An increase of tropospheric ozone at global scale has been reported by recent studies [55,56]. The obtained tropospheric ozone trend over Irene station is consistent with results reported in previous studies focused on Southern Africa [10,57]. An increase in free-tropospheric ozone is reported by [10] based on ozonesonde records from the early 1990s to 2008 at Irene, while [57] studied the tropospheric ozone climatology over Irene and observed an upward trend from 1998 to 2013 defined by an increase in ozone tropospheric columns from 55 to 65.6 DU in spring and from 32 to 55 DU in summer. Based on the previous works and the present study, one can conclude that the tropospheric ozone column has increased since the early 1990 up to now (2017). However, as explained by many authors, the sources are different from one season to another. Ozone increase in summer and spring seasons is linked to an increase of photochemical activities such biomass burning (agriculture reason) and greenhouse gases (urban-industry reason) implemented in Africa [4,5,57], while the winter sources can be associated not only to power-generating plants and domestic biofuel, but also by long-range transport of growing pollution in the Southern Hemisphere at an intercontinental scale [10].

The present study consisted of splitting the TCO into a tropospheric and a stratospheric column, with the aim to estimate trends in both time series. A negative change in SPC from 1998 to 2017 is observed. However, the sources of this downward trend are not well documented. Recent works highlighted negative change of stratospheric ozone at global scale [56,58,59] and argued that this negative change does not reveal an efficiency of Montreal protocol. They associate this ozone decrease mainly to dynamical variability on long- and short-time scales.

5.3.2. Mann–Kendall Trend Analysis

The trend estimation was also calculated using Mann–Kendall trend analysis method in this study. This method offers an advantage over other techniques because it provides a non-parametric test that does not require data to be normally distributed, and it also is not dependent on the magnitude of data. To approximate the time periods of the beginning of a significant trend and trend change detection points, the sequential version of the Mann–Kendall trend test was applied. This was used side-by-side with the Theil–Sen method [60,61]. The Theil–Sen method is based on the use of the Theil–Sen function, which is described in the Openair manual as documented by Carslaw [62].

Table 4 shows the summary of the Mann–Kendall test for total, stratospheric, and tropospheric columns of ozone, respectively. In terms of the TCO time series, the Mann–Kendall z-score and p-value seem to indicate that there is a significant upwards trend in data. On the other hand, the SPC data indicates a negative and significant downwards trend, while the TCO indicates a negative but insignificant trend. It is also noteworthy to mention that significant trends are those that have a z-score value of less than −1.96 and greater than +1.96, while an insignificant trend is considered to be the trend that has a z-score value between ±1.96. Moreover, the z-score statistic results presented here agree with the trend results estimated by the Trend-Run model.

Figure 10 shows the sequential Mann–Kendall statistics of progressive (Prog) u(t) and retrograde (Retr) u’(t) time series (left panel), and the Theil–Sen plots (right panel) for Irene TCO, SPC, and TPC time series. The TCO (Figure 10a) seems to have a change of trend after year 2005 where the forward/progressive trends seem to indicate a downwards trend that started to be significant in the year 2007. However, this downwards and significant trend is observed to revert back to a negative and insignificant trend where the time series is observed to sit in the region between 0 and −1.96 values. The Theil–Sen plot for TCO is shown in Figure 10b, where the 228 points were used for Theil–Sen trend estimation. On the Figure 10b,d,f, the symbols shown next to the Theil–Sen trend estimate are associated to how statistically significant the trend is. The symbols which are used depend on the value of p-value where for p < 0.001 = ***, p < 0.01 = **, p < 0.05 = * and p < 0.1 = +. The three-stars symbol indicates a highly statistically significant trend estimate while the plus sign represents the trend estimate that is not statistically significant. The Theil–Sen result shown in Figure 10b indicates that TCO experiences a significant decline of −0.11% per year over Irene for the study period. The superimposed blue line indicates the obtained linear trend, and the dashed orange lines indicate the 95% confidence interval. The SPC SQ-MK (Figure 10c) also indicates a decline of ozone concentration that begins at the start of the time series and becomes significant by the late year 2004. In 2004, a change in trend was detected as well. This was indicated by the point where the forward and the backwards trend cross each other. The significant downwards trend observed in the SPC time series is also supported by the Theil–Sen plot in Figure 10d. The Theil–Sen plot in Figure 10b indicates a significant decline of ozone concentration in the stratosphere by −0.17% per year over Irene. On the other hand, both the SQ-MK and Theil–Sen methods indicate an upwards trend of ozone concentration in the troposphere (Figure 10e,f). In terms of the change detection point, the TPC SQ-MK diagram in Figure 10e indicates that the change point of the trend was during the year 2002. The Theil–Sen test method in Figure 10f indicates a significant increase of ozone concentration by 0.36% per year over Irene.

As it can be seen in Table 3 and Table 4, the results of the trend estimates by the Trend-Run model and the Theil–Sen method are broadly similar and consistent. Indeed, both methods result in a decrease in the stratospheric column (−0.56% and −1.7% per decade, respectively) and an increase in the tropospheric column (+2.37% and +3.6%, per decade, respectively). It should also be noted that the SQ-MK method offers the advantage to detect trend turning-points over time. Moreover, the results presented here indicate that the slowing down of the total ozone decline is somewhat due to the contribution of the tropospheric ozone concentration. Thus, for accurate assessment of the ozone variability and trend in the atmosphere, it is always better to divide the total column of ozone in stratospheric and tropospheric columns.

6. Conclusions

The work presented in this paper investigated the variability and trends of tropospheric and stratospheric ozone over Irene, South Africa using ground-based and satellite observations between 1998 and 2017. The ozone data included TCO, tropospheric partial column, and derived stratospheric ozone column. The comparison between ground-based observations of ozone and satellite observations showed a good agreement between the two datasets in terms of correlation and relative difference. This allowed for ozone-times to be constructed, as data gaps in the ground-based data could be filled with satellite data.

The constructed ozone time series was used to investigate the forcings, which have a significant influence on the variability of tropospheric and stratospheric ozone by using the wavelet analysis as well as the Trend-Run model to determine the contribution of each forcing. Furthermore, the Mann–Kendall test was used to detect significant trends in the ozone time series. The combined contribution of the AO, SAO, QBO, ENSO, and 11-year solar cycle was 69% and 79% for stratospheric and tropospheric ozone, respectively. It was evident that the annual oscillations have the highest contribution towards ozone variability, contributing approximately 78% and 35% in the troposphere and stratosphere, respectively. In the stratosphere, the QBO response was very important and reached 8% and was out of phase with the temporal evolution in both stratosphere and troposphere. Both ENSO and SAO were in phase with the time evolution of ozone. The contribution of ENSO was most evident in the troposphere while the contribution of SAO was higher in the stratosphere. The signature of the solar cycle is well observed in the time series of stratospheric ozone and where its contribution to the total variability was very high (18.43%) compared to the troposphere (1.01%). Here, its response was negative, which means that the increasing phase of sunspot index is associated with a decrease of tropospheric ozone levels.

The Mann–Kendall test confirmed the Trend-Run model results and indicated an upward trend of tropospheric column of ozone and a downward trend of stratospheric column of ozone during the study period. The inter-annual increase of tropospheric column of ozone has become systematic with a rate of about 2.37% per decade. However, the trend of the stratospheric ozone is not monotonic: the negative trend of ozone is observed to be very significant from 2004 to 2012 and then from 2014 towards the end of the study period. These results suggest that the positive tropospheric ozone trend has resulted in a slowed decline of total ozone.

Author Contributions

Conceptualization, H.B. and A.M.T.; methodology and software, H.B., A.M.T., and N.M., validation and data curation, G.C., H.B., A.M.T., original draft preparation and writing, H.B., A.M.T., and N.M.; review and editing, D.J.d.P., N.M., V.S., G.C., and N.B.; project administration and funding acquisition, H.B. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded jointly by the CNRS (Centre National de la Recherche Scientifique) and the NRF (National Research Foundation) in the framework of the IRP ARSAIO and by the South Africa/France PROTEA Program (project No. 42470VA), and by Université de la Réunion through the OMNCG (Observatoire des Milieux Naturels et des Changements Globaux) federation of the OSU-Reunion (Observatoire des Sciences de l’Univers—La Réunion).

Acknowledgments

This work is undertaken in the framework of the French South-African International Research Group IRP-ARSAIO (International Research Project—Atmospheric Research in Southern Africa and the Indian Ocean) supported by the NRF and CNRS and by the Protea program (project No. 42470VA). The authors are thankful to all the Irene Weather Observers, under the supervision of Peter Mabasa and Misokuhle Bulayane for the many daily total ozone observations over a period of three decades.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Geographical location of the Irene site (red symbol), a subtropical site in the southern hemisphere.

Figure 2. (a) Monthly average values of total columns of ozone (TCO) measured by Dobson (red), Total Ozone Mapping Spectrometer (TOMS) (dotted red line), and Ozone Monitoring Instrument (OMI) (dotted blue line) instruments over Irene station. (b) Time-evolution of TCO constructed combining Dobson and satellite data.

Figure 3. Vertical distribution of ozone concentration recorded on 15 February 2017 by balloon-sondes (blue line) and Microwave Limb Sounder (MLS) (black line) satellite over Irene.

Figure 4. Monthly mean climatological temperature (a) and ozone-concentration (b) time-height cross-sections as derived from ozonesonde observations at Irene from 1998 to 2017. Cold-point tropopause (CPT) and lapse-rate tropopause (LRT) monthly variations are superimposed on Figure 4a with white dashed and continuous red lines, respectively.

Figure 5. (a) The time-evolution of monthly average values of tropospheric partial column of ozone (TPC) recorded from ozonesonde (1998–2007, 2012–2017) (see blue line). OMI/MLS TPC values are shown with a black line (2004–2017) and (b) relative difference obtained during the period where radiosonde and satellite measurements are achieved.

Figure 6. (a) The time-evolution of total ozone adjusted using satellite data on ozonesonde observation (in blue line) and (b) time-evolution of TPC as obtained from a combination between radiosonde and adjusted satellite data.

Figure 7. Time-evolution of stratospheric partial column of ozone obtained by subtracting the tropospheric column of ozone from the total column of ozone.

Figure 8. Wavelet and global power spectra obtained by applying the Morlet wavelet to (a) total column ozone (TCO), (b) stratospheric partial column (SPC), and (c) tropospheric partial column (TPC) time series. The black contour encloses regions of greater than the 95% confidence level and the U-shaped curve indicates the cone of influence.

Figure 9. Time evolution of monthly ozone values: panel (a) and panel (b) represent stratospheric and tropospheric partial columns, respectively. Observational data are represented by blue line and simulated data are presented in black. The straight red lines illustrate the linear trend obtained by the Trend-Run model.

Figure 10. The inter-annual variation of (left panel) sequential statistics values of forward/progressive (Prog) u(t) (solid red line) and backwards/retrograde u’(t) (black solid line) obtained by the sequential Mann–Kendall (SQ-MK) test for TCO (a), SPC (c), and TPC (e). The right-hand side panel is the Theil–Sen plot for the TCO (b), SPC (d), and TPC (f) time series.

Table 1

Monthly mean climatological heights of LRT as derived from balloon-sonde temperature observations at Irene and during the study period (1998–2017), and number of profiles used.

	JAN	FEB	MAR	APR	MAY	JUN	JUL	AUG	SEP	OCT	NOV	DEC
LRT (km)	15.97	16.19	15.98	15.72	14.70	14.84	14.82	14.64	14.36	15.28	15.31	15.54
N. profiles	25	20	30	30	28	28	27	26	24	32	26	26

Table 2

Regression coefficients and percentage of contributions of forcings as derived from the Trend-Run model applied to tropospheric and stratospheric columns of ozone at Irene site.

Forcing	Stratospheric Column of Ozone (SPC)		Tropospheric Column of Ozone (TPC)
	Reg. Coeff.	% Contribution ± 1σ	Reg. Coeff.	% Contribution ± 1σ
Annual (AO)	7.20	35.38 ± 1.29	32.12	72.82 ± 22.82
Semi-annual (SAO)	1.84	3.94 ± 0.33	5.54	2.79 ± 1.09
QBO	−2.28	8.23 ± 0.62	−0.31	0.27 ± 0.04
ENSO	1.66	0.42 ± 0.14	1.40	2.22 ± 0.34
11-years Solar	7.69	18.43 ± 0.63	−1.07	1.01 ± 0.80
Total contributions		69.37		79.11
Residual term		30.63		20.89

Table 3

Tropospheric and stratospheric ozone trend values as derived by the Trend-Run model at Irene station.

Ozone Trend Estimates ± 1σ	O3 Stratospheric Column	O3 Tropospheric Column
DU/decade	−1.27 ± 0.21	0.76 ± 0.05
%/decade	−0.56 ± 0.01	2.37 ± 0.17

Table 4

Mann–Kendall statistics summery.

Data	z-Score	p-Value	Trend Estimate
Total Column of Ozone	2.28	0.002	−0.11% per year
Stratospheric Column of Ozone	−2.75	0.005	−0.17% per year
Tropospheric Column of Ozone	−1.22	0.222	+0.36% per year