- Email: [email protected]
Detecting variability changes in Arctic total ozone column
ARTICLE IN PRESS
Journal of Atmospheric and Solar-Terrestrial Physics 68 (2006) 1383–1395 www.elsevier.com/locate/jastp
Detecting variability changes in Arctic total ozone column Franc- ois Borchia,, P. Naveaub, P. Keckhuta, A. Hauchecornea a
Service d’Ae´ronomie, IPSL CNRS, Verrie`res le Buisson, France Laboratoire des Sciences du Climat et l’Environnement, IPSL CNRS, Saclay, France
b
Received 9 January 2006; received in revised form 17 April 2006; accepted 2 May 2006 Available online 7 July 2006
Abstract To better assess inter-annual variance changes in ozone time series, we study and propose a statistical procedure based on recent advances in wavelet multi-resolution analysis. This approach, novel to the field of ozone analysis, has the advantages of detecting significant changes in variance and of characterizing the distribution of these changes, including not only their timings but also their strengths. As a test case for our method, we study total ozone column time series from the TOMS (Total Ozone Mapping Spectrometer) Version 7 in the Arctic region over the period 1981–1992. This statistical procedure allows us to analyse the well-known winter–spring transition and the inter-annual variability behaviour. For this specific example, two distinct types of winter–spring transition can be clearly identified for the ozone series. As expected from dynamical studies of the vortex evolution, the first type is constituted of the years for which there exists a very abrupt and temporally localized ozone variability change. More precisely, we detect a significant peak in ozone variance at two multi-resolution wavelet levels. The second type is characterized by smooth ozone variance changes over time and space. Moreover, the Arctic vortex breakdown timing seems to occur, on average, 1 month earlier than our detected changes in ozone variability. These results suggest a potential link between these two phenomena. r 2006 Elsevier Ltd. All rights reserved. Keywords: Stratosphere; Total ozone column; Arctic vortex; Change-point models; Wavelet analysis
1. Introduction The mean change of stratospheric ozone has been subject to many studies (WMO, 2003), but few statistical investigations on its short-term variance have been published because detecting rapid temporal variance changes in stratospheric ozone represents a series of statistical challenges. Classical Corresponding author. Tel.: +33 1 64 47 43 47; fax: +33 1 69 20 29 99. E-mail address: [email protected] (F. Borchi).
1364-6826/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.jastp.2006.05.011
linear multi-regression techniques cannot be directly implemented for these stratospheric studies. This is due to the complexity of the different spatiotemporal scales contained in such stratospheric ozone observations and the non-stationary nature of these time series. Indeed, variability of polar stratospheric ozone depends on various internal (ozone chemistry, atmospheric circulation, temperature changes, etc.) and external (solar activity, volcano eruptions, etc.) processes and complex feedbacks (Schnadt et al., 2002). Particularly, polar ozone transport is strongly associated with
ARTICLE IN PRESS 1384
F. Borchi et al. / Journal of Atmospheric and Solar-Terrestrial Physics 68 (2006) 1383–1395
planetary waves activity and the strength and deformation of the vortex (Fusco and Salby, 1999; Marchand et al., 2003). Multi-resolution analysis is becoming more popular in atmospheric sciences. Whitcher and his collaborators wrote a series of articles on multiresolution analysis with some applications to geophysics (Whitcher et al., 2000; Whitcher and Jensen, 2000; Whitcher et al., 2002). With respect to the analysis of ozone time series, Echer (2004) recently performed a wavelet decomposition to total ozone column that was observed by the Nimbus-7 Total Ozone Mapping Spectrometer (TOMS) instrument from 1979 to 1992 and detected a series of different time scales: high frequency oscillations (2–4 months), semi-annual variation (6 months), annual variation, quasi-biennial oscillation (QBO), El Nino-Southern Oscillation (ENSO) and solar cycle. Although the scope of our work is closely related to this past study in terms of data and methods, our research differs in two fundamental ways. First, our study does not focus on decomposing total ozone column at different scales per se but instead our goal is to assess inter-annual variances changes to detect breakpoint timings in stratospheric polar ozone measured at different locations over 12 years. In other words, the multiresolution analysis helps us but it is only a first step. We propose and implement other statistical procedures to go much further than displaying, comparing and commenting wavelet decompositions. This remark leads to the second difference with Echer’s study. By taking advantage of test statistics, we are able to detect significant variance changes within the small temporal scales obtained from multi-resolution decomposition. To our knowledge, such a statistical detection procedure of breakpoint timings is novel with respect to ozone time series analysis and it provides a more complete description of ozone variation. Before closing this introduction, it is important to recall that many authors (Coy et al., 1997; Shindell et al., 1998; Randel and Wu, 1999; Waugh and Randel, 1999; Waugh et al., 1999; Zhou et al., 2000) have studied the springtime persistence of the Arctic vortex that have an impact on ozone depletion. Our objective is different here: we aim to detect the interannual variance changes in ozone time series, but not the vortex evolution. To clarify this point, we remind the reader that the polar vortex breaking can either be characterized by high potential vorticity(PV) with extremes in zonal winds U within the
region of steep PV gradients at the vortex edge (Butchard and Remsberg, 1986), or by having its size shrinks to less than 1% of the Earth surface area (Zhou et al., 2000). Alternatively, zonal wind U and PV can be used together. Waugh et al. (1999) compared these three diagnostics (U alone, PV alone, both) when determining the breakdown timing of both the Arctic and Antarctic vortices over the last 40 years. He concluded that the three methods could provide similar results about vortex persistence. In comparison, we will neither look at PV nor at wind data. Instead, we solely focus on the total ozone column time series from TOMS. Working with ozone, instead of PV or wind data, shed a new light on the analysis of breakpoint timings because they are different for total ozone column and the vortex. Hence, our work should be viewed as a complement to previous dynamical studies based on vortex persistence, which validates our methodology. This paper is organized as follows. In Section 2, we present the TOMS data set under study and the multi-resolution approach is briefly explained. In Section 3.1, we detail the statistical challenges associated with the detection of abrupt ozone changes. A wavelet-based analysis is undertaken to solve this detection problem. Section 3.2 presents a summary of breakpoint distributions for each year. It is followed by a binary classification for the 12 years. In Section 3.3, an inter-annual breakpoint evolution is performed to determine the differences between winter and spring ozone variabilities. We conclude in Section 4. 2. Statistical analysis 2.1. Description of the ozone time series To illustrate our statistical approach, we choose to analyse data sets from the TOMS version 7 (McPeters, 1996). This instrument provides a large coverage of ozone abundance measurements over a long continuous period with a good temporal sampling. This allows us to investigate the shortterm variability on an inter-annual basis. The data coverage during the Nimbus 7 TOMS lifetime spans from November 1, 1978 to May 6, 1993. The individual TOMS measurements have been averaged into grid cells covering 11 of latitude by 1.251 of longitude from 901S to 901N, from 1801W to 1801E. While total ozone column integrates different layers (each one having its own variability), it is
ARTICLE IN PRESS F. Borchi et al. / Journal of Atmospheric and Solar-Terrestrial Physics 68 (2006) 1383–1395
strongly associated with stratospheric levels around 20 km where ozone has a maximum. Total ozone column is inferred from the differential absorption of scattered sunlight in the ultraviolet range. TOMS directly measures the ultraviolet sunlight scattered by the Earth’s atmosphere. This method does not allow for obtaining ozone data during periods
without sunlight. This constraint implies that most of the winter season in polar regions are not analysed in our study. More precisely, we have selected a spatial grid of six latitudes equally spaced from 64.51N to 74.51N and of 18 longitudes from 1801W to 1801E. This grid selection allows us to investigate ozone variability at high latitudes. Log-diff Transform
450
0.2
400
0.1 log(o(t)/o(t-1))
Ozone
Long=20.625 Lat=68.5
350
0.0
300
-0.1
250
-0.2 0
1385
100 200 300 Time in days (1990)
0
100 200 300 Time in days (1990) Log-diff Transform
Long=20.625 Lat=68.5
0.2 500
0.1 log(o(t)/o(t-1))
Ozone
450
400
350
0.0
-0.1 300
-0.2
250 0
100 200 300 Time in days (1987)
0
100 200 300 Time in days (1987)
Fig. 1. TOMS daily ozone time series in 1990 (top row) and in 1987 (lower row) at 68.51N and longitude 20.6251E. The second column represents the transformed original time series after taking the logarithmic difference defined by Eq. (1). The dotted vertical line corresponds to the annual breaking point timing found by our wavelet analysis (see Figs. 2 and 3 for more the determination of such timings).
ARTICLE IN PRESS 1386
F. Borchi et al. / Journal of Atmospheric and Solar-Terrestrial Physics 68 (2006) 1383–1395
The upper-left panel of Fig. 1 represents the daily ozone evolution in 1990 at a specific location, latitude 68.51N and longitude 20.6251E, near Kiruna in Sweden. In this panel, it is clear that an abrupt ozone change in variability occurs around day 120, at the end of April (the vertical dashed line). This example corresponds to the typical evolution of ozone, a rapid increase associated with a strong variability during the winter followed by a slow decay with lower variability during the spring. Hence, we call such a variance change as a ‘‘breakpoint timing’’, i.e. it is the date when the strong winter variability in the ozone time series switches abruptly to the smaller spring variability at a given location and for a given year. But, such behaviour in ozone variance changes does not occur for each year. For example, the lower-left panel of Fig. 1 indicates that the variability changes in 1987 took place at a smoother pace than 1990. Before presenting the details of our statistical analysis, we would like to make a few comments. First, we would like to emphasize that our definition of breakpoint timings, neither corresponds to yearly ozone maxima (see the left panels of Fig. 1), nor to the ones obtained with PV area diagnostic (less than 1% of the Earth’s surface area) in Zhou et al. (2000). Second, our strategy is to associate a breakpoint date for every particular location of our spatial grid, whenever such an association is possible. Consequently, this approach allows us to estimate the distribution of breakpoint timings for a given year. This probabilistic scheme could be seen as restrictive with respect to our main goal, we do not derive a unique date for a given year and an entire region. However due to the complexity of ozone variability in such region, we prefer to assume that each location is influenced by the complex phenomena that bring a strong stochastic component, i.e. spatially similar but different breaking dates. Under this assumption, a probabilistic approach is not only possible but necessary. It may not even be possible to associate breakpoint timings for locations with weak ozone variability for a given year. For example in 1987, we will see that the breakpoint dates for many locations are not always clearly temporally defined. This point evokes the important question of how to determine the level of confidence of all breakpoint timings. This issue will be carefully addressed in this paper. For each location we will test the significance of all breakpoint timings derived from our statistical procedure.
2.2. Removing the seasonal variation One possible option to detect rapid changes in daily ozone time series is to follow the next two steps. First, the seasonal variation would be removed, either by a parametric regression (e.g. polynomial fit) or by a filtering procedure (low pass filter). Second, the variability contained in the residuals would be studied. Although of interest, this two-step approach contains a few drawbacks. A parametric method to model the seasonality implies extra parameters and additional assumptions. A filtering procedure may discard some important information for a range of scales. To remove the seasonal variation, we instead apply the following transform to our daily time series: O ðtÞ ¼ log OðtÞ log Oðt 1Þ, for each year and location.
ð1Þ
In the above formula, O(t) represents the TOMS total ozone column for a given day, say t. The right panels of Fig. 1 show the result of such a logarithmic transformation. This transformation has been used extensively in finance (Gencay et al., 2002). From Fig. 1, one can see that the seasonal variation has been successfully removed. Besides its simplicity of implementation, this log-difference has the advantage of removing any smooth (i.e. polynomial) trend. Others procedures could have also been applied. For the remainder of this paper, we will work with the transformed data instead of the raw ozone time series. 2.3. Multi-resolution analysis Wavelets are fundamental building block functions, analogous to the trigonometric sine and cosine functions. The Fourier transform extracts details from the signal frequency, but all information about the location of a particular frequency within the signal is lost. In comparison, the multiresolution analysis known as the ‘‘zoom-in, zoomout’’ property makes wavelets particularly appealing for this study, because they are localized in time and the signals are examined using widely varying levels of focus. In this article, we work with a nondecimated wavelet transform, which gives a much better temporal resolution at coarser scales than ordinary wavelet transform. Below we provide a brief summary of the non-decimated wavelet transform. Books by Chui (1992), Daubechies (1992) and
ARTICLE IN PRESS F. Borchi et al. / Journal of Atmospheric and Solar-Terrestrial Physics 68 (2006) 1383–1395
Mallat (1989) may be consulted for a more in-depth discussion about wavelets. For a given location during a given year, the transformed ozone time series is denoted as O ¼ ðO ðt1 Þ; . . . ; O ðtn ÞÞ with ti ¼ i=n, i ¼ 1; . . . ; n, and n ¼ 2J . The Discrete Wavelet Transform (DWT) of the vector O* is simply a matrix product d ¼
W O , where d is an n 1 vector of discrete wavelet coefficients indexed by two integers, djk with j ¼ 1; . . . ; J, and W is an orthogonal n n matrix associated with the wavelet basis. In particular, the coefficient djk provides local information about the behaviour of the time series O*(t) concentrating on effects of scale around 2j and near time k 2j.
Multiresolution wavelet analysis 0.2
d1
0.1
0.0
-0.1
0.10
d2
0.05
0.00
-0.05
-0.10 0.06 0.04
d3
0.02 0.00 -0.02 -0.04 -0.06 0.02
d4
0.01 0.00 -0.01 -0.02 -0.03 50
100
1387
150 200 Time in days (1990)
250
300
Fig. 2. Multi-resolution analysis of the transformed ozone time series plotted in the right-up corner of Fig. 1.
ARTICLE IN PRESS 1388
F. Borchi et al. / Journal of Atmospheric and Solar-Terrestrial Physics 68 (2006) 1383–1395
This property explains why wavelets are said to be localized in time. One particular problem with DWT is that, unlike the discrete Fourier transform, it is not translation invariant. In addition, the DWT has poor time resolution as j becomes larger. A finer sampling rate at all levels is needed to detect variance changes at all levels. The non-decimated wavelet transform (NWT) solves these problems. It is defined as the set of all DWT’s formed from the n possible shifts of the data by amounts i/n, i ¼ 1; . . . ; n. Thus, unlike the DWT that has 2J–j+1 coefficients on the j resolution level, there are n equally spaced wavelet coefficients in the NWT for each resolution level. We denote these coefficients djt with tA{t1,y,tn}. Consequently the NWT becomes translation invariant. Due to its finer sampling rate at all levels; the NWT provides a better exploratory tool for analyzing changes in the scale behaviour of the signal under study. The advantages of the NWT over the DWT in time series analysis are demonstrated in Nason et al. (2000). To illustrate how the NWT decomposes signals, we plot in Fig. 2 the first four NWT resolution levels in the time domain of the transformed TOMS ozone time series in 1990 at a specific location (68.51W and 20.6251E). In Fig. 2, we see that the first level captures the most rapid changes from the original time series. The second, third and fourth levels represent increasing larger scales. The first decomposition level is particularly of interest to identify breakpoint timings. In comparison, the level 4 does not bring any information about small scale variabilities. The vertical dashed line in Fig. 2 corresponds to our estimated breakpoint date (the question of how this date is found will be addressed in the next section). This capacity of decomposing at different independent (orthogonal) scales also explains our a posteriori choice of working with a wavelet analysis. Finally, we would like to give some general remarks about the mathematical aspects of our multi-resolution analysis. A very few missing values in our time series were replaced by using a robust regression. Having regularly spaced observations with no missing values greatly simplifies the implementation of wavelet schemes without changing the overall quality of our statistical analysis. Although there is no physical reason that the length of our ozone time series n is a power of 2, we will assume n ¼ 2J . This greatly facilitates the mathematical exposition of our wavelet analysis. In practice, all wavelet procedures can be adapted by padding the data in order to work with the
necessary diadic length. The ‘‘Daubechies 4’’ wavelets were used in all our computations (Daubechies, 1992). 3. Variance changes in ozone time series 3.1. Locating the changes in variability in total ozone column The key component for detecting breakpoint timings is the normalized cumulative sum of squares test statistic (Brown et al., 1975; Hsu, 1977; Inclan and Tiao, 1994; Gencay et al., 2002). To define such a statistic, we first need to introduce the following ratio for each level of resolution j: 0 1,0 1 k n1 X X PðjÞ ¼ @ (2) d2 A @ d 2 A, j;t
k
t¼Lj 1
j;t
t¼Lj 1
for k ¼ Lj 1; . . . ; n 2 and where Lj is the length of the wavelet filter. The quantity PkðjÞ is a ratio of two sums of squares, the numerator being the sum before the instant k and the denominator corresponding to the total sum over all the time periods. If the wavelet coefficients of the time series under study have a constant variance at the jth level, one expects that PkðjÞ it will increase regularly. But, if a rapid change of variance occurs, then itPðjÞ k has to change accordingly. We use this property to define the test statistic D, the normalized cumulative sum of squares: k Lj þ 2 ðjÞ ðjÞ D ¼ max max Pk ; k N Lj k Lj þ 1 max PðjÞ . ð3Þ k k N 1 Whitcher et al. (2002) showed that D(j) can accurately locate a change of variance at level j. The value of D(j) improves as the ratio of the wavelet variance at level j before and after the breakpoint moves away from unity. Hence, the estimated location of the breakpoint timing is the date at which D(j) is achieved. From Fig. 2, it was clear that the first level provides the finest sampling rate and the variance change is the most striking at this resolution. Consequently the first resolution level is chosen to be the one from which the breakpoint timing is defined, i.e. it is the date at which D(1) is achieved. To illustrate this detection procedure, the temporal evolution of the statistic D(j) is displayed for levels 1, 2 and 3 in Fig. 3, the solid line
ARTICLE IN PRESS F. Borchi et al. / Journal of Atmospheric and Solar-Terrestrial Physics 68 (2006) 1383–1395
corresponds to the year 1990 and the dashed line to 1987. For the year 1990, the time series used to computed D(j) are the ones shown in Fig. 2. In the first panel (j ¼ 1) of Fig. 3, we observe the typical ‘‘triangle shape’’ associated with a strong change of variance, the solid line first increases, reaches a maximum and then steadily decreases. The breakpoint is simply defined as the date of the peak, i.e. the vertical solid line for 1990 and the vertical dashed line for 1987. In 1990, the breakpoint was
1389
well localized in time at a particular location (68.51W and 20.6251E), the statistic D(j) for j ¼ 2, 3 behaves similarly to level 1. In particular, the peaks at all three levels occur around the same date (vertical solid line). In comparison, the year 1987 corresponds to a very different situation. First, the intensity of D(j) is much weaker than in 1990, specially for j ¼ 3. Second, peaks for 1987 are not clearly defined; there are certainly two different peaks for j ¼ 2. Third, the typical ‘‘triangle’’ shape
Variability for level 1
0.4
0.3
0.2
0.1
0.0 0
100
200 Time in days
300
0
100
200 Time in days
300
0
100
200 Time in days
300
Variability for level 2
0.4
0.3
0.2
0.1
0.0
Variability for level 3
0.4
0.3
0.2
0.1
0.0
Fig. 3. Comparison between the statistics D defined by Eq. (2) and computed for each of the first three resolution levels of the transformed ozone time series. The solid line and dotted line corresponds to the years 1990 and 1987, respectively. The chosen location for both years is 68.51N and longitude 20.6251E. The top, middle and lower panels, corresponds to the first, second and third resolution levels, respectively.
ARTICLE IN PRESS F. Borchi et al. / Journal of Atmospheric and Solar-Terrestrial Physics 68 (2006) 1383–1395
1390
of D(j) observed in 1990 is not as well defined in 1987, specially for j ¼ 1, 3. All these factors indicate that the total ozone column variability in 1987 was not very well localized in time at this specific location. As mentioned before, this begs the question of how strong D(1) has to be in order to keep or reject a breakpoint date obtained from this statistical procedure. To answer this question, we use the following approximation for large samples (Inclan and Tiao, 1994; Whitcher et al., 2002): ðjÞ
1 p ¼ PðD 4xp Þ2
1 X
In this case, all of these 5 years have a fairly flat and uniform density distribution and no strong mode can be found. To confirm this clustering, we propose another way of presenting the two groups of years by showing the scatter plot between the inter-quartile range (IQR) of the significant breakpoint timings and the following statistic for each year: X
S1 ¼
½ðDð1Þ xp ÞIðDð1Þ 4xp Þ
all locations l
ð1Þ expðl
2
N j x2p Þ,
(4)
X
l¼1
where N j ¼ ½n2j has to be a large number, i.e. for the first two resolution levels. This approximation allows us to compute the 95% significant level, i.e. xp with p ¼ 0:05 in Eq. (3). One may ask how the statistic D(j) varies in space for a given year. In Fig. 4, the x-axis represents the dates at which D(j) is achieved for all locations and the y-axis corresponds to the value of D(j). As before, we focus our attention on the two years 1987 (left panels) and 1990 (right panels). The statistic D(j) was computed for the three levels (j ¼ 1 for the top panels, j ¼ 2 for the middle panels and j ¼ 3 for the bottom panels). In addition, the horizontal dashed lines on the first two levels indicate the 5% critical levels for j ¼ 1, 2, the points above these thresholds can be considered significant breakpoint timings. Comparison of 1987–1990 shows a strong difference in the intensity of D(j). In 1990, almost all locations are above the critical levels for j ¼ 1, 2, i.e. a significant change in variance occurs for most of the locations. The opposite situation occurs in 1987; very few locations can be associated with significant breakpoints. It is natural to wonder if other years behave like 1987, 1990 or differently. We address this question in the next section. 3.2. Classification of year We first computed D(1) for all 12 years and for all locations. Then, a density distribution of the dates at which D(1) is achieved was estimated and plotted in Fig. 5 for each year. From this graph, we identify a first group of seven years (1981, 1982, 1985, 1986, 1988, 1989 and 1990) that have a common strong modal density distribution and break-point dates are displayed in Table 1. A second group of years (1983, 1984, 1987, 1991 and 1992) is also identified.
ð1Þ
IðD 4xp Þ
! ð5Þ
all locations
where xp is found by using Eq. (3) and I(A) represents the indicator function, equal to 1 when A is true and 0 otherwise. In Fig. 6, we display the variable S1 defined by Eq. (5) versus the IQR of the breakpoint timing. We recall that the statistic IQR is the difference between the second quartile (the first 25% cutting point) and the third quartile (75%). This statistic measures the variability among breakpoint timings for each year. Fig. 6 provides an even clearer dichotomy between the two groups by showing that the variability through IQR is weaker for ‘‘modal’’ years than for ‘‘flat’’ years and is also linked to the number and intensity of significant locations via the yearly statistic S1. This bi-modal behaviour is important. It means that the total ozone column should be examined according to the type of year under study; otherwise the variability changes will not be correctly taken into account in inter-annual ozone variability studies. The classification of stratospheric variations in two groups is not new, and some authors (Labitzke, 1982) have already reported two preferred internal modes or type of circulation. One corresponds to a strong well-isolated vortex with cold polar temperatures. In contrast, the other mode exhibits a perturbed vortex with warmer mean temperatures induced by successive major stratospheric warmings that occur throughout the winter. It is important to note that we obtain via our statistical method, exactly the same two clusters of years than Labitzke (1982) who took advantage of dynamical parameters and their horizontal patterns. So, our approach, which is dedicated to ozone variance changes, not only confirms past results, but it also indicates a potential link between ozone variability and vortex persistence.
ARTICLE IN PRESS F. Borchi et al. / Journal of Atmospheric and Solar-Terrestrial Physics 68 (2006) 1383–1395 1990 Level 1 (all locations)
0.5
0.5
0.4
0.4
Variability of level 1
Variability of level 1
1987 Level 1 (all locations)
0.3 0.2 0.1 0.0
0.3 0.2 0.1 0.0
80
100 120 140 160 180 Breakpoint timing in Days
200
80
100 120 140 160 180 Breakpoint timing in Days
200
1990 Level 2
0.5
0.5
0.4
0.4
Variability of level 2
Variability of level 2
1987 Level 2
0.3 0.2 0.1 0.0
0.3 0.2 0.1 0.0
80
100 120 140 160 180 Breakpoint timing in Days
200
80
1987 Level 3
100 120 140 160 180 Breakpoint timing in Days
200
1990 Level 3
0.5
0.5
0.4
0.4
Variability of level 3
Variability of level 3
1391
0.3 0.2 0.1 0.0
0.3 0.2 0.1 0.0
80
100 120 140 160 180 Breakpoint timing in Days
200
80
100 120 140 160 180 Breakpoint timing in Days
200
Fig. 4. Amplitude of the statistic D for level 1 (y-axis) and its timing (x-axis) for all the locations. The right column corresponds to the year 1990, and the left one to 1987. The rows correspond to the first, second and third resolution levels. The dotted horizontal line indicates the 95% significant level. The 5% critical level for level 3 is above 0.5.
In addition, seven density peaks for the years associated with a clear ozone variance change can be seen in Fig. 5 and their dates of break-points are displayed in Table 1. Those timings occur, on average, 1 month later than the vortex breaking dates given by Waugh et al. (1999). This indicates another potential link between ozone variability and
vortex persistence. Thus, ozone variability in the lower stratosphere is principally due to the combined effect of air mass advection and large horizontal gradients. In this respect, the arctic vortex formation and its subsidence lead to much larger horizontal gradients during winter than in the summer period. When the vortex breaks in few
ARTICLE IN PRESS F. Borchi et al. / Journal of Atmospheric and Solar-Terrestrial Physics 68 (2006) 1383–1395
1392
1981
1982 Density level 1
0.015
0.000
0.015
0.000 80 100
140 Days
180
1984
80 100
0.015
0.015
180
1987
140 Days
180
80 100
1988
0.015
0.000
0.000 140 Days
180
1990
140 Days
180
80 100
1991
0.000
0.015
180
0.015
0.000
0.000 140 Days
180
0.030 Density level 1
Density level 1
0.015
140 Days 1992
0.030
80 100
0.015
0.000 80 100
0.030
180
1989 Density level 1
Density level 1
0.015
140 Days
0.030
0.030
80 100
0.015
0.000 80 100
0.030
180
1986
0.000 140 Days
140 Days
0.030 Density level 1
Density level 1
Density level 1
180
1985
0.000
Density level 1
140 Days
0.030
80 100
0.015
0.000 80 100
0.030
Density level 1
1983 0.030
0.030 Density level 1
Density level 1
0.030
80 100
140 Days
180
80 100
140 Days
180
Fig. 5. Yearly evolution of the breakpoint timing density distribution from 1981 to 1992.
Table 1 Break-point dates for the group constituted of the years for which there exists a very abrupt and temporally localized ozone variability change Year
Break-points date
1981 1982 1985 1986 1988 1989 1990
19th May 7th April 18th May 8th April 23rd April 11th April 30th April
days, ozone may not have enough time to instantaneously mix with surrounding air and to reach a chemical equilibrium with O2. In this situation, some dynamical and homogeneous structures can persist several weeks (Marchand et al., 2003). Our study indicates a delay of about 1 month. However, this is a preliminary investigation and a detailed analysis of this phenomenon is beyond the scope of this paper. It is also interesting to notice that the years 1983, 1984, 1987 and 1992 (years with fairly flat and uniform density distribution and no strong mode)
ARTICLE IN PRESS F. Borchi et al. / Journal of Atmospheric and Solar-Terrestrial Physics 68 (2006) 1383–1395
84 70
92
Breakpoint timing IQR
60
50
83
87 91
40
85 30
88 86
81
82 89 20
90 0.03
0.04
0.05
0.06 0.07 S1
0.08
0.09
1393
Although these periods do not correspond to a classical season, we label them generically ‘‘winter’’ (before the breakpoint) and ‘‘spring’’ (after the breakpoint). Then we compute the winter and spring standard deviations for each location. The results of such computations are presented by the box plots (first quartile, median and third quartile) in Fig. 7. Only the breakpoint timings with a significant level of significance were incorporated in the box plot computations. The grey box plots correspond to winter and the white ones to spring. Consequently, the winter box plots have to be above the spring ones, the variability during the winter is known to be much higher than in the spring. Note that the years for which there is an overlap between winter and spring box plots (1987, 1991, 1992) belong to the first group of our classification derived in Section 3.2. This result confirms that the potential influence of the vortex breakdown on ozone for these years is not very well localized in time. Therefore, the distinction between the two seasons is more difficult to make in term of temporal ozone variability.
0.10
Fig. 6. Classification between the two groups of year; x-axis: the statistic S1 defined by Eq. (5), y-axis: the inter-quartile range (IQR) of the significant breakpoint timings over the 12 years. Note the IQR is defined as the difference between the first quartile (25%) and the third quartile (75%).
are related to explosive volcanic activity. The atmospheric consequences of El-Chichon (March– April 1982, Mexico) were mostly seen during 1983 and 1984. The signal from Nevado del Ruiz (November 1985, Colombia) may have had a climate influence in 1987. The well-known Pinatubo eruption (June 1991, Philippines) and the Cerro Hudson eruption (August 1991, Chile) had a climate effect in 1992. It may be of interest to explore such potential links between our clustering and volcanic activity in future work. 3.3. Winter and spring variabilities Most of Section 3.1 was dedicated to the detection of significant breakpoint timings in total ozone column time series. In this section, we now use these breakpoint dates to better assess temporal ozone variability. To reach this goal, each time series with a significant breakpoint is divided into two parts, before and after the breakpoint.
4. Conclusions In order to improve our knowledge about ozone variability at small temporal scales, we applied an innovative multi-resolution statistical procedure to TOMS total ozone column over the period 1981–1992. We derived the yearly distribution of significant breakpoint timings in stratospheric polar ozone. From those timings, the well-known winter– spring transition for each year was clearly identified in our data. In particular, we were able to differentiate two types of years. The first one is characterized by very abrupt ozone changes in variability. In contrast, the second type presents slow ozone variance changes over time. This yearly dichotomy could be linked to the two types of vortex persistence behaviour: either a strong well isolated vortex with cold polar temperatures or a perturbed vortex with warmer mean temperatures induced by successive major stratospheric warmings that occurred throughout the winter. Our clustering into two classes of years is in complete agreement with previous dynamical studies on the same topic, thus validating our statistical methodology. Also, the presence of abrupt changes in the ozone time series clearly indicates that future research on its inter-annual variability should take into account,
ARTICLE IN PRESS 1394
F. Borchi et al. / Journal of Atmospheric and Solar-Terrestrial Physics 68 (2006) 1383–1395
Fig. 7. Yearly evolution of the standard deviations before and after the breakpoint timings. The grey box plots correspond to standard deviations computed before breakpoint timings, the white ones after breakpoint timings. Only the locations with significant breakpoint timings for a given year were used in the computation of the standard deviation. The boxes are drawn with widths proportional to the square-roots of the number of observations in the group. The circles correspond to outliers.
dynamically and statistically, those breakpoint timings. While our study shows that Arctic vortex breakdown timings occurred, on average, 1 month earlier than our detected changes in ozone variability, the mechanisms connecting ozone variability and vortex persistence in the NH are complex and not yet well understood. This means that more research is necessary to discriminate the respective factors (e.g. solar irradiance, vortex activities) that could explain NH temporal ozone variability. Acknowledgements We wish to thank members of the NASA TOMS Nimbus Experiment and Information Processing Teams for providing the version 7 TOMS ozone data. References Brown, R.G., Durbin, J., Evans, J.M., 1975. Techniques for testing the constancy of regression relationships over time. Journal of the Royal Statistical Society B 37, 149–163.
Butchard, N., Remsberg, E.E., 1986. The area of the stratospheric polar vortex as a diagnostic for trace transport on an isentropic surface. Journal of Atmospheric Science 43, 1319–1339. Chui, C., 1992. An Introduction to Wavelets. Academic Press, New York. Coy, L.E., Nash, E.R., Newman, P.A., 1997. Meteorology of the polar vortex, spring 1997. Geophysical Research Letters 24, 2693–2696. Daubechies, I., 1992. Ten Lectures on Wavelets. SIAM, Philadelphia. Echer, E., 2004. Multi-resolution analysis of global total ozone column during 1979–1992 Nimbus-7 TOMS period. Annales Geophysicae 22, 1487–1493. Fusco, A., Salby, M.L., 1999. Interannual variations of total ozone and their relationship to variations of planetary wave activity. Journal of Climate 12, 1619–1620. Gencay, R., Selcuk, F., Whitcher, B., 2002. An Introduction to Wavelets and Other Filtering Methods in Finance and Economics. Academic Press, San Diego. Hsu, D.-A., 1977. Tests for variance shift at an unknown time point. Applied Statistics 26, 279–284. Inclan, C., Tiao, G.C., 1994. Use of cumulative sums of squares for retrospective detection of changes of variance. Journal of the American Statistical Association 89, 913–923. Labitzke, K., 1982. On the interannual variability of the middle stratosphere during the northern winters. Journal of Meteorilogical Society of Japan 60, 124–139.
ARTICLE IN PRESS F. Borchi et al. / Journal of Atmospheric and Solar-Terrestrial Physics 68 (2006) 1383–1395 Mallat, S., 1989. A theory for multiresolution signal decomposition: the wavelet representation. IEEE Transactions on Pattern Analysis and Machanical Intelligence 11, 674–693. Marchand, M., Godin, S., Hauchecorne, A., Lefe`vre, F., Bekki, S., Chipperfield, M., 2003. Influence of polar ozone loss on northern midlatitude regions estimated by a high-resolution chemistry transport model during winter 1999–2000. Journal Geophysical Research 108 (D5), 8326. McPeters, R., 1996. Nimbus 7 total ozone mapping spectrometer (TOMS) data products user’s guide. NASA Reference Publ 1384, 67pp. Nason, G., von Sachs, R., Kroisandt, G., 2000. Wavelet processes and adaptive estimation of the evolutionary wavelet spectrum. Journal of Research Statistical Society B 62, 271–292. Randel, W.J., Wu, F., 1999. Cooling of the Arctic and Antarctic polar stratosphere due to ozone depletion. Journal of Climate 12, 1467–1479. Schnadt, C., Dameris, M., Ponater, M., Hein, R., Grewe, V., Steil, B., 2002. Interaction of atmospheric chemistry and climate and its impact on stratospheric ozone. Climate Dynamics 18, 501–517. Shindell, D.T., Rind, D., Lonergan, P., 1998. Increased polar stratospheric ozone losses and delayed eventual recovery
1395
owing to increasing greenhouse-gas concentrations. Nature 392, 589–592. Waugh, D.W., Randel, W.J., 1999. Climatology of Arctic and Antarctic polar vortices using elliptical diagnostics. Journal of Atmospheric Science 56, 1594–1613. Waugh, D.W., Randel, W.J., Pawson, S., Newman, P.A., Nash, E.R., 1999. Persistence of the lower stratospheric polar vortices. Journal of Geophysical Research 104 (D22), 27191–27201. Whitcher, B., Jensen, M.J., 2000. Wavelet estimation of a local long memory parameter. Exploration Geophysics 31, 94–103. Whitcher, B., Guttorp, P., Percival, D.B., 2000. Wavelet analysis of covariance with application to atmospheric time series. Journal of Geophysical Research 105 (D11), 14941–14962. Whitcher, B., Byers, S.D., Guttorp, P., Percival, D.B., 2002. Testing for homogeneity of variance in time series: long memory, wavelets and the Nile river. Water Resources Research 38 (5), 1029. Wmo, A., 2003. Scientific assessment of ozone depletion: 2002. Global Ozone Research and Monitoring Project Rep 47, 498pp. Zhou, S., Gelman, M.E., Miller, A., McCormack, J.P., 2000. An inter-hemisphere comparison of the persistent stratospheric polar vortex. Geophysical Research Letters 27 (8), 1123–1126.
Journal of Atmospheric and Solar-Terrestrial Physics 68 (2006) 1383–1395 www.elsevier.com/locate/jastp
Detecting variability changes in Arctic total ozone column Franc- ois Borchia,, P. Naveaub, P. Keckhuta, A. Hauchecornea a
Service d’Ae´ronomie, IPSL CNRS, Verrie`res le Buisson, France Laboratoire des Sciences du Climat et l’Environnement, IPSL CNRS, Saclay, France
b
Received 9 January 2006; received in revised form 17 April 2006; accepted 2 May 2006 Available online 7 July 2006
Abstract To better assess inter-annual variance changes in ozone time series, we study and propose a statistical procedure based on recent advances in wavelet multi-resolution analysis. This approach, novel to the field of ozone analysis, has the advantages of detecting significant changes in variance and of characterizing the distribution of these changes, including not only their timings but also their strengths. As a test case for our method, we study total ozone column time series from the TOMS (Total Ozone Mapping Spectrometer) Version 7 in the Arctic region over the period 1981–1992. This statistical procedure allows us to analyse the well-known winter–spring transition and the inter-annual variability behaviour. For this specific example, two distinct types of winter–spring transition can be clearly identified for the ozone series. As expected from dynamical studies of the vortex evolution, the first type is constituted of the years for which there exists a very abrupt and temporally localized ozone variability change. More precisely, we detect a significant peak in ozone variance at two multi-resolution wavelet levels. The second type is characterized by smooth ozone variance changes over time and space. Moreover, the Arctic vortex breakdown timing seems to occur, on average, 1 month earlier than our detected changes in ozone variability. These results suggest a potential link between these two phenomena. r 2006 Elsevier Ltd. All rights reserved. Keywords: Stratosphere; Total ozone column; Arctic vortex; Change-point models; Wavelet analysis
1. Introduction The mean change of stratospheric ozone has been subject to many studies (WMO, 2003), but few statistical investigations on its short-term variance have been published because detecting rapid temporal variance changes in stratospheric ozone represents a series of statistical challenges. Classical Corresponding author. Tel.: +33 1 64 47 43 47; fax: +33 1 69 20 29 99. E-mail address: [email protected] (F. Borchi).
1364-6826/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.jastp.2006.05.011
linear multi-regression techniques cannot be directly implemented for these stratospheric studies. This is due to the complexity of the different spatiotemporal scales contained in such stratospheric ozone observations and the non-stationary nature of these time series. Indeed, variability of polar stratospheric ozone depends on various internal (ozone chemistry, atmospheric circulation, temperature changes, etc.) and external (solar activity, volcano eruptions, etc.) processes and complex feedbacks (Schnadt et al., 2002). Particularly, polar ozone transport is strongly associated with
ARTICLE IN PRESS 1384
F. Borchi et al. / Journal of Atmospheric and Solar-Terrestrial Physics 68 (2006) 1383–1395
planetary waves activity and the strength and deformation of the vortex (Fusco and Salby, 1999; Marchand et al., 2003). Multi-resolution analysis is becoming more popular in atmospheric sciences. Whitcher and his collaborators wrote a series of articles on multiresolution analysis with some applications to geophysics (Whitcher et al., 2000; Whitcher and Jensen, 2000; Whitcher et al., 2002). With respect to the analysis of ozone time series, Echer (2004) recently performed a wavelet decomposition to total ozone column that was observed by the Nimbus-7 Total Ozone Mapping Spectrometer (TOMS) instrument from 1979 to 1992 and detected a series of different time scales: high frequency oscillations (2–4 months), semi-annual variation (6 months), annual variation, quasi-biennial oscillation (QBO), El Nino-Southern Oscillation (ENSO) and solar cycle. Although the scope of our work is closely related to this past study in terms of data and methods, our research differs in two fundamental ways. First, our study does not focus on decomposing total ozone column at different scales per se but instead our goal is to assess inter-annual variances changes to detect breakpoint timings in stratospheric polar ozone measured at different locations over 12 years. In other words, the multiresolution analysis helps us but it is only a first step. We propose and implement other statistical procedures to go much further than displaying, comparing and commenting wavelet decompositions. This remark leads to the second difference with Echer’s study. By taking advantage of test statistics, we are able to detect significant variance changes within the small temporal scales obtained from multi-resolution decomposition. To our knowledge, such a statistical detection procedure of breakpoint timings is novel with respect to ozone time series analysis and it provides a more complete description of ozone variation. Before closing this introduction, it is important to recall that many authors (Coy et al., 1997; Shindell et al., 1998; Randel and Wu, 1999; Waugh and Randel, 1999; Waugh et al., 1999; Zhou et al., 2000) have studied the springtime persistence of the Arctic vortex that have an impact on ozone depletion. Our objective is different here: we aim to detect the interannual variance changes in ozone time series, but not the vortex evolution. To clarify this point, we remind the reader that the polar vortex breaking can either be characterized by high potential vorticity(PV) with extremes in zonal winds U within the
region of steep PV gradients at the vortex edge (Butchard and Remsberg, 1986), or by having its size shrinks to less than 1% of the Earth surface area (Zhou et al., 2000). Alternatively, zonal wind U and PV can be used together. Waugh et al. (1999) compared these three diagnostics (U alone, PV alone, both) when determining the breakdown timing of both the Arctic and Antarctic vortices over the last 40 years. He concluded that the three methods could provide similar results about vortex persistence. In comparison, we will neither look at PV nor at wind data. Instead, we solely focus on the total ozone column time series from TOMS. Working with ozone, instead of PV or wind data, shed a new light on the analysis of breakpoint timings because they are different for total ozone column and the vortex. Hence, our work should be viewed as a complement to previous dynamical studies based on vortex persistence, which validates our methodology. This paper is organized as follows. In Section 2, we present the TOMS data set under study and the multi-resolution approach is briefly explained. In Section 3.1, we detail the statistical challenges associated with the detection of abrupt ozone changes. A wavelet-based analysis is undertaken to solve this detection problem. Section 3.2 presents a summary of breakpoint distributions for each year. It is followed by a binary classification for the 12 years. In Section 3.3, an inter-annual breakpoint evolution is performed to determine the differences between winter and spring ozone variabilities. We conclude in Section 4. 2. Statistical analysis 2.1. Description of the ozone time series To illustrate our statistical approach, we choose to analyse data sets from the TOMS version 7 (McPeters, 1996). This instrument provides a large coverage of ozone abundance measurements over a long continuous period with a good temporal sampling. This allows us to investigate the shortterm variability on an inter-annual basis. The data coverage during the Nimbus 7 TOMS lifetime spans from November 1, 1978 to May 6, 1993. The individual TOMS measurements have been averaged into grid cells covering 11 of latitude by 1.251 of longitude from 901S to 901N, from 1801W to 1801E. While total ozone column integrates different layers (each one having its own variability), it is
ARTICLE IN PRESS F. Borchi et al. / Journal of Atmospheric and Solar-Terrestrial Physics 68 (2006) 1383–1395
strongly associated with stratospheric levels around 20 km where ozone has a maximum. Total ozone column is inferred from the differential absorption of scattered sunlight in the ultraviolet range. TOMS directly measures the ultraviolet sunlight scattered by the Earth’s atmosphere. This method does not allow for obtaining ozone data during periods
without sunlight. This constraint implies that most of the winter season in polar regions are not analysed in our study. More precisely, we have selected a spatial grid of six latitudes equally spaced from 64.51N to 74.51N and of 18 longitudes from 1801W to 1801E. This grid selection allows us to investigate ozone variability at high latitudes. Log-diff Transform
450
0.2
400
0.1 log(o(t)/o(t-1))
Ozone
Long=20.625 Lat=68.5
350
0.0
300
-0.1
250
-0.2 0
1385
100 200 300 Time in days (1990)
0
100 200 300 Time in days (1990) Log-diff Transform
Long=20.625 Lat=68.5
0.2 500
0.1 log(o(t)/o(t-1))
Ozone
450
400
350
0.0
-0.1 300
-0.2
250 0
100 200 300 Time in days (1987)
0
100 200 300 Time in days (1987)
Fig. 1. TOMS daily ozone time series in 1990 (top row) and in 1987 (lower row) at 68.51N and longitude 20.6251E. The second column represents the transformed original time series after taking the logarithmic difference defined by Eq. (1). The dotted vertical line corresponds to the annual breaking point timing found by our wavelet analysis (see Figs. 2 and 3 for more the determination of such timings).
ARTICLE IN PRESS 1386
F. Borchi et al. / Journal of Atmospheric and Solar-Terrestrial Physics 68 (2006) 1383–1395
The upper-left panel of Fig. 1 represents the daily ozone evolution in 1990 at a specific location, latitude 68.51N and longitude 20.6251E, near Kiruna in Sweden. In this panel, it is clear that an abrupt ozone change in variability occurs around day 120, at the end of April (the vertical dashed line). This example corresponds to the typical evolution of ozone, a rapid increase associated with a strong variability during the winter followed by a slow decay with lower variability during the spring. Hence, we call such a variance change as a ‘‘breakpoint timing’’, i.e. it is the date when the strong winter variability in the ozone time series switches abruptly to the smaller spring variability at a given location and for a given year. But, such behaviour in ozone variance changes does not occur for each year. For example, the lower-left panel of Fig. 1 indicates that the variability changes in 1987 took place at a smoother pace than 1990. Before presenting the details of our statistical analysis, we would like to make a few comments. First, we would like to emphasize that our definition of breakpoint timings, neither corresponds to yearly ozone maxima (see the left panels of Fig. 1), nor to the ones obtained with PV area diagnostic (less than 1% of the Earth’s surface area) in Zhou et al. (2000). Second, our strategy is to associate a breakpoint date for every particular location of our spatial grid, whenever such an association is possible. Consequently, this approach allows us to estimate the distribution of breakpoint timings for a given year. This probabilistic scheme could be seen as restrictive with respect to our main goal, we do not derive a unique date for a given year and an entire region. However due to the complexity of ozone variability in such region, we prefer to assume that each location is influenced by the complex phenomena that bring a strong stochastic component, i.e. spatially similar but different breaking dates. Under this assumption, a probabilistic approach is not only possible but necessary. It may not even be possible to associate breakpoint timings for locations with weak ozone variability for a given year. For example in 1987, we will see that the breakpoint dates for many locations are not always clearly temporally defined. This point evokes the important question of how to determine the level of confidence of all breakpoint timings. This issue will be carefully addressed in this paper. For each location we will test the significance of all breakpoint timings derived from our statistical procedure.
2.2. Removing the seasonal variation One possible option to detect rapid changes in daily ozone time series is to follow the next two steps. First, the seasonal variation would be removed, either by a parametric regression (e.g. polynomial fit) or by a filtering procedure (low pass filter). Second, the variability contained in the residuals would be studied. Although of interest, this two-step approach contains a few drawbacks. A parametric method to model the seasonality implies extra parameters and additional assumptions. A filtering procedure may discard some important information for a range of scales. To remove the seasonal variation, we instead apply the following transform to our daily time series: O ðtÞ ¼ log OðtÞ log Oðt 1Þ, for each year and location.
ð1Þ
In the above formula, O(t) represents the TOMS total ozone column for a given day, say t. The right panels of Fig. 1 show the result of such a logarithmic transformation. This transformation has been used extensively in finance (Gencay et al., 2002). From Fig. 1, one can see that the seasonal variation has been successfully removed. Besides its simplicity of implementation, this log-difference has the advantage of removing any smooth (i.e. polynomial) trend. Others procedures could have also been applied. For the remainder of this paper, we will work with the transformed data instead of the raw ozone time series. 2.3. Multi-resolution analysis Wavelets are fundamental building block functions, analogous to the trigonometric sine and cosine functions. The Fourier transform extracts details from the signal frequency, but all information about the location of a particular frequency within the signal is lost. In comparison, the multiresolution analysis known as the ‘‘zoom-in, zoomout’’ property makes wavelets particularly appealing for this study, because they are localized in time and the signals are examined using widely varying levels of focus. In this article, we work with a nondecimated wavelet transform, which gives a much better temporal resolution at coarser scales than ordinary wavelet transform. Below we provide a brief summary of the non-decimated wavelet transform. Books by Chui (1992), Daubechies (1992) and
ARTICLE IN PRESS F. Borchi et al. / Journal of Atmospheric and Solar-Terrestrial Physics 68 (2006) 1383–1395
Mallat (1989) may be consulted for a more in-depth discussion about wavelets. For a given location during a given year, the transformed ozone time series is denoted as O ¼ ðO ðt1 Þ; . . . ; O ðtn ÞÞ with ti ¼ i=n, i ¼ 1; . . . ; n, and n ¼ 2J . The Discrete Wavelet Transform (DWT) of the vector O* is simply a matrix product d ¼
W O , where d is an n 1 vector of discrete wavelet coefficients indexed by two integers, djk with j ¼ 1; . . . ; J, and W is an orthogonal n n matrix associated with the wavelet basis. In particular, the coefficient djk provides local information about the behaviour of the time series O*(t) concentrating on effects of scale around 2j and near time k 2j.
Multiresolution wavelet analysis 0.2
d1
0.1
0.0
-0.1
0.10
d2
0.05
0.00
-0.05
-0.10 0.06 0.04
d3
0.02 0.00 -0.02 -0.04 -0.06 0.02
d4
0.01 0.00 -0.01 -0.02 -0.03 50
100
1387
150 200 Time in days (1990)
250
300
Fig. 2. Multi-resolution analysis of the transformed ozone time series plotted in the right-up corner of Fig. 1.
ARTICLE IN PRESS 1388
F. Borchi et al. / Journal of Atmospheric and Solar-Terrestrial Physics 68 (2006) 1383–1395
This property explains why wavelets are said to be localized in time. One particular problem with DWT is that, unlike the discrete Fourier transform, it is not translation invariant. In addition, the DWT has poor time resolution as j becomes larger. A finer sampling rate at all levels is needed to detect variance changes at all levels. The non-decimated wavelet transform (NWT) solves these problems. It is defined as the set of all DWT’s formed from the n possible shifts of the data by amounts i/n, i ¼ 1; . . . ; n. Thus, unlike the DWT that has 2J–j+1 coefficients on the j resolution level, there are n equally spaced wavelet coefficients in the NWT for each resolution level. We denote these coefficients djt with tA{t1,y,tn}. Consequently the NWT becomes translation invariant. Due to its finer sampling rate at all levels; the NWT provides a better exploratory tool for analyzing changes in the scale behaviour of the signal under study. The advantages of the NWT over the DWT in time series analysis are demonstrated in Nason et al. (2000). To illustrate how the NWT decomposes signals, we plot in Fig. 2 the first four NWT resolution levels in the time domain of the transformed TOMS ozone time series in 1990 at a specific location (68.51W and 20.6251E). In Fig. 2, we see that the first level captures the most rapid changes from the original time series. The second, third and fourth levels represent increasing larger scales. The first decomposition level is particularly of interest to identify breakpoint timings. In comparison, the level 4 does not bring any information about small scale variabilities. The vertical dashed line in Fig. 2 corresponds to our estimated breakpoint date (the question of how this date is found will be addressed in the next section). This capacity of decomposing at different independent (orthogonal) scales also explains our a posteriori choice of working with a wavelet analysis. Finally, we would like to give some general remarks about the mathematical aspects of our multi-resolution analysis. A very few missing values in our time series were replaced by using a robust regression. Having regularly spaced observations with no missing values greatly simplifies the implementation of wavelet schemes without changing the overall quality of our statistical analysis. Although there is no physical reason that the length of our ozone time series n is a power of 2, we will assume n ¼ 2J . This greatly facilitates the mathematical exposition of our wavelet analysis. In practice, all wavelet procedures can be adapted by padding the data in order to work with the
necessary diadic length. The ‘‘Daubechies 4’’ wavelets were used in all our computations (Daubechies, 1992). 3. Variance changes in ozone time series 3.1. Locating the changes in variability in total ozone column The key component for detecting breakpoint timings is the normalized cumulative sum of squares test statistic (Brown et al., 1975; Hsu, 1977; Inclan and Tiao, 1994; Gencay et al., 2002). To define such a statistic, we first need to introduce the following ratio for each level of resolution j: 0 1,0 1 k n1 X X PðjÞ ¼ @ (2) d2 A @ d 2 A, j;t
k
t¼Lj 1
j;t
t¼Lj 1
for k ¼ Lj 1; . . . ; n 2 and where Lj is the length of the wavelet filter. The quantity PkðjÞ is a ratio of two sums of squares, the numerator being the sum before the instant k and the denominator corresponding to the total sum over all the time periods. If the wavelet coefficients of the time series under study have a constant variance at the jth level, one expects that PkðjÞ it will increase regularly. But, if a rapid change of variance occurs, then itPðjÞ k has to change accordingly. We use this property to define the test statistic D, the normalized cumulative sum of squares: k Lj þ 2 ðjÞ ðjÞ D ¼ max max Pk ; k N Lj k Lj þ 1 max PðjÞ . ð3Þ k k N 1 Whitcher et al. (2002) showed that D(j) can accurately locate a change of variance at level j. The value of D(j) improves as the ratio of the wavelet variance at level j before and after the breakpoint moves away from unity. Hence, the estimated location of the breakpoint timing is the date at which D(j) is achieved. From Fig. 2, it was clear that the first level provides the finest sampling rate and the variance change is the most striking at this resolution. Consequently the first resolution level is chosen to be the one from which the breakpoint timing is defined, i.e. it is the date at which D(1) is achieved. To illustrate this detection procedure, the temporal evolution of the statistic D(j) is displayed for levels 1, 2 and 3 in Fig. 3, the solid line
ARTICLE IN PRESS F. Borchi et al. / Journal of Atmospheric and Solar-Terrestrial Physics 68 (2006) 1383–1395
corresponds to the year 1990 and the dashed line to 1987. For the year 1990, the time series used to computed D(j) are the ones shown in Fig. 2. In the first panel (j ¼ 1) of Fig. 3, we observe the typical ‘‘triangle shape’’ associated with a strong change of variance, the solid line first increases, reaches a maximum and then steadily decreases. The breakpoint is simply defined as the date of the peak, i.e. the vertical solid line for 1990 and the vertical dashed line for 1987. In 1990, the breakpoint was
1389
well localized in time at a particular location (68.51W and 20.6251E), the statistic D(j) for j ¼ 2, 3 behaves similarly to level 1. In particular, the peaks at all three levels occur around the same date (vertical solid line). In comparison, the year 1987 corresponds to a very different situation. First, the intensity of D(j) is much weaker than in 1990, specially for j ¼ 3. Second, peaks for 1987 are not clearly defined; there are certainly two different peaks for j ¼ 2. Third, the typical ‘‘triangle’’ shape
Variability for level 1
0.4
0.3
0.2
0.1
0.0 0
100
200 Time in days
300
0
100
200 Time in days
300
0
100
200 Time in days
300
Variability for level 2
0.4
0.3
0.2
0.1
0.0
Variability for level 3
0.4
0.3
0.2
0.1
0.0
Fig. 3. Comparison between the statistics D defined by Eq. (2) and computed for each of the first three resolution levels of the transformed ozone time series. The solid line and dotted line corresponds to the years 1990 and 1987, respectively. The chosen location for both years is 68.51N and longitude 20.6251E. The top, middle and lower panels, corresponds to the first, second and third resolution levels, respectively.
ARTICLE IN PRESS F. Borchi et al. / Journal of Atmospheric and Solar-Terrestrial Physics 68 (2006) 1383–1395
1390
of D(j) observed in 1990 is not as well defined in 1987, specially for j ¼ 1, 3. All these factors indicate that the total ozone column variability in 1987 was not very well localized in time at this specific location. As mentioned before, this begs the question of how strong D(1) has to be in order to keep or reject a breakpoint date obtained from this statistical procedure. To answer this question, we use the following approximation for large samples (Inclan and Tiao, 1994; Whitcher et al., 2002): ðjÞ
1 p ¼ PðD 4xp Þ2
1 X
In this case, all of these 5 years have a fairly flat and uniform density distribution and no strong mode can be found. To confirm this clustering, we propose another way of presenting the two groups of years by showing the scatter plot between the inter-quartile range (IQR) of the significant breakpoint timings and the following statistic for each year: X
S1 ¼
½ðDð1Þ xp ÞIðDð1Þ 4xp Þ
all locations l
ð1Þ expðl
2
N j x2p Þ,
(4)
X
l¼1
where N j ¼ ½n2j has to be a large number, i.e. for the first two resolution levels. This approximation allows us to compute the 95% significant level, i.e. xp with p ¼ 0:05 in Eq. (3). One may ask how the statistic D(j) varies in space for a given year. In Fig. 4, the x-axis represents the dates at which D(j) is achieved for all locations and the y-axis corresponds to the value of D(j). As before, we focus our attention on the two years 1987 (left panels) and 1990 (right panels). The statistic D(j) was computed for the three levels (j ¼ 1 for the top panels, j ¼ 2 for the middle panels and j ¼ 3 for the bottom panels). In addition, the horizontal dashed lines on the first two levels indicate the 5% critical levels for j ¼ 1, 2, the points above these thresholds can be considered significant breakpoint timings. Comparison of 1987–1990 shows a strong difference in the intensity of D(j). In 1990, almost all locations are above the critical levels for j ¼ 1, 2, i.e. a significant change in variance occurs for most of the locations. The opposite situation occurs in 1987; very few locations can be associated with significant breakpoints. It is natural to wonder if other years behave like 1987, 1990 or differently. We address this question in the next section. 3.2. Classification of year We first computed D(1) for all 12 years and for all locations. Then, a density distribution of the dates at which D(1) is achieved was estimated and plotted in Fig. 5 for each year. From this graph, we identify a first group of seven years (1981, 1982, 1985, 1986, 1988, 1989 and 1990) that have a common strong modal density distribution and break-point dates are displayed in Table 1. A second group of years (1983, 1984, 1987, 1991 and 1992) is also identified.
ð1Þ
IðD 4xp Þ
! ð5Þ
all locations
where xp is found by using Eq. (3) and I(A) represents the indicator function, equal to 1 when A is true and 0 otherwise. In Fig. 6, we display the variable S1 defined by Eq. (5) versus the IQR of the breakpoint timing. We recall that the statistic IQR is the difference between the second quartile (the first 25% cutting point) and the third quartile (75%). This statistic measures the variability among breakpoint timings for each year. Fig. 6 provides an even clearer dichotomy between the two groups by showing that the variability through IQR is weaker for ‘‘modal’’ years than for ‘‘flat’’ years and is also linked to the number and intensity of significant locations via the yearly statistic S1. This bi-modal behaviour is important. It means that the total ozone column should be examined according to the type of year under study; otherwise the variability changes will not be correctly taken into account in inter-annual ozone variability studies. The classification of stratospheric variations in two groups is not new, and some authors (Labitzke, 1982) have already reported two preferred internal modes or type of circulation. One corresponds to a strong well-isolated vortex with cold polar temperatures. In contrast, the other mode exhibits a perturbed vortex with warmer mean temperatures induced by successive major stratospheric warmings that occur throughout the winter. It is important to note that we obtain via our statistical method, exactly the same two clusters of years than Labitzke (1982) who took advantage of dynamical parameters and their horizontal patterns. So, our approach, which is dedicated to ozone variance changes, not only confirms past results, but it also indicates a potential link between ozone variability and vortex persistence.
ARTICLE IN PRESS F. Borchi et al. / Journal of Atmospheric and Solar-Terrestrial Physics 68 (2006) 1383–1395 1990 Level 1 (all locations)
0.5
0.5
0.4
0.4
Variability of level 1
Variability of level 1
1987 Level 1 (all locations)
0.3 0.2 0.1 0.0
0.3 0.2 0.1 0.0
80
100 120 140 160 180 Breakpoint timing in Days
200
80
100 120 140 160 180 Breakpoint timing in Days
200
1990 Level 2
0.5
0.5
0.4
0.4
Variability of level 2
Variability of level 2
1987 Level 2
0.3 0.2 0.1 0.0
0.3 0.2 0.1 0.0
80
100 120 140 160 180 Breakpoint timing in Days
200
80
1987 Level 3
100 120 140 160 180 Breakpoint timing in Days
200
1990 Level 3
0.5
0.5
0.4
0.4
Variability of level 3
Variability of level 3
1391
0.3 0.2 0.1 0.0
0.3 0.2 0.1 0.0
80
100 120 140 160 180 Breakpoint timing in Days
200
80
100 120 140 160 180 Breakpoint timing in Days
200
Fig. 4. Amplitude of the statistic D for level 1 (y-axis) and its timing (x-axis) for all the locations. The right column corresponds to the year 1990, and the left one to 1987. The rows correspond to the first, second and third resolution levels. The dotted horizontal line indicates the 95% significant level. The 5% critical level for level 3 is above 0.5.
In addition, seven density peaks for the years associated with a clear ozone variance change can be seen in Fig. 5 and their dates of break-points are displayed in Table 1. Those timings occur, on average, 1 month later than the vortex breaking dates given by Waugh et al. (1999). This indicates another potential link between ozone variability and
vortex persistence. Thus, ozone variability in the lower stratosphere is principally due to the combined effect of air mass advection and large horizontal gradients. In this respect, the arctic vortex formation and its subsidence lead to much larger horizontal gradients during winter than in the summer period. When the vortex breaks in few
ARTICLE IN PRESS F. Borchi et al. / Journal of Atmospheric and Solar-Terrestrial Physics 68 (2006) 1383–1395
1392
1981
1982 Density level 1
0.015
0.000
0.015
0.000 80 100
140 Days
180
1984
80 100
0.015
0.015
180
1987
140 Days
180
80 100
1988
0.015
0.000
0.000 140 Days
180
1990
140 Days
180
80 100
1991
0.000
0.015
180
0.015
0.000
0.000 140 Days
180
0.030 Density level 1
Density level 1
0.015
140 Days 1992
0.030
80 100
0.015
0.000 80 100
0.030
180
1989 Density level 1
Density level 1
0.015
140 Days
0.030
0.030
80 100
0.015
0.000 80 100
0.030
180
1986
0.000 140 Days
140 Days
0.030 Density level 1
Density level 1
Density level 1
180
1985
0.000
Density level 1
140 Days
0.030
80 100
0.015
0.000 80 100
0.030
Density level 1
1983 0.030
0.030 Density level 1
Density level 1
0.030
80 100
140 Days
180
80 100
140 Days
180
Fig. 5. Yearly evolution of the breakpoint timing density distribution from 1981 to 1992.
Table 1 Break-point dates for the group constituted of the years for which there exists a very abrupt and temporally localized ozone variability change Year
Break-points date
1981 1982 1985 1986 1988 1989 1990
19th May 7th April 18th May 8th April 23rd April 11th April 30th April
days, ozone may not have enough time to instantaneously mix with surrounding air and to reach a chemical equilibrium with O2. In this situation, some dynamical and homogeneous structures can persist several weeks (Marchand et al., 2003). Our study indicates a delay of about 1 month. However, this is a preliminary investigation and a detailed analysis of this phenomenon is beyond the scope of this paper. It is also interesting to notice that the years 1983, 1984, 1987 and 1992 (years with fairly flat and uniform density distribution and no strong mode)
ARTICLE IN PRESS F. Borchi et al. / Journal of Atmospheric and Solar-Terrestrial Physics 68 (2006) 1383–1395
84 70
92
Breakpoint timing IQR
60
50
83
87 91
40
85 30
88 86
81
82 89 20
90 0.03
0.04
0.05
0.06 0.07 S1
0.08
0.09
1393
Although these periods do not correspond to a classical season, we label them generically ‘‘winter’’ (before the breakpoint) and ‘‘spring’’ (after the breakpoint). Then we compute the winter and spring standard deviations for each location. The results of such computations are presented by the box plots (first quartile, median and third quartile) in Fig. 7. Only the breakpoint timings with a significant level of significance were incorporated in the box plot computations. The grey box plots correspond to winter and the white ones to spring. Consequently, the winter box plots have to be above the spring ones, the variability during the winter is known to be much higher than in the spring. Note that the years for which there is an overlap between winter and spring box plots (1987, 1991, 1992) belong to the first group of our classification derived in Section 3.2. This result confirms that the potential influence of the vortex breakdown on ozone for these years is not very well localized in time. Therefore, the distinction between the two seasons is more difficult to make in term of temporal ozone variability.
0.10
Fig. 6. Classification between the two groups of year; x-axis: the statistic S1 defined by Eq. (5), y-axis: the inter-quartile range (IQR) of the significant breakpoint timings over the 12 years. Note the IQR is defined as the difference between the first quartile (25%) and the third quartile (75%).
are related to explosive volcanic activity. The atmospheric consequences of El-Chichon (March– April 1982, Mexico) were mostly seen during 1983 and 1984. The signal from Nevado del Ruiz (November 1985, Colombia) may have had a climate influence in 1987. The well-known Pinatubo eruption (June 1991, Philippines) and the Cerro Hudson eruption (August 1991, Chile) had a climate effect in 1992. It may be of interest to explore such potential links between our clustering and volcanic activity in future work. 3.3. Winter and spring variabilities Most of Section 3.1 was dedicated to the detection of significant breakpoint timings in total ozone column time series. In this section, we now use these breakpoint dates to better assess temporal ozone variability. To reach this goal, each time series with a significant breakpoint is divided into two parts, before and after the breakpoint.
4. Conclusions In order to improve our knowledge about ozone variability at small temporal scales, we applied an innovative multi-resolution statistical procedure to TOMS total ozone column over the period 1981–1992. We derived the yearly distribution of significant breakpoint timings in stratospheric polar ozone. From those timings, the well-known winter– spring transition for each year was clearly identified in our data. In particular, we were able to differentiate two types of years. The first one is characterized by very abrupt ozone changes in variability. In contrast, the second type presents slow ozone variance changes over time. This yearly dichotomy could be linked to the two types of vortex persistence behaviour: either a strong well isolated vortex with cold polar temperatures or a perturbed vortex with warmer mean temperatures induced by successive major stratospheric warmings that occurred throughout the winter. Our clustering into two classes of years is in complete agreement with previous dynamical studies on the same topic, thus validating our statistical methodology. Also, the presence of abrupt changes in the ozone time series clearly indicates that future research on its inter-annual variability should take into account,
ARTICLE IN PRESS 1394
F. Borchi et al. / Journal of Atmospheric and Solar-Terrestrial Physics 68 (2006) 1383–1395
Fig. 7. Yearly evolution of the standard deviations before and after the breakpoint timings. The grey box plots correspond to standard deviations computed before breakpoint timings, the white ones after breakpoint timings. Only the locations with significant breakpoint timings for a given year were used in the computation of the standard deviation. The boxes are drawn with widths proportional to the square-roots of the number of observations in the group. The circles correspond to outliers.
dynamically and statistically, those breakpoint timings. While our study shows that Arctic vortex breakdown timings occurred, on average, 1 month earlier than our detected changes in ozone variability, the mechanisms connecting ozone variability and vortex persistence in the NH are complex and not yet well understood. This means that more research is necessary to discriminate the respective factors (e.g. solar irradiance, vortex activities) that could explain NH temporal ozone variability. Acknowledgements We wish to thank members of the NASA TOMS Nimbus Experiment and Information Processing Teams for providing the version 7 TOMS ozone data. References Brown, R.G., Durbin, J., Evans, J.M., 1975. Techniques for testing the constancy of regression relationships over time. Journal of the Royal Statistical Society B 37, 149–163.
Butchard, N., Remsberg, E.E., 1986. The area of the stratospheric polar vortex as a diagnostic for trace transport on an isentropic surface. Journal of Atmospheric Science 43, 1319–1339. Chui, C., 1992. An Introduction to Wavelets. Academic Press, New York. Coy, L.E., Nash, E.R., Newman, P.A., 1997. Meteorology of the polar vortex, spring 1997. Geophysical Research Letters 24, 2693–2696. Daubechies, I., 1992. Ten Lectures on Wavelets. SIAM, Philadelphia. Echer, E., 2004. Multi-resolution analysis of global total ozone column during 1979–1992 Nimbus-7 TOMS period. Annales Geophysicae 22, 1487–1493. Fusco, A., Salby, M.L., 1999. Interannual variations of total ozone and their relationship to variations of planetary wave activity. Journal of Climate 12, 1619–1620. Gencay, R., Selcuk, F., Whitcher, B., 2002. An Introduction to Wavelets and Other Filtering Methods in Finance and Economics. Academic Press, San Diego. Hsu, D.-A., 1977. Tests for variance shift at an unknown time point. Applied Statistics 26, 279–284. Inclan, C., Tiao, G.C., 1994. Use of cumulative sums of squares for retrospective detection of changes of variance. Journal of the American Statistical Association 89, 913–923. Labitzke, K., 1982. On the interannual variability of the middle stratosphere during the northern winters. Journal of Meteorilogical Society of Japan 60, 124–139.
ARTICLE IN PRESS F. Borchi et al. / Journal of Atmospheric and Solar-Terrestrial Physics 68 (2006) 1383–1395 Mallat, S., 1989. A theory for multiresolution signal decomposition: the wavelet representation. IEEE Transactions on Pattern Analysis and Machanical Intelligence 11, 674–693. Marchand, M., Godin, S., Hauchecorne, A., Lefe`vre, F., Bekki, S., Chipperfield, M., 2003. Influence of polar ozone loss on northern midlatitude regions estimated by a high-resolution chemistry transport model during winter 1999–2000. Journal Geophysical Research 108 (D5), 8326. McPeters, R., 1996. Nimbus 7 total ozone mapping spectrometer (TOMS) data products user’s guide. NASA Reference Publ 1384, 67pp. Nason, G., von Sachs, R., Kroisandt, G., 2000. Wavelet processes and adaptive estimation of the evolutionary wavelet spectrum. Journal of Research Statistical Society B 62, 271–292. Randel, W.J., Wu, F., 1999. Cooling of the Arctic and Antarctic polar stratosphere due to ozone depletion. Journal of Climate 12, 1467–1479. Schnadt, C., Dameris, M., Ponater, M., Hein, R., Grewe, V., Steil, B., 2002. Interaction of atmospheric chemistry and climate and its impact on stratospheric ozone. Climate Dynamics 18, 501–517. Shindell, D.T., Rind, D., Lonergan, P., 1998. Increased polar stratospheric ozone losses and delayed eventual recovery
1395
owing to increasing greenhouse-gas concentrations. Nature 392, 589–592. Waugh, D.W., Randel, W.J., 1999. Climatology of Arctic and Antarctic polar vortices using elliptical diagnostics. Journal of Atmospheric Science 56, 1594–1613. Waugh, D.W., Randel, W.J., Pawson, S., Newman, P.A., Nash, E.R., 1999. Persistence of the lower stratospheric polar vortices. Journal of Geophysical Research 104 (D22), 27191–27201. Whitcher, B., Jensen, M.J., 2000. Wavelet estimation of a local long memory parameter. Exploration Geophysics 31, 94–103. Whitcher, B., Guttorp, P., Percival, D.B., 2000. Wavelet analysis of covariance with application to atmospheric time series. Journal of Geophysical Research 105 (D11), 14941–14962. Whitcher, B., Byers, S.D., Guttorp, P., Percival, D.B., 2002. Testing for homogeneity of variance in time series: long memory, wavelets and the Nile river. Water Resources Research 38 (5), 1029. Wmo, A., 2003. Scientific assessment of ozone depletion: 2002. Global Ozone Research and Monitoring Project Rep 47, 498pp. Zhou, S., Gelman, M.E., Miller, A., McCormack, J.P., 2000. An inter-hemisphere comparison of the persistent stratospheric polar vortex. Geophysical Research Letters 27 (8), 1123–1126.