Climate of the Past Tropical Pacific spatial trend patterns in observed sea level : internal variability and / or anthropogenic signature ?

In this study we focus on the sea level trend pattern observed by satellite altimetry in the tropical Pacific over the 1993–2009 time span (i.e. 17 yr). Our objective is to investigate whether this 17-yr-long trend pattern was different before the altimetry era, what was its spatio-temporal variability and what have been its main drivers. We try to discriminate the respective roles of the internal variability of the climate system and of external forcing factors, in particular anthropogenic emissions (greenhouse gases and aerosols). On the basis of a 2-D past sea level reconstruction over 1950– 2009 (based on a combination of observations and ocean modelling) and multi-century control runs (i.e. with constant, preindustrial external forcing) from eight coupled climate models, we have investigated how the observed 17-yr sea level trend pattern evolved during the last decades and centuries, and try to estimate the characteristic time scales of its variability. For that purpose, we have computed sea level trend patterns over successive 17-yr windows (i.e. the length of the altimetry record), both for the 60-yr long reconstructed sea level and the model runs. We find that the 2-D sea level reconstruction shows spatial trend patterns similar to the one observed during the altimetry era. The pattern appears to have fluctuated with time with a characteristic time scale of the order of 25–30 yr. The same behaviour is found in multicentennial control runs of the coupled climate models. A similar analysis is performed with 20th century coupled climate model runs with complete external forcing (i.e. solar plus volcanic variability and changes in anthropogenic forcing). Results suggest that in the tropical Pacific, sea level trend fluctuations are dominated by the internal variability of the ocean–atmosphere coupled system. While our analysis cannot rule out any influence of anthropogenic forcing, it concludes that the latter effect in that particular region is still hardly detectable.


Introduction
Long term sea level rise is a critical issue of the global climate change because of its potential negative impact on many coastal regions of the world . For this reason, it has been extensively studied in recent years. Since 1993, sea level is accurately monitored by satellite altimetry (using the Topex/Poseidon, Jason-1, Jason-2, ERS-1/2 and Envisat satellite missions) with high accuracy and global coverage. Recent studies based on these observations showed that sea level is rising at a global mean rate of 3.3 mm yr −1 since 1993 (e.g. Ablain et al., 2009;Cazenave and Llovel, 2010;Nerem et al., 2010). They showed as well that sea level does not rise uniformly but displays strong regional variability (see Fig. 1a). To highlight this regional variability, the uniform trend (global mean) of 3.3 mm yr −1 has been removed from Fig. 1a. In some regions such as the western Pacific, the North Atlantic around Greenland or the south-eastern Indian ocean, sea level rose up to 4 times faster than the global mean over 1993-2009. Meanwhile, other regions such as the eastern Pacific or the north-western Indian ocean show lower rates of sea level rise. These large deviations from the global mean trend suggest that in different parts of the world, low lying lands are not facing the same risk of sea level rise. Hence, when assessing the potential impacts of sea level rise, it is of primary importance to understand the time variability of observed regional sea level trend patterns and the causes which drive them.
A number of previous studies have shown that sea level trend patterns over the altimetry era mainly result from ocean temperature and salinity changes (e.g. Bindoff et al., 2007). This was evidenced by the comparison between altimetrybased and steric trend patterns deduced from in-situ hydrographic measurements (Ishii and Kimoto, 2009;Levitus, 2005;Levitus et al., 2009;Lombard et al., 2005a, b) and ocean general circulation models (OGCMs) outputs (Wunsch et al., 2007;Kohl and Stammer, 2008;Carton and Giese, 2008;Lombard et al., 2009). Analyses of in situ ocean temperature measurements showed that the thermosteric component is the most important contribution to the observed sea level regional variability. OGCM runs with or without data assimilation have confirmed that point. However salinity changes may also play some role at regional scale, partly compensating temperature effects, as shown by Wunsch et al. (2007) and confirmed by other studies (e.g. Kohl and Stammer, 2008;Lombard et al., 2009).
Other phenomena may also contribute to the regional variability in rates of sea level change. This is the case of gravitational and deformational effects of the solid Earth in response to mass redistributions of the last deglaciation and ongoing land ice melting Milne and Mitrovica, 2008;Milne et al., 2009;Mitrovica et al., 2001Mitrovica et al., , 2009). However, these effects are currently small and have not yet been detected in the altimetry-based sea level patterns.
Observations have shown that thermosteric spatial patterns are not stationary but fluctuate in time and space in response to driving mechanisms such as the ENSO (El Niño-Southern Oscillation), the NAO (North Atlantic Oscillation) and the PDO (Pacific Decadal Oscillation) (Levitus, 2005;Lombard et al., 2005a;Di Lorenzo et al., 2010;Lozier et al., 2010). Thus regional sea level trend patterns observed over the altimetry era were likely different prior to 1993.
The spatial and temporal variability of steric sea level change is tightly linked to complex ocean dynamics. It results from the redistribution of heat and fresh water both horizontally and vertically through sea-air fluxes and changes in the ocean circulation. Several studies have identified surface wind stress as the main driving mechanism of circulation-based heat and salt redistribution over the past few decades (e.g. Kohl and Stammer, 2008), in particular in the Indo-Pacific region (Timmermann et al., 2010;Han et al., 2010). But according to Kohl and Stammer (2008), surface fluxes and particularly buoyancy fluxes may have played an increasing role during the past 2 decades.
This important role of heat and fresh water redistribution in steric sea level trend patterns was previously noticed by Wunsch et al. (2007). These authors argued that given the long memory time of the ocean, observed patterns reflect an integration of the present forcing with internal changes and forcing that occurred in the past. This is another argument suggesting that the sea level trend patterns observed by satellite altimetry over the last 17 yr (the altimetry era) are not steady.
The purpose of the present study is to address two important scientific questions related to the sea level regional variability: (1) if sea level trend patterns are not stationary through time, how have they evolved during the last decades and what are their characteristic time scales? (2) What are the factors that drive them: are they mainly due to internal variability of the climate system or do they already reflect external forcing factors, in particular anthropogenic forcing? In this study we focus on the sea level spatial trend pattern of the tropical Pacific. To try to answer the above questions concerning this region, we use different observational data sets and climate model outputs. For the sea level observations we use altimetry data since 1993 and a new version of a past sea level reconstruction over 1950. For the climate model outputs, we use runs of eight coupled global climate models (CGCM here after) from the Coupled Model Intercomparison Project 3 database (hereafter CMIP3). These datasets are presented hereafter in more details.
Using the altimetry data set, we compute the observed spatial trend pattern over 1993-2009. Using the reconstructed 2-D sea level fields since 1950, we compute sea level trend patterns over successive 17-yr windows (the length of the altimetry data set) for the past 6 decades. The objective is to identify the dominant modes of variability and the characteristic life time of the 17-yr-long spatial trend pattern over the last 60 yr. Then we compute the sea level trend patterns over successive 17-yr windows from the CGCM multi-centennial control run outputs. These runs with constant, preindustrial external forcing give us an estimation of the modelled 17-yr sea level trend pattern produced by the internal variability of the climate system and its evolution with time. We check whether the tropical Pacific sea level trend pattern resembles the observed one or not. A similar analysis is done with the CGCM climate model runs over the 20th century. The latter runs include human-induced changes in atmospheric greenhouse gases and aerosols (as discussed later in this paper, some of these models do not take into account natural forcing such as solar or volcanic variability). Results are compared to those of the control runs. The issue is indeed to detect or not the signature of natural and anthropogenic forcing factors on the sea level trend patterns of the tropical Pacific.

Satellite altimetry sea level data (1993-2009)
We used altimetry-based 2-D sea level fields from AVISO (http://www.aviso.oceanobs.com/en/data/products/ sea-surface-height-products/global/index.html). The data consisted of gridded sea level anomalies at weekly intervals on a 1/4 • regular grid, from January 1993 to December 2009. We used the DT-MSLA "Ref" series computed at CLS (Collecte Localisation Satellite) by combining several altimeter missions, namely: Topex/Poseidon, Jason 1 and 2, Envisat and ERS 1 and 2. It is a global, homogenous, intercalibrated dataset based on global crossover adjustment (Le Traon and Ogor, 1998) using Topex/Poseidon and then Jason 1 as reference missions. It is corrected for the long wavelength orbit errors , ocean tides, and wet/dry troposphere and ionosphere (see Ablain et al., 2009 for more details). The inverted barometer (IB) correction has also been applied in order to minimize aliasing effects (Volkov et al., 2007) through the MOG2D barotropic model correction that includes the dynamic ocean response to short period (<20 day) atmospheric wind and pressure forcing and the static IB correction at periods above 20 day (see Carrere and Lyard, 2003 for details).

Two-dimensional past sea level reconstruction (1950-2009)
To determine the sea level trend pattern variability over the last decades (i.e. prior to the altimetry era), we updated the previous past sea level reconstruction developed by Llovel et al. (2009)  . Let us briefly summarize Llovel et al. (2009)'s reconstruction. The method, based on the reduced optimal interpolation described by Kaplan et al. (2000), consists of interpolating long tide gauge records with a time varying linear combination of spatial modes of a 2-D sea level field. The method has 2 steps. In the first step, an Empirical Orthogonal Function (EOF) decomposition (Preisendorfer, 1988;Toumazou and Cretaux, 2001) of 2-D sea level fields is performed (based on outputs of the OPA/NEMO ocean circulation model). This decomposition allows separating the spatially well-resolved signal of the model into spatial modes (EOFs) and their related principal components (PCs). The second step consists of computing new PCs over a longer period  covered by the selection of the considered 99 long tide gauge records. This is done through a least-squares optimal procedure that minimizes the differences between the reconstructed field and the tide gauge records at the tide gauge locations. Compared to the previous sea level reconstruction of Church et al. (2004), the originality of Llovel et al. (2009)'s reconstruction is the use of long-term sea level patterns (EOFs) deduced from a 44-yr long run of the OPA/NEMO ocean model, instead of the shorter altimetry record. This, in principle, allows better capturing of the decadal variability of the spatial trend patterns (see Llovel et al., 2009 for more details).
In the present study, we use a new version of the reconstruction based on three modifications. First, we followed Christiansen et al. (2010) and made use of a covariance matrix error and of the so-called EOF0 (i.e. the EOF mode corresponding to the geographically averaged but time-variable sea level, processed separately from the other EOF modes as in Church et al., 2004's reconstruction). This drastically improves the accuracy of the reconstructed trends over the reconstructed period. Second, on the basis of the lastest data available from the Permanent Service For Mean Sea Level (PSMSL: http://www.psmsl.org), we updated the tide gauge records used by Llovel et al. (2009)  the computation of the EOFs patterns, instead of using the OPA/NEMO model, which has a coarse resolution of 2 • on average, we preferred to use the ORCA025-B83 run of the DRAKKAR/NEMO model which has a higher resolution of 1/4 • . Indeed, Penduff et al. (2010) showed that ocean models with higher resolution bring substantial improvements in the representation of the sea level spatial variability, in particular at interannual time scales. The ORCA025-B83 run is based on the free surface ocean circulation model NEMO version 2.3 (Madec, 2008). This simulation is very close to the simulation ORCA025-G70, analysed and compared to satellite altimetry data by Lombard et al. (2009). It does not assimilate any observational data (e.g. satellite altimetry or in situ data) as in ORCA025-G70. It is forced by the realistic hybrid surface forcing "DRAKKAR forcing set 4.1" described in details by Brodeau et al. (2010). This forcing is based on the CORE dataset assembled by W. Large (Large and Yeager, 2004), the ECMWF (European Center for Medium-Range Weather Forecast) ERA 40 reanalysis (Uppala et al., 2005), and the ECMWF operational analyses for recent years. This new reconstruction has been validated at global scale  and in the tropical Pacific region  by comparison with independent tide gauge records not used in the reconstruction process. Figure 1b shows the reconstructed sea level trend patterns over 1950-2009 (uniform trend removed). As expected, regional sea level trend patterns reconstructed over the last 60 yr differ largely from those observed during the last 17 yr (see Fig. 1a for comparison).

Coupled Climate Model (CGCM) runs
Concerning the CGCM simulations, we analyzed both multicentennial control runs and runs covering the 20th century starting in the mid 19th century and ending in the 2000s (named picntrl and 20c3m runs respectively in the CMIP3 nomenclature).
The control runs and the 20th century runs differ by the external forcing. For the 20th century runs, external forcing includes changes in greenhouse gases, tropospheric and stratospheric ozone, anthropogenic stratospheric sulfates, black and organic carbon, volcanic aerosols, solar irradiance and the distribution of land cover. In the control runs, all external forcing variables are kept constant at their (preindustrial) 1860 values. These picntrl runs are intended to provide an estimate of the internal variability of the climate system. The 20c3m runs provide an estimate of the 20th century climate. They use observed, time-varying external forcing from around 1860 to 2000.
Among all CGCMs available in CMIP3 we have selected models that provide both sea surface temperature and sea level variables in their control run and 20th century run outputs. Models with less than 300 yr of control run were discarded. This led us to use a subset of 8 models: GFDL cm2.1, CNRM cm3, GISS model er, IAP fgoals g1.0, IPSL cm4, MIROC 3.2 medres, NCAR ccsm3 0, and UKMO hadcm3. The vertical and horizontal resolutions for the atmospheric and oceanic modules of each model are gathered in Table 1 along with the references in which the models are described in details. Note that for each of the selected models, one control run and one 20th century run at least were analyzed. When several 20th century runs were available (see Table 1), they were all analyzed one by one. Note as well that the selected models differ in their 20c3m run external forcing. For example, the IPSL 20c3m run does not include solar and volcanic variability. The IAP 20c3m does not include the volcanic variability. All other models include both solar and volcanic variability (see Table 1). But we decided to keep the IPSL and the IAP models in our selected subset because only few models provide the outputs necessary for this study. The final discussion will show that this does not impact our conclusions.
The sea surface height (SSH) fields given in the model outputs are incomplete and can not be directly compared to the observations in terms of global mean. Indeed the models do not contain the global mean steric sea level signal because they use the Boussinesq assumption that enforces the total ocean volume to remain constant. On the other hand, the models correctly simulate the regional sea level changes because the Boussinesq assumption has no impact on the latter (see Greatbatch, 1994). Hence the SSH output variable is well adapted to the regional analysis performed in this study.

Results
All computations were done on the basis of the monthly time series from altimetry, the 2-D past reconstruction and the CGCMs. In this study, we are interested in the regional variations of the sea level trends at inter-annual to multi-decadal time scales. Hence, two processings have been applied to the sea level time series prior to our analysis. (1) The global mean sea level trend was removed from the sea level time series because this study focuses on the regional variability around the global mean (this also removes any internal drift of the CGCM runs at the same time), (2) the time series were averaged on a yearly basis and filtered for the multicentennial signals to focus on the inter-annual to the multidecadal time scales. These low frequency signals were filtered out with a high-pass filter with a 80-yr cutoff. The high-pass filter built here is a fast Fourier transform convolution with a Hamming window cutting at 1/80 yr −1 (Brigham, 1974). The choice of the 80-yr filter cutoff enables us to keep signals with periods shorter than 70 yr, in which we are most interested (see further), and to ensure a reliable (somewhat conservative) filtering of the multi-centennial signals.

Observed 17-yr trend patterns over the tropical Pacific since 1950
Over the 17 yr of altimetry era, sea level trends show characteristic patterns (see Fig. 1a). The most prominent feature is the strong east-west dipole in the tropical Pacific region, positive in the western part and negative in the eastern part. This pattern has been persistent for several years now. It was already observed over the first 10 yr of the altimetry period (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003) by Cazenave and Nerem (2004) (see their Fig. 7). With the 2-D past sea level reconstruction, we can gain some insights on how the tropical Pacific sea level trend pattern evolved over longer time scales (here, up to 60 yr). We computed spatial trend patterns of the reconstructed sea level over successive 17-yr windows (i.e. the length of the altimetry data set). This gives as an output a set of 43 17-yr-long trend maps (starting in 1950 and shifting by one year the 17-yr time span, this provides 43 trend maps). Figure 2 shows 3 of these reconstructed 17-yr sea level trend patterns over three time spans: 1993-2009, 1976-1992 and 1959-1975. They all exhibit a strong ENSO-like dipole pattern. As expected, we note that the 1993-2009 reconstructed sea level pattern is very close both in shape and amplitude to the observed one. The 1959-1975 trend pattern is somehow different in shape but still exhibits a strong ENSO-like pattern with slightly lower amplitude than over 1993-2009. This indicates that a trend pattern similar to the presently observed one already existed in the late 1960s and the early 1970s. On the other hand, the 1976-1992 pattern is opposite to the 1993-2009 one almost everywhere. This suggests that the tropical Pacific trend pattern seems to have fluctuated with time and was opposite to the present one in the late 1970s and the 1980s.
To investigate this further, we performed an EOF decomposition of the set of 43 reconstructed 17-yr trend maps. On the basis of these EOFs, we computed the first rotated EOF (see von Storch and Zwiers, 1999) to obtain the EOF that maximizes the explained variance among the 43 17-yr-long trend maps. This EOF is presented in Fig. 3. It accounts for 37 % of the total variance. The second rotated EOF explains 25 % of the total variance and the third one 7 %. The leading EOF spatial pattern is very similar to the altimetry spatial pattern (Fig. 1a). This shows that the spatial pattern observed by satellite altimetry during the last 17 yr turns out to be the most frequently observed pattern among all the 17yr sea level trend maps of the last 60 yr (in terms of explained variance). Its PC indicates that this pattern has fluctuated with time (see Fig. 3). It was opposite to its current value in the 1970s and early 1980s and went back to values similar to what we observe now in the late 1960s (with a lower amplitude), as already suggested by Fig. 2. These fluctuations follow the low frequency variations of the extended NINO3 index (a proxy of El Niño; Kaplan et al., 1998). This index is the average of the sea surface temperatures (SST) over the NINO3 region (150 • W-90 • W, 5 • S-5 • N). Here it has been filtered with a 10-yr running mean (see Fig. 3). The correlation coefficient between the two curves shown in Fig. 3 is 0.63, with a significance level (SL) >99 %. This result suggests that the reconstructed sea level trend pattern fluctuation reflects a natural mode of the climate system: a low frequency modulation of ENSO.
Another way to analyse the reconstructed 17-yr-long sea level trend fluctuations of the tropical Pacific is to look at two regions defined by boxes A and B (see Fig. 2). Box A is located in the western Pacific (15 • N-15 • S by 120 • E-200 • E) Fig. 3. First rotated EOF of the set of 17-yr trend maps computed with the reconstruction over the tropical Pacific. It explains 37 % of the total variance. The PC is the black curve. The NINO3 index from Kaplan et al. (1998), filtered with a 10-yr running mean and detrended, is superimposed in blue. Their correlation coefficient is 0.63 (SL > 99 %). and box B in the eastern Pacific (15 • N-15 • S by 200 • E-280 • E). We computed the mean sea level (global mean sea level trend removed) in each box and compared it to the extended NINO3 index. The resulting curves are presented in Fig. 4a. The mean sea level in box A (western Pacific) covaries (but with opposite sign) with the mean sea level in box B (eastern Pacific). This confirms that the tropical Pacific sea level behaves as an east-west dipole that fluctuates with time. This fluctuation closely follows the ENSO mode of variability represented by the NINO3 index: the correlation between NINO3 index and eastern Pacific mean sea level is 0.73 (SL > 99 %).
We also computed the mean sea level trend in box B over successive 17-yr windows. Its time amplitude is displayed in Fig. 4b along with the NINO3 index time series smoothed with a 10-yr running mean. As suggested by the rotated EOF, this analysis confirms that the 17-yr trends in the tropical Pacific fluctuate with time following some low frequencies of the ENSO mode.
In order to isolate the low frequencies of the extended NINO3 index connected with the 17-yr trends fluctuations, we computed the NINO3 index power spectrum in Fig. 4c.  (Kaplan et al., 1998) power spectrum (in black) and the best fit AR2 process in grey as the null hypothesis. The grey dash lines indicate the 95 % confidence level around the null hypothesis The power spectrum of the best fit random process is also shown. It gives an estimation of what would be the power spectrum of a random process with a variability similar to the NINO3 index one (in terms of mean, variance and power spectrum). We added its 95 % confidence interval in Fig. 4c (it is the area of the power spectrum that is covered by 95 % of the randomly generated series). All frequency bands in which the NINO3 index power rises above the 95 % confidence level reveal the presence of a robust, deterministic ENSO frequency (at the 95 % confidence level), while the signal contained inside the confidence interval can be considered as undistinguishable from random fluctuations (null hypothesis).
The choice of an appropriate random process is a key issue. It has been done here in 2 steps. In the first step we chose the random process among simple random linear processes: the autoregressive models (AR) (this is the classical approach for climate records, see von Storch and Zwiers, 1999). Here we considered an AR process of order 2 because, according to the partial autocorrelation function of the NINO3 index, it appears to be indistinguishable from zero for lag greater than 2 (AR processes of a given order n are known to have a partial autocorrelation function that is indistinguishable from zero for lag greater than n, see von Storch and Zwiers, 1999). Moreover we computed an objective test for AR order (based on the Akaike information criterion, see von Storch and Zwiers, 1999) which indicates that an AR2 should best fit the data as well. In a second step, given the order of the AR process, we computed the parameters that best fit the NINO3 index using a least squares procedure. Figure 4c confirms that the AR2 process spectrum fits the NINO3 index spectrum well. Thus the extended NINO3 index is well represented by a linearly damped oscillator driven by white noise (AR2 model). It shows a range of preferred time scales centred around 4.3 yr (see the black arrow in Fig. 4c). Some spectral peaks significantly differ from the AR2 model. They can be found in the interannual band around the 3.7 yr and 5.8 yr periods and in the inter-decadal band around the 13.2 yr period. But none of these peaks can account for the low frequency modulation of ENSO identified earlier (with periods between 20 yr and 30 yr, see Fig. 4b). The NINO3 index shows actually some power in the 20 yr to 30 yr waveband (see Fig. 4c) but it remains indistinguishable from a random fluctuation either because it is of random origin or because the short length of the NINO3 record does not allow a good estimation of its power spectrum.
The time span covered by the sea level reconstruction is relatively short, covering only 60 yr. Nevertheless, during this time span, the 17-yr long spatial trend pattern observed by satellite altimetry in the tropical Pacific fluctuated with time, revealing periods during which sea level rise accelerated or decelerated. Presently, the sea level trend pattern is similar to what it was in the 1960s. This long-term fluctuation seems to be connected to some multi-decadal frequency of ENSO variability in the 20 yr to 30 yr waveband. But given the short length of the NINO3 index record, these low frequency ENSO fluctuations are not significant at the 95 % confidence level.

PIcntrl runs
The same strategy was applied to analyse the trend patterns of the tropical Pacific from CGCM control runs. We considered 17-yr windows to compute sea level trends from the multi-centennial models outputs. The first rotated EOF and the averaged 17-yr trends over boxes A and B were then derived from the resulting set of 17-yr trend maps following the www.clim-past.net/8/787/2012/ Clim. Past, 8, 787-802, 2012 method described in Sect. 4.1 for the reconstructed 2-D sea level fields. In Fig. 5 the spatial patterns of the first rotated EOFs for each CGCM control run are presented, while Fig. 6 shows the power spectra of their respective PCs. Note that, as for the rotated EOF in Fig. 3, the maximum of the PC has been normalised to 1 and the same colour scale has been used.
Looking at Fig. 5, we note that all CGCM control runs show spatial patterns with a magnitude comparable to reconstructed (Fig. 3) and observed fields (Fig. 1a). The magnitudes range from −12 mm yr −1 to +12 mm yr −1 . Furthermore, 6 models out of 8 (GFDL, CNRM, GISS, IAP, NCAR and UKMO) show a clear east-west dipole in the tropical Pacific as in the observations. But, in general, patterns differ in shape. In 4 models (GFDL, CNRM, NCAR, and UKMO), patterns are fairly similar to each other and exhibit some common features with the observed ones. They are dominated by a strong positive signal south-east of Papua-New Guinea and a more modest positive anomaly north of Papua-New Guinea that are seen in the observations as well (see Figs. 1a and 3 for comparison). For the NCAR and the UKMO runs, the positive trend anomalies extend slightly too far eastward in the equatorial Pacific. They exhibit as well negative trend anomalies east of the Philippines. The latter are also seen in the observations but are actually centred farther northward (around 20 • North).
The power spectra of the first PCs are shown in Fig. 6. Note that the IAP run is an exception: it is the only one showing a significant variance of its PC in the inter-annual waveband. Except for this run, all PCs concentrate their variance into a few peaks in the multi-decadal waveband. This indicates that in the CGCM control runs, the tropical Pacific 17-yr trends appear to fluctuate with time at multi-decadal time scales as in the observations. The GFDL and CNRM runs present a unique peak centred at 28 yr indicating that their spatial trend pattern fluctuate at time scales of 25-33 yr, as in the observations. The GISS, NCAR and UKMO runs show instead three peaks centred at 20 yr, 28 yr and 40 yr (up to 50 yr in the UKMO run).
To investigate this further, we analysed the 17-yr sea level trends averaged over boxes A and B for each model. As for the reconstructed sea level, the two time series exhibit a strong anti-correlation (not shown) indicating that in the control runs as well, 17-yr sea level trends in the tropical Pacific behave as an east-west dipole. Furthermore, the box B signal correlates well with the NINO3 index smoothed with a 10-yr running mean: the correlation coefficients are between 0.58 (for the CNRM cm3 control run) and 0.82 (for the IAP control run) (SL > 99 %).
We computed the power spectra of the NINO3 index and the box B 17-yr sea level trends of each control run (see Fig. 7). The long time period covered by the control runs allows capturing the multi-decadal variability of both signals, unlike with the reconstruction and the satellite altimetry observations. NINO3 indices of the control runs are shown in Fig. 7 together with their respective autoregressive (AR2) null hypothesis and 95 % confidence interval. The AR2 random processes have been chosen using the same procedure as described in Sect. 4.1. As for the observations, AR2 processes appear to properly fit the NINO3 index spectra. This indicates that for each CGCM control run also, the linear damped oscillator models (AR2) reproduce fairly well the NINO3 index in terms of power spectrum (for the NCAR run, the monotonic spectrum suggests that an AR1 process may have been sufficient). Only the GFDL, GISS and UKMO runs appear to have their range of preferred time scales centred at 4 yr as in the observations (see the black arrow on Fig. 4c) and a fairly good distribution of their variance between inter-annual and multi-decadal time scales. The GISS run shows a too low total variance. The CNRM, IAP, IPSL and NCAR runs have their range of preferred time scales centred around 2-3 yr instead of 4 yr, and show too large variance in the inter-annual timescales compared to the observations. The MIROC run shows too much variance in the multi-decadal time scales. It is interesting to note that most of the CGCM runs (except the IPSL and the IAP ones) show several peaks significant at the 95 % level of confidence in the multi-decadal waveband. But their central periods differ. In particular, the GFDL, CNRM, GISS, MIROC and UKMO runs agree on the presence of a significant peak at periods around 20-30 yr.
Spectra of the 17-yr sea level trends averaged over box B are shown in Fig. 7 along with their null hypothesis. The null hypothesis here is the assertion that the 17-yr box B sea level trends are indistinguishable from the 17-yr trends of the AR2 process that best fit the sea level in box B. This choice was driven by the same reasons as in Sect. 4.1. We also added the 95 % confidence in Fig. 7.
The spectra shown in Fig. 7 are not significant in the interannual waveband. At these time scales, the power spectra of the 17-yr trends' time series are actually dominated by the auto-correlation coming from the overlapping of the 17yr windows. For the multi-decadal waveband, we note that 6 CGCM runs (GFDL, CNRM, IAP, MIROC, NCAR and UKMO) show some significant peaks at the 95 % confidence level. They confirm that 17-yr sea level trends oscillate significantly at multi-decadal time scales in the tropical Pacific. Furthermore, except for the IAP run, each significant peak is associated with a significant peak of its respective NINO3 index spectrum (see Fig. 7), suggesting that for the majority of the control runs, the significant low frequency fluctuation of the 17-yr sea level trends at multi decadal time scales follows an ENSO-related low frequency modulation. This is confirmed by the squared coherence function of both signals (black curves on Fig. 7). This function, computed using Welch's overlapped averaged periodogram method (Rabiner and Gold, 1975), gives the coherence between two signals (value between 0 and 100 %). Figure 7 shows that for 5 models (GFDL, CNRM, MIROC, NCAR  and UKMO), the 17-yr signal in box B is very coherent with the NINO3 index at low frequency. In particular, at each frequency where a significant peak of both the 17-yr sea level trend and its respective NINO3 index can be found (see the grey vertical bars in Fig. 7), the coherence between both signals is very high: more than 70 % for the GFDL, CNRM, MIROC and UKMO models and 60 % for the NCAR model. This confirms a fairly strong relationship between both signals at these low frequencies.
In summary, CGCM control runs show results fairly similar to what was suggested by the reconstructed and the observed ones. Indeed, the tropical Pacific trends computed over successive 17-yr windows also show significant low frequency modulations in 6 out of 8 CGCM control runs. In 5 of them, the fluctuations appear to be tightly linked to significant ENSO low frequency modulation. Furthermore, 4 of these models (GFDL, CNRM, NCAR and UKMO) exhibit almost the same 17-yr spatial trend patterns and they display common features with the reconstruction and the satellite altimetry observations.
There are still some discrepancies between the CGCMs control runs and the reconstruction or the observations. The spatial trend pattern of the first rotated EOF of the GFDL, CNRM, NCAR and UKMO control runs differ from the reconstructed and the observed ones north of 15 • N (Fig. 3). The power spectra of the observed NINO3 index and the control runs NINO3 index (compare Figs. 4c and 7) differ as well. The observed NINO3 index does not show any significant multi-decadal variability at time scales superior to 15 yr while the NINO3 indices from control runs do. But discrepancies with the observations are not necessarily significant because the time period covered by the observations is much shorter than for the control runs. In effect, when trend patterns and power spectra are computed from observations they are likely dominated by high frequency or random features that would be smoothed out in the longer control runs. In Fig. 7. Power spectra of the box B sea level trends computed over successive 17-yr windows and NINO3 index from the coupled climate model control runs. Top panels: power spectra of box B 17-yr sea level trends (solid red line). They are shown together with the 17-yr trend of the best fit AR2 null hypothesis (solid black line) and the 95 % confidence levels(dashed grey line). Middle panels: power spectra of the NINO3 index (solid blue line), shown together with the best fit AR2 null hypothesis (solid black line) and the 95 % confidence levels(dashed grey line). Bottom Panels: magnitude squared coherence between the box B 17-yr sea level trends and their respective NINO3 index. long control runs we expect to see lower frequency features that would be missed by the shorter observation datasets. We verified this point in the GFDL control run. We split the 500yr GFDL control run into 8 independent 60-yr long samples (the time length covered by the reconstruction). For each of these, we computed the first rotated EOF from the 17-yr trend maps as done for the reconstruction in Sect. 4.1. Among the 8 resulting rotated EOFs, 6 were very similar to the rotated EOF computed over the whole control run time span (Fig. 5), while 2 were similar to the observed one (Fig. 3). In Fig. 8 we present one of them: the first rotated EOF computed from the 4th sample (years 220 to 280) of the GFDL control run. It explains 49 % of the total variance. To be fully consistent with the observations, Fig. 8 also shows the power spectrum of the NINO3 index computed over a 155 yr time period, i.e. the time period covered by the NINO3 index (it has been computed over the years 125-280 of the GFDL control run). Comparing Figs. 8 and 3, we note that the model strikingly resembles the observations in the 4th sample case. The spatial pattern of the first rotated EOF is very close to the observed one and the PC follows as well the NINO3 index with a quasi periodicity between 20-30 yr. The NINO3 index is also in better agreement with the data when computed over a 155 yr time period: the multi-decadal peak observed in Fig. 7 around 28-yr appears now under the 95 % confidence level as for the observations (Fig. 4c). This indicates that a 17-yr sea level trend variability similar to the reconstructed one can be found among the 60 yr long samples of a control run.

20c3m runs
In this section we consider the CGCMs 20th century runs (20c3m) to check whether any differences with the control runs can be found or not. Our approach is to consider that the best estimation of the internal modes of variability of the climate system are provided by the control runs and to use them as references against which we test the 20th century runs. The choice of the control runs as reference is motivated by the fact that control runs are multi-centennial runs unperturbed by any changes in the external forcing. Hence, they provide an estimation of the modes of variability of a preindustrial, unperturbed, steady climate. The point is to see whether 20th century runs make any significant differences with respect to a steady state climate in terms of 17-yr sea level trends (the null hypotheses here is the assertion that the 20th century runs are indistinguishable from control runs in terms of 17-yr sea level trend patterns and variability).
As in the previous analysis, we computed for each CGCM 20c3m run sea level trend maps over successive 17-yr windows and we performed two comparisons with the control Fig. 8. First rotated EOF of the set of 17-yr trend maps computed with the GFDL control run between years 220 and 280. It explains 43 % of the total variance. Top panel: spatial pattern of the EOF Middle panel: the PC is the black curve. The NINO3 index, filtered with a 10-yr running mean and detrended, is superimposed in blue. Bottom panel: NINO3 index power spectrum (in blue) and the best fit AR2 process in black as the null hypothesis. The grey dashed lines indicate the 95 % confidence level around the null hypothesis. runs: one with the 17-yr spatial trend patterns and one with the 17-yr averaged trend over box B.
For the first comparison, we considered the 17-yr spatial trend patterns of the rotated EOF of the control runs as the references. Then we projected the set of 17-yr trend maps computed from each 20c3m run outputs on these patterns (each 20c3m run outputs was projected on its respective control run spatial pattern). The resulting PCs show how the reference spatial patterns from the control runs fluctuate through the simulated 20th century climates. Their power spectra are plotted in Fig. 9. When several 20c3m runs were available for a given model, each PC was also plotted on the same graph. We added as well the power spectra of the control run leading PC (previously presented in Fig. 6). For a consistent comparison with the 20c3m power spectra, they were computed over subsets of 140 yr (length of the 20c3m runs) instead of the whole multi-centennial control run. We used a chunk spectral estimator (see von Storch and Zwiers, 1999) to compute them. This method consists in dividing the time series into a number of chunks of equal length and to estimate the spectrum by averaging the spectra obtained over each subset. This allows also estimating the 95 % confidence level. We added them in Fig. 9. But, note that here, only very few independent subsets (2 to 4) of 140 yr could be computed out of the 350 or 500 yr long control runs. So the 95 % confidence level plotted in Fig. 9 can not be considered as fully reliable but only indicative.
It is interesting to note that in Fig. 9, all CGCM 20c3m runs show some low frequency modulation of the spatial pattern as in the reconstruction and the control runs. Among models with complete external forcing (i.e. total anthropogenic forcing plus solar and volcanic variability, see Sect. 3.3 and Table 1), only 1 model (GISS) out of 6 shows significant differences in the leading PC between its control run and its 20c3m run. Three models (GFDL, MIROC and NCAR) show instead differences that are not significant. The 2 models left (CNRM and UKMO) actually show differences that reach the 95 % confidence level. But these 2 models only have 350 yr of control run. This enabled us to compute only 2 chunks to deduce the 95 % confidence level which makes their estimation particularly unreliable. So we assume that, for these 2 models, the differences are not significant. Concerning the IPSL model, which only takes into account the anthropogenic forcing (see Sect. 3.3 and Table 1), significant differences can be seen in the leading PC between the control run and the 20c3m run. But for the IAP model which includes anthropogenic forcing and solar variability (but no volcanic forcing), no significant differences are observed. Finally there are a majority of models (6 out of 8) which do not show significant difference in 17-yr sea level trend patterns between their control run and their 20th century runs. In particular, 5 out of the 6 models which include complete external forcing do not show significant differences.
In the second comparison, we considered the 17-yr sea level trends averaged over box B in the 20c3m runs. Their power spectra are plotted in Fig. 10 along with the power spectra of the 17-yr averaged trends in box B computed from the control runs. The latter power spectra (from the control runs) were already presented in Fig. 7. Here the same power spectra were computed but over subsets of 140 yr (length of the 20c3m runs) instead of the whole multi-centennial control run in order to be consistent with the 20c3m power spectra. We used a chunk spectral estimator to do so. There are too few independent subsets of 140 yr among the   10. Power spectra of the box B sea level trend computed over successive 17-yr windows and NINO3 index from the coupled climate models. The red curves indicate the power spectra from the 20c3m runs. The orange curves indicate the power spectra from their respective control runs and the black curves their respective best fit AR2 process. The grey dashed lines indicate the 95 % confidence level computed from the control runs best fit AR2 processes. multi-centennial control runs to get a reliable estimation of a confidence interval on the control run spectra. So, instead of using the control run spectra as a reference (null hypothesis), we prefered to use their best fit AR2 process computed earlier and shown in Fig. 7. This gives more reliable estimations of the confidence intervals. The AR2 process power spectra and their 95 % confidence level are plotted in Fig. 10. Among models with complete external forcing, 2 models out of 6 (CNRM and UKMO) show a significant difference between the power spectra computed from their control run and their 20c3m run. 3 models (GFDL, GISS and NCAR) show both some 20c3m runs that differ significantly from their control run and some which do not. For each of these 3 models, only one run among all available 20c3m runs shows significant differences (i.e. 1 out of 4 available for the GFDL model, 1 out of 8 for the GISS model and 1 out of 2 for the NCAR model). Finally one model (MIROC) does not show any significant differences between its 20c3m run and its control run.
In this comparison, the IPSL (with anthropogenic forcing only) does not show significant differences between its control run and its 20c3m run, while the IAP model (whose 20c3m run contains greenhouse gas emissions and solar variability) shows some. Finally, among the 19 20c3m runs available from our dataset (see Table 1), 13 of them do not show any significant differences to their respective control runs (this ratio increases to 12 out of 17 when considering only 20c3m runs with complete external forcing). Only 6 20c3m runs show some differences. In total, only 5 20c3m runs show differences if we consider only those with complete external forcing. Nevertheless, it is interesting to note that for each of these 5 runs, the power spectra reveal that the 17-yr trends fluctuations in the 20c3m runs present a peak at higher frequency with more variance than in their respective control run.
Previous studies have shown that the sea level trend pattern observed in the tropical Pacific through 17-yr long precise altimetry observations (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)) is largely of thermal origin (Ishii and Kimoto, 2009;Levitus, 2005;Levitus et al., 2009;Lombard et al., 2005a, b), the thermosteric sea level trends being themselves driven by surface wind stress (Carton et al., 2005;Kohl and Stammer, 2008;Timmermann et al., 2010). In this region, observed and thermosteric sea level trends are tightly linked to the ENSO mode of variability known to occur on a broad range of time scales, from interannual to multi-decadal (Knutson and Manabe, 1998;Lau and Weng, 1999;Vimont et al., 2002;Vimont, 2005). So the 17-yr trends observed by satellite altimetry in the tropical Pacific are expected to fluctuate on these time scales. Nevertheless, only a few studies have suggested that sea level trends observed prior to the altimetry era changed with time (Wunsch et al., 2007;Khol and Stammer, 2008). Here we show that 17-yr sea level trends are un-steady and fluctuated in the past with a characteristic time scale of around 25 yr (see Sect. 4.1). On the basis of a past sea level reconstruction , we find that the tropical Pacific trend patterns over successive 17-yr windows fluctuate in time and show some periods during which sea level rise accelerates or decelerates (or equivalently, trend patterns display increasing or decreasing intensity). The trend pattern behaves as a fluctuating east-west dipole following a low frequency modulation of ENSO. The relatively short time span of the 60-yr long reconstruction makes it difficult to precisely determine the characteristic time-scales of the pattern fluctuations but some value around 25 yr is suggested by the observations.
The CGCM control runs with constant, preindustrial external forcing, show fairly similar results. Indeed, 4 out of 8 CGCM control runs (GFDL, CNRM, NCAR, and UKMO) show significant 17-yr sea level trend fluctuations tightly linked to significant ENSO low frequency modulations as well. They display the same 17-yr spatial trend pattern which differs slightly from the reconstructed one. We have shown with the GFDL control run that this is probably due to the different time periods covered by the reconstruction and the control runs. So, the internal variability of the climate system, simulated here by the CGCM control runs, seems to be well able to explain most of the sea level trend pattern fluctuations of the tropical Pacific seen in the reconstruction and observed by altimetry. Note nevertheless that the different control runs do not always agree on the characteristic periods of these fluctuations. The GFDL and CNRM control runs exhibit periodicities in the range 25-33 yr close to what is suggested by the reconstruction. But the NCAR and MIROC control runs on the one hand and the UKMO on the other hand show periodicities in the range 18-22 yr and in the range 50-60 yr, respectively.
The CGCM 20th century runs with external forcing (including anthropogenic forcing) show similar sea level trend behaviours for the tropical Pacific as in the control runs. Actually, a majority of the 20c3m runs which include a complete external forcing (12 out of 17) do not show any significant differences to their respective control run, either in terms of temporal or spatial structures. Consequently, because the 20c3m and control runs provide similar results, we conclude that the internal variability of the climate system is still the dominant contributor to the fluctuations of the observed 17yr spatial trend pattern in the tropical Pacific. In effect, our analysis does not detect any clear signature of external forcing, whether of anthropogenic origin or of natural origin (solar and volcanic variability). In other words, over the short altimetry record (17 yr), the amplitude of the noise represented here by the internal climate variability is so strong in the tropical Pacific that it prevents us from detecting the signal of anthropogenic forcing on the regional variability of the sea level trends in this region. These conclusions are also true when considering windows shorter than 17-yr for both data and model outputs. For example, using 10-yr and 15-yr-long windows led to similar results (with the same dominant low frequency variability of the tropical Pacific trend pattern).
Nevertheless, a minority of 20th century runs (see Sect. 4.2) seem to suggest that the impact of the external forcing could be possibly seen in the characteristic frequency band of the 17-yr trend pattern fluctuations. Several 20th century runs indeed show higher frequency oscillations for the 17-yr trend pattern than their respective control run (as explained earlier). But this remains very unclear. To get a clearer picture, we would need both more 20th century runs (to perform statistics on how many runs support this assertion) and longer runs with external forcing to compute more accurate power spectra. These runs should be available within the CMIP5 project. This will be the subject of future investigations.