## Abstract

Two sea surface temperature (SST) time series, the Extended Reconstructed SST version 3 (ERSST.v3) and the Hadley Centre Sea Ice and Sea Surface Temperature dataset (HadISST), are used to investigate SST multidecadal variability in the Mediterranean Sea and to explore possible connections with other regions of the global ocean. The consistency between these two time series and the original International Comprehensive Ocean–Atmosphere Dataset version 2.5 (ICOADS 2.5) over the Mediterranean Sea is investigated, evaluating differences from monthly to multidecadal scales. From annual to longer time scales, the two time series consistently describe the same trends and multidecadal oscillations and agree with Mediterranean ICOADS SSTs. At monthly time scales the two time series are less consistent with each other because of the evident annual cycle that characterizes their difference.

The subsequent analysis of the Mediterranean annual SST time series, based on lagged-correlation analysis, multitaper method (MTM), and singular spectral analysis (SSA), revealed the presence of a significant oscillation with a period of about 70 yr, very close to that of the Atlantic multidecadal oscillation (AMO). An extension of the analysis to other World Ocean regions confirmed that the predominance of this multidecadal signal with respect to longer period trends is a unique feature of the Mediterranean and North Atlantic Ocean, where it reaches its maximum at subpolar latitudes. Signatures of multidecadal oscillations are also found in the global SST time series after removing centennial and longer-term components.

The analysis also reveals that Mediterranean SST and North Atlantic indices are significantly correlated and coherent for periods longer than about 40 yr. For time scales in the range 40–55 yr the coherence between the Mediterranean and subpolar gyre temperatures is higher than the coherence between the Mediterranean SST and North Atlantic Oscillation (NAO) or AMO. Finally, the results of the analysis are discussed in the light of possible climate mechanisms that can couple the Mediterranean Sea with the North Atlantic and the Global Ocean.

## 1. Introduction

Instrumental records of increasing duration and spatial coverage document substantial ocean variability, on time scales ranging from months to decades. Some of this variability can be easily related to known forcing mechanisms, but in many cases the relationship between observed variability and forcing is less clear. It has been noted, for example, that the climatic system displays variability on many time scales, even in the case of a constant external forcing, that is, in the absence of changes in insolation or atmospheric composition (Dijkstra and Ghil 2005; Pisacane et al. 2006; Knight et al. 2005; Grist et al. 2009).

The global atmospheric and oceanic circulation has a number of preferred patterns of variability conventionally measured by specific climate indices. All these have expressions in surface climate fields that can have a profound regional influence, modulating the location and strength of the local fluxes of heat, moisture, and momentum. Classical examples of these “climate modes” are the North Atlantic Oscillation (NAO) (Hurrell 1995; Trenberth et al. 2007; Knight et al. 2006) and the Atlantic multidecadal oscillation (AMO) (Enfield et al. 2001). The NAO is a leading pattern of weather and climate variability over the Northern Hemisphere, which contributes to the redistribution of air masses between the Arctic and the subtropical Atlantic, influencing weather conditions and SST over the Atlantic and adjacent seas. AMO is measured by the detrended SST anomalies averaged over the North Atlantic from 0° to 70°N (Enfield et al. 2001). Other definitions of AMO have been obtained by subtracting the global mean SST to remove the large-scale global signal that is associated with global processes (Trenberth and Shea 2006) or by using empirical orthogonal functions (EOFs) (Parker et al. 2007). AMO is characterized by a cycle of about 70 years and has been the object of increasing interest during recent years. There are two main hypotheses about its origin. The first relates this multidecadal cycle to the internal variability of the thermohaline circulation (THC) and to the intrinsic nonlinear behavior of the climate system (Huck et al. 1999; Delworth and Mann 2000; te Raa and Dijkstra 2002; Knight et al. 2005), while the second, based on experimental and numerical evidence, brings into play the role of the free oscillation motion of the atmosphere–ocean–ice coupled system (Jungclaus et al. 2005; Dong and Sutton 2005; Dima and Lohmann 2007). Recent numerical experiments (Grosfeld et al. 2008) indicate the Atlantic as the main driver for the AMO with a clear coherence between surface air temperature and sea level pressure (SLP), exerting influence on the Euro–Atlantic sector. The analysis of the surface temperatures performed by Schlesinger and Ramankutty (1994) suggested that the multidecadal oscillation arises from predictable internal variability of the ocean–atmosphere system.

Analyses of coupled models and reconstructed datasets over the global ocean already pointed out the influence of the AMO over the climate of Euro–Atlantic sector (Sutton and Hodson 2005; Knight et al. 2006; Hodson et al. 2009). Sutton and Hodson (2005) documented for the first time the influence of AMO on variations in the summertime climate of both North America and western Europe. Mohino et al. (2011) found that the decadal Sahel rainfall variability is connected to AMO and related positive anomalies over the Mediterranean Sea. None of these works studies in detail the Mediterranean SST multidecadal variability even if the Mediterranean Sea is often included in their global analysis (e.g., Fig. 1 of Sutton and Hodson 2005).

The main objective of this paper is to place the Mediterranean Sea at the center of the debate by investigating in detail the multidecadal variability of the Mediterranean SST field and its relation with the global ocean.

This objective has been pursued using the most recent reconstructed (and non reconstructed) SST datasets (Rayner et al. 2003; Rayner et al. 2006; Smith et al. 2008).

Recent critical analyses of SST datasets have highlighted the need for fully validated climatic data records to draw conclusions on trends or low-frequency oscillations (Folland and Parker 1995; Thompson et al. 2008). We first analyzed the robustness and reliability of the longest instrumental SST time series of the Mediterranean Sea now available (section 2). Then connections between Mediterranean SST and AMO are demonstrated using the singular spectral analysis (SSA) and the multitaper method (MTM) (section 3). Moreover, the Mediterranean multidecadal oscillations are compared with other ocean regions to investigate their relevance and correspondences with global-scale ocean variability (sections 4 and 5). In the last section, the results of the analysis are discussed in the light of possible climate mechanisms that can couple the Mediterranean Sea with the North Atlantic and global ocean. Some speculative but intriguing hypotheses on feedbacks between Mediterranean outflow and the North Atlantic THC based on modeling (Artale et al. 2002; Artale et al. 2006; Calmanti et al. 2006) and observations (Reid 1979; Lozier and Stewart 2008) are also proposed.

## 2. Data and methods

### a. The SST datasets

The variability of the SST field over the last 150 yr has been analyzed using both reconstructed and nonreconstructed SST datasets. One is the International Comprehensive Ocean–Atmosphere Dataset (ICOADS) SST version 2.5, which includes all available global surface marine data from 1800 to September 2010 (see Worley et al. 2005 for ICOADS version 2.1 and Woodruff et al. 2011 for ICOADS version 2.5). This dataset consists of monthly summary statistics of several quality-checked variables at 1° resolution starting from 1960 and at 2° resolution starting from 1800. The dataset contains data voids where in situ observations were not enough to compute significant monthly statistics.

The Extended Reconstructed SST (ERSST) reconstruction uses the ICOADS SST and is based on a separate analysis of the low- and high-frequency variations. The former are analyzed using a simple average and smoothing of the relatively sparse data, while the analysis of the interannual and shorter-period variations is based on a fit of observed high-frequency SST anomalies to a set of empirical orthogonal teleconnections (Van den Dool et al. 2000) derived from the more recent spatially complete SST analysis (Smith and Reynolds 2003; Smith and Reynolds 2004). In this work we used the updated version 3 (ERSST.v3; Smith et al. 2008) of the dataset that also includes the adjustment for data before 1941, taking into account differences in measurement practices before that year (Smith and Reynolds 2002).

The third dataset we use is the Hadley Centre Sea Ice and Sea Surface Temperature dataset (HadISST), which is based on the Met Office Marine Data Bank (MDB) but also includes ICOADS SSTs. In this dataset, sea surface temperature gaps are reconstructed using a two-stage reduced-space optimal interpolation (OI) procedure, followed by superposition of quality-improved gridded observations onto the reconstructions to restore local detail (further details are given in Rayner et al. 2003). HadISST also addresses the problem of the bias in historic SSTs because of changes in SST sampling methods in the 1940s and earlier, using a bias adjustment method proposed by Folland and Parker (1995). ERSST and HadISST use bias-adjusted Advanced Very High Resolution Radiometer (AVHRR) SST beginning in 1985 and 1982, respectively.

### b. Intercomparison of SST datasets over the Mediterranean Sea

The annual sea surface temperature anomaly (SSTA) during the last 100–150 yr can be more easily estimated using reconstructed SST datasets, such as ERSST and HadISST, which are the most extensively used. However, globally complete, that is, interpolated, datasets must be used with care (see, e.g., Rayner et al. 2003), particularly in specific regions or periods of time where data are too sparse. Since for most of the time span there are no independent sea “truth” data to validate the reconstructed SST, the only way to evaluate the consistency and quality of the reconstructed time series over particular regions is to intercompare information extracted from different datasets. Here, we have compared the ERSST and HadISST datasets in the Mediterranean region, both on annual and monthly scales, to evaluate differences in long-term trends or oscillations and seasonal components. Focusing on the Mediterranean Sea, we notice a good availability of ICOADS (2 × 2 degree resolution) data for most of the analyzed period (Fig. 1). In fact, excluding the periods of the two World Wars, during the last century ICOADS observations are available in about 80% of Mediterranean grid boxes or more. This makes the Mediterranean Sea one of the best sampled areas of the World Ocean.

This abundance of data suggests the possibility of estimating the time evolution of the mean sea surface temperature of the Mediterranean Sea directly from the ICOADS SSTs. So the Mediterranean monthly mean time series can be computed by simple space average of all valid ICOADS SST anomalies with respect to a mean SST derived from the entire ICOADS record.

As suggested by several authors (Bottomley et al. 1990; Folland and Parker 1995), ICOADS SST prior to 1942 has been adjusted using corrections provided by Smith and Reynolds (2002) before any processing step. Smith and Reynolds (2002) showed a general consistency between their corrections and those suggested by Folland and Parker (1995) even though some more marked difference can be observed in winter at high latitudes where the Smith and Reynolds correction is stronger (Smith and Reynolds 2003).

Starting from monthly ICOADS 2.5 data, yearly averages were computed only for those grid points where 12 monthly values were available for the specific year under examination. In years in which the yearly SST was estimated in less than 15% of the Mediterranean Sea grid points, the sea surface temperature was not computed, resulting in data gaps in the annual time series. The 15% limit was empirically fixed after several trials. Figure 2 shows the percentage of grid boxes with yearly averages available in the Mediterranean Sea for each year. Figure 3 compares the three time series for the yearly averages.

The ICOADS SSTs, which are not interpolated, follow quite closely those from the two other series, with differences reaching 0.2°C between 1880 and 1920. The three time series basically describe the same history both in terms of trends and of multidecadal oscillations and can consequently be used interchangeably to study the long-term variation of the Mediterranean SST field. Moreover, the agreement between the two interpolated SST series (ERSST.v3, HadISST) and the noninterpolated one suggests that damping or other possible distortions of the interpolation are not dominant in the Mediterranean basin when yearly averages are used.

As shown in Fig. 4a, on a monthly scale the differences between the two analyses are larger. The differences exhibit a marked annual cycle, suggesting that a form of compensation between seasons may yield the better agreement found for the yearly averages. Similar differences are found between ERSST.v3 and ICOADS (Fig. 4b) and HadISST and ICOADS (Fig. 4c). Figures 4a and 4c show an increased amplitude of the annual cycle before 1941, when instrumental bias corrections have been applied. A portion of the difference in the amplitude of the annual cycle can possibly be attributed to the different bias adjustment used for the three datasets. Figure 4b, where ERSST and ICOADS are compared, supports this hypothesis since the difference between these two datasets, which use the same bias correction, is smaller and there is no clear difference between the periods preceding and following 1941.

This intercomparison does not allow us to identify which of the series is responsible for this annual cycle, but we can still use two additional SST series covering the last three decades to further investigate the problem: the Optimally Interpolated Pathfinder SST (OI_Path) over the Mediterranean Sea (Marullo et al. 2007) and the Reynolds OIv2 maps based on AVHRR and in situ data (Reynolds et al. 2007) (OIPath and OIv2 hereafter).

The validation performed by Marullo et al. (2007) showed that OIPath is able to reproduce in situ measurements with a mean bias of less than 0.1 K and RMSE of about 0.5 K and that errors do not drift from year to year (and drift very little seasonally) or with the percent interpolation error. The absence of any major annual cycle in the difference between in situ data and OIPath SSTs makes this dataset the ideal reference to investigate seasonal fluctuation of errors in the reconstructed SST fields.

Reynolds OIv2 are expected to have similar or even better performances being forced to be unbiased with respect to in situ data. Figures 5a and 5b show the differences between the three SST time series under investigation and OIv2 and OIPath, respectively, during the satellite era. Also, in this case a clear annual oscillation, with amplitude larger than 0.5°C, is present. If satellite estimates are free from seasonally oscillating errors as suggested by Marullo et al. (2007), the annual cycle must be attributed to all the reconstructed datasets. This hypothesis is further confirmed by Fig. 5c where the monthly differences between in situ (obtained by the Coriolis dataset, http://www.coriolis.eu.org/Data-Services-Products/View-Download/Data-selection) and OIPath satellite data show that these estimations are free from any annual cycle distortion.

Rayner et al. (2003) have already noted that in the Northern Hemisphere summer, in areas that are well sampled by in situ data, HadISST is systematically warmer than OIv2 SSTs with differences that reach about 0.7°C in some locations and suggested that this is symptomatic of a general difference between HadISST and COADS or OIv2. Figures 5a and 5b confirm that, at least in the Mediterranean Sea, a similar bias also occurs between ERSST and OIv2 or OIPath. Some more inference can be achieved by averaging the S1 = (ERSST − OIPath), S2 = (HadISST − OIPath), and S3 = (all in situ − OIPath) monthly records of Figs. 5b and 5c to obtain a single mean monthly cycle for each of the three time series. The differences S1 − S3 and S2 − S3 will approximately represent the mean monthly bias of ERSST and HadISST reconstruction with respect to in situ observations (Fig. 6).

While taking into account the limitations of this estimate due to the uneven distribution of the in situ data and to intrinsic differences related to point measurements with respect to averages over large boxes or satellite pixels, Fig. 6 reinforces the Rayner et al. (2003) conclusions about the summer HadISST bias and confirms that the same bias also affects ERSST estimates. The mean annual difference between HadISST and in situ and between ERSST and in situ are −0.03°C and −0.13°C respectively, that is, very close to zero.

In conclusion, we can state that the three time series are interchangeable and anyone of them can be used, depending upon convenience, for climatic studies, provided that they are yearly averaged, to compensate the annual cycle observed in Figs. 4 and 5. On the other hand, the monthly time series should be used with more care.

### c. Spectral analysis methods

We have computed yearly SST averages from ERSST and HadISST for the Mediterranean Sea, for the Global Ocean, and for eight Global Ocean subregions shown in Fig. 10 (shown later) . Then, annual time series of anomalies with respect to the 1971–2000 mean have been computed and used to investigate the variability of the SST field. Three different methods have been used. First, an autocorrelation analysis was applied to the Mediterranean Sea and to the global ocean ERSST and HadISST anomalies time series, with time lags between 0 and 100 yr, to highlight possible multidecadal oscillations. Then, the time series were also analyzed using SSA and the MTM. The SSA technique is a nonparametric spectral estimation method to extract information from short and noisy time series and thus provides insight into the unknown or only partially known dynamics of the underlying system that generated the series (Broomhead and King 1986; Fraedrich 1986; Vautard and Ghil 1989). This method is based on embedding a time series *X*(*t*): *t* = 1, *N* in a vector space of dimension *M*. In practice, SSA proceeds by diagonalizing the *M* × *M* lag-covariance matrix of *X*(*t*) to obtain spectral information on the time series. The entries *c _{ij}* of this matrix, which depend only on the lag |

*i − j*|, are

The *M* eigenvectors *E _{k}* of the lag-covariance matrix are the temporal EOFs. The eigenvalues

*λ*of account for the partial variance in the direction

_{k}*E*, and the sum of the eigenvalues gives the total variance of the original time series

_{k}*X*(

*t*). This is analogous to the classical principal component analysis (PCA) method where the lag-covariance matrix is substituted by the covariance matrix.

SSA can be combined with advanced spectral analysis methods such as MTM to refine the interpretation of oscillatory behavior. In this sense, SSA and MTM can be considered complementary. The MTM is a technique developed by Thomson (1982) to estimate the power spectrum of a stationary ergodic finite-variance random process, which is useful for the analysis of relatively short time series whose spectrum may contain both broadband and line components. The tapering is a solution to approach the leakage problem of finite time series. Classical methods use a unique data taper multiplying a portion of the data at the beginning and at the end of the series by an appropriate-weighting function. This implies that portions of the time series are excluded from the analysis as a trade-off from reducing spectral leakage in the frequency domain. MTM uses a small set of tapers rather than a single taper. The data are multiplied by orthogonal tapers constructed to minimize the spectral leakage because of the finite length of the time series. These optimal tapers belong to a family of functions known as “Discrete Prolate Spheroidal Sequences” (DPSS) obtained as eigenvectors of a suitable Rayleigh–Ritz minimization problem (Slepian 1978). The algorithm to calculate the prolate–Slepian tapers is described by Lees and Park (1995). A detailed description of the MTM is given by Ghil et al. (2002).

To determine whether broadband or line components are present, a reshaped spectrum can be determined in which the contributions from harmonic signals are removed (Thomson 1982; Ghil et al. 2002). The comparison between the reshaped, in which contributions from pure harmonic signals detected by an *F*-variance ratio test with a minimum significance threshold have been removed, and unreshaped MTM gives the opportunity to distinguish signals that are best approximated by harmonic or narrowband oscillations from amplitude and phase modulated and possibly intermittent oscillations.

The use of the MTM implies the choice of the number of tapers (*K*) and the integer bandwidth parameter (*p*). The choice of the bandwidth *2pf _{R}* [where

*f*is the Rayleigh frequency

_{R}*1*/(

*N*Δ

*T*) with

*N*being the number of samples and Δ

*T*the sampling interval], that is, the choice of

*p*and

*K*, is the classical trade-off between spectral resolution and the stability or variance reduction properties of the spectral estimate.

## 3. Analysis of the Mediterranean multidecadal SST variability

Figure 3 shows that Mediterranean Sea SST time series exhibit warm phases (positive anomalies respect to the 1971–2000 average) during 1860–80, 1925–70, and 1985–today, while cold phases occurred during 1880–1925, 1970–85, and presumably before 1860. This suggests that an oscillation with a period of about 70 yr, similar to that of the AMO, is also present in the Mediterranean Sea. The occurrence of such an oscillation is confirmed by the lagged autocorrelation plot of Fig. 7, where a clear maximum for a time lag of 70 yr is obtained using both the ERSST and the HadISST temperatures.

To get further physical insight and information on the whole spectrum, we then use the SSA and the MTM methods.

We have applied the SSA method to the detrended Mediterranean ERSST (the longest time series). To obtain the signal-to-noise (S/N) separation we computed and plotted the eigenvalue spectrum of the detrended time series (Fig. 8a). The detrending was applied to remove any possible deterministic trend leaving eigenvalues 1 and 2 together to explain the multidecadal oscillation. The first two eigenvalues form a pair that contributes about 40% of the total variance (Fig. 8b) and are separated by a very steep slope from the noise, which is characterized by much lower values and a mild slope. The distance of these two modes above the noise indicates that they are deterministic rather then stochastic (Vautard and Ghil 1989; Vautard et al. 1992). The first two leading EOFs have approximately the same amplitude (considering the range of error bars) and are in quadrature (Fig. 8c). Finally, it is worth noting that with only the first two EOFs we capture the main low-frequency variability of this SST time series (Fig. 8d). This fact could be interpreted as the occurrence of a ghost limit cycle related to a physical oscillation of the dynamical system that has generated our SST time series (Ghil et al. 2002).

The low-frequency variability of the Mediterranean Sea was also investigated applying the multitaper method to the detrended yearly SST derived from the ERSST.v3 dataset. Following Ghil et al. (2002), we set *p* = 2 and *K* = 3 (less or equal to 2*p*-1) as a compromise between the need to resolve decadal scales and the benefit of multiple spectral degrees of freedom. We also detrended the ERSST.v3 time annual series over the ocean areas under investigation to better distinguish between low-frequency bands and the trend.

Figure 9 shows the comparison between the harmonic analysis test (Fig. 9a) and the adaptively weighted estimation (Fig. 9b) against red noise background for the ERSST.v3 series of the Mediterranean Sea. The *F*-test criterion for harmonic signals (Fig. 9a) yields 6 peaks at the 99% confidence level and 10 peaks at the 95% confidence level that can eventually be associated to some harmonic, phase-coherent oscillation. Some of the peaks (7.3, 4.3, and 3.2 yr) that exceed the 99% confidence level in the harmonic test plot (Fig. 9a) are associated with very weak power in the spectrum (Fig. 9b) while others (73, 6.3, and 2.8 yr) correspond to significant spectral bands in the spectrum.

The harmonic test is also useful in the reshaping procedure. Figure 9b shows five peaks that satisfy the harmonic detection test at the 90% level and are also significant relative to the red noise in the spectrum. They include a low-frequency band that peaks at 73 yr and four higher frequency that peak at 6.3, 3.9, 2.8, and 2.2 yr. The higher-frequency peaks are very close to the preferred scale of variability of the ENSO’s quasi-biennial and low-frequency modes (Rasmusson et al. 1990; Jiang et al. 1995; Yiou et al. 2000) while the lower frequency is close to that of the AMO. The comparison between the reshaped and the unreshaped spectrum in Fig. 9b suggests that the four higher-frequency peaks (two significantly well above the 99% confidence level and two over the 90% confidence level) are likely to contain a harmonic phase-coherent oscillation with, at least, a 90% level of confidence. The same is not true for the low-frequency band at 73 yr where the reshaped and the unreshaped spectra are identical. Even if this could be interpreted either as the consequence of the shortness of the SST time series or as the signature of a phase and amplitude modulation and/or an intermittently oscillatory behavior, the multidecadal band is a significant feature of the spectrum. This inference is also supported by the SSA analysis that suggests the physical correspondence of the multidecadal oscillation with some geophysical phenomenon (see Fig. 8b and related comments).

However, some concern still remains about the possibility to truly determine the existence of a 73-yr periodicity in a 155-yr time series. Previous authors (e.g., Wittenberg 2009) have verified that century-long time-scale records were not sufficient to resolve the temporal variability of ENSO, implying that a different 150-yr record may not contain the same oscillation. Having no possibility to extend back in time the instrumental data records, one can still investigate whether coupled GCM simulations of several hundreds of years present plausible scenarios where such or similar multidecadal oscillations are present. To increase the size of our Mediterranean data record we have analyzed, using the same MTM spectral tool, the third climate configuration of the Met Office Unified Model (HadCM3; Met Office 2006) control run, which is a coupled atmosphere–ocean circulation model that has been run for 1000 yr. This analysis has evidenced, in both harmonic signal *F* test and MTM spectrum, the presence of a multidecadal peak centered at about 100 yr well above the 99% significance level. In this case, the Rayleigh frequency resolution for a harmonic signal resulting from the *F* test is ±*f _{R}* = ±1000

^{−1}cpy (cycles per year) implying that this peak is between 90 and 110 yr. At the same frequency, the MTM spectrum shows the presence of a narrowband (±

*pf*bandwidth) that is likely to contain this harmonic phase-coherent oscillation with a 99% confidence limit. The presence of similar centennial fluctuations of the North Atlantic Ocean in HadCM3 has been already reported by previous authors (Vellinga and Wu 2004; Frankcombe et al. 2010).

_{R}In the case of the shorter ERSST Mediterranean time series, the broader low-frequency band contained the 73-yr peak of the harmonic analysis test but with Rayleigh frequency resolution of = ±155^{−1} cpy, implying that the peak was between 50 and 138 yr. Consequently, the two peaks at about 100 and 70 yr, found in the HadCM3 and in the ERSST time series, respectively, cannot be considered significantly separated, given the length of the two time series. This result, although not definitive proof, is consistent with the hypothesis that a different 155-yr record should also contain a similar 70-yr peak.

## 4. Mediterranean Sea versus global ocean

The Mediterranean multidecadal oscillation observed in Fig. 3 and confirmed by the subsequent time series analysis (section 3) is not evident in the global SSTs time series. The SSA and MTM analyses (not shown) of the global SST time series performed using the same parameters as for the Mediterranean Sea do not reveal any significant pair of EOFs corresponding to the 70-yr oscillation. In this case, the first two temporal EOFs, which contribute about 68% of the total variance, suggest the presence of a longer than 70-yr-period oscillation. The third mode, which only accounts for 9% of the variance, is associated to a period of 68 yr but does not form any pair with the next mode that contributes to the total variance only with 5%. As suggested by Schlesinger and Ramankutty (1994), the absence of a dominant multidecadal oscillation can be the effect of long-term components not completely removed applying a simple linear detrending of the time series before the SSA analysis. In fact, their analysis of the Intergovernmental Panel on Climate Change (IPCC) global surface temperature, detrended for the global-mean response simulated by a simple climate–ocean model, revealed the existence of a 65–70-yr oscillation. To evaluate the importance of the multidecadal signal in the global SST time series, we then applied the same procedure used by Dima and Lohmann (2007) for the analysis of Fram Strait sea ice reconstruction. Long-term components (centennial and longer), resulting from the SSA of the global SST time series, have been identified and the corresponding reconstructed time series have been subtracted from the original one to produce a residual time component. Appling the SSA to the residual time series, the first two eigenvalues, corresponding to a period of 56 yr, appear in quadrature, and explain about 45% of the total variance of this residual time series, showing that the multidecadal signal was masked by the longer-period trends. The difference between the 56-yr oscillation period found in our analysis and the 65–70-yr period found by Schlesinger and Ramankutty (1994) can be either due to the different length of the two time series (155 and 135 yr, respectively) or to the difference between the temperature observations themselves, ERSST being limited to sea surface temperatures while the IPCC data used by Schlesinger and Ramankutty (1994) also include land surface observations. Schlesinger and Ramankutty (1994) concluded that their multidecadal signal in the global time series is the statistical result of a regional long-time-scale oscillation that exclusively occurs in the Northern Hemisphere. Taking into account the results of their IPCC time series analysis, it makes sense to investigate whether the 70-yr multidecadal signal is still recognizable in the ERSST dataset only in the Northern Hemisphere. A first screening can be done by applying the lagged autocorrelation procedure used for the Mediterranean Sea to other regions of the oceans. We have divided the oceans in eight boxes (subregions) and produced an autocorrelation plot for each of them. Figure 10 suggests that the North Atlantic Ocean is the only other area of the World Ocean that clearly exhibits the same oscillation with period of 70 yr observed in the Mediterranean Sea. Other ocean regions could also have multidecadal SST oscillations even if they appear much less pronounced in the lagged autocorrelation analysis. The SSA analysis applied to the 8 regions of Fig. 10 (not shown) confirmed that the 70-yr AMO-like oscillation is the dominant mode of the North Atlantic, while other subregions of the Northern Hemisphere reveal multidecadal oscillations only after removal of the dominant centennial and longer components. This is in agreement with the results of Fig. 2 of Schlesinger and Ramankutty (1994). Our SSA analysis of SST residuals in the North Pacific sectors (boxes 7 and 3 in Fig. 10) confirms the presence of the AMO-like oscillation observed by Zhang and Delworth (2007). Our analysis indicates that this oscillation has a period of 55–60 yr and explains 30%–35% of the variance of the residuals time series.

Since the North Atlantic box covers a wide latitudinal range going from tropical to subpolar regions, it is interesting to investigate the latitudinal distribution of the low-frequency peak, so to quantify the relative role of the polar and of the southern regions in modulating the multidecadal oscillation. To do so, we have applied the MTM analysis to latitudinal bands of 20° in the North Atlantic, which move with steps of 2° from 10°N to 60°N using ERSST data. The latitudinal behavior of the amplitude of the low-frequency peak resulting from harmonic analysis test of each time series is shown in Fig. 11.

Peaks that exceed the 99% confidence limit centered between 64 and 73 yr are observed only for latitudinal bands north of 30°N with a maximum at 50°N. Note that the difference between the three periods of Fig. 11 is not significant considering the Rayleigh frequency resolution ±*f _{R}* (

*f*= 1/

_{R}*NΔT*) that results from the annual time series length. All the peaks indicated in Fig. 11 correspond to significant bands (over 99% confidence limit) in the power spectrum so, even if the harmonic test alone does not definitively prove the existence of a particular oscillation, the simultaneous existence of harmonic signals and spectrum bands that satisfy a 99% confidence limit is a strong indication that such an oscillation exists. In addition, the qualitative analysis of the time series plots (Fig. 10) further supports such an assumption.

The latitudinal variations of the amplitude of the harmonic signal demonstrate that multidecadal oscillations observed in the residual of the global SST time series are essentially due to signal of the northern regions of the North Atlantic, where the AMO signal is dominant. The results of Fig. 11 are in agreement with the AMO spatial pattern resulting from the EOF analysis of Parker et al. (2007) (see their Fig. 6) and with the hypothesis of Frankcombe et al. (2010) of the Arctic origin of the AMO pattern.

## 5. Correlation of SST with NAOI

The existence of a multidecadal global oscillation, controlled by the Northern Hemisphere regions, poses the problem of its internal or external origin. As discussed by Schlesinger and Ramankutty (1994), there are three possible sources of this oscillation: an external oscillatory forcing, such as a variation in the solar constant; a random forcing, such as white noise; and an internal oscillation of the atmosphere–ocean system. They argued that, owing to the uneven spatial distribution of the oscillation, the first two must be excluded because the first should have produced a global response and this did not happen, as demonstrated by Schlesinger and Ramankutty (1994); the second, the white noise forcing, should have produced an oceanwide response with some kind of signature also in the SST field and this did not happen, as demonstrated by their and our analysis.

Based on these arguments, they suggested that the most probable cause of multidecadal oscillations could be an internal oscillation of the atmosphere–ocean system. If this is true we should be able to find some kind of correlation between multidecadal SST variations and some key climatic index, representative of the coupled ocean–atmosphere system in the Northern Hemisphere, also in the Mediterranean Sea. NAO is one of these indices.

Actually there are a variety of NAO indices available from several climate centers and various authors. Most of them are derived either from sea level pressure anomaly differences between two stations or from the principal components time series of the leading sea level pressure EOF. The station-based NAO indices have the advantage of going back to the midnineteenth century even if, due to the movement of the NAO centers through the year, they are only able to capture the NAO variability for a part of the year. Figure 12 reports winter NAO indices provided by three different sources: the University Corporation for Atmospheric Research (UCAR) NAO (Hurrell 1995), Climate Research Unit (CRU) (Hurrell 1995; Jones et al. 1997), and NAO index (NAOI) (Li and Wang 2003). The three indices show many similarities and some differences. More specifically, for frequencies around 10 yr or higher there is a complete correspondence between the three series, while lower frequencies have different amplitudes. In particular, the NAOI exhibits a higher-amplitude multidecadal oscillation that vaguely resembles the AMO one.

Li and Wang (2003) defined their NAOI as the difference in the normalized SLP regionally zonal averaged over the North Atlantic sector from 80°W to 30°E and between 35° and 65°N. They demonstrated that NAOI is able to describe well the large-scale NAO circulation features and suggested that this index can be used to characterize the NAO for all seasons. The following analysis is based on NAOI.

The correlation between NAOI and SST was computed for each grid point over the North Atlantic and Mediterranean Sea for each season and year using HadISST (Fig. 13). In Fig. 13 the SST and NAOI correlation pattern shows the typical tripole structure (Cayan 1992a,b; Sutton et al. 2000) with alternate significant anticorrelation in the subpolar region and in the tropical Atlantic Ocean and positive correlation in between. This corresponds to the typical SST pattern that accompanies the positive (negative) NAO regime with low (high) temperatures in the higher-latitude regions and south of 20°N and higher (lower) SSTs between 25° and 45°N (Flatau et al. 2003). A significant negative correlation is also observed in the eastern Mediterranean Sea but only in winter, while few grid points passed the test in the western Mediterranean Sea and the Adriatic Sea during autumn. Similar results are obtained using ERSST with relatively smoother maps (not shown). The only noteworthy difference was observed in autumn when the northern pole of the correlation maps did not pass the 90% confidence-level test.

On the basis of these results we focused our analysis on the winter [January–March (JFM)] correlation between NAOI and SST in the subpolar gyre (SPG), in the subtropical gyre (STG), and the Mediterranean Sea. The SPG and the STG SSTs have been computed by averaging the SST in all the grid points within the northern and southern centers of the Atlantic tripole, respectively. Figure 14 is based on the HadISST time series and shows the time variability of the winter SST in the three selected areas together with the AMO and the NAOI (note that the NAOI has been plotted with the sign changed to facilitate the interpretation). Similar calculations made considering the cold half of the year [December–April (DJFMA)] rather than JFM give very almost identical results (not shown).

As expected, the NAOI correlates (negatively) quite well with the three time series and with AMO for most of the period under investigation. All the relative minima and maxima of the SPG SST series correspond to equivalent NAOI peaks (with the exception of the 1925 NAOI peak) even though, after 1940, NAOI seems to lead SPG SST by 1–2 yr. Either NAOI and AMO or any of the three SSTs show evidence of a multidecadal oscillation with period of about half of the length of the time series.

At annual scales similar correlations also exist, as qualitatively suggested by the autocorrelation analysis of the SST time series (Fig. 10). A quantitative measure of the functional coupling between Mediterranean Sea SST and other North Atlantic time series is the estimate of the magnitude squared coherence (MSC). The MSC is a normalized cross-spectral density function and measures the strength of association and relative linearity between stationary processes on a scale from 0 to 1 (Wang et al. 2004). Given two functions *x*(*t*) and *y*(*t*) the MSC is defined as

where *P _{xy}*(

*f*) is the complex cross-spectral density, and

*P*(

_{xx}*f*) and

*P*(

_{yy}*f*) are the auto-spectral densities at frequency

*f.*The confidence limit of the estimated coherence at a given level

*α*(95% is

*α*= 0.05) has been calculated as , where

*K*is the number of tapers used for the MTM (Emery and Thomson 2001).

The MSC between the Mediterranean SST time series and every other time series of the North Atlantic Ocean, that is, NAOI, AMO, subpolar gyre region SST, and subtropical gyre region SST, is shown in Fig. 15. All the Atlantic indices are significantly coherent with the Mediterranean time series (confidence limit over 90%) for periods longer than about 40 yr. Moreover, a second coherence peak around 18 yr is present. Other higher-frequency peaks, exceeding the 95% confidence limit, are also observed for coherences between Mediterranean SST and the Atlantic indices.

In the multidecadal-frequency range, the coherence between NAO and Mediterranean SST never exceed the 99% confidence limit, oscillating around 95%. The AMO coherence with Mediterranean SST exceeds the 95% confidence limit for periods longer than ~85 yr and the 99% confidence limit for periods longer than 100 yr. A similar signal is visible also in the in the STG–Med plot even if less pronounced. The index showing the most pronounced peak of coherence with the Mediterranean at multidecadal time scales is the SPG, which passes the 99% confidence limit in the range between 40 and 55 yr. The two time series are also coherent at 95% confidence limit for longer periods. This means that low-frequency Mediterranean oscillations are more coherent with SPG SST than with AMO and NAO.

## 6. Conclusions and discussion

In this paper we analyzed the most extensive SST datasets for the Mediterranean Sea currently available to investigate low-frequency (multidecadal) oscillations and explored possible connections with other regions of the global ocean.

Before examining the time series we evaluated the robustness and reliability of the two reconstructed SST datasets (HadISST and ERSST.v3) over the Mediterranean Sea, finding that damping or other potential distortions due to interpolation are not dominant when yearly averages are used. On the contrary, larger differences between the two interpolated datasets have been found at monthly scales (±0.4°C). These differences exhibit a marked annual cycle, more pronounced prior to 1941, when different SST bias adjustments were applied. The same annual cycle was found comparing the three analyzed datasets (HadISST, ERSST.v3, and ICOADS 2.5) with recent satellite estimates for the period 1985–2005 validated against in situ observations. The analysis reveals that both the reconstructed datasets overestimate SST during summer by 0.3°–0.5°C and underestimate the temperature in winter, resulting in an annual cycle of ~0.8°C (Fig. 6). This cycle is obviously removed by the yearly average; therefore, the annual bias with respect to observations is of −0.03° and −0.13°C for HadISST and ERSST.v3, respectively. Consequently, the three time series are interchangeable and any of them can be used, depending upon convenience, for climatic studies when averaged over the year, while the monthly reconstructed time series should be used with more care.

Following the guidelines proposed by Ghil et al. (2002), we have applied the multitaper method and the singular spectral analysis to the yearly SST derived from the ERSST.v3 dataset, to investigate the low-frequency variability of the Mediterranean Sea. The harmonic analysis shows a low-frequency peak near 70 yr with confidence limits at 99% corresponding to a significant spectral band in the power spectrum and four significant high-frequency peaks (6.3, 3.9, 2.8, and 2.2 yr) that qualify as harmonics in the reshaping procedure (Fig. 9). The higher-frequency peaks are linked to the scale of variability of the ENSO quasi-biennial and low-frequency modes (Rasmusson et al. 1990; Jiang et al. 1995; Yiou et al. 2000). The lower-frequency band around 70 yr (very close to the AMO period) does not pass the reshaping test. Our spectral analysis suggests that the four higher-frequency peaks are likely to contain a harmonic phase-coherent oscillation with, at least, a 90% level of confidence; however, the same is not true for the low-frequency band, where the reshaped and the unreshaped spectra are identical.

However, the combined results of the MTM and SSA spectral analysis clearly indicate that the 70-yr low-frequency band is a significant feature of the Mediterranean spectrum. The analysis of the HadCM3 1000-yr control run confirmed the existence of a Mediterranean multidecadal oscillation and supports the results of the shorter observational time series analyses.

This Mediterranean multidecadal oscillation is not evident in the global SST time series, as also confirmed by our SSA and MTM analysis of the global SST, which did not display any multidecadal periodicity but suggested the presence of a longer than 70-yr-period oscillation. A multidecadal oscillation, corresponding to a period of 56 yr, appears only after removal of the long-term components (centennial and longer), showing that the multidecadal signal was masked by the longer-period trends. This result is in agreement with Schlesinger and Ramankutty (1994), who found a similar signal in the residual of the IPCC temperature dataset that includes surface temperature over land. Schlesinger and Ramankutty (1994) concluded that their multidecadal signal in the global time series is the statistical result of a regional long-time-scale oscillation that exclusively occurs in the Northern Hemisphere.

The analysis of the SST field over other regions of the global ocean confirmed that the North Atlantic Ocean is the only other area of the World Ocean that clearly exhibits the same oscillation observed in the Mediterranean (Fig. 10). Other ocean regions could also have multidecadal SST oscillations, even if they appear much less pronounced in the lagged autocorrelation analysis and visible in the SSA analysis only when applied to SST residuals. The latitudinal behavior of the amplitude of the low-frequency peak resulting from the harmonic analysis test in the North Atlantic Ocean shows that this multidecadal signal has maximum at subpolar latitudes (Fig. 11). This result is in agreement with the AMO spatial pattern (Parker et al. 2007) and with the hypothesis of Frankcombe et al. (2010) of the Arctic origin of the AMO pattern. This latitudinal variation of the amplitude of the harmonic signal demonstrates that multidecadal oscillations observed in the residual of the global SST time series are essentially due to the signal of the northern regions of the Atlantic, thus confirming the results of Schlesinger and Ramankutty (1994).

A quantitative measure of the functional coupling between Mediterranean Sea SST and other North Atlantic time series has been established by magnitude squared coherence analysis (Fig. 15), where winter Mediterranean SST is highly coherent (confidence limit over 90%) with the Atlantic indices for periods longer than about 40 yr. However, the coherence between winter Mediterranean SST and NAO oscillates around the 90%–95% confidence limit, while coherence with AMO exceeds the 95%–99% confidence limit for periods longer than 85–100 yr. Only SPG shows a pronounced peak of coherence with the Mediterranean at multidecadal time scales; the coherence passes the 99% confidence limit in the range between 40 and 55 yr and the 95% confidence limit for longer periods. These results show that the winter Mediterranean SST oscillations are more coherent with SPG than with NAO or AMO.

Since the Mediterranean and North Atlantic Oscillations are coherent at multidecadal scale, one should consider the existence of a physical mechanisms that connect the two areas. Dima and Lohmann (2007) have recently proposed a deterministic mechanism for AMO based on atmosphere–ocean–ice interaction, synthesized in Fig. 8 of their paper, which relates the multidecadal SST variability to the THC strengthening. They argued that the atmospheric response to the positive SST anomalies in the North Atlantic is represented by a SLP low over the SST pattern, including Eurasia and therefore the Mediterranean region. Afterward, the effect is transferred to the North Pacific through atmospheric teleconnections producing a SLP high, which determines a maximum pressure gradient over Fram Strait, increasing Arctic sea ice and freshwater export with a consequent THC reduction. This sequence is cyclically repeated producing the 70-yr oscillation.

The results of our data analysis, which identify the maximum of the 70-yr signal at subpolar latitudes (Fig. 11) and a limited signature of this oscillation in the Northwest Pacific, is consistent with this schema. However, the presence of a similar oscillation over the Mediterranean Sea is out of the Dima and Lohmann (2007) schema. The possible contribution of the Mediterranean outflow of saltier waters through Gibraltar strait in modulating the North Atlantic THC, suggested by several authors (Artale et al. 2006; Calmanti et al. 2006; Lozier and Stewart 2008), even if relatively small (e.g., Rahmstorf 1998), can be considered as a possible way to explain the correspondence of North Atlantic and Mediterranean multidecadal oscillation. Moreover, Lozier and Stewart (2008), in an analysis of historical hydrographic data in the Northeastern Atlantic, suggest a connection between the northward penetration of the Mediterranean Overflow Water (MOW) and the location of the subpolar front, leaving open the discussion of the possible effect of water mass transformation in the subpolar regions.

Recently, Frankcombe et al. (2010) proposed that the multidecadal variability of the North Atlantic is dominated by two main time scales (20–30 and 50–70 yr). The first has an ocean internal origin caused by the variability of the Atlantic meridional overturning circulation, while the 50–70-yr variability is related to atmospheric forcing and exchange process between Atlantic and Arctic Ocean. Under this new schema, the hypothesis of an atmospheric origin also for the Mediterranean 70-yr signal is reinforced, also taking into account the absence of lag between the AMO and Mediterranean SST.

Our analysis of the Mediterranean SST variability does not allow a definitive answer to the question about whether the forcing of the observed multidecadal signal has an atmospheric origin or whether there is some contribution due to the Mediterranean THC. Coupled ocean–atmosphere models could contribute in investigating the origin of the Mediterranean multidecadal oscillation and in separating the concurring contributions of the atmosphere and the Mediterranean thermohaline circulation. However, given the importance of the water mass exchanges through the strait of Gibraltar, it is very important to include the strait explicitly in the coupled model with a physical connection between the Mediterranean Sea and the Atlantic Ocean. The results of Sannino et al. (2009) clearly indicate the importance of explicitly including the Gibraltar connection also in coupled atmosphere–ocean models to investigate multidecadal fluctuations in the Mediterranean Sea and to correctly separate atmospheric and internal ocean contributions to such oscillations.

## Acknowledgments

First of all we thank the three anonymous reviewers and the editor for the critical review of the manuscript and for the helpful suggestions.

ERSST.v3 data provided by NOAA Satellite and Information Service–National Environmental Satellite, Data, and Information Servicee (NESDIS). (http://lwf.ncdc.noaa.gov/oa/climate/research/sst/ersstv3.php).

HadISST data provided by Met Office, Hadley Centre, British Atmospheric Data Centre (http://badc.nerc.ac.uk/data/hadisst/). ICOADS data provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado, from their Web site at http://www.esrl.noaa.gov/psd/. We thank Dr. Tom Smith for providing global bias adjustment data for prior to 1941 SSTs and Drs. Robert H. Evans, Viva Banzon, and Roberto Iacono for stimulating discussions.