Comparative Spectral Analysis and Correlation Properties of Observed and Simulated Total Column Ozone Records

We present a statistical analysis of total column ozone records obtained from satellite measurements and from two global climate chemistry models on global scale. Firstly, a spectral weight analysis is performed where the relative strength of semiannual, annual and quasi-biennial oscillations are determined with respect to the integrated power spectra. The comparison reveals some anomalies in the model representations at each spectral component. The tails of the spectra demonstrate that both models underestimate high frequency (daily) ozone variability, which might have a complex origin, since several dynamical processes affect short time changes of the ozone level at a given location. Secondly, detrended fluctuation analysis is exploited to analyze two-point correlations of anomaly time series. Both models reproduce the characteristic geographic dependence of correlation strength over the overlapping area with empirical observations (latitude band between 60°S and 60°N). The values of precise correlation exponents are hard to obtain over regions where quasi-biennial oscillations or other strong nonstationarities (ozone hole) are present. In spite of all the numerical difficulties, significant long range correlations are detected for total ozone over all geographic locations.


Introduction
The European project RECONCILE has comprehensively addressed several key questions in the context of polar ozone depletion, with the objective to quantify the rates of some of the most relevant, yet still uncertain physical and chemical processes [1].The project used an extended toolkit of laboratory experiments, two field missions in the Arctic winter 2009/10 by the research aircraft M55-Geophysica (12 flights of 57 hours), a Match ozonesonde campaign in the winter 2010/11, together with microphysical and chemical transport modeling.One of the main goals was to construct an improved scheme of ozone chemistry for global chemistry climate model (CCM).Since CCMs provide the only tool for long-time projections, their validation and the evaluation of performance are crucial aspects of predicting the time evolution of the stratospheric ozone layer on global scale [2][3][4].
CCM activities within RECONCILE are carried out with the LMDz-Reprobus model coupling the Reprobus stratospheric chemistry package with the LMDz general circulation model [5][6][7].The performance of the model was tested against observations and other CCM models [6,[8][9][10][11].The most extended comparison project was the CCM Validation Activity in two runs (CCMVal-1, CCMVal-2), where 18 leading global CCMs are investigated by several observationally-based performance metrics in order to quantify the ability of models in reproducing the dynamics of stratospheric ozone [12].Various metrics based on climatological mean values, seasonal cycle, variability amplitudes and patterns, etc. are evaluated and compared [12].Here we focus on the spectral properties and temporal correlations of daily total column ozone (TO) time series on the global scale.Empirical data are obtained from the database compiled and maintained by Bodeker Scientific (http://www.bodekerscientific.com)from satellite and ground based observations [13][14][15].Two CCMs are used for the comparison: the LMDz-Reprobus mentioned above, and a quite similar model built up at the Meteorological Research Institute, Tsukuba, Japan [16,17].The latter is chosen from the CCMVal-2 database because its setup is very close to LMDz-Reprobus (see the next section), however it is among the 6 models from the 18 producing dynamically the quasi-biennial oscillations (QBO) in the equatorial band [12].This is known to be a difficult task, because it requires a high spatial and temporal resolution to reproduce a realistic spectrum (temporal and spatial) of upward propagating waves in the tropics [16,18,19].The main question we pose is whether the presence of a QBO in a CCM has a global impact on the correlation properties.A similarly important point is that we intend to produce a reference analysis of two CCMs in order to detect any global effect of an improved ozone chemistry and microphysics to be developed from the findings of the RECONCILE project.
In Section 2, we provide a concise description of empirical and simulated datasets and the methodology of spectral and correlation analysis.Section 3 is devoted to the results of the statistical tests.Global maps of the spectral weight of semiannual, annual and QBO components of the time series are presented and compared for empirical and CCM data.Two-point correlations are investigated by detrended fluctuation analysis (DFA).The signature of long-range correlations (power-law behavior) is detected globally with different strength at different geographic locations.Zonal mean exponents indicate that the reproduction of correlation properties is satisfactory for both models.In Section 4 we discuss the CCM results for the polar regions, where continuous empirical data are missing.Both CCMs predict an increased "memory" (stronger correlations), which seems to be in accordance with our present knowledge on the role of polar vortices.

Data and Methods
Semi-empirical total column ozone time series are obtained from Bodeker Scientific covering the time interval 1979-2011 with daily temporal, and 1.0 • × 1.25 • (lat/lon) spatial resolution.Missing satellite data are complemented by spatial linear interpolation [20,21].(For this reason, in this work we use the term "empirical data" instead of "measurements".)Satellite instruments use backscattered solar radiation to infer total column ozone, therefore measurements cannot be made during polar darkness.This fact limits the spatial coverage between 60 • S and 60 • N latitudes.Figure 1(a) illustrates a time series at an equatorial location north to Papua New Guinea (black line).It is easy to see that an annual cycle is missing, instead the variability is dominated by quasi-biennial oscillations.The two global chemistry climate models selected for comparisons are the LMDz-Reprobus and MRI as mentioned in the Introduction.Both have similar dynamical cores based on the primitive equations with a finite volume horizontal discretization on lat-lon grid [12].Horizontal resolutions are 2.5 • × 3.75 • (lat/lon) for the LMDz, and 2.8 • × 2.8 • (T42) for the MRI model, respectively.The vertical grids are terrain following hybrid pressure levels in both models with somewhat different resolutions: 50 levels for LMDz, and 68 for MRI.The higher vertical resolution in MRI and weaker horizontal diffusivity in the middle atmosphere than in the troposphere facilitate a dynamical representation of QBO.
The CCMVal-2 project defined a transient run from 1960 (with a 10-year spin-up period) to 2006.Total column ozone (TO) data were obtained at daily temporal resolution by numerically integrating ozone mixing ratios at the individual pressure levels.All forcings in the simulations were taken from observation and included the most important anthropogenic and natural sources (changes in trace gases, solar variability, major volcanic eruptions, QBO, etc.) [12].
The chemistry module Reprobus calculates the chemical evolution of 55 species using 160 gas-phase and 6 heterogeneous reactions, and the sedimentation of polar stratospheric cloud particles is also taken into account.MRI comprises a stratospheric chemistry including 34 long-lived species of 7 families and 15 short-lived species with 79 gas phase reactions and 34 photo-dissociations.
Figure 1 demonstrates a comparison with the satellite data at two "problematic" locations at the equator and at the edge of Antarctica.While QBO is present in the former case (Figure 1(a)), in the latter case (Figure 1(b)) the missing observations make it more difficult to accurately reproduce total column ozone values.Note that the MRI model outputs exhibit a slight overall upward shift of the mean TO levels as it is clearly seen in Figure 1.
Next we provide a short overview of the statistical methodology developed in [22], where long range correlations are already detected for ground measurements [23] and satellite TO data (TOMS Nimbus-7) [22].Here the analysis covers a definitely longer time interval than in [22] and it is extended for global CCM time series.
In order the characterize the relative strength of the different components, the spectral weight of a given peak is determined by calculating the area under the peak and normalizing by the total area under the whole spectrum.The estimates are somewhat inaccurate for the wide QBO peaks because the record lengths are limited and the resolution in the Fourier space strongly decays at low frequencies (see Figure 2(a)).

Spectral Weight Analysis
Power spectra at each grid point are determined by standard Fourier transformation for the three datasets described above.Visual inspection at several locations revealed that three main peaks characterize the TO time series depending on the geographic location: semiannual, annual and QBO components near the equator (see Figure 2 It is quite illuminating to inspect the high frequency tails of the power spectra shown in Figure 3 for the same geographic locations.A comparison with the empirical data (black lines) reveals that the high frequency variability is underestimated by both CCMs, especially in the equatorial band.(The same is visible already in Figure 1(a)).This might be the consequence of the limited spatial resolution of the global models producing rather smooth meteorological fields, where the representation of subgrid scale processes (e.g., small scale turbulent mixing, breaking internal waves, etc.) is restrained.However, this is a difficult question, because several collateral processes contribute to the short time changes of ozone levels.Note that the period on the horizontal axes is given in units of day.The same labels, coloring and notations are used.Eleven-point running mean smoothing is applied to suppress high frequency noise.Note that a direct quantitative cross-comparison of the empirical and model spectra cannot be easily performed.The standard method of simple normalization does not work because of the different shapes (see Figure 3).When a given spectrum exhibits a faster decay at the high frequency tail, normalization gives too large relative weight for the low frequency part.Also matching characteristic peaks (e.g., the annual oscillation) is problematic, because there is no guarantee that a model accurately reproduces a particular component of a signal and fails only at other frequencies.

Detrended Fluctuation Analysis
The method of detrended fluctuation analysis (DFA) has been developed for quantifying two-point correlations in non-stationary signals, because the autocorrelation function (or power spectrum) can erroneously indicate long range correlations in the absence of strict stationarity.Since DFA has become a standard tool in time series analysis [24][25][26], we give a very concise summary of the procedure.
Apparent periodic components in a signal spoil the results [27,28], therefore the first step is to obtain an anomaly time series a(t) (t = 1...N ) by computing the difference between an instantaneous TO value and the climatic mean for the given calendar day (obtained from as many years as available).Note that this procedure effectively removes stationary semiannual and annual oscillations from a signal, however it cannot wipe out a global trend or quasi-periodic components such as QBO.DFA regards an anomaly time series as successive increments of a random walk process, thus the next step is the "reconstruction" of the "trajectory" by numerical integration: y(τ ) = τ t=1 a(t), τ = 1...N .This signal is divided into equal segments (windows) of length n, and local trends are determined by a polynomial fit of order p, separately for each segment.The root mean square fluctuation value in a segment is obtained from the variance around the local trend.The DFAp fluctuation function F (n) is calculated as an average of local fluctuation values for a given window size n.This calculation can remove a weak polynomial trend of order p − 1 from the original signal, therefore it provides a more robust tool to quantify correlations than the direct determination of the autocorrelation function or power spectrum.
The presence of overall trends can be tested by comparing the DFAp fluctuation functions at increasing p values.In the case of ozone anomaly time series, we found that the shape of the curves does not change when p > 2, therefore we fixed the order of local polynomial detrending as p = 3. (A needlessly high order slows down the computations and results in higher inaccuracies, because the magnitude of fluctuation functions F p (n) monotonously decreases with an increasing p.)There is no methodological consensus around the maximum window size.Various ranges have been used in the literature, and values between n max = N/10 and N/4 are the most common choices.Detailed tests with computed model series suggest that even higher n max values can bear some useful information when long range correlations are expected [29], therefore we determined the DFA3 fluctuation function for the whole length N of each record.Finally we note that the curves are smoothed by using a sliding window technique, where the window of a given length n is shifted by 1 day along the time axis.
Figure 4 illustrates on double logarithmic scale the DFA3 fluctuation functions (p = 3) for three typical locations: at the equator (Figure 4(a)), somewhat away from it (Figure 4(b)), and at a midlatitude site (Figure 4(c)).Long range correlations are present in the signal when the fluctuation function follows a power law F (n) ∼ n α (straight line in a log-log plot) with exponent values α > 1/2.Clean power law behavior is rarely obtained for empirical data, especially for such complex systems as the atmosphere [22,23,28,30].The reason is that many different processes of different characteristic time scales affect the instantaneous value of any atmospheric parameter, therefore nonstationarities, various trends, cycles etc. "contaminate" the time series.For example, the apparent kink visible in Figure 4(a) for the empirical (black) and MRI (blue) data denoted by a red arrow is a clear consequence of QBO which is not present in the LMDz time series (red).As we mentioned, it is difficult to filter out quasi-periodic components.We tested one possible method, i.e., the Wiener filter for the QBO frequency band, however we decided not to apply it in the global analysis, because it affects long-range correlations in geographic locations where the QBO is not present at all.Instead, we determined a power law fit in two ranges for each F (n) curve: one is for the whole range of n (dashed lines in Figure 4), and the other for a restricted range between 60 days and 15 years (thin solid lines in Figure 4(b,c)).We are fully aware of the fact that the fluctuation functions deviate from a clean power law especially at the terminals of the curves, therefore we use the fitted exponent values as an approximate parameter of correlations.There is no way to visually check thousands of curves in order to determine an optimal fitting range for each of them.

Spectral Peak Maps
We determined the spectral weights of semiannual, annual and QBO peaks for each power spectrum by the method described in Subsection 2.1.Figure 5 displays the geographic distribution of these weights for the empirical and the two simulated datasets.Semiannual oscillations are better reproduced by the LMDz-Reprobus model around the equator, however the area of strongest spectral weights is shifted towards the East by about 70 degrees (Figure 5, first column).This component is very weak in the MRI signals, however both models indicate its appearance over the Antarctic where empirical data is not available.A remarkable small detail is the visible local maximum over the Tibetan Plateau in both models (albeit weak in MRI again), which is missing from the empirical data.
The annual oscillations are much better represented by the LMDz-Reprobus model (Figure 5, second column).The reproduction of the patterns at high latitude values over both hemispheres provides a strong confidence that the annual oscillations are properly simulated over the polar regions.The numerical values are somewhat higher in the LMDz model than in the empirical data, which is the consequence of too small variability at high frequencies (observed already in Figure 3) decreasing the overall area under each spectrum.The annual component is practically missing in the MRI signals around the equator and over the Arctic; even an area of maximal spectral weights appears over North America, which is missing from the empirical data.The missing QBO component from the LMDz-Reprobus time series is well-known and the reasons are understood [31] (Figure 5, third column).The reproduction is more satisfying in the MRI model, however note the different color coded ranges for the empirical and model data: QBO is present but much weaker than in the measured records.

DFA Exponent Maps
Third order DFA3 fluctuation functions are determined for the three datasets after removing the climatic mean values (see Section 2.2).Although the curves are only approximately straight lines in the log-log plots, we performed power law fits in order to provide a simple characterization on a global scale.The exponent values are definitely higher than 1/2 everywhere, indicating the presence of long memory processes in the dynamics and chemistry of ozone formation and decay.
The geographic distribution of the exponent values are color coded and displayed in Figure 6.The left column shows the results for whole-range fits, while exponents by restricted-range fits (60 days-15 years) are plotted in the right column.In general, both models seem to properly reproduce the fluctuation functions exhibited by the empirical data.The weak patterns around latitudes −60 • and 60 • appearing in the empirical data (Figure 6, top panels) are also present in the maps of model calculations, thus stronger correlations over the polar regions seem to be plausible.This observation is consistent also with our knowledge on stratospheric dynamics: meridional mixing is weaker across the walls of polar vortices, thus longer "memory" is expected in the internal regions.The stronger correlations (higher exponent values) around the equator were found already in [22] for shorter satellite records.The most probable reason for the weaker correlations in the midlatitude bands is the two-way coupling between the troposphere and stratosphere [32], e.g., the stronger dynamic activity in the troposphere in these geographic ranges excite more energetic internal gravity waves entering the stratosphere and inducing more intense mixing.4(a), red arrow): the apparent exponent values are systematically higher than in the case of whole range fits.Since the LMDz model does not have QBO, such effect is missing (red curves in Figure 6).The changes of exponent values away from the equator are due to the tail effects, both for small and large segment sizes n.Deviations from a pure power law behavior are also detected by the numerical values of the standard error of fitted slope coefficient.This quantity is color coded and plotted in Figure 8 both for the whole and restricted range fits.The most prominent feature for the empirical data (Figure 8(top)) is a band of somewhat higher fitting error around the equator.This is explained by the presence of QBO, which cannot be easily filtered out, and it is manifested as a marked kink on DFA3 fluctuation curves (cf. Figure 4(a)).This band is much weaker for the MRI simulations (Figure 8(bottom)), which is consistent with the results of spectral analysis ((Figure 5(bottom right)) indicating a QBO component of a suppressed amplitude compared with empirical data.Anomalously strong fitting errors appear above the Antarctic for both sets of CCM simulations (Figure 8 middle and bottom).Note that a precursor pattern is also visible in the empirical data close to −60 • (lat), between −150 • and −50 • (long), which is definitely enhanced toward the south pole.An explanation is given in Figure 9, where DFA3 fluctuation functions and power spectra are plotted for two locations over the Antarctic.The main problem is caused by the extreme nonstationary behavior of the ozone time series in this region associated with the ozone hole [33,34].Figure 1(b) clearly demonstrates a strong overall trend, which can be wiped out by a polynomial fit without any difficulty.However, besides the drop of mean ozone levels, the amplitude of the annual cycle is almost doubled since the 1980s, and this kind of variability change cannot be filtered out by computing climatic mean values.Figure 9(c) shows that a residual annual (and even semiannual) periodicity remains in the signal, which leads to the appearance of the kinks in Figure 9(a,b).Besides [22] and the present work, the most comprehensive analyses of total ozone correlations were performed in [35,36] on global scale.In the latter studies, monthly mean data were evaluated at rather low spatial resolution (5 • × 10 • , lat-long).The period of 29 years provided short time series of length of 348 points.Each time series was fitted by a combined statistical model of 21 free parameters consisting of the mean, the seasonal cycle, the QBO component, the solar cycle, a long-term trend, and residuals (noise).Apart from side effects of possible over-fitting, e.g., Figure 3 in [36] demonstrates fully inconsistent results depending on the method of estimating correlation exponents.The detected behavior changes from uncorrelated noise to extremely strong long range correlations at the same geographic latitude.We also found that numerical values are sensitive to the way of filtering and range of fitting (see Figure 7), however our results reflect consistency both for the three datasets and in a geographic sense: it is quite improbable that correlation properties change hectically from one latitude band to another.

Conclusions
Based on the presented analysis, we can summarize the results as follows: • Both the LMDz and MRI models provide an adequate reproduction of empirical data in the period of comparison (see Figure 1), however the latter seems to slightly overestimate total column ozone values.• Spectral analysis reveals that both models underestimate high frequency (daily) variability (see Figure 3), MRI performs better from this point of view.Such deficiency might have a rather complex origin, because several dynamical processes (advection, diffusion and mixing) contribute to short time fluctuations.• Spectral weight analysis indicates some anomalies in the dynamical core of both CCMs.The geographic locations of semiannual oscillations are shifted in the models.MRI is known to reproduce the QBO mode along the equator (albeit it is definitely weaker than in the empirical records), however it underperforms with respect of the annual component.• The presence of quasi-biennial oscillations decreases the accuracy of DFA analysis (see Figure 4(a)).Band pass Fourier filtering can remove such spectral component (see [22]), however one should be careful, because the DFA method is based on residual (anomaly) analysis.Specifically, the standard Wiener filtering consists of two steps: firstly, the Fourier transform of a given time series is produced and the amplitudes are set to be zero or to some assumed (practically interpolated) base line values in the required frequency ranges, and secondly, an inverse Fourier transform produces the filtered time series.When this method is applied on a global scale, obviously the filtered components disappear everywhere, also from signals where QBO oscillations are not present.Care should be taken especially when the target is to detect long range correlations, because such procedure might distort the original properties of the time series.Too strong filtering can yield (geographically) inconsistent results as probably demonstrated in [35,36].
Note that the low model values of high frequency variabilities are reflected by the fact that the initial parts of fluctuation functions remain systematically below the curve for empirical data (Figure 4).• The DFA3 fluctuation functions are adequately reproduced by both CCMs (see Figures 6 and 7).
We do not claim that the curves represent pure power law behavior (see Figure 4); still the overall trends are close to linear on log-log scale, indicating the presence of long range correlations.
(The appearance of clean power laws is not even expected, because many processes of various characteristic time scales contribute to the instantaneous value of any atmospheric variable at a given location.)• Since the reproduction of correlation properties in the models is quite satisfactory in the overlapping geographic latitude band between −60 • and 60 • , we can have a strong confidence in the prediction of both models that correlations are stronger over both polar regions.This observations is also supported by considering the role of polar vortices in the meridional transport processes.• Classical methods of time series analysis have a limited applicability over the polar regions, especially over the Antarctic.Such strong nonstationary signals as represented by the total ozone level during the decades of ozone hole make any residual (anomaly) analysis very difficult.Standard methods to remove semiannual and annual oscillations fail when the amplitudes change strongly (see Figure 9), however an over-filtering might yield information loss about correlation properties.

Figure 1 .
Figure 1.(a) Empirical and simulated total column ozone (TO) time series for an equatorial location over the Pacific Ocean in Dobson units (DU).Black line belongs to the Bodeker Scientific record (label: Bd), red and blue denote CCM simulations by the LMDz-Reprobus (LMDz) and MRI models; (b) Model time series over the Antarctic (see label).Note the different vertical scales.

Figure 2 .
Figure 2. Non-normalized power spectra for three different geographic locations (see legends) as a function of period.Note the logarithmic horizontal scales.Black lines with circles belong to the empirical Bodeker Scientific records (label: Bd), red and blue denote CCM simulations by the LMDz-Reprobus (LMDz) and MRI models.(a) Equatorial location (Gabon, Africa); (b) South Pacific Ocean; (c) Québec, Canada.
(a)), a marked semiannual component and missing QBO away from the equator at some places (Figure 2(b)), and a clean annual periodicity in the midlatitude bands, especially on the Northern hemisphere (Figure 2(c)).

Figure 3 .
Figure 3.The high frequency tails of the spectra in Figure 2 on double logarithmic scales.Note that the period on the horizontal axes is given in units of day.The same labels, coloring and notations are used.Eleven-point running mean smoothing is applied to suppress high frequency noise.

Figure 4 .
Figure 4. Logarithm of DFA3 fluctuation function log 10 [F (n)] as a function of logarithmic window size log 10 n for three different geographic locations (see legends).Black lines belong to the empirical Bodeker Scientific records (label: Bd), red and blue denote CCM simulations by the LMDz-Reprobus (LMDz) and MRI models.(a) Gulf of Guinea, near the equator; (b) Indian Ocean, close to Christmas Island; (c) Guo Dao, Xinjiang, China.Thin dashed lines indicate whole range fits for the LMDz data, thin straight lines in (b) and (c) denote the same for the restricted range (60 days-15 years), see text.

Figure 5 .
Figure 5. Geographic distribution of the spectral peak weights for the empirical Bodeker Scientific records (label: Bd) and CCM simulations by the LMDz-Reprobus (label: LMDz) and MRI models (label: MRI).The columns from left to right display the weights of semiannual, annual and QBO peaks, respectively.Note the different color-scale ranges for the QBO maps.

Figure 6 .Figure 7
Figure 6.Geographic distribution of the DFA3 exponent values for the empirical Bodeker Scientific records (label: Bd) and CCM simulations by the LMDz-Reprobus (label: LMDz) and MRI models (label: MRI).The left column displays fitted values for the whole range of available data, the right column shows the results for the restricted range fits from 60 days to 15 years (from 1.78 to 3.74 on the log 10 axis), see Figure 4. DFA3 (whole range fit) DFA3 (60 days − 15 years fit)

Figure 7 .
Figure 7. Zonal mean values of the DFA3 exponents for the empirical Bodeker Scientific records (label: Bd) and CCM simulations by the LMDz-Reprobus (label: LMDz) and MRI models (label: MRI) as a function of geographic latitude.(a) Fitted values for the whole range of available data; and (b) restricted range fits from 60 days to 15 years (from 1.78 to 3.74 on the log 10 axis), see Figure 4. Error bars represent the variability of exponent values along a circle of latitude.

Figure 8 .
Figure 8. Geographic distribution of the standard fitting error of DFA3 exponent values for the empirical Bodeker Scientific records (label: Bd) and CCM simulations by the LMDz-Reprobus (label: LMDz) and MRI models (label: MRI).The left column displays fitted values for the whole range of available data, the right column shows the results for the restricted range fits from 60 days to 15 years (from 1.78 to 3.74 on the log 10 axis), see Figure 4.Note the different color scales at the two columns.

Figure 9 .
Figure 9. (a) and (b): Logarithm of DFA3 fluctuation function log 10 [F (n)] as a function of logarithmic window size log 10 n for two different geographic locations (see legends) over the Antarctic.Red and blue denote CCM simulations by the LMDz-Reprobus (LMDz) and MRI models, empirical data are not available; Red arrow on (b) indicates the characteristic kink appearing everywhere for locations over the Antarctic; (c) Power density spectra for the total column ozone anomaly records (annual cycle removed) simulated by the LMDz (red) and MRI (blue) models over an Antarctic location (see legend).