variability is less sensitive to model resolution. This reaffirms that in our simulations, mesoscale eddies and nonlinear processes are likely to play a major role in producing intrinsic variability, and suggests that they are weakly involved in the SLA response to direct atmospheric forcing (which is probably linear at first order). The chaotic behavior of this eddy-driven intrinsic component may explain why temporal correlations between AVISO and fully forced simulations decrease as resolution increases, both in the present study and in Penduff P10.

Our results have several implications regarding the interannual variability simulated by ocean models at various resolutions. As seen in the $T - 2^\circ \rightarrow T - 1/4^\circ$ transition, resolving eddies 1) brings the distribution and magnitude of the total interannual variability very close to observational estimates, 2) does not modify much the large-scale distribution and magnitude of the forced variability but 3) strongly enhances the intrinsic component. Since the $T - 1/4^\circ$ simulation compares much better than $T - 2^\circ$ with AVISO observations, it is likely that the $1/4^\circ$ diagnostics presented here have some relevance for the real ocean: a substantial part of the real interannual SLA variability is probably generated by the ocean under the sole influence of mean and seasonal atmospheric forcing.

We also showed that the extratropical variability simulated by our Intergovernmental Panel on Climate Change (IPCC)-class coarse-resolution model driven by the same interannual forcing amounts to only a fraction of the observed variability, mostly because the intrinsic variability generated by laminar dynamics (and parameterized eddy fluxes) is very weak. In other words, laminar ocean models are much less influenced by chaotic variability than eddying models. At least within regions with a substantial intrinsic variability, the predictability of interannual variability (i.e., in agreement with the phase of observed events) may therefore be degraded as the resolution increases in ocean-only high-resolution forecasts (like in our $T - 1/4^\circ$ hindcast). This advocates for ensemble high-resolution ocean forecasts.

Berloff and McWilliams (1999) and others have shown in idealized setups that increasing Reynolds numbers (via, e.g., smaller viscosities and/or finer resolutions) tend to enhance the magnitude and chaotic character of the intrinsic variability; is there evidence of such a sensitivity within the existing range of eddying OGCM simulations? “Eddy resolving” (e.g., $1/10^\circ$) OGCMs generally produce more realistic mean states and more mesoscale activity than their “eddy admitting” (e.g., $1/4^\circ$) counterparts (e.g., Hurlburt and Hogan 2000; Smith

et al. 2000; Treguier et al. 2005). At interannual time scales in the Kuroshio region, Taguchi et al. (2007)'s $1/10^\circ$ global simulations and our $1/4^\circ$ counterparts do not differ much regarding the magnitude of their intrinsic variabilities, total variabilities, and ratios between both. This regional resemblance suggests that these $1/4^\circ$ and $1/10^\circ$ models capture comparable features of total and intrinsic interannual variabilities in this particular region, although this agreement cannot be assessed at global scale (Taguchi et al.'s study is regional). It is also likely that the detailed structure of the intrinsic variability depends on the model considered. Investigating these dependencies to model parameters and resolution remains to be done at regional and global scales.

The absence of air-sea coupling allowed us to separate forced and intrinsic ocean variabilities. Brightly shaded regions in the top panel in Fig. 8 locate many areas where the low-frequency SLA variability is mostly intrinsic. One may expect that in such regions, low-frequency intrinsic SLA fluctuations have signatures in SST that could then force low-frequency atmospheric fluctuations. At least in these regions, low-frequency atmospheric variability might therefore be driven by slowly varying SST changes of oceanic origin, which are not resolved in coarse-resolution ocean models. Western boundary current extensions (Berloff and McWilliams 1999) and large areas of the Southern Ocean are good candidates for such a scenario. Identifying the implications of our results for the fully coupled climate system is left for future studies.

Like the forced ocean variability (Sverdrup-type responses, Rossby and Kelvin waves, etc.), the intrinsic variability has subsurface, multivariate signatures [see idealized studies by, e.g., Spall (1996); Hogg and Blundell (2006)] that we ignored, mostly because observations with the adequate space-time span are restricted to the surface. Although certain important climate indices like the meridional overturning circulation (MOC) (Bingham and Hughes 2009), the transport of major currents (Kelly 1991) or overflows (Hoyer and Quadfasel 2001) have a signature in SLA, we did not investigate their intrinsic variations since we lack adequate observations. Dedicated observational and modeling studies are thus needed to monitor and understand the intrinsic variability of subsurface fields and climate-relevant quantities. As also stated by Berloff and McWilliams (1999) in their conclusion, the (intrinsic) low-frequency variability is poorly known in the real ocean; quantitative estimates about its manifestations and magnitude in the real ocean requires more “realistic” simulations than those used in process studies. The present study may be seen as such an attempt to provide process-oriented studies with
a more realistic context, in close connection with real
observations.

*Acknowledgments.* Suggestions made by Pavel Berloff and an anonymous reviewer led to significant improvement of this paper. This work is a contribution of the Drakkar project, which is supported by the Centre National d'Etudes Spatiales (CNES) through the Ocean Surface Topography Science Team (OST/ST), by the Centre National de la Recherche Scientifique (CNRS), the Institut National des Sciences de l'Univers (INSU), and the Groupe Mission Mercator Coriolis (GMMC). Partial support of the European Commission under Grant Agreement FP7-SPACE-2007-1-CT-218812-MYOCEAN is also gratefully acknowledged. Computations presented in this study were performed at Institut du Développement et des Ressources en Informatique Scientifique (IDRIS) and the Centre Informatique National de l'Enseignement Supérieur (CINES). The T − 1/4° (ORCA025-MJM01) simulation was performed at CINES as part of the Grands Challenges GENCI/CINES 2008. Most of this research was done during TP's stay at the Department of Oceanography at FSU, now named the Earth, Ocean and Atmospheric Science Department. TP wishes to acknowledge a number of interesting conversations with Peter Rhines, Michael Ghil, Greg Holloway, Pavel Berloff, and Julie Deshayes.

APPENDIX A

Periods

We aim to illustrate from SLA fields some relation-
ships between the mesoscale activity (HF band) and the
intrinsic low-frequency variability (LF band) in the *I* −
1/4° simulation. This requires the availability of LF and
HF time series simultaneously, which only occurs after
300 yr of spinup (only monthly *T* − 1/4° outputs are
available before year 300). We verified that the de-
trended interannual variabilities σ<sup>*I*−1/4°</sup><sub>LF</sub> computed over
23 successive 12-yr periods vary moderately throughout
the last 276 yr of *I* − 1/4°: the standard deviation of this
23-member time series remains below 10% of its mean,
peaking at 20% at 30°N. In other words, comparing
detrended SLA variabilities in *I* − 1/4° and *T* − 1/4° over
12-yr periods chosen after different spinup durations
(301 and 35 yr, respectively) is consistent, and also al-
lows a comparison between HF and LF variabilities in
*I* − 1/4°. In summary, *T* − 1/4° statistics are computed over
1993–2004; they are compared in section 3 to AVISO
statistics computed over the same period and in section 4
to *I* − 1/4° statistics computed over years 302–13.

APPENDIX B

Detrending

Our main aim is to compare observed, total, intrinsic,
and forced LF variabilities of SLA at low and high fre-
quency over the same range of available time scales, that
is, over the 12-yr period that is common to AVISO and
interannual runs. Linear trends were removed from all
LF collocated (observed and simulated) local SLA time
series before computing and comparing statistics (stan-
dard deviations, correlations). Indeed, some components
of local trends may differ within our seven-member SLA
dataset for three main reasons. First, part of these trends
might simply come from the truncation of interdecadal
oscillations, which are present in AVISO, but are likely
affected by the limited duration of model spinups. Sec-
ond, the potentially intrinsic (i.e., chaotic) component of
interdecadal variabilities may induce SLA fluctuations
with different phases in AVISO and each simulation,
hence yield arbitrary trends in individual 12-yr time
series and contaminate statistics. Third, the 12-yr
truncation of long-term, forced climate trends present
in AVISO may be altered by spurious SLA drifts in
forced simulations because freshwater fluxes are
uncertain (Béranger et al. 2006). In summary, we re-
moved the linear trends from all LF time series be-
cause their origins in simulations and in observations
may be different, hence cannot be compared consis-
tently. Beyond these reasons, removing trends also
focuses our analysis on a range of time scales that is
“resolved” in every 12-yr time series. Note, however,
that this removal had only a small quantitative impact
on our results.<sup>B1</sup>

REFERENCES

Barnier, B., and Coauthors, 2006: Impact of partial steps and momentum advection schemes in a global circulation model at eddy-permitting resolution. *Ocean Dyn.*, **56**, 543–567, doi:10.1007/s10236-006-0082-1.

Béranger, K., B. Barnier, S. Gulev, and M. Crépon, 2006: Comparing 20 years of precipitation estimates from different sources over the world ocean. *Ocean Dyn.*, **56**, 104–138.

Berloff, P. S., and J. C. McWilliams, 1999: Large-scale, low-frequency variability in wind-driven ocean gyres. *J. Phys. Oceanogr.*, **29**, 1925–1949.

<footnote>B1 P10 chose not to detrend SLA time series because their four
simulations were forced identically and were compared to AVISO
after the same period of spinup. This choice explains the small
differences found between the dashed lines in Fig. 5 and their
counterpart (black and magenta lines) in P10’s Fig. 6.</footnote>
—, W. Dewar, S. Kravtsov, and J. C. McWilliams, 2007: Ocean eddy dynamics in a coupled ocean–atmosphere model. *J. Phys. Oceanogr.*, **37**, 1103–1121.

Biastoch, A., C. Böning, and J. Lutjeharms, 2008: Agulhas leakage dynamics affects decadal variability in Atlantic overturning circulation. *Nature*, **456**, 489–492.

Bingham, R. J., and C. W. Hughes, 2009: Signature of the Atlantic meridional overturning circulation in sea level along the east coast of North America. *Geophys. Res. Lett.*, **36**, L02603, doi:10.1029/2008GL036215.

Brodeau, L., B. Barnier, T. Penduff, A.-M. Treguier, and S. Gulev, 2010: An ERA40-based atmospheric forcing for global ocean circulation models. *Ocean Modell.*, **31**, 88–104.

Cabanes, C., T. Huck, and A. Colin de Verdière, 2006: Contributions of wind forcing and surface heating to interannual sea level variations in the Atlantic Ocean. *J. Phys. Oceanogr.*, **36**, 1739–1750.

Cessi, P., and S. Louazel, 2001: Decadal oceanic response to stochastic wind forcing. *J. Phys. Oceanogr.*, **31**, 3020–3029.

Chang, K.-I., M. Ghil, K. Ide, and C.-C. A. Lai, 2001: Transition to aperiodic variability in a wind-driven double-gyre circulation model. *J. Phys. Oceanogr.*, **31**, 1260–1286.

de Ruijter, W., 1982: Asymptotic analysis of the Agulhas and Brazil Current systems. *J. Phys. Oceanogr.*, **12**, 361–373.

Dewar, W., 2001: On ocean dynamics in midlatitude climate. *J. Climate*, **14**, 4380–4397.

—, 2003: Nonlinear midlatitude ocean adjustment. *J. Phys. Oceanogr.*, **33**, 1057–1082.

Dijkstra, H. A., and M. Ghil, 2005: Low-frequency variability of the large-scale ocean circulation: A dynamical systems approach. *Rev. Geophys.*, **43**, RG3002, doi:10.1029/2002RG000122.

Eden, C., and H. Dietze, 2009: Effects of mesoscale eddy/wind interactions on biological new production and eddy kinetic energy. *J. Geophys. Res.*, **114**, C05023, doi:10.1029/2008JC005129.

Fyfe, J. C., and O. A. Saenko, 2006: Simulated changes in the extratropical Southern Hemisphere winds and currents. *Geophys. Res. Lett.*, **33**, L06601, doi:10.1029/2005GL025332.

Gershunov, A., and T. Barnett, 1998: Interdecadal modulation of ENSO teleconnections. *Bull. Amer. Meteor. Soc.*, **79**, 2715–2726.

Hall, N., B. Barnier, T. Penduff, and J.-M. Molines, 2004: Interannual variation of Gulf Stream heat transport in a high resolution model forced by reanalysis data. *Climate Dyn.*, **23**, 341–351.

Hazeleger, W., and S. Drijfhout, 2000: A model study on internally generated variability in subtropical mode water formation. *J. Geophys. Res.*, **105**, 13 965–13 979.

Hogg, A. M., and J. Blundell, 2006: Interdecadal variability of the Southern Ocean. *J. Phys. Oceanogr.*, **36**, 1626–1645.

Hoyer, J. L., and D. Quadfasel, 2001: Detection of deep overflows with satellite altimetry. *Geophys. Res. Lett.*, **28**, 1611–1614.

Hurlburt, H. E., and P. J. Hogan, 2000: Impact of 1/8° to 1/64° resolution on Gulf Stream model-data comparisons in basin-scale subtropical Atlantic Ocean models. *Dyn. Atmos. Oceans*, **32**, 283–330.

Jiang, S., F.-F. Jin, and M. Ghil, 1995: Multiple equilibria, periodic, and aperiodic solutions in a wind-driven, double-gyre, shallow-water model. *J. Phys. Oceanogr.*, **25**, 764–786.

Kelly, K. A., 1991: The meandering Gulf Stream as seen by the Geosat altimeter: Surface transport, position, and velocity

variance from 73° to 46°W. *J. Geophys. Res.*, **96**, 16 712–16 738.

LaCasce, J. H., and P. E. Isachsen, 2010: The linear models of the ACC. *Prog. Oceanogr.*, **84**, 139–157.

Madec, G., and the NEMO Team, 2008: NEMO ocean engine. Institut Pierre-Simon Laplace Note du Pôle de Modélisation 27.

Orsi, A. H., T. Whitworth III, and W. D. Nowlin Jr., 1995: On the meridional extent and fronts of the Antarctic Circumpolar Current. *Deep-Sea Res.*, **42**, 641–673.

Pacanowski, R., 1987: Effect of equatorial currents on surface stress. *J. Phys. Oceanogr.*, **17**, 833–838.

Penduff, T., B. Barnier, W. K. Dewar, and J. J. O'Brien, 2004: Dynamical response of the oceanic eddy field to the North Atlantic Oscillation: A model-data comparison. *J. Phys. Oceanogr.*, **34**, 2615–2629.

—, M. Juza, L. Brodeau, G. C. Smith, B. Barnier, J.-M. Molines, A.-M. Treguier, and G. Madec, 2010: Impact of global ocean model resolution on sea-level variability with emphasis on interannual time scales. *Ocean Sci.*, **6**, 269–284.

Philander, S. G. H., 1990: *El Niño, La Niña and the Southern Oscillation*. Academic Press, 289 pp.

Pierini, S., 2006: A Kuroshio Extension system model study: Decadal chaotic self-sustained oscillations. *J. Phys. Oceanogr.*, **36**, 1605–1625.

Qiu, B., 2003: Kuroshio Extension variability and forcing of the Pacific decadal oscillations: Responses and potential feedback. *J. Phys. Oceanogr.*, **33**, 2465–2482.

Ridgway, K. R., and J. R. Dunn, 2007: Observational evidence for a Southern Hemisphere oceanic supergyre. *Geophys. Res. Lett.*, **34**, L13612, doi:10.1029/2007GL030392.

Rintoul, S. R., C. W. Hughes, and D. Olbers, 2001: The Antarctic Circumpolar Current system. *Ocean Circulation and Climate: Observing and Modelling the Global Ocean*, G. Siedler, J. Church, and J. Gould, Eds., International Geophysics Series, Vol. 77, Academic Press, 271–302.

Saji, N., B. Goswami, P. Vinayachandran, and T. Yamagata, 1999: A dipole in the tropical Indian Ocean. *Nature*, **401**, 360–363.

Sallée, J. B., K. Speer, and R. Morrow, 2008: Southern Ocean fronts and their variability to climate modes. *J. Climate*, **21**, 3020–3039.

Schott, F. A., and J. P. McCreary Jr., 2001: The monsoon circulation of the Indian Ocean. *Prog. Oceanogr.*, **51**, 1–123.

Scott, R. B., and B. K. Arbic, 2007: Spectral energy fluxes in geostrophic turbulence: Implications for ocean energetics. *J. Phys. Oceanogr.*, **37**, 673–688.

Simonnet, E., M. Ghil, K. Ido, and R. Temam, 2003: Low-frequency variability in shallow-water models of the wind-driven ocean circulation. Part II: Time-dependent solutions. *J. Phys. Oceanogr.*, **33**, 729–752.

—, —, and H. A. Dijkstra, 2005: Homoclinic bifurcations in the quasi-geostrophic double-gyre circulation. *J. Mar. Res.*, **63**, 931–956.

Smith, R., M. Maltrud, F. Bryan, and M. Hecht, 2000: Numerical simulation of the North Atlantic Ocean at 1/10°. *J. Phys. Oceanogr.*, **30**, 1532–1561.

Smith, T. M., 2000: Tropical Pacific sea level variations (1948–98). *J. Climate*, **13**, 2757–2769.

Spall, M., 1996: Dynamics of the Gulf Stream/deep western boundary current crossover. Part II: Low-frequency internal oscillations. *J. Phys. Oceanogr.*, **26**, 2169–2182.

Taguchi, B., S.-P. Xie, N. Schneider, M. Nonaka, H. Sasaki, and Y. Sasai, 2007: Decadal variability of the Kuroshio Extension:
Observations and an eddy-resolving model hindcast. *J. Climate*, **20**, 2357–2377.

Treguier, A.-M., S. Theetten, E. Chassignet, T. Penduff, R. Smith, L. Talley, J. O. Beismann, and C. Böning, 2005: The North Atlantic subpolar gyre in four high-resolution models. *J. Phys. Oceanogr.*, **35**, 757–774.

von Storch, H., and F. W. Zwiers, 1999: *Statistical Analysis in Climate Research*. Cambridge University Press, 484 pp.

Webb, D. J., and B. A. de Cuevas, 2003: The region of large sea surface height variability in the southeast Pacific Ocean. *J. Phys. Oceanogr.*, **33**, 1044–1056.

Wirth, A., J. Willebrand, and F. Schott, 2002: Variability of the Great Whirl from observations and models. *Deep-Sea Res. II*, **49**, 1279–1295.

Zebiak, S., 1993: Air–sea interaction in the equatorial Atlantic region. *J. Climate*, **6**, 1567–1568.
Intrinsic Variability of Sea Level from Global 1/12° Ocean Simulations: Spatiotemporal Scales

GUILLAUME SÉRAZIN

CNRS, LGGE (UMR5183), Grenoble, and, Sciences de l'Univers au CERFACS, CERFACS/CNRS, URA1857,
Université Paul Sabatier, Toulouse, France

THIERRY PENDUFF, SANDY GRÉGORIO, BERNARD BARNIER, AND JEAN-MARC MOLINES

CNRS, Université Grenoble Alpes, LGGE (UMR5183), Grenoble, France

LAURENT TERRAY

Sciences de l'Univers au CERFACS, CERFACS/CNRS, URA1857, Toulouse, France

(Manuscript received 28 July 2014, in final form 3 December 2014)

ABSTRACT

In high-resolution ocean general circulation models (OGCMs), as in process-oriented models, a substantial amount of interannual to decadal variability is generated spontaneously by oceanic nonlinearities: that is, without any variability in the atmospheric forcing at these time scales. The authors investigate the temporal and spatial scales at which this intrinsic oceanic variability has the strongest imprints on sea level anomalies (SLAs) using a 1/12° global OGCM, by comparing a “hindcast” driven by the full range of atmospheric time scales with its counterpart forced by a repeated climatological atmospheric seasonal cycle. Outputs from both simulations are compared within distinct frequency–wavenumber bins. The fully forced hindcast is shown to reproduce the observed distribution and magnitude of low-frequency SLA variability very accurately. The small-scale (L < 6°) SLA variance is, at all time scales, barely sensitive to atmospheric variability and is almost entirely of intrinsic origin. The high-frequency (mesoscale) part and the low-frequency part of this small-scale variability have almost identical geographical distributions, supporting the hypothesis of a nonlinear temporal inverse cascade spontaneously transferring kinetic energy from high to low frequencies. The large-scale (L < 12°) low-frequency variability is mostly related to the atmospheric variability over most of the global ocean, but it is shown to remain largely intrinsic in three eddy-active regions: the Gulf Stream, Kuroshio, and Antarctic Circumpolar Current (ACC). Compared to its 1/4° predecessor, the authors’ 1/12° OGCM is shown to yield a stronger intrinsic SLA variability, at both mesoscale and low frequencies.

**1. Introduction**

The atmospheric and oceanic general circulations are barotropically and baroclinically (Eady 1949; Charney 1947) unstable and spontaneously generate geostrophic turbulence. In the ocean, this mesoscale turbulence emerges at the scale of *O*(10–100) km and *O*(10–100) days and in turn strongly interacts with the general

circulation (e.g., Holland 1978). Mesoscale turbulence is also a well-known, strong manifestation of intrinsic ocean variability (i.e., it emerges without any atmospheric variability) at relatively small time and space scales. Mesoscale turbulence is now (partly) resolved in many ocean general circulation model (OGCM) global simulations, as well as in several recent coupled climate-oriented simulations. Bryan (2013) (and references therein) report various examples of how increases in OGCM resolution have substantially improved the fidelity of ocean simulations, and our understanding of dynamical processes underlying the multiscale oceanic variability. Beyond the inclusion of mesoscale phenomena into the resolved spectrum, the transition from laminar (1°–2°) to eddying 1/4° and finer) OGCM

Denotes Open Access content.

Corresponding author address: Guillaume Sérazin, Laboratoire de Glaciologie et Géophysique de l'Environnement LGGE/CNRS, BP 96, 38402 Saint-Martin d'Hères, Grenoble, France.
E-mail: guillaume.serazin@legi.grenoble-inp.fr

DOI: 10.1175/JCLI-D-14-00554.1

© 2015 American Meteorological Society
resolutions also yields substantial modifications in the simulated variability at larger space and time scales. Comparing a series of global ocean–sea ice “hindcasts” at increasing resolution with observations, Penduff et al. (2010) showed how the transition from laminar to eddy-permitting regimes brings the interannual variability of the ocean sea level much closer to altimeter measurements, in terms of both spatial distribution and of magnitude (increase of low-frequency variance larger than 60% over half of the global ocean).

Increasing resolution yields an increase in Reynolds number, which was also shown to favor the spontaneous emergence of low-frequency intrinsic variability (LFIV) in the ocean. Quasigeostrophic and shallow-water model experiments implemented on flat-bottom rectangular geometries have demonstrated that an interannual to decadal LFIV may appear under constant (or seasonal) forcing devoid of low frequencies and affect various components of the general circulation, in particular within the main unstable currents: strength and trajectory of surface and subsurface currents, intensity of recirculation gyres, thickness and volume of mode water pools, etc. (Spall 1996; Hazeleger and Drijfhout 2000; Simonnet and Dijkstra 2002; Dewar 2003; Pierini 2006; Berloff et al. 2007; Quattrocchi et al. 2012).

OGCM simulations tend to confirm these idealized predictions and illustrate various imprints of the LFIV in realistic contexts: path of oceanic jets (Taguchi et al. 2007; Thompson and Richards 2011; Douglass et al. 2012), subtropical mode waters (Douglass et al. 2013), eastern boundary circulation patterns (Combes and Di Lorenzo 2007), and larger-scale circulation features such as the Atlantic overturning circulation (Thomas and Zhai 2013; Grêgorio et al. 2014, manuscript submitted to J. Phys. Oceanogr., hereafter GPSH) or wide areas of the Southern Ocean (O’Kane et al. 2013). Penduff et al. (2011, hereafter P11) provided evidence that the strong increase in sea level anomaly (SLA) interannual variance, reported by Penduff et al. (2010) when switching from 2° to 1/4° resolution, may be related to the strong LFIV that emerges only in the eddying regime. This intrinsic component is indeed spontaneously produced by the same 1/4° model without any interannual forcing.

Our experimental strategy remains similar to P11's: a global OGCM is first driven by an atmospheric forcing containing a broad range of time scales to produce a realistic hindcast and then by a yearly repeated climatological atmospheric seasonal cycle to isolate the intrinsic variability. The present study further investigates the contribution of intrinsic processes to the SLA variability in the global ocean and extends P11's work in three ways. First, we highlight the spectral distribution of atmospherically forced and intrinsic SLA

variabilities within high- and low-frequency bands, at small and large spatial scales. Second, we are now using global simulations at 1/12° resolution in complement to P11's 1/4° simulations, providing us with the opportunity to assess the impact of model resolution on intrinsic variability. Finally, the present analysis is performed over a 42-yr period instead of 12 yr.

Section 2 presents the model configurations and numerical experiments used in this study. The simulated low-frequency SLA variability is then compared to altimeter observations. Postprocessing and filtering methods are described in section 3. In section 4, we compare the spatial distributions of low-frequency SLA variabilities between both 1/12° experiments at small (<6°) and large (>12°) spatial scales. Section 5 presents the sensitivity of intrinsic variability to the increase in resolution from 1/4° to 1/12°. Conclusions are given in section 6.

## 2. Numerical simulations and assessment

### a. Numerical simulations

As done earlier by Taguchi et al. (2007), P11, and O’Kane et al. (2013), our approach consists in comparing two eddying global OGCM simulations with two different atmospheric forcings. Our first 1/12° simulation (the so-called *T* experiment) is driven by the full range of atmospheric time scales between 1958 and 2012. It is intended to simulate the total ocean variability, which combines the intrinsic and atmospherically forced components. Our second 1/12° simulation (the *I* experiment), driven during 85 yr by a repeated climatological atmospheric seasonal cycle, isolates the low-frequency intrinsic variability that emerges without any low-frequency forcing. Note that the 1/4° resolution climatological simulation described and studied in P11 has been reprocessed and will be compared to the 1/12° resolution simulation in section 5. This 1/4° *I* experiment was conducted over 327 yr.

All the simulations were performed in the framework of the Drakkar project<fnref n="1" /> using the Nucleus for European Modeling of the Ocean (NEMO; Madec 2008) ocean/sea ice numerical model. The 1/12° simulations use the ORCA12 configuration with NEMO version 3.4, the 1/4° simulation uses the ORCA025 configuration with NEMO2.3. The three simulations<fnref n="2" /> share the same 46-level vertical discretization, a partial cell representation

<footnote n="1">¹ http://www.drakkar-ocean.eu</footnote>

<footnote n="2">² The 1/12° *I* and *T* experiments are referred to as ORCA12-GJM02 and ORCA12-MJM88 in the Drakkar database, respectively, and the 1/4° *I* experiment is referred to as ORCA025-MJM01.</footnote>
of topography, and a momentum advection scheme that conserves energy and entrophy (Barnier et al. 2006; Penduff et al. 2007; Le Sommer et al. 2009); a total variance diminishing (TVD) tracer advection scheme; an isopycnal Laplacian tracer diffusion operator; a vertical mixing scheme based on the TKE turbulent closure model (Blanke and Delecluse 1993); and a convective adjustment scheme based on enhanced vertical mixing in case of static instability.

The full and climatological forcing functions used to drive the 1/12° simulations come from the same original forcing dataset, which is referred to as the Drakkar forcing set (DFS4.4; Dussin and Barnier 2013); this 1958–2012 dataset is based on satellite observations (monthly precipitations and daily radiative heat fluxes) and on ERA-40 before 31 December 2001 and ERA-Interim afterward (6-hourly 10-m air temperature, humidity, and winds).<fnref n="3" /> The method we used to derive the climatological forcing from the full forcing is described in P11, who demonstrated (see also GPSH) that both forcing functions yield very similar mean states. Our 1/12° simulations are described in detail in Molines et al. (2014). Note that the seasonally forced 1/12° simulation has been analyzed in Treguier et al. (2014) and Deshayes et al. (2013) and that our pair of 1/12° simulations is also used by GPSH. The 327-yr 1/4° climatological simulation to which we compare the 1/12° I experiment in section 5 is presented in detail in P11. It was reprocessed exactly as its 1/12° counterpart (see section 3). More information about model configurations and solutions may be found in the aforementioned papers.

## b. Model assessment

Archiving, Validation, and Interpretation of Satellite Oceanographic Data (AVISO) provides us with 20 yr of global altimeter observations on a 1/4° × 1/4° regular Cartesian grid.<fnref n="4" /> Our fully forced, 1/12° global hindcast is first compared to this reference over the period July 1994–July 2011 in terms of low-frequency SLA standard deviation. We use monthly-mean, nonlinearly detrended, low-pass filtered (periods longer than 18 months) SLA fields to compute low-frequency SLA standard deviations

from both AVISO and the 1/12° T experiment; the post-processing method is detailed in the next section. Figure 1 demonstrates the capability of the 1/12° model to simulate the SLA low-frequency variability field with a remarkable accuracy. The location and intensity of the low-frequency variability in western boundary currents and their extension are well reproduced. Equatorial low-frequency activity is found in the same locations with comparable intensities in the observed and simulated oceans. The simulated path of the Antarctic Circumpolar Current variability maximum is close to that derived from altimetric observations. The 1/12° hindcast therefore yields a very realistic low-frequency variability field; it may be fruitfully analyzed, jointly with its climatologically forced counterpart, to characterize the intrinsic part of the low-frequency variability measured by altimetry.

## 3. Postprocessing and scale separation

### a. Separating temporal scales

The first years of model spinup are discarded from our analyses. Comparisons between the 1/12° simulations are performed from successive monthly sea surface height (SSH) averages over the same duration (47 yr): that is, over the period 1966–2012 from the T experiment and years 16–62 from the I experiment. The postprocessing of simulated SSH fields starts with the removal of global spatial averages from each monthly map (see Greatbatch 1994; P11). The resulting time series $\eta(t)$ may then be split as

$$ \eta(t) = \bar{\eta} + \tilde{\eta}(t) + \eta^s(t) + \eta^{LF}(t) + \eta^{HF}(t) \quad (1) $$

at every grid point, where $\bar{\eta}$ denotes the time average, $\tilde{\eta}$ is the long-term trend, $\eta^s$ is the mean annual cycle, $\eta^{LF}$ is the low-frequency (LF; interannual and slower) variability, and $\eta^{HF}$ is the high-frequency (HF; subannual) variability. Our study is focused on the intrinsic high- and low-frequency variabilities; time averages are removed first.

* Long-term trends $\tilde{\eta}(t)$ denote very low-frequency signals that need to be removed from raw time series to avoid biases in variance analyses. These trends may be of two kinds. The first kind is related to model adjustment toward its final equilibrium state, which may take a few centuries. This signal is unavoidable in OGCM simulations because of uncertainties that include those in atmospheric forcing, initial states, and physical and numerical approximations. However, it is small compared to the actual (well resolved) temporal variability. The second kind of trend may appear in finite-length time series when signals with periods longer than the time series itself are present. Both kinds of trends may yield a linear or nonlinear $\tilde{\eta}(t)$

<footnote n="3">³ DFS4.4 differs from DFS4 (Brodeau et al. 2010) only after 31 December 2001, where ERA-Interim is used instead of the ECMWF analysis for all forcing variables: wind vector, air temperature and humidity, downward shortwave and longwave radiation, total precipitation, and snowfall. Corrections described in Brodeau et al. (2010) ensure a smooth transition in 2001–02. ERA-Interim variables, whose native resolution is 0.7° and 3 h, were projected at the ERA-40 resolution (1.125° and 6 h) to build DFS4.4. Note that DFS4 was also compared with the more widely used CORE.v2 forcing in the latter reference.</footnote>

<footnote n="4">⁴ http://www.aviso.altimetry.fr/duacs/</footnote>
<figure><img src="image_57.png" /><figcaption>FIG. 1. Low-frequency ($T > 18$ months) SLA standard deviations computed over the period July 1994–July 2011 from (top) the 1/4° AVISO dataset and (bottom) the 1/12° T experiment.</figcaption></figure>

evolution over the length of the time series. Our detrending approach is based on a nonlinear, second-order local regression method (LOESS; Cleveland and Devlin 1988), which high-pass filters model time series with a 20-yr cutoff period. This method preserves the length of original time series without adverse edge effects (Fig. 2). Our cutoff period ensures the removal of both kinds of trends, and confines our analyses and results on time scales between 2 months and 20 yr in both simulations.

• The mean annual cycle $\eta^s(t)$ is then computed from SLA time series and subtracted to obtain deseasonalized, detrended SLA time series.

• The low- and high-frequency [$\eta^{LF}(t)$ and $\eta^{HF}(t)$, respectively] variabilities are then isolated from the time series resulting from the former steps. Following P11, our low-frequency signals include time scales longer than 18 months. Unlike these authors, however, our high-frequency signals include time scales ranging between 2 and 18 months, with no mean annual cycle. Low- and high-frequency components are separated using a Lanczos temporal filter (Duchon 1979), chosen for

its efficiency and its ability to provide clean signals with small Gibbs oscillations. Note that this is a linear filtering method and it has edge effects because of the fixed size of the convolution kernel (relative to the filter order). Thus, 3 yr on both extremities of the time series are lost in this computation, yielding low- and high-frequency 41-yr time series from both 1/12° runs: years 1969–2009 in the *T* experiment and years 19–59 in the *I* experiment.

b. Separating spatial scales

The main purpose of this paper is to document how intrinsic processes contribute to the total SLA variability in various wavenumber–frequency ($\kappa-\omega$) classes. High- and low-frequency datasets are further split into three ranges of spatial scales by applying a two-dimensional, isotropic Lanczos low-pass filter twice with different cutoff lengths. We shall focus in the following on structures shorter than about 6° [labeled as small scales (SS)] and structures larger than about 12° [large scales (LS)]. Structures with sizes between 6° and
<figure><img src="image_58.png" /><figcaption>FIG. 2. Example of the nonlinear trend calculated (thick line) from a 42-yr time series of SLA at (52.5°N, 179°W) in the ACC (thin line) obtained in the *T* experiment.</figcaption></figure>

12° are also produced during the filtering process but are not investigated here. Motivations leading to 6° and 12° are similar to those given in Penduff et al. (2010) and Taguchi et al. (2007).

In practice, the spatial cutoffs are defined in terms of model grid points rather than distances: hence, the approximate scales mentioned above. The spatial filter is adapted to process SLA fields in the vicinity of land points (islands and coasts): the kernel coefficients used in the convolution are weighted to take into account only wet points. To limit computational costs, the size of the convolution kernel was kept the same in both spatial filtering processes: this size is close to four 6° wavelengths and two 12° wavelengths. The 6° filter is therefore more selective than the 12° filter.

In summary, the postprocessing chain described in section 3 decomposes the (detrended) total and intrinsic variabilities within six main $κ-ω$ bands (as summarized in Fig. 3). Figure 4 illustrates three steps of this processing in November 1997: detrended monthly SLA fields after temporal low-pass filtering (top panels) and after further splitting into small and large scales (middle and bottom panels). Note that the high-frequency small-scale (HFSS) variability will also be referred to as mesoscale because geostrophic turbulence is found in this spectral domain.

## 4. Spatial scales of intrinsic SLA variability

We focus first on the LF component of SLA variabilities in both model experiments, taking all spatial scales into account. Figure 5 (top) compares in zonal

<figure><img src="image_59.png" /><figcaption>FIG. 3. Filter cutoffs in the wavenumber-frequency ($κ-ω$) space. The hatched horizontal bar at 12 months indicates that the mean annual cycle harmonics were removed from high-frequency fields.</figcaption></figure>

average their meridional distributions. The total (fully forced) SLA low-frequency standard deviation remains within 3–4 cm over most of the 60°S–60°N latitudinal range, with absolute maxima at midlatitudes and minima poleward of 60° in both hemispheres. Its intrinsic counterpart exhibits a marked contrast between low-latitude and midlatitude areas. Both the intrinsic SLA low-frequency variability and its relative contribution $σ_I^2/σ_T^2$ reach their absolute minima in the intertropical band (less than 1 cm and 10%, respectively, between 10°S and 10°N) and their absolute maxima at midlatitudes in both hemispheres. Poleward of 60°N the intrinsic low-frequency variance reaches its secondary minima; most of the SLA variance there is therefore related to the low-frequency atmospheric variability in the *T* experiment.

### a. Small-scale intrinsic variability

The resemblance between the blue lines in the top and middle panels of Fig. 5 illustrates the importance of small-scale features in the distribution of intrinsic low-frequency variability; this section is focused on this low-frequency small-scale (LFSS) intrinsic variability. Low-frequency small-scale sea level anomalies have very similar distributions and amplitudes with or without atmospheric variability (middle panels in Figs. 4 and 5). The LFSS SLA variability is therefore largely intrinsic: that is, spontaneously produced by the eddying ocean over most regions (bottom panel in Fig. 6).<fnref n="5" /> This result is consistent with Taguchi et al. (2007), who showed that the structure of the

<footnote n="5">⁵The few exceptions to this concerns regions where LFSS variance is very weak ($σ < 0.5$ cm): that is, in shallow areas and equatorial basins.</footnote>
<figure><img src="image_60.png" /><figcaption>FIG. 4. SLA global maps on November 1997 from the (left) *T* and (right) *I* experiments in the low-frequency range: (top) containing all spatial scales, (middle) restricted to small spatial scales (*L* < 6°), and (bottom) restricted to large-scale signals only (*L* > 12°).</figcaption></figure>

interannually varying Kuroshio has small (meridional) scales and is shaped by purely oceanic processes (i.e., identical with and without low-frequency forcing). A similar conclusion was drawn in the Southern Ocean by O’Kane et al. (2013), who showed that small-scale low-frequency modulations of the Antarctic Circumpolar Current (ACC) are captured in both fully forced and climatological simulations. Intrinsic, interannual fluctuations of sea level can largely be explained by mesoscale processes in the Gulf of Alaska as well (Okkonen et al. 2001; Combes and Di Lorenzo 2007). In other words, the bottom panel in Fig. 6 suggests that the conclusions of these regional studies may well be valid over most of the global ocean.

As recalled in the introduction, mesoscale eddies emerge through hydrodynamic instability processes at small space and time scales. These motions are captured in our HFSS spectral box, and their (intrinsic) standard deviation in the *I* experiment is shown in the top panel of Fig. 7. This map is very close to its counterpart in the *T* experiment (not shown), confirming that most HFSS (mesoscale) activity is generated through intrinsic

processes and is barely sensitive to low-frequency atmospheric variability.

Interestingly, large values of LFSS intrinsic variability are found where mesoscale eddies are strong (Fig. 6): that is, in the ACC, East Australian Current, Brazil Current/Malvinas Confluence region, Agulhas Current, Kuroshio, and Gulf Stream/North Atlantic Current systems. Indeed, LFSS and HFSS intrinsic variability maps (Figs. 6 and 7) are very similar: the bottom panels in Fig. 7 show that SLA standard deviations in both spectral ranges have similar meridional distributions (left) and are highly correlated in space. These results are consistent with the spontaneous emergence of mesoscale eddies through (intrinsic) instability processes in both simulations, with only a weak sensitivity to atmospheric variability. The marked coincidence between HFSS and LFSS intrinsic variability maps suggests that the energy of these mesoscale motions may spontaneously cascade toward longer time scales and feed the LFSS spectral component. Such a temporal inverse cascade has indeed been diagnosed recently from idealized and realistic simulations (Arbic et al. 2012, 2014), albeit
<figure><img src="image_61.png" /><figcaption>FIG. 5. Zonal average of total (red) and intrinsic (blue) low-frequency standard deviations: (top) all spatial scales, (middle) small scales, and (bottom) large scales. Corresponding intrinsic to total variance ratios are shown (right y axis in %) in green in each panel.</figcaption></figure>

over a narrower range of time scales. Whether this non-linear process may actually transfer mesoscale kinetic energy up to interannual and longer time scales remains to be assessed and is currently being investigated.

b. *Large-scale low-frequency intrinsic variability*

The bottom panels in Fig. 4 indicate that strong low-frequency large-scale SLAs, which are hindcasted in the *T* experiment (strong El Niño event, Indian Ocean dipole,

and Pacific decadal oscillation positive phases in November 1997), are not spontaneously generated by the eddy-ing ocean; they require direct forcing by the atmosphere (or full air-sea coupling) to exist. The zonally averaged LFLS SLA variability is indeed much stronger with the full forcing over most of the global ocean (Fig. 5), especially in the intertropical band, where SLA intensities reach their global maximum in the *T* experiment and their global minimum in the *I* experiment.
<figure><img src="image_62.png" /><figcaption>FIG. 6. Low-frequency small-scale (top) total and (middle) intrinsic standard deviations. (bottom) Intrinsic to total variance ratio in the same range of scales.</figcaption></figure>

In western boundary current extensions and in the ACC, however, our 1/12° model spontaneously generates low-frequency large-scale SLA anomalies with small but nonnegligible zonally averaged intensities (bottom-right panel in Fig. 4 and bottom panel in Fig. 5). LFLS

intrinsic variability is strongest in eddy-active regions (middle panel in Fig. 8): that is, where small-scale standard deviations were proven largest (HFSS and LFSS bands). The geographical correspondence between these three spectral maxima further suggests that
<figure><img src="image_63.png" /><figcaption>FIG. 7. (top) High-frequency small-scale intrinsic variability. (bottom left) Normalized distributions of intrinsic variances (zonal averages) as a function of latitude at LFSS (red) and HFSS (green). (bottom right) The spatial correlation of these latter two fields as a function of latitude.</figcaption></figure>

intrinsic, eddy-driven inverse cascade processes may be transferring mesoscale (HFSS) kinetic energy toward larger temporal (LFSS) and spatial (LFLS) scales (e.g., Berloff et al. 2007; Arbic et al. 2014).

The contribution of intrinsic processes to the fully forced, low-frequency large-scale SLA variance is therefore substantial in several regions (bottom panel in Fig. 8): it exceeds 30% over wide regions of the Northern Hemisphere and reaches 40% in zonal average around 50°S, with local maxima near 100% in Drake Passage, the Brazil Current/Malvinas Confluence, and the Agulhas Current.

**5. Influence of model resolution**

P11 presented an analysis of the low-frequency intrinsic variability from a seasonally forced simulation at 1/4°. This eddy-permitting simulation has been reprocessed in the same manner as was our 1/12° simulations in order to assess the influence of model resolution on intrinsic variability. The comparison is based on SLA standard deviations computed over 36 yr from both simulations.

**a. Increase in mesoscale variability**

Increasing resolution from 1/4° to 1/12° is known to strongly enhance the strength of simulated mesoscale

activity in NEMO as well in other models (e.g., Hurlburt et al. 2009). Indeed, the zonally averaged HFSS SLA variability, which is largely intrinsic as mentioned earlier, markedly increases as resolution is refined at midlatitudes and in the ACC belt (top panel in Fig. 9). The top panel in Fig. 10 shows that this increase in intrinsic, mesoscale SLA standard deviation reaches up to 2 cm in most regions poleward of 10°–15°. Increased resolution, however, yields less mesoscale activity in the vicinity of the Kuroshio and Gulf Stream extensions and the ACC. Most of these local decreases are actually due to the displacement of energetic fronts and associated recirculation gyres and eddy fields, toward a more realistic mean state in the 1/12° model.

**b. Increase in low-frequency intrinsic variability**

Increased resolution does not only enhance mesoscale turbulence. The bottom panel in Fig. 9 shows a concomitant increase in low-frequency intrinsic SLA standard deviation at all spatial scales. This increase occurs in regions where mesoscale activity is enhanced by the resolution increase, not only in zonal average but also at regional and frontal scales (Figs. 9 and 10). This tight geographical correspondence between changes in
<figure><img src="image_64.png" /><figcaption>FIG. 8. As in Fig. 6, but at LFLS scales.</figcaption></figure>

mesoscale and LF intrinsic variability was also striking at 1/12° (Fig. 7) and further supports the idea that LF intrinsic variability may be directly fed by mesoscale activity through nonlinear temporal inverse cascade processes, both within and farther away from eddy-active regions.

Interestingly, GPSH are showing from the same simulations that in contrast with SLA, most of the intrinsic LF variability diagnosed at 1/12° is captured at 1/4° with respect to the Atlantic meridional overturning stream-function (AMOC). This shows that a 1/4° ocean model that somewhat underestimates local intrinsic variability
<figure><img src="image_65.png" /><figcaption>FIG. 9. Zonally averaged intrinsic standard deviation of SLA (cm) at 1/4° (solid) and 1/12° (dashed) at (top) HFSS and (bottom) LF for all scales.</figcaption></figure>

imprints on SLA may yet capture the main processes generating intrinsic variability of the important climate-related, large-scale oceanic index of the AMOC. Possible reasons for this remain to be identified.

## 6. Conclusions

The intrinsic variability that eddying ocean models spontaneously generate has been analyzed through its imprint on sea level anomalies at various scales. The NEMO-based 1/12° ORCA12 global OGCM configuration was driven by the full range of atmospheric time scales ("hindcast" simulation) and then by a repeated atmospheric annual cycle (climatological simulation). Resulting SLA fields have been deseasonalized, filtered, and compared in three spatiotemporal ranges: the high-frequency small-scale (HFSS) range ($T < 18$ months, $L < 6°$), which includes mesoscale turbulence;

the low-frequency small-scale (LFSS) range ($T > 18$ months, $L < 6°$); and the low-frequency large-scale (LFLS) range ($T > 18$ months, $L > 12$). The hindcast was shown to reproduce the distribution and magnitude of low-frequency ($T > 18$ months) variability observed by satellite altimetry with a high degree of accuracy.

* Just as for mesoscale activity, most of the LFSS SLA variability is spontaneously generated by the eddying ocean with only weak sensitivity to the atmospheric variability. Regions of strong mesoscale variability closely coincide with regions where LFSS variability is large. This coincidence suggests that mesoscale energy may spontaneously cascade toward longer time scales through nonlinear processes, as recently suggested by Arbic et al. (2014).

* Most of the LFLS SLA variability is forced by the atmosphere in our simulations, in particular at low latitudes and in eastern basins. In eddy-active regions, however, intrinsic processes generate an important part of this large-scale low-frequency SLA variance: from 30%–50% in western boundary current systems up to 90–100% where the South Atlantic connect with surrounding basins.

* Decreasing the model resolution from 1/12° to 1/4° yields a decrease of both mesoscale and low-frequency intrinsic SLA variance, along with displacements of certain main currents. The 1/4° low-frequency intrinsic variability distribution, however, remains comparable with its 1/12° counterpart, suggesting that generation processes are reasonably well captured in the eddy-permitting regime.

Many features of the intrinsic low-frequency variability remain unknown. Our results could be complemented by two- or three-dimensional, possibly multivariate, characterizations of intrinsic modes of variability in various regions [e.g., as done in the Southern Ocean by O'Kane et al. (2013)]. Of particular interest would be a precise description of how intrinsic variability imprints sea surface temperature or upper-ocean heat content at low frequencies, given the sensitivity of the atmospheric variability to these thermal constraints at long time scales (Gulev et al. 2013; Brachet et al. 2012). Current investigations are indeed highlighting substantial SST imprints where intrinsic variability imprints SLA.

The mechanisms that generate intrinsic low-frequency variability in complex OGCM simulations also remain under debate. Based on 5-yr SLA time series, simulated by a global eddying OGCM driven by a full forcing, Arbic et al. (2014) showed that nonlinear eddy-eddy interactions can spontaneously transfer kinetic energy from mesoscale eddies, which are largely chaotic, to longer time scales. In other words, this
<figure><img src="image_66.png" /><figcaption>FIG. 10. Difference ($\sigma_I^{1/12} - \sigma_I^{1/4}$; cm) of intrinsic standard deviation of SLA between 1/12° and 1/4° simulations, at (top) HFSS and (bottom) LF. The 1/12° field has been subsampled on the 1/4° grid.</figcaption></figure>

temporal inverse cascade process might well transfer the mesoscale intrinsic chaotic variability toward lower frequencies. This analysis is currently being extended to much longer (multidecadal) SLA time series provided by our seasonally forced 1/4° and 1/12° I experiments, in order to assess whether the same mechanism may actually generate intrinsic variability up to decadal (or possibly longer) time scales. Note that very little is known today about the oceanic intrinsic variability at very long time scales (i.e., at multidecadal or longer periods). GPSH have recently shown from the full 327-yr 1/4° I experiment that the AMOC spontaneously fluctuates up to such time scales; assessing the possible imprint of these processes on SLA, other variables, and climate indices remains to be done.

We have focused on the scales of the low-frequency intrinsic variability that spontaneously emerges without low-frequency forcing. The behavior of intrinsic variability modes in fully forced simulations (with low-frequency forcing)—and more generally the actual

constraint exerted by the atmospheric forcing on regions where this intrinsic variability is strong—is another important (and complex) issue. The Kuroshio system is already known to exhibit intrinsic modes of variability, which might be triggered by the low-frequency atmospheric forcing when present (Taguchi et al. 2007). This phenomenon has been studied in the framework of dynamical system theory with process-oriented idealized models; results show that stochastic forcings (Pierini 2011) or slowly varying forcings (Pierini 2014) could excite and lead to the resonance of intrinsic modes that are not spontaneously generated under a constant wind forcing. Whether this process is at work and dominates in complex OGCM simulations is still unknown and would deserve further investigations that we leave for the future; multiple high-resolution OGCM simulations, such as those performed in idealized contexts, or ensemble simulations would help address this question.

Finally, our ocean/sea ice simulations were forced by prescribed atmospheres. This decoupled approach,
inspired from several idealized studies (e.g., reviewed in Dijkstra and Ghil 2005), allows the isolation and investigation of the intrinsic variability itself, ignoring its interaction with the atmosphere. Our results may help interpret how the transition from laminar to eddying oceans affects the low-frequency variability of the fully coupled climate system. We argue that a better understanding of intrinsic SLA variability from such eddying simulations will also help interpreting the present 20-yr altimetric observed dataset, as the models tend to get more accurate and closer to these data.

*Acknowledgments.* The authors acknowledge the constructive comments made by three anonymous reviewers, which led to a significant improvement of this paper. This work is a contribution to CHAOCEAN, OCCIPUT, and MyOcean2 projects. It benefited from the Drakkar international coordination network (GDRI) established between the Centre National de la Recherche Scientifique (CNRS), the National Oceanography Centre in Southampton (NOCS), GEOMAR in Kiel, and IFREMER. For this work, Drakkar also benefited from a grant from the Groupe Mission Mercator Coriolis (GMMC) through the LEFE program of the Institut National des Sciences de l'Univers (INSU). CHAOCEAN is supported by the Centre National d'Études Spatiales (CNES) through the Ocean Surface Topography Science Team (OST/ST). OCCIPUT is supported by the Agence Nationale de la Recherche (ANR) through Contract ANR-13-BS06-0007-01. The research leading to these results has received funding from the European Community's 160 Seventh Framework Programme FP7/2007-2013 under Grant Agreement 283367 (MyOcean2). The computations presented in this study were performed at the Centre Informatique National de l'Enseignement Supérieur (CINES) under the allocation made by GENCI x2013010727. The ORCA12-GJM02 simulation was performed as part of the Grands Challenges GENCI/CINES 2013. The altimeter products were produced by SSALTO/DUACS and distributed by AVISO, with support from CNES. GS is supported by CNES and Région Midi-Pyrénées, TP is supported by CNRS, SG by MyOcean2, and LT is supported by CERFACS.

## REFERENCES

Arbic, B. K., R. B. Scott, G. R. Flierl, A. J. Morten, J. G. Richman, and J. F. Shriver, 2012: Nonlinear cascades of surface oceanic geostrophic kinetic energy in the frequency domain. *J. Phys. Oceanogr.*, **42**, 1577–1600, doi:10.1175/JPO-D-11-0151.1.

—, M. Müller, J. G. Richman, J. F. Shriver, A. J. Morten, R. B. Scott, G. Sérazin, and T. Penduff, 2014: Geostrophic turbulence

in the frequency–wavenumber domain: Eddy-driven low-frequency variability. *J. Phys. Oceanogr.*, **44**, 2050–2069, doi:10.1175/JPO-D-13-054.1.

Barnier, B., and Coauthors, 2006: Impact of partial steps and momentum advection schemes in a global ocean circulation model at eddy-permitting resolution. *Ocean Dyn.*, **56** (5–6), 543–567, doi:10.1007/s10236-006-0082-1.

Berloff, P. S., A. M. Hogg, and W. Dewar, 2007: The turbulent oscillator: A mechanism of low-frequency variability of the wind-driven ocean gyres. *J. Phys. Oceanogr.*, **37**, 2363–2386, doi:10.1175/JPO3118.1.

Blanke, B., and P. Delecluse, 1993: Variability of the tropical Atlantic Ocean simulated by a general circulation model with two different mixed-layer physics. *J. Phys. Oceanogr.*, **23**, 1363–1388., doi:10.1175/1520-0485(1993)023<1363:VOTTAO>2.0.CO;2.

Brachet, S., F. Codron, Y. Feliks, M. Ghil, H. Le Treut, and E. Simonnet, 2012: Atmospheric circulations induced by a midlatitude SST front: A GCM study. *J. Climate*, **25**, 1847–1853, doi:10.1175/JCLI-D-11-00329.1.

Brodeau, L., B. Barnier, A.-M. Treguier, T. Penduff, and S. Gulev, 2010: An ERA40-based atmospheric forcing for global ocean circulation models. *Ocean Modell.*, **31** (3–4), 88–104, doi:10.1016/j.ocemod.2009.10.005.

Bryan, F. O., 2013: Introduction: Ocean modeling—Eddy or not. *Ocean Modeling in an Eddying Regime*, Geophys. Monogr., Vol. 177, Amer. Geophys. Union, 1–3.

Charney, J. G., 1947: The dynamics of long waves in a baroclinic westerly current. *J. Meteor.*, **4**, 136–162, doi:10.1175/1520-0469(1947)004<0136:TDOLWI>2.0.CO;2.

Cleveland, W. S., and S. J. Devlin, 1988: Locally weighted regression: An approach to regression analysis by local fitting. *J. Amer. Stat. Assoc.*, **83**, 596–610, doi:10.1080/01621459.1988.10478639.

Combes, V., and E. Di Lorenzo, 2007: Intrinsic and forced interannual variability of the Gulf of Alaska mesoscale circulation. *Prog. Oceanogr.*, **75**, 266–286, doi:10.1016/j.pocean.2007.08.011.

Deshayes, J., and Coauthors, 2013: Oceanic hindcast simulations at high resolution suggest that the Atlantic MOC is bistable. *Geophys. Res. Lett.*, **40**, 3069–3073, doi:10.1002/grl.50534.

Dewar, W. K., 2003: Nonlinear midlatitude ocean adjustment. *J. Phys. Oceanogr.*, **33**, 1057–1082, doi:10.1175/1520-0485(2003)033<1057: NMOA>2.0.CO;2.

Dijkstra, H. A., and M. Ghil, 2005: Low-frequency variability of the large-scale ocean circulation: A dynamical systems approach. *Rev. Geophys.*, **43**, RG3002, doi:10.1029/2002RG000122.

Douglass, E. M., S. R. Jayne, F. O. Bryan, S. Peacock, and M. Maltrud, 2012: Kuroshio pathways in a climatologically forced model. *J. Oceanogr.*, **68**, 625–639, doi:10.1007/s10872-012-0123-y.

—, Y.-O. Kwon, and S. R. Jayne, 2013: A comparison of North Pacific and North Atlantic subtropical mode waters in a climatologically-forced model. *Deep-Sea Res. II*, **91**, 139–151, doi:10.1016/j.dsr2.2013.02.023.

Duchon, C., 1979: Lanczos filtering in one and two dimensions. *J. Appl. Meteor.*, **18**, 1016–1022, doi:10.1175/1520-0450(1979)018<1016:LFIOAT>2.0.CO;2.

Dussin, R., and B. Barnier, 2013: The making of DFS 5.1. Drakkar Project Rep., 40 pp. [Available online at http://www.drakkar-ocean.eu/publications/reports/dfs5-1-report.]

Eady, E. T., 1949: Long waves and cyclone waves. *Tellus*, **1**, 33–52, doi:10.1111/j.2153-3490.1949.tb01265.x.
Greatbatch, R. J., 1994: A note on the representation of steric sea level in models that conserve volume rather than mass. *J. Geophys. Res.*, **99** (C6), 12 767–12 771, doi:10.1029/94JC00847.

Gulev, S. K., M. Latif, N. Keenlyside, W. Park, and K. P. Koltermann, 2013: North Atlantic Ocean control on surface heat flux on multidecadal timescales. *Nature*, **499**, 464–467, doi:10.1038/nature12268.

Hazeleger, W., and S. S. Drijfhout, 2000: A model study on internally generated variability in subtropical mode water formation. *J. Geophys. Res.*, **105** (C6), 13 965, doi:10.1029/2000JC900041.

Holland, W. R., 1978: The role of mesoscale eddies in the general circulation of the ocean–numerical experiments using a wind-driven quasi-geostrophic model. *J. Phys. Oceanogr.*, **8**, 363–392, doi:10.1175/1520-0485(1978)008<0363: TROMEI>2.0.CO;2.

Hurlburt, H. E., and Coauthors, 2009: High-resolution global and basin-scale ocean analyses and forecasts. *Oceanography*, **22**, 110–127, doi:10.5670/oceanog.2009.70.

Le Sommer, J., T. Penduff, S. Theetten, G. Madec, and B. Barnier, 2009: How momentum advection schemes influence current-topography interactions at eddy permitting resolution. *Ocean Modell.*, **29**, 1–14, doi:10.1016/j.ocemod.2008.11.007.

Madec, G., 2008: NEMO ocean engine. Institut Pierre-Simon Laplace (IPSL) Note du Pole de Modélisation 27, 217 pp.

Molines, J.-M., B. Barnier, T. Penduff, A. M. Treguier, and J. Le Sommer, 2014: ORCA12.L46 climatological and interannual simulations forced with DFS4.4: GJM02 and MJM88. Drakkar Group Experiment Rep. GDRI-DRAKKAR-2014-03-19, 50 pp. [Available online at http://www.drakkar-ocean.eu/publications/reports/orca12_reference_experiments_2014.]

O’Kane, T. J., R. J. Matear, M. A. Chamberlain, J. S. Risbey, B. M. Sloyan, and I. Horenko, 2013: Decadal variability in an OGCM Southern ocean: Intrinsic modes, forced modes and metastable states. *Ocean Modell.*, **69**, 1–21, doi:10.1016/j.ocemod.2013.04.009.

Okkonen, S. R., G. A. Jacobs, E. Joseph Metzger, H. E. Hurlburt, and J. F. Shriver, 2001: Mesoscale variability in the boundary currents of the Alaska gyre. *Cont. Shelf Res.*, **21** (11–12), 1219–1236, doi:10.1016/S0278-4343(00)00085-6.

Penduff, T., J. Le Sommer, B. Barnier, A.-M. Treguier, J.-M. Molines, and G. Madec, 2007: Influence of numerical schemes on current-topography interactions in 1/4° global ocean simulations. *Ocean Sci.*, **3**, 509–524, doi:10.5194/os-3-509-2007.

—, M. Juza, L. Brodeau, G. C. Smith, B. Barnier, J.-M. Molines, A.-M. Treguier, and G. Madec, 2010: Impact of global ocean model resolution on sea-level variability with emphasis on interannual time scales. *Ocean Sci.*, **6**, 269–284, doi:10.5194/os-6-269-2010.

—, —, B. Barnier, J. Zika, W. K. Dewar, A.-M. Treguier, J.-M. Molines, and N. Audiffren, 2011: Sea level expression of intrinsic and forced ocean variabilities at interannual time scales. *J. Climate*, **24**, 5652–5670, doi:10.1175/JCLI-D-11-00077.1.

Pierini, S., 2006: A Kuroshio Extension system model study: Decadal chaotic self-sustained oscillations. *J. Phys. Oceanogr.*, **36**, 1605–1625, doi:10.1175/JPO2931.1.

—, 2011: Low-frequency variability, coherence resonance, and phase selection in a low-order model of the wind-driven ocean circulation. *J. Phys. Oceanogr.*, **41**, 1585–1604, doi:10.1175/JPO-D-10-05018.1.

—, 2014: Kuroshio Extension bimodality and the North Pacific Oscillation: A case of intrinsic variability paced by external forcing. *J. Climate*, **27**, 448–454, doi:10.1175/JCLI-D-13-00306.1.

Quattrocci, G., S. Pierini, and H. A. Dijkstra, 2012: Intrinsic low-frequency variability of the Gulf Stream. *Nonlinear Processes Geophys.*, **19**, 155–164, doi:10.5194/npg-19-155-2012.

Simonnet, E., and H. A. Dijkstra, 2002: Spontaneous generation of low-frequency modes of variability in the wind-driven ocean circulation. *J. Phys. Oceanogr.*, **32**, 1747–1762, doi:10.1175/1520-0485(2002)032<1747:SGOLFM>2.0.CO;2.

Spall, M. A., 1996: Dynamics of the Gulf Stream/deep western boundary current crossover. Part II: Low-frequency internal oscillations. *J. Phys. Oceanogr.*, **26**, 2169–2182, doi:10.1175/1520-0485(1996)026<2169:DOTGSW>2.0.CO;2.

Taguchi, B., S.-P. Xie, N. Schneider, M. Nonaka, H. Sasaki, and Y. Sasai, 2007: Decadal variability of the Kuroshio Extension: Observations and an eddy-resolving model hindcast. *J. Climate*, **20**, 2357–2377, doi:10.1175/JCLI4142.1.

Thomas, M. D., and X. Zhai, 2013: Eddy-induced variability of the meridional overturning circulation in a model of the North Atlantic. *Geophys. Res. Lett.*, **40**, 2742–2747, doi:10.1002/grl.50532.

Thompson, A. F., and K. J. Richards, 2011: Low frequency variability of Southern Ocean jets. *J. Geophys. Res.*, **116**, C09022, doi:10.1029/2010JC006749.

Treguier, A. M., and Coauthors, 2014: Meridional transport of salt in the global ocean from an eddy-resolving model. *Ocean Sci.*, **10**, 243–255, doi:10.5194/os-10-243-2014.
<figure><img src="image_67.png" /></figure>

Intrinsic variability of the Atlantic meridional overturning circulation at
interannual-to-multidecadal timescales

Sandy Grégorio* and Thierry Penduff

CNRS, Univ. Grenoble Alpes, LGGE (UMR5183), F-38041 Grenoble, France.

Guillaume Sérazin

CNRS, LGGE (UMR5183), F-38041 Grenoble, France.

Sciences de l'Univers au CERFACS, CERFACS/CNRS, URA1857, Univ. Paul Sabatier, Toulouse,
France.

Jean-Marc Molines and Bernard Barnier

CNRS, Univ. Grenoble Alpes, LGGE (UMR5183), F-38041 Grenoble, France.

Joël Hirschi

National Oceanography Centre, Southampton, SO14 3ZH, UK

*Corresponding author address: Laboratoire de Glaciologie et Géophysique de l'Environnement,

LGGE/CNRS, BP 96, F-38402 Saint-Martin d'Hères cedex, France.

E-mail: sandy Gregorio@legi.grenoble-inp.fr
ABSTRACT
The low-frequency variability of the Atlantic Meridional Overturning Circulation (AMOC) is investigated from 2°, 1/4°, and 1/12° global ocean/sea-ice simulations, with a specific focus on its internally-generated (i.e. "intrinsic") component. A 327-year climatological 1/4° simulation, driven by a repeated seasonal cycle (i.e. a forcing devoid of interannual timescales), is shown to spontaneously generate a significant fraction *R* of the interannual-to-decadal AMOC variance obtained in a 50-year "fully-forced" hindcast (with reanalyzed atmospheric forcing including interannual timescales). This intrinsic variance fraction *R* slightly depends on whether AMOCs are computed in geopotential or density coordinates, and on the period considered in the climatological simulation, but the following features are quite robust when mesoscale eddies are simulated (at both 1/4° and 1/12° resolutions): *R* barely exceeds 5-10% in the subpolar gyre, but reaches 30-50% at 34°S, up to 20-40% near 25°N and 40-60% near the Gulf Stream. About 25% of the meridional heat transport interannual variability is attributed to intrinsic processes at 34°S and near the Gulf Stream. Fourier and wavelet spectra, built from the 327-year 1/4° climatological simulation, further indicate that spectral peaks of intrinsic AMOC variability [i] are found at specific frequencies ranging from interannual to multidecadal, [ii] often extend over the whole meridional scale of gyres, [iii] stochastically change throughout these 327 years, and [iv] sometimes match the spectral peaks found in the fully-forced hindcast in the North Atlantic. Intrinsic AMOC variability is also detected at multidecadal timescales, with a marked meridional coherence between 35°S and 25°N (15-30 year periods), and throughout the whole basin (50-90-year periods).
## 1. Introduction

The Atlantic meridional overturning circulation (AMOC) consists of two main cells that extend across the whole latitudinal range of the basin. In the upper cell, warm surface waters are advected northward: they progressively lose their heat to the atmosphere in the northern basin, become denser and sink at middepth at high latitudes to feed the southward return flow of the North Atlantic Deep Water (NADW). In the bottom cell, Antarctic Bottom Water (AABW) enters the basin from the south along the sea floor, raises upward as it mixes with lighter water, and finally re-circulates and exits the basin southward within the NADW mid-depth flow. The AMOC upper cell's northward limb transports a substantial amount of heat northward, which reaches its maximum (about 1.3 PW) at 24°N-26°N; this advective heat flux accounts for about 25% of the total (atmosphere and ocean) poleward heat transport at those latitudes (Hall and Bryden 1982; Johns et al. 2011). The convergence of this heat transport strongly influences the air-sea heat fluxes in the Gulf Stream region; changes in the AMOC northward limb transport may therefore have substantial climatic impacts (Dong et al. 2007).

Observational programmes, developed in the last decade, are monitoring the AMOC strength, structure, and variability: the RAPID-MOCHA program (Hirschi et al. 2003; Baehr et al. 2004; Cunningham et al. 2007; Kanzow et al. 2007; McCarthy et al. 2012) started in 2004 and continuously monitors the AMOC at the 26.5°N latitude; the MOVE program (Send et al. 2002) started in 2000 and is focused on 16°N; the OVIDE programme (Lherminier et al. 2007) has been monitoring a section from Greenland to Portugal every other year since 2002. Besides, the mechanisms of AMOC variability have also been studied numerically (e.g. Delworth et al. 1993; Bentsen et al. 2004) throughout the whole basin and on much longer time scales, either from ocean simula-
tions driven by atmospheric reanalyses (in laminar or eddying regimes) or from coupled ocean-atmospheric simulations (mostly in the laminar regime).

Most studies to date about the AMOC variability in the eddying regime have been concerned with the response of the AMOC to a prescribed atmospheric variability. Using Ocean General Circulation Models (OGCMs), Biastoch et al. (2008) and Hirschi et al. (2013) have documented another aspect of the interannual-to-decadal variability of the AMOC computed along geopotential levels. Both studies revealed that when mesoscale eddies are (even partially) resolved, non-linearities become significant and the phase of AMOC timeseries is no longer locked to the phase of the atmospheric forcing: it becomes partly chaotic, up to decadal timescales. Using a laminar global OGCM with an increase in resolution in the Agulhas region driven by an atmospheric reanalysis, Biastoch et al. (2008) showed that resolving mesoscale eddies locally strongly modifies the phase of the decadal variability of the AMOC northward limb transport (approximated by the AMOC timeseries at 1000 m); this phase perturbation was maximum in the Agulhas region but did propagate as far as the northern hemisphere subtropics within a few months. Hirschi et al. (2013) then compared AMOC northward limb transport timeseries obtained from two eddy-permitting (1/4°) global OGCM hindcasts with different initial conditions but identical surface forcing. The chaotic interannual-to-decadal AMOC variability component, revealed by substantial differences between individual timeseries, was shown to account for a significant proportion of the total AMOC interannual variance found in individual hincasts: around 15-35% throughout the South Atlantic, 10% between the equator and 25°N, 15-20% between 25 and 47°N, and about 5% further north <fnref n="1" />. Both Biastoch et al. (2008) and Hirschi et al. (2013) thus identified substantial chaotic AMOC variability, some of which may be generated in the South Atlantic and propagate to the North.

<footnote n="1">¹Hirschi et al. (2013)'s Fig. 14c presents this proportion in terms of AMOC standard deviations at low frequencies.</footnote>
The present study uses a different strategy to further investigate how the low-frequency (LF, i.e. at timescales longer than 1 year) variability of the AMOC upper cell northward limb transport is impacted by intrinsic oceanic processes. As done in Penduff et al. (2011), we mainly compare two global 1/4° OGCM simulations: one 50-year hindcast forced by the full range of observed atmospheric timescales, and one 327-year run forced each year by the same climatological atmospheric annual cycle. The "fully-forced" run is comparable to Hirschi et al. (2013)'s hindcasts and yields the so-called total AMOC variability over the last decades. The seasonally-forced run isolates the AMOC LF variability that is intrinsic, i.e. which is generated by the ocean without any forcing at these timescales. Using this or similar approaches, Penduff et al. (2011) and Thomas and Zhai (2013) showed the substantial imprint of interannual-to-decadal intrinsic variability on sea-surface height and AMOC, respectively. Penduff et al. (2011) also demonstrated that switching to the laminar regime (2°) almost suppresses sea-level intrinsic variability.

The simulation strategies mentioned above yield complementary views on the LF variability that spontaneously emerges in the eddying ocean. Both Biastoch et al. (2008) and Hirschi et al. (2013) used a pair of fully-forced hindcasts, one of which was perturbed (through a local increase of resolution, and through perturbed initial conditions, respectively). Differences between both hindcasts were used by these authors to estimate the imprint of intrinsic variability under full forcing, yet with two limitations: timescales of interest could not extend beyond one decade or so since atmospheric reanalyses have limited duration, and since the estimate of intrinsic variability intensity (derived from 2-member ensemble dispersions) probably has a sizeable inaccuracy. Unlike these authors, Thomas and Zhai (2013) used atmospheric forcing fields devoid of any low-frequency, hence isolating in the North Atlantic the purely intrinsic (i.e. disentangled from its atmospherically-forced counterpart) LF ocean variability over, still relatively short, 15-year simulations. Note that Penduff et al. (2011) also used short 12-year datasets to study the sea-level
intrinsic variability. In the present study, we follow Penduff et al. (2011)'s and Thomas and Zhai (2013)'s climatological forcing strategy to study the purely intrinsic LF variability, but now over a much longer, 327-year global integration: interannual-to-decadal intrinsic variability is first estimated from statistics based on 9 successive 32-year periods; its multi-decadal and intermittent behavior is then analyzed from 309-year timeseries.

Our first objective here is to quantify from the 327-year seasonally-forced simulation the intensity of the intrinsic interannual-to-decadal AMOC variability, its stability over successive 32-year periods, and its relative contribution to the total AMOC variance simulated with full forcing. Our second objective is to discuss from these results the potential contribution of intrinsic variability on existing and future AMOC monitoring data (i.e. RAPID, OVIDE, and OSNAP, SAMOC, respectively; see Schiermeier 2013) and on the meridional heat transport (MHT) at these key latitudes. Our third objective is to examine the sensitivity of these results to model resolution, ranging from 2° (laminar regime) through 1/4° (eddy-permitting) up to 1/12° (eddy-resolving). We will finally describe the temporal behavior of the intrinsic AMOC LF variability throughout the whole seasonally-forced run, including timescales longer than one decade.

Our numerical simulations and post-processing techniques are presented in Section 2. The intensity of intrinsic and total AMOC variabilities will be compared over ~ 1 – 10 year timescales throughout the basin (Section 3), first at 1/4° resolution, then focussing on resolution impacts. The longer-term behavior of the intrinsic variability, its dominant timescales and meridional coherence of its multidecadal timescales will then be studied from the long 327-year seasonally-forced 1/4° simulation (Section 4). Conclusions are given in Section 5.
130

## 2. Model simulations and processing

131

### a. Model simulations

132

133

134

135

136

Each pair of simulations used in this study includes a fully-forced "F" run and a seasonally-forced "S" run that only differ in their forcing. F simulations are intended to mimick the evolution of the observed ocean, and S simulations isolate the intrinsic part of the low-frequency variability. Each pair corresponds to a specific horizontal resolution (2°, 1/4°, and 1/12°) and associated subgrid-scale parameterizations.

137

138

139

140

141

142

143

144

These 6 simulations were performed in the framework of the DRAKKAR project using the Nu-cleus for European Modelling of the Ocean (NEMO; Madec 2008) ocean/sea-ice numerical model with the same 46-level vertical discretization. The 6 simulations share: a partial cell representation of topography and a momentum advection scheme that conserves energy and enstrophy (Barnier et al. 2006; Penduff et al. 2007; Le Sommer et al. 2009); a total variance diminishing (TVD) tracer advection scheme; an isopycnal laplacian tracer diffusion operator; a vertical mixing scheme based on the TKE turbulent closure model (Blanke and Delecluse 1993); and a convective adjustment scheme based on enhanced vertical mixing in case of static instability.

145

146

147

148

149

150

Table 1 summarizes the main differences between these simulations. Note that the 2° and 1/4° configurations and simulations used here are described in detail in Penduff et al. (2010, 2011). The 1/12° model configurations are described in detail in Molines et al. (2014); the seasonally-forced 1/12° simulation has been analyzed in Treguier et al. (2014) and Deshayes et al. (2013). More information about model configurations and solutions may be found in the aforementioned papers as well as in Sérazin et al. (2014).

151

152

Each pair of simulations comprises a 50-year fully-forced "F" run driven by the full range of atmospheric timescales available in original datasets (6-hourly forcing from 1958 to 2007), and
a seasonally-forced "S" run driven by its corresponding climatological annual cycle; the forcing
153
of both simulations within a pair is based on the same original forcing dataset (see Table 1); the
154
method we used to generate the climatological annual cycle from the full forcing is described in
155
Penduff et al. (2011). The original forcing dataset used for the 2° and 1/4° pairs is referred to as
156
DFS4. The original forcing dataset DFS4.4 used for our 1/12° pair is described in detail in Dussin
157
and Barnier (2013) and differs from DFS4 only after December 31st 2001 2.
158

b. AMOC computation

The first 18 years are discarded from all simulations: this removes from AMOC timeseries
160
the imprint of initial dynamical adjustments, and ensures that their long-term adjustment trends
161
are well captured (and removed) by the LOESS high-pass filtering process. Total (fully-forced)
162
AMOC LF statistics are thus computed from F runs over the 32-year 1976-2007 period. Intrinsic
163
(seasonally-forced) AMOC LF statistics are computed from various segments of S-run timeseries,
164
depending on the questions addressed and available data: either from the whole S run (i.e. years
165
19-327 at 1/4°), or from successive 32-year segments: two segments (years 19-50 and 51-82) at
166
1/12°, and nine segments (19-50, 51-82, ..., 275-306) at 1/4° and 2°, that will be combined into a
167
pseudo-ensemble for some statistics.
168

169
AMOCs are computed in both geopotential and potential density referenced to 2000 m, σ₂,
170
coordinates from monthly model outputs. These monthly values are then averaged over 12 months
171
to yield yearly AMOCₓ and AMOCₐ₂ estimates (also see Section 2e). The impact of the model
172
resolution will only be estimated in terms of AMOCₓ since the absence of monthly temperature and
173
salinity (hence density) files prevented us from computing AMOCₐ₂ in certain simulations. All

<footnote n="2">²Use of ERA-interim instead of ERA40 for all forcing variables: wind vector, air temperature and humidity, downward short and long wave radiation, total precipitation and snowfall. Corrections described in Brodeau et al. (2010) ensure a smooth transition in 2001-2002. ERA-interim variables, whose native resolution is 0.7° and 3 hours, were projected at ERA40 resolution (1.125° and 6 hours).</footnote>
other analyses, however, will take advantage of AMOC$_{\sigma_2}$ estimates. In accordance with Zhang (2010), our 1/4° AMOC$_{\sigma_2}$ captures the diapycnal water mass transformation (DWT). At high latitudes in particular (i.e. a region where steep isopycnals are common) there can therefore be substantial differences between AMOC$_{\sigma_2}$ and AMOC$_z$ as the latter represents sinking and not DWT. Similar values are typically found for both AMOC$_{\sigma_2}$ and AMOC$_z$ in low to mid-latitudes where isopycnals are flatter. We shall follow Zhang (2010) on a second point: the transports of the AMOC upper cell northward limb are estimated every year at each latitude as the maximum of AMOC streamfunctions in depth and density (for AMOC$_z$ and AMOC$_{\sigma_2}$, respectively).

### c. Comparison of mean AMOC$_{\sigma_2}$ in 1/4° runs

Fig. 1a-b shows that our F and S 1/4° simulations yield very similar time-averaged AMOC$_{\sigma_2}$ streamfunctions, in particular in the 36-37 kg m⁻³ density range. In this range, the magnitude (~20.5 Sv) and location (55°N, 36.7 kg m⁻³) of both overturning maxima are almost indistinguishable. The meridional distributions of upper cell time-mean intensities (Fig. 1c) differ by less than 2% throughout the whole basin, except between 66 and 70°N where this difference reaches 10%. In other words, the differences in atmospheric forcings barely affect the time-mean AMOC$_{\sigma_2}$ structure and intensity over 50 years, hence ensuring that differences in variabilities obtained in F and S runs are not due to different mean states.

### d. Non-linear detrending

Any OGCM progressively adjusts towards its own "climate", hence explaining non-linear trends in certain timeseries. The blue line in Fig. 2 shows that this adjustment yields a 8% decrease in AMOC$_{\sigma_2}$ at 26.5°N throughout the 327-year 1/4° S run, and progressively tends to level off.
The long-term trends found in $AMOC_z$ and $AMOC_{\sigma_2}$ timeseries in each simulation were removed at all latitudes by fitting local regression models (LOESS: Cleveland et al. 1996; Cleveland and Loader 1996), which acts as a non-linear high-pass temporal filter with negligible end-point effects. In order to remove signals with periods comparable to (or longer than) the length of timeseries, i.e. to ensure proper estimates of variances and spectra, this cutoff period was set to 11 and 110 years for the processing of 32-year and 309-year timeseries, respectively (compare red and blue lines in Fig. 2). This filtering process thus retains timescales between 2 and 11 years in 32-year segments, and between 2 and 110 years in 309-year segments. It is applied on all AMOC timeseries prior to statistical computations (standard deviation (std), spectra, etc).

### e. Impact of temporal sampling on $AMOC_{\sigma_2}$ variability

As mentioned earlier, yearly $AMOC_{\sigma_2}$ were computed from model fields archived at the monthly, "nominal" frequency. Here, we assess the possible contribution of density-velocity covariances at periods shorter than 2 months on the interannual $AMOC_{\sigma_2}$ variability. A "reference" $AMOC_{\sigma_2}$ timeseries was built from 5-day model outputs, available over the final 27 years of the 1/4° S run, and is compared to its "nominal" counterpart. Fig. 3 shows that the nominal $AMOC_{\sigma_2}$ std is similar to its reference counterpart at most latitudes. A ~ 15% difference, however, is found between the Equator and 10°N. Additional investigations (not shown) showed that Tropical Instability Waves, whose typical periods are shorter than 2 months, generate a slowly varying shallow overturning cell in this latitude band; filtering this process in the AMOC computation with monthly model fields induces this 15% difference. Over most of the basin however, our "nominal" dataset provides an estimate of the interannual $AMOC_{\sigma_2}$ std to within a few percent.
f. Intrinsic variability envelopes 216

After detrending AMOC timeseries from all simulations at all latitudes, we computed the total std's ($s_T$) using the 32-year (1976-2007) AMOC timeseries from F runs, and the intrinsic std's ($s_I$) from the nine successive 32-year timeseries available in S runs (only two at 1/12°); let us note that $s_T$, $s_I$ and all derived quantities depend on latitude. Fig. 2 indeed suggests that $s_I$ varies at 26°N during the 1/4° S run 3. At 1/4° and 2° resolutions, our best estimate of $s_I$ is given by the average $\overline{s_I}$ of these nine $s_I$'s; the uncertainty of $\overline{s_I}$ is simply estimated by the standard deviation $s(s_I)$.

The contribution of intrinsic processes to the total AMOC variance will be defined as $R = s_I^2/s_T^2$, and expressed in percent; the best estimate for this quantity is simply $\overline{R} = \overline{s_I^2}/s_T^2$. The statistics introduced so far may also yield an accuracy estimate for R. Let us assumed that the accuracies of $s_T$ and $s_I$ are equal, i.e. $s(s_T) = s(s_I)$, and that a 9-member $s_T$ timeseries obtained from fully-forced passes would be uncorrelated with our 9-member $s_I$ timeseries. This would be consistent with Hirschi et al. (2013)'s paradigm: variations in successive estimates of total AMOC std (with full forcing) are due to intrinsic variability, which is chaotic in nature and thus has a random phase. In that case, the uncertainty (first order approximation) for the ratio R simply reads:

$$s(R) = 2 \frac{s(s_I)}{\overline{s_I}} \overline{R} \sqrt{1 + \overline{R}}. \qquad (1)$$

The intrinsic std's and R estimates will thus be presented in terms of envelopes at 2° and 1/4°, i.e. will refer to $\overline{s_I} \pm s(s_I)$ and to $\overline{R} \pm s(R)$, respectively. Regarding the 1/12° simulations, both available 32-year estimates of $s_I$ and R are shown for AMOC$_z$ as a function of latitude in Fig. 7.

<footnote n="3">³Note that changes in the 9 $s_I$ estimates are decorrelated from the long-term decrease of the original AMOC itself (Fig. 2); these changes are actually due to the presence of intrinsic variability at timescales longer than those accessible from 32-year segments (i.e. 11 year due to the LOESS detrending).</footnote>
234 MHT yearly timeseries where computed at each latitude from the monthly outputs of both 1/4° simulations, and non-linearly detrended as AMOC timeseries. We selected the same 32-year periods as for AMOC, and calculated (one) fully-forced and (nine successive) intrinsic MHT inter-annual variance estimates. As their equivalent R for AMOC, nine M ratios were calculated from these MHT variances, yielding the envelope shown in Fig. 5e, and commented in the next section.
235
236
237
238

3. Total and intrinsic AMOC variability at interannual-to-decadal timescales

The total and intrinsic AMOC variabilities are now compared from the 32-year timeseries de-
rived from (detrended) S and F simulations, i.e. over timescales ranging between 2 and 11 years.
240
241
242
The impact of intrinsic variability on the MHT is also examined over these 32-year periods.

a. Eddy-permitting regime (1/4°): AMOC<sub>σ2</sub>

Let us first compare the variability of the two-dimensional AMOC<sub>σ2</sub> streamfunctions in the 1/4°
simulations F and S (Fig. 4). The largest AMOC variability in the F run (panel a), found north
of 40°N, reaches its maxima along the northward and southward limbs of this subpolar upper cell
(indicated by the white contours), suggesting that the total variability affects the density of the
water masses carried by both limbs. The total variability of the transport of the upper AMOC cell,
i.e. which may be seen in the figures around 36.7 kg m<sup>-3</sup> where the AMOC itself is maximum,
is smaller and only exceeds 1 Sv between 43 and 53°N. The subpolar intrinsic variability (panel
b) exhibits a similar structure with maxima in the southward and northward limbs (density fluc-
tuations) and weak variability of the upper cell transport (near the AMOC maximum). Subpolar
intrinsic variability is thus much weaker than its fully-forced counterpart: the ratio R between
intrinsic and total subpolar AMOC variances remains below 10% (panel c).
244
245
246
247
248
249
250
251
252
253
254
One of the main peaks south of 40°N in fully-forced and intrinsic variabilities is found around the latitude of the Gulf Stream separation (35°N) across the whole density range, including the isopycnal where the AMOC itself is maximum: this indicates that the transport of the upper cell varies in both runs at these latitudes. Fig. 4c shows that about half of the total interannual variability is accounted for by the intrinsic component there. The wide range of densities associated with this AMOC variability peak suggests that the associated circulation fluctuations have a barotropic character, as indeed expected within strong recirculations observed off the US coast (e.g. Bryden et al. 2009). On the other hand, our two 1/4° simulations, like many eddy-permitting model solutions, exhibit an unrealistic, almost standing eddy near Cape Hatteras (see e.g. Barnier et al. 2006); whether this half-intrinsic AMOC variability peak is robust and realistic is therefore questionable.

A weaker intrinsic variability peak is found near the RAPID array (25°N) with a substantial contribution (40%) to the fully-forced AMOC variance. This maximum does not seem as suspicious as that mentioned above since no obvious discrepancy is found in the model solutions; the robustness of this peak at 1/12° will be assessed in the following.

The largest fraction of intrinsic variability is found at the southern end of the domain throughout most density classes. This peak is likely to be linked to the spontaneous interannual variability in the Agulhas leakage (e.g. Biastoch et al. 2008). The large intrinsic fraction (R) seen in light waters there suggest that in the upper South Atlantic the variability of the meridional circulation can only partly be predicted from surface forcing. The large surface intrinsic contribution is confined to latitudes south of 25°S, suggesting that the variability of the horizontal surface circulation can largely be predicted from the surface forcing further north. In deeper layers associated with the AMOC upper cell (its southward limb), the intrinsic contribution remains quite large (30-50%) between 25°N and the Equator, which is similar to the results found for AMOC<sub>z</sub> by Hirschi et al. (2013).
From now onwards we shall restrict our variability analyses to timeseries of $AMOC_{\sigma_2}$ maxima (in density) at each latitude and each year. The timeseries correspond to the evolution of the strength of the northward limb of the upper AMOC cell described in the introduction. Fig. 5a shows that the total $AMOC_{\sigma_2}$ std in the fully-forced 1/4° simulation globally increases northward with typical values of 0.6-0.8 Sv. It reaches its maxima in the North Atlantic Current between 45 and 50°N (~ 1.13 Sv) and in the subtropics between 10 and 20°N (~ 0.9 Sv). Its absolute minimum (~ 0.35 Sv) lies poleward of 65°N <fnref n="4" />, with 3 relative minima (~ 0.6 Sv) within 23-33°N, around the equator and south of 27°S.

Fig. 5b shows the envelope of intrinsic $AMOC_{\sigma_2}$ std, derived from the 9 $s_I$ estimates computed over the successive 32-year periods: the plain and dashed lines respectively represent $\bar{s}_I$, $\bar{s}_I - s(s_I)$ and $\bar{s}_I + s(s_I)$ as defined in Section 2f. The dispersion $s(s_I)$ remains close to 10-15% of its average $\bar{s}_I$ throughout the basin, showing a rather uniform meridional distribution of $AMOC_{\sigma_2}$ std dispersion among the nine 32-year segments of the S 1/4° run. Switching from full to seasonal forcing (i.e. from the F to S run) reduces the interannual AMOC variance, but does not suppress it, especially at the basin's southern end where the intrinsic std reaches 0.4 Sv.

The contribution $R$ of intrinsic processes to the total $AMOC_{\sigma_2}$ variance in the 1/4° simulation is shown in Fig. 5c as its envelope defined in Section 2f. Roughly 15-30% of the total variance $s_T^2$ happen to be due to intrinsic processes over the Atlantic tropics and subtropics, with a marked increase of $R$ from north to south. Three additional features are visible in this panel and might be relevant for the interpretation of observational MOC variability estimates. Interestingly, and despite substantial differences in experimental strategies and model setups, the 3 following features are also visible in Hirschi et al. (2013).

<footnote n="4">⁴A slight local maximum in $\sigma_T$ may however be seen along the 36.7 kg m⁻³ isopycnal in Fig. 4a near this latitude: this difference is explained by temporal changes in isopycnals along which $AMOC_{\sigma_2}$ reaches its successive maxima, which are in turn used in Fig. 5a.</footnote>
* R reaches its absolute maximum (50±15%) around 30°S. This latitude sits a few degrees north of the SAMOC monitoring array, which is currently being deployed at 34.5°S, and where R remains substantial (40±10%). There, a substantial part of the interannual AMOC variability may not be directly related to that of the atmosphere and purely due to oceanic processes.

* Between 24 and 26°N, the total and intrinsic AMOC$_{σ_2}$ variabilities reach a local minimum and maximum, respectively: R reaches its secondary maximum (40±10%) at 25°N, and R = 30 ± 10% at the latitude of the RAPID array. Whether this result is robust to other model parameters remains to be checked, but these results suggest that interannual AMOC variability measured at 26.5°N may be affected by ocean-driven processes.

* As expected from the opposite basin-scale gradients of total and intrinsic std's, intrinsic processes reach their minimum contribution (~5%) between 45 and 60°N. This subpolar area almost entirely encompasses the paths of OVIDE tracks along which several observational AMOC$_{σ_1}$ estimates have been derived since 2002 (Treguier et al. 2006; Mercier et al. 2013). The forthcoming OSNAP monitoring array will also be deployed within this gyre. Despite the misalignment between non-zonal OVIDE tracks and our quasi-zonal model grid there, and despite the difference in density reference levels, our results suggest that in this area almost all the interannual-to-decadal AMOC variance is related with the atmosphere with a minor contribution of oceanic processes. Interestingly, R at 1/4° reaches a localized maximum around 65°N (near the Grenland-Iceland-Scotland ridge), where the total variability $s_T^2$ reaches a minimum. A local R maximum is found right above this ridge in our 1/12° experiments as well (R ~ 23 %, Fig. 6). This topographic feature may therefore host oceanic processes in our 1/4° and 1/12° simulations that locally increase R, i.e. which moderate the strong connection between the atmosphere and interannual AMOC variance.
The interannual intrinsic variabilities of $AMOC_{\sigma_2}$ and MHT are significantly correlated across a wide range of latitudes over the 9 successive 32-year periods of the 1/4° S run (Fig. 5d): 50-74 % of the MHT and $AMOC_{\sigma_2}$ intrinsic variances are shared between 34 and 5°S, 18 and 43°N, and 45 and 54°N. The intrinsic MHT and $AMOC_{\sigma_2}$ variabilities have comparable meridional distributions (not shown) at these latitude ranges. The ratio *M* between intrinsic and total MHT variances (its envelope is shown in Fig. 5e) is substantial in the South Atlantic, as was noted for *R*: *M* reaches its absolute maximum (~26±7%) near the SAMOC latitude around 30°S. Between 15°S and 23°N, intrinsic processes contribute to less than 10% of the MHT variance, due to a marked maximum of the fully-forced MHT variability (not shown). *M* reaches ~12±3% at RAPID (26.5°N) and ~23±6% in the Gulf Stream extension (37-42°N). As for $AMOC_{\sigma_2}$, the contribution of intrinsic processes to the total MHT variance reaches its absolute minimum north of 45°N where OVIDE and OSNAP measurements have and will be made. The interannual variability of the northward heat transport is therefore mostly controlled by the atmosphere in the subpolar gyre, but is partly due to intrinsic processes at mid latitudes.

### b. Laminar, eddy-permitting and eddy-resolving regimes: $AMOC_z$

As explained earlier, we switch from $AMOC_{\sigma_2}$ to $AMOC_z$ within this section. Let us first investigate how this switch affects overturning variabilities by comparing Fig. 5 and 7 in terms of total std's (panels a) at 1/4° (black lines). The total standard deviations of both quantities at 1/4° are very close to each other south of 40°N, a result similar to that reported by Zhang (2010, her Fig. 2b) in the laminar regime. Poleward of 40°N however, std's are systematically larger in terms of $AMOC_{\sigma_2}$ than of $AMOC_z$. This reflects the fact that $AMOC_{\sigma_2}$ variabilities integrate both the zonally-integrated vertically-accumulated transport and the horizontal gyre circulation, and that the contribution of this latter component increases toward subpolar latitudes (e.g. Marsh et al.
347 These authors also reported stronger AMOC$_z$ values at subpolar latitudes in coarse reso-
lution simulations, as most of the meridional transports occur in the overturning part in laminar
348 regimes but are more evenly split between gyre and overturning parts in eddying regimes. This
349 resolution impact is confirmed by our results in terms of interannual variability (stronger at 2°, red
350 line in Fig. 7a).
351

We now investigate how the model resolution influences intrinsic variability. Fig. 7b demon-
352 strates that switching from the eddy-permitting (1/4°) to the laminar regime (2°) almost entirely
353 suppresses interannual intrinsic variability of AMOCz, as it does for sea-level (Penduff et al. 2011).
354 The basin's intrinsic variability minimum is found north of 40°N both at 1/4° and 1/12°, but is
355 stronger in the latter case. At least in this subpolar region, the intensity of interannual intrin-
356 sic variability therefore seems to increase with the ability of the model to (even partly) resolve
357 mesoscale eddies. Fig. 5a indicates that this does not imply that the total AMOC std increases
358 under full forcing when switching from 1/4° to 1/12° resolution. Interestingly, the 1/4°-to-1/12°
359 resolution increase barely changes the general distribution and magnitude of this intrinsic std be-
360 tween 35°S and 40°N: over most of this latitude range the two estimates of intrinsic std at 1/12°
361 remain within (or very close to) the 1/4° envelope. Fig. 7b shows that in both eddying simula-
362 tions, the strongest imprint of LF intrinsic variability on AMOCz occurs around the latitude of the
363 Gulf Stream; substantial and relatively uniform σl levels are then found between this latitude and
364 35°S. This meridional distribution is in contrast with the imprint of LF intrinsic variability on sea-
365 surface height (SSH) at both 1/4° and 1/12° (see Penduff et al. 2011; Sérazin et al. 2014), which
366 exhibits very distinct peaks in the highly turbulent Gulf Stream and Agulhas Retroflexion regions,
367 with almost no SSH imprint in between. Both the horizontal and overturning circulations are thus
368 impacted by LF intrinsic variability in these two eddy-active regions. Its relatively uniform im-
369 print on AMOC between these latitudes is likely related with the marked meridional coherence of
370
371 AMOC intrinsic variability anomalies across this latitude range at interannual-to-decadal periods (as discussed in Section 4d for longer timescales).

373 Fig. 7c summarizes the impact of resolution on *R*, the intrinsic-to-total interannual variance ratio. The interannual variability of the forcing explains nearly all interannual *AMOC*<sub>*z*</sub> variance at 2° throughout the basin and remains largely dominant (up to 95%) in the Subpolar Gyre when 375 eddies are present (1/4° and 1/12°). In both eddying simulations however, the contribution of 376 nonlinear intrinsic processes in the total AMOC variance tends to increase toward the south, and 377 reaches about one third around 35°N up to 40-50% in the Agulhas region. Our comments about 378 the *R* extrema found at OVIDE and SAMOC in terms of *AMOC*<sub>*σ*₂</sub> (previous section) thus remain 379 valid in terms of *AMOC*<sub>*z*</sub>, but *R* is smaller at RAPID (20±5%) than in terms of *AMOC*<sub>*σ*₂</sub>: intrinsic 380 variability here is partly associated with meridional flows occuring at different densities but same 381 depths. 382

383 Thomas and Zhai (2013) used a different model (MITgcm), resolution (1/10°) and experimental 384 strategy, but comparing our Fig. 7c with the black line in their Fig. 4 is instructive. First, their 385 estimate of the *R* ratio in the eddying regime compares quite well with ours in the 20-45°N latitude 386 range (and with our 1/12° results at subpolar latitudes). Second, this comparison strongly suggests 387 that most of the intrinsic interannual variability in the northern tropics (0-20°N) comes from the 388 Southern basin. Indeed, their use of a closed southern boundary at 14°S yields a negligible intrinsic 389 interannual variability at low latitudes, while it remains substantial and homogeneous across 34°S- 390 30°N in our global configurations as well as in Hirschi et al. (2013)'s. This result is also consistent 391 with Biastoch et al. (2008) who showed how intrinsic AMOC variability generated in the turbulent 392 Agulhas region can influence the tropical north Atlantic.
393

Interestingly, our 1/4° model proves able to capture the essential processes generating intrinsic variability since 1/12° and 1/4° results do not differ much from each other, at least south of the subpolar gyre. Most of our results there are also robust to the use of $AMOC_z$ or $AMOC_{\sigma_2}$.

394

395

396

## 4. Long-term $AMOC_{\sigma_2}$ intrinsic variability at 1/4°

397

We now take advantage of the 327-year seasonally-forced 1/4° simulation to describe the long-term evolution of the intrinsic $AMOC_{\sigma_2}$ variability. We first highlight its general features through a time-frequency analysis at 4 selected latitudes (a); we then focus on the evolution of interannual-to-decadal intrinsic variability (b) and its possible relationship with the atmospheric variability (c), and on the time-latitude distribution of the multidecadal intrinsic variability throughout the basin (d).

398

399

400

401

402

403

### a. Time-frequency analysis of the intrinsic variability

404

405

406

407

408

409

410

411

Detrended 309-year timeseries of intrinsic $AMOC_{\sigma_2}$ anomalies are presented at four latitudes in Fig. 8, along with corresponding wavelet power spectra. Morlet wavelets (see Torrence and Compo 1998) were chosen, spectra were normalized by $AMOC_{\sigma_2}$ LF variances, and confidence levels, shown as contours, were obtained by fitting lag-1 autoregressive red noise processes. The four chosen latitudes correspond to the observational arrays mentioned earlier, and/or to intrinsic $AMOC_{\sigma_2}$ extrema visible in Fig. 5: both R maxima found at SAMOC and RAPID (panels a and b), the R minimum around OVIDE (c) and the R maximum above the Greenland-Iceland-Scotland ridge (d).

412

413

The main feature visible in these wavelet spectra (and at other latitudes as well) is the temporal irregularity of the frequency content throughout the 1/4° S run. Significant spectral peaks appear and disappear during the run; these occurrences are not particularly confined at the begin-

414
415
416
417
418

ning or the end of the simulation, suggesting that the long-term adjustment of the model state
does not adversely contaminate these results. The wavelet transforms also show that this intrinsic
AMOC<sub>σ₂</sub> variability concerns a wide frequency range, with significant power at interannual-to-
decadal timescales throughout the Atlantic.

419

b. *Interannual-to-decadal intrinsic variability*

420
421
422
423
424
425
426
427
428
429

AMOC$_{\sigma_2}$ Fourier spectra are computed at each latitude over the nine successive 32-year seg-
ments of the 1/4° S run (left panels in Fig. 9). Confidence levels are estimated as in Torrence and
Compo (1998) by modeling background spectra as red noise (lag-1 autoregressive process). These
panels reveal generic features of the intrinsic variability that persist throughout the simulation: at
low and mid latitudes, significant spectral peaks are visible most of the time, extending over wide
(up to 30°N) latitudinal ranges at well-defined frequencies. These "ridges" are particularly clear
south of 25°N: several of them happen to span the whole South Atlantic at about 3, 5, and 8-year
periods. Comparable structures also appear in the northern subtropics, with somewhat shorter
meridional extents though. North of about 40°N, intrinsic variability has a reduced contribution to
the AMOC variance; its frequency peaks are more confined in latitude and fluctuate more in time.

430
431
432
433
434
435
436
437

These results are summarized in color in the lower right panel of Fig. 9, which shows the prob-
ability that intrinsic variability reaches a significant spectral peak during these 9 periods. Ev-
erywhere south of about 25°N, intrinsic *AMOC*<sub>σ₂</sub> variability often occurs (30-50% of the 309-
year period) at 8 and 4.5 year periods; additional occurrences are detected around 3 years in the
South Atlantic, suggesting that future SAMOC measurements may include intrinsic AMOC sig-
nals around these 3 interannual peaks. Between 25 and 40°N, significant spectral peaks are mostly
found around 5-6 years. More detailed investigations in our 1/4° and 1/12° S runs confirm that
AMOC at the RAPID array (26.5°N) has intrinsic variability peaks near 4-5 and 10 years. This
438
439
440
441

variability is in line with Rossby-wave-like intrinsic modes that continuously cross the basin be-
tween 15 and 35°N (see Penduff et al. 2011, their section 4.a.2) in both simulations, and whose
spatio-temporal features are very similar to those found by Hazeleger and Drijfhout (2000) in
idealized simulations.

442
443
444
445
446
447

Significant peaks occur at more diverse periods further north; however, intrinsic variability
seems to preferably peak around 3.5 and 8 years poleward of 45°N, with slightly longer peri-
ods (5 and 10 years) at the northern boundary of the basin. In summary, the amplitude, latitudinal
extent, and dominant frequencies of intrinsic *AMOC*<sub>σ₂</sub> variability vary in a substantial, complex
and intermittent way during the integration, as expected from the irregular shape of the wavelets
(Fig. 8).

448
449
450
451
452
453
454
455

The impact of model resolution on interannual intrinsic variability spectra is difficult to assess accurately. *AMOC*<sub>*z*</sub> intrinsic variability spectra at 1/12° are indeed available over two 32-year periods only (not shown), and are quite different from each other (as in the 1/4° case, see Fig. 9).

Although the intensity of intrinsic variability is similar over these successive periods (Fig. 7b), no robust estimate of the 1/12° intrinsic variability spectral peaks may be derived from these two realizations. However, meridionally-coherent spectral peaks appear between about 35°S and 35°N at both resolutions and over many 32-year periods of the S runs, suggesting that this feature of intrinsic variability is robust.

456

c. *Possible impacts of the atmospheric forcing*

457
458
459
460

The spectrum of AMOC<sub>σ₂</sub> variability was also computed over 1976-2007 from the fully-forced
1/4° run; its significant peaks are superimposed in white in the lower right panel of Fig. 9 in order
to assess potential correspondences between AMOC preferred periods, with and without inter-
annual forcing. Interestingly, both simulations exhibit common spectral peaks in a few regions:
around 5-6 years at the northern limit of the basin (68-70°N), as well as around 30-40°N: Gulf Stream surface velocity observations (Ezer et al. 2013) confirm the actual existence of this latter peak, which might therefore be of intrinsic origin. Correspondence between seasonally- and fully-forced AMOC spectral peaks is less clear but may also be present in the subpolar gyre, around 8 years between 45 and 60°N; this peak could therefore be of intrinsic origin but also excited by the North Atlantic Oscillation (NAO; see Eden and Willebrand 2001). In this region, we may hypothesize that the interannual forcing excites and paces intrinsic variability modes that spontaneously emerge in the S run. A similar paradigm was proposed by Pierini (2014) in the Kuroshio, but assessing its relevance in the present case would require further investigations.

At lower latitudes, white contours reveal that switching on the interannual forcing does not affect the meridional alignment of spectral peaks. This feature was clearly noticed in the S run, strongly suggesting that the meridional coherence of the forced (and observed) interannual-to-decadal AMOC variability is inherent to oceanic dynamics. These intrinsic and fully-forced spectral peaks, however, are found at different periods at mid and low latitudes. Whether these correspondences and differences in AMOC variability peaks are realistic and robust is unclear: an ensemble of simulations (e.g. different resolutions, initial conditions, parameters) would be necessary to assess whether AMOC variability periods are biased in our S run and may better match their fully-forced counterparts, or if the fully-forced periods would be the same in different conditions. On the other hand, frequency shifts due to the forcing are known to exits in nonlinear systems like the eddying ocean (e.g. Kartashova 2010), and such effects might also explain these differences. These assessments are left for future studies.
d. *Multidecadal intrinsic variability*

Wavelet power spectra (Fig. 8) have shown that intrinsic *AMOC*<sub>σ₂</sub> variability in the 1/4° S run can also be significant at multidecadal periods: 15-30 years at 64°N and 30°S, up to 50-90 years at 54.3°N and 30°S (panels a, c, d). Multidecadal peaks are intermittently visible at other latitudes as well but they do not appear as significant in Fig. 8, perhaps masked by more energetic, simultaneous shorter-period signals. In order to highlight the meridional-temporal structure of these multidecadal signals, intrinsic *AMOC*<sub>σ₂</sub> timeseries were bandpass-filtered at each latitude using a temporal Lanczos filter in the 15-30 year and the 50-90 year bands (see Fig. 10); these plots reveal that these intermittent multidecadal oscillations have a substantial coherence over large meridional extents.

The 15-30 year intrinsic AMOC anomalies are strongest at the southern end of the domain. They extend over the whole south Atlantic (e.g. years 50-80 in the top panel), often reach the RAPID latitude without much attenuation (see at 26°N, middle of the top panel), with occasional extensions further north. In contrast, the 50-90 year intrinsic AMOC anomalies tend to be strongest in the 30-60°N latitudinal band. Hovmoeller diagrams in both frequency bands show that most AMOC intrinsic anomalies propagate in latitude, but at phase speeds that happen to vary in space and time, and with occasional changes in direction. South of 30°N for instance, most anomalies propagate northward, with an acceleration of phasespeeds (phase lines are closer to the vertical) after year 210 in the 15-30 year band, or even a tendency for southward propagation before year 130 in the 50-90 year band. Despite these features, which further illustrate the non-stationary character of intrinsic variability, most AMOC signals in this 30°S-30°N latitude band propagate northward in both frequency bands with phase speeds lying within the range 3-30 cm s<sup>-1</sup>, as reported by Biastoch et al. (2008) at shorter (decadal) timescales.
The coherence of AMOC anomalies in the 15-30 year band is somewhat blurred by meridionally-confined disturbances in the 30-50°N latitudinal range, but remains much clearer in the 50-90 year band. Interestingly, most intrinsic AMOC anomalies seem to propagate southward from 50°N to 30°N in this intergyre region, in particular in the 50-90 year band (at about 3 cm s⁻¹), and in the 15-30 year band between years 120 and 260. North of 50°N in the subpolar gyre, meridional propagation of AMOC anomalies remains dominantly northward in both frequency bands (also in the 6-15 and 30-50 year bands, not shown); these phase speeds are smaller (0.5-1.5 cm s⁻¹) at these high latitudes, and may be related with the poleward advection of anomalous isopycnal transports across the eastern North Atlantic. A number of intrinsic AMOC anomalies thus have an approximate C-shaped pattern north of 30°N at these multidecadal scales (as well as in the 30-50 year period range, not shown), and seem to propagate away from the 50°N latitude line. Whether this feature is also found in other OGCM simulations needs to be verified. Simplified model simulations would then help identify the possible source of these signals, which approximately lies at the latitude where the North Atlantic Current crosses the basin. More generally, explaining the nonlinear dynamics that are likely to generate and propagate these superimposed, non-stationary and complex multidecadal signals lies beyond the scope of this study, and is left for future process-oriented investigations.

# 5. Conclusions and discussion

Three global ocean/sea-ice model simulations, forced at the surface by a repeated annual cycle, have been analyzed in terms of AMOC intrinsic variability at interannual-to-multidecadal timescales. Intrinsic variability emerges spontaneously at these timescales when mesoscale eddies are (even partly) resolved, i.e. in our 1/4° and 1/12° resolution models. At 2° resolution, the same model does not produce any low-frequency spontaneously. Equatorward of 40°N, the
intrinsic AMOC low-frequency variance $s_T^2$ does not significantly increase when switching from 1/4° to 1/12° resolution, suggesting that its dynamical sources are properly simulated at low- and mid-latitudes in the eddy-permitting regime.

The same 3 model configurations were also driven over the last decades by a forcing based on atmospheric reanalyses which now include interannual and longer timescales. These fully-forced hindcasts yield estimates at 3 resolutions of the "total" interannual-to-decadal AMOC variability, which was compared to its intrinsic counterpart. At 1/4° and 1/12°, about one-third to one-half of the fully-forced interannual-to-decadal AMOC variance ($s_T^2$) is accounted for by intrinsic processes around 25°N and 30°S, i.e. close to the RAPID and SAMOC in-situ arrays that are presently (or will soon be) monitoring the AMOC variability. This suggests that only a part of the interannual AMOC variance at these two latitudes may be explained by atmospheric causes. In the subpolar North Atlantic, on the other hand, our results suggest that overturning observational estimates, such as those derived from the existing OVIDE (and future OSNAP) measurements, are strongly related to the atmospheric interannual variability since the intrinsic component only accounts for 5-10% of the total AMOC variance. About a quarter of the meridional heat transport interannual variability happens to be driven by intrinsic processes near the Gulf Stream and at the southern opening of the basin.

Wavelet power spectra revealed that the AMOC intrinsic variability exhibits an intermittent, stochastic character throughout the 327-year 1/4° seasonally-forced simulation. At mid and low latitudes, in the 2-11 year range, marked spectral peaks spread over wide latitude ranges (as in idealized model studies, e.g. Hazeleger and Drijfhout 2000), but the intensity and frequency of these spectral peaks slowly vary in time, apparently in a random way. The cumulated time over which the most recurrent intrinsic, latitudinally-extended spectral peaks are present and significant amounts to about half of the 327 years. In the North Atlantic, some recurrent intrinsic variability peaks
552

553

554

overlap the peaks of the AMOC driven by full forcing, suggesting that interannual atmospheric variability modes, such as the North Atlantic Oscillation, may excite certain intrinsic variability modes in the Gulf Stream or the subpolar gyre.

555

556

557

558

559

560

561

562

Intrinsic AMOC variability extends beyond decadal timescales. Irregular oscillations are found at the longest periods accessible from our dataset, extending all the way from 35°S to 25°N in the 15-30 year band, and throughout the whole Atlantic basin in the 50-90 year band. In both ranges of periods, multidecadal AMOC anomalies tend to propagate northward between 34°S and about 30°N, but seem to propagate northward and southward away from the North Atlantic Current system. Further investigations are needed to identify and quantify the multivariate imprints of intrinsic variability, the dynamical processes involved in its generation, structure, propagation and non-stationarity, and its dependency to model formulation.

563

564

565

566

567

568

569

570

571

572

573

574

The analysis of our climatological simulations, complemented by Sérazin et al. (2014) in terms of sea level intrinsic variability, raise additional open questions. The appearance of a substantial intrinsic variability when oceanic eddies are resolved suggests that the next generation of coupled models with eddying ocean components will exhibit additional (or perhaps stronger) low-frequency variability modes. Intensified or additional stochasticity might also be expected in future coupled simulations given the intermittency spontaneously generated by the model during 327 years. These results also echo conclusions by Hirschi et al. (2013) and Thomas and Zhai (2013): the low-frequency AMOC variability is not fully deterministic but stochastic in the eddying regime. These findings claim for the use of probabilistic approaches and ensemble modelling strategies in climate-oriented oceanographic studies as proposed e.g. in Penduff et al. (2014), potentially using 1/4° models which proved able to capture most of the intrinsic AMOC variance obtained at 1/12°.
Acknowledgments. The authors wish to thank Jean-Michel Brankart, Anne-Marie Treguier, Lau-
rent Terray, and George Nurser for interesting discussions and suggestions. This work is a contri-
bution to MyOcean2, CHAOCEAN and OCCIPUT projects. It benefited from the DRAKKAR in-
ternational coordination network (GDRI) established between the Centre National de la Recherche
Scientifique (CNRS), the National Oceanography Centre in Southampton (NOCS), GEOMAR
in Kiel and IFREMER. The research leading to these results has received funding from the
European Community's Seventh Framework Programme FP7/2007-2013 under grant agreement
n°283367 (MyOcean2). CHAOCEAN is supported by the Centre National d'Études Spatiales
(CNES) through the Ocean Surface Topography Science Team (OST/ST). OCCIPUT is supported
by the Agence Nationale de la Recherche (ANR) through contract ANR-13-BS06-0007-01. The
computations presented in this study were performed at the Centre Informatique National de
l'Enseignement Supérieur (CINES) and at the Institut du Développement et des Ressources en In-
formatique Scientifique (IDRIS) under allocations made by GENCI. The ORCA025-MJM01 and
ORCA12-GJM02 simulations were performed as part of the Grands Challenges GENCI/CINES
2008 and 2013, respectively. SG is supported by MyOcean2, GS is supported by CNES and
Région Midi-Pyrénées, and TP by CNRS. We also thank two anonymous reviewers for their con-
structive comments that helped us improve the paper.

References
---

Baehr, J., J. Hirschi, J.-O. Beismann, and J. Marotzke, 2004: Monitoring the meridional overturn-
ing circulation in the north atlantic: A model-based array design study. *J. Mar. Res.*, **62** (3),
283–312, doi:10.1357/0022240041446191.

Barnier, B., and Coauthors, 2006: Impact of partial steps and momentum advection schemes in
a global ocean circulation model at eddy-permitting resolution. *Ocean Dynam.*, **56** (5-6), 543–
598

567, doi:10.1007/s10236-006-0082-1.

599

Bentsen, M., H. Drange, T. Furevik, and T. Zhou, 2004: Simulated variability of the atlantic meridional overturning circulation. *Climate Dynam.*, **22** (6-7), 701–720, doi:10.1007/s00382-004-0397-x.

600

601

602

Biastoch, A., C. W. Böning, and J. R. E. Lutjeharms, 2008: Agulhas leakage dynamics affects decadal variability in atlantic overturning circulation. *Nature*, **456** (7221), 489–492, doi:10.1038/nature07426.

603

604

605

Blanke, B., and P. Delecluse, 1993: Variability of the tropical atlantic ocean simulated by a general circulation model with two different mixed-layer physics. *J. Phys. Oceanogr.*, **23** (7), 1363–1388, doi:10.1175/1520-0485(1993)023<1363:VOTTAO>2.0.CO;2.

606

607

608

Brodeau, L., B. Barnier, A.-M. Treguier, T. Penduff, and S. Gulev, 2010: An era40-based atmospheric forcing for global ocean circulation models. *Ocean Modell.*, **31** (3-4), 88–104, doi: 10.1016/j.ocemod.2009.10.005.

609

610

611

Bryden, H. L., A. Mujahid, S. A. Cunningham, and T. Kanzow, 2009: Adjustment of the basin-scale circulation at 26°N to variations in gulf stream, deep western boundary current and ekman transports as observed by the rapid array. *Ocean Sci.*, **5** (4), 421–433, doi: 10.5194/os-5-421-2009.

612

613

614

615

Cleveland, W. S., E. Grosse, and M.-J. Shyu, 1996: A package of c and fortran routines for fitting local regression models: Loess user's manual. Technical report, Bell Labs.

616

617

Cleveland, W. S., and C. Loader, 1996: Smoothing by local regression: Principles and methods. *Statistical Theory and Computational Aspects of Smoothing*, W. Härdle, and M. G. Schimek, Eds., Contributions to Statistics, Physica-Verlag HD, 10–49.

618

619
620

Cravatte, S., G. Madec, T. Izumo, C. Menks, and A. Bozec, 2007: Progress in the 3-d circulation of the eastern equatorial pacific in a climate ocean model. *Ocean Modell.*, **17** (1), 28-48, doi: 10.1016/j.ocemod.2006.11.003.

621

622

623

Cunningham, S. A., and Coauthors, 2007: Temporal variability of the atlantic meridional over-
turning circulation at 26.5°n. *Science*, **317** (5840), 935-938, doi:10.1126/science.1141304.

624

625

Delworth, T., S. Manabe, and R. J. Stouffer, 1993: Interdecadal variations of the thermohaline circulation in a coupled ocean-atmosphere model. *J. Climate*, **6** (11), 1993–2011, doi:10.1175/1520-0442(1993)006<1993:IVOTTC>2.0.CO;2.

626

627

628

Deshayes, J., and Coauthors, 2013: Oceanic hindcast simulations at high resolution suggest that the atlantic moc is bistable. *Geophys. Res. Lett.*, **40** (12), 3069–3073, doi:10.1002/grl.50534.

629

630

Dong, S., S. L. Hautala, and K. A. Kelly, 2007: Interannual variations in upper-ocean heat content and heat transport convergence in the western north atlantic. *J. Phys. Oceanogr.*, **37** (11), 2682–2697, doi:10.1175/2007JPO3645.1.

631

632

633

Dufour, C. O., J. Le Sommer, J. D. Zika, M. Gehlen, J. C. Orr, P. Mathiot, and B. Barnier, 2012: Standing and transient eddies in the response of the southern ocean meridional overturning to the southern annular mode. *J. Climate*, **25** (20), 6958–6974, doi:10.1175/JCLI-D-11-00309.1.

634

635

636

Dussin, R., and B. Barnier, 2013: The making of DFS 5.1. Drakkar project report, LGGE, Grenoble, France, 1–40 pp. URL http://www.drakkar-ocean.eu/forcing-the-ocean.

637

638

Eden, C., and J. Willebrand, 2001: Mechanism of interannual to decadal variability of the north atlantic circulation. *J. Climate*, **14** (10), 2266–2280, doi:10.1175/1520-0442(2001)014(2266:MOITDV)2.0.CO;2.

639

640
641

642

643

Ezer, T., L. P. Atkinson, W. B. Corlett, and J. L. Blanco, 2013: Gulf stream's induced sea level rise and variability along the u.s. mid-atlantic coast. *J. Geophys. Res. Oceans*, **118** (2), 685-697, doi:10.1002/jgrc.20091.

644

645

Hall, M. M., and H. L. Bryden, 1982: Direct estimates and mechanisms of ocean heat transport. *Deep Sea Res.*, **29** (3), 339-359, doi:10.1016/0198-0149(82)90099-1.

646

647

648

Hazeleger, W., and S. S. Drijfhout, 2000: A model study on internally generated variability in
subtropical mode water formation. *J. Geophys. Res.*, **105 (C6)**, 13965–13979, doi:10.1029/2000JC900041.

649

650

651

Hirschi, J., J. Baehr, J. Marotzke, J. Stark, S. Cunningham, and J.-O. Beismann, 2003: A monitoring design for the atlantic meridional overturning circulation. *Geophys. Res. Lett.*, **30** (7), doi:10.1029/2002GL016776.

652

653

654

Hirschi, J. J.-M., A. T. Blaker, B. Sinha, A. Coward, B. de Cuevas, S. Alderson, and G. Madec, 2013: Chaotic variability of the meridional overturning circulation on subannual to interannual timescales. *Ocean Sci.*, **9** (5), 805–823, doi:10.5194/os-9-805-2013.

655

656

Johns, W. E., and Coauthors, 2011: Continuous, array-based estimates of atlantic ocean heat transport at 26.5°N. *J. Climate*, **24** (10), 2429-2449, doi:10.1175/2010JCLI3997.1.

657

658

Kanzow, T., and Coauthors, 2007: Observed flow compensation associated with the moc at 26.5°n
in the atlantic. *Science*, **317** (5840), doi:10.1126/science.1141293.

659

660

Kartashova, E., 2010: *Nonlinear Resonance Analysis: Theory, Computation, Applications*. Cambridge University Press, Cambridge.
661

Le Sommer, J., T. Penduff, S. Theetten, G. Madec, and B. Barnier, 2009: How momentum advection schemes influence current-topography interactions at eddy permitting resolution. *Ocean Modell.*, **29** (1), 1–14, doi:10.1016/j.ocemod.2008.11.007.

662

663

664

Lherminier, P., H. Mercier, C. Gourcuff, M. Alvarez, S. Bacon, and C. Kermabon, 2007: Transports across the 2002 greenland-portugal ovide section and comparison with 1997. *J. Geophys. Res.*, **112** (C7), C07 003, doi:10.1029/2006JC003716.

665

666

667

Madec, G., 2008: NEMO ocean engine. Note du Pole de modlisation 27, Institut Pierre-Simon Laplace (IPSL), France.

668

669

Marsh, R., B. A. de Cuevas, A. C. Coward, J. Jacquin, J. J.-M. Hirschi, Y. Aksenov, A. J. G. Nurser, and S. A. Josey, 2009: Recent changes in the north atlantic circulation simulated with eddy-permitting and eddy-resolving ocean models. *Ocean Modell.*, **28** (4), 226–239, doi:10.1016/j.ocemod.2009.02.007.

670

671

672

673

McCarthy, G., and Coauthors, 2012: Observed interannual variability of the atlantic meridional overturning circulation at 26.5°N. *Geophys. Res. Lett.*, **39** (19), L19 609, doi:10.1029/2012GL052933.

674

675

676

Mercier, H., and Coauthors, 2013: Variability of the meridional overturning circulation at the greenland-portugal ovide section from 1993 to 2010. *Prog. Oceanogr.*, doi:10.1016/j.pocean. 2013.11.001.

677

678

679

Molines, J.-M., B. Barnier, T. Penduff, A.-M. Treguier, and J. Le Sommer, 2014: Orca12.146 climatological and interannual simulations forced with dfs4.4: Gjm02 and mjm88. Drakkar project report, LGGE, 1–50 pp. URL http://www.drakkar-ocean.eu/global-models/orca12.

680

681
682

Penduff, T., M. Juza, B. Barnier, J. Zika, W. K. Dewar, A.-M. Treguier, J.-M. Molines, and N. Au-diffren, 2011: Sea level expression of intrinsic and forced ocean variabilities at interannual time scales. *J. Climate*, **24 (21)**, 5652–5670, doi:10.1175/JCLI-D-11-00077.1.

683

684

685

Penduff, T., M. Juza, L. Brodeau, G. C. Smith, B. Barnier, J.-M. Molines, A.-M. Treguier, and G. Madec, 2010: Impact of global ocean model resolution on sea-level variability with emphasis on interannual time scales. *Ocean Sci.*, **6 (1)**, 269–284, doi:10.5194/os-6-269-2010.

686

687

688

Penduff, T., J. Le Sommer, B. Barnier, A.-M. Treguier, J. Molines, and G. Madec, 2007: Influence of numerical schemes on current-topography interactions in 1/4° global ocean simulations. *Ocean Sci.*, **3 (4)**, 509–524.

689

690

691

Penduff, T., and Coauthors, 2014: Ensembles of eddying ocean simulations for climate. CLIVAR Exchanges, Special Issue on High Resolution Ocean Climate Modelling, **19 (65)**, 26–29.

692

693

Pierini, S., 2014: Kuroshio extension bimodality and the north pacific oscillation: A case of intrinsic variability paced by external forcing. *J. Climate*, **27 (1)**, 448–454, doi:10.1175/JCLI-D-13-00306.1.

694

695

696

Schiermeier, Q., 2013: Oceans under surveillance. *Nature*, **497**, 167–168, doi:10.1038/497167a.

697

Send, U., T. Kanzow, W. Zenk, and M. Rhein, 2002: Monitoring the atlantic meridional overturning circulation at 16°N. CLIVAR Exchanges, Special issue on CLIVAR Atlantic, **7 (25)**, 31–33.

698

699

Sérazin, G., T. Penduff, S. Grégorio, B. Barnier, J.-M. Molines, and L. Terray, 2014: Intrinsic variability of sea-level from global 1/12° ocean simulations: spatio-temporal scales. *J. Climate*, doi:http://dx.doi.org/10.1175/JCLI-D-14-00554.1.

700

701
702

Thomas, M. D., and X. Zhai, 2013: Eddy-induced variability of the meridional overturning circulation in a model of the north atlantic. *Geophys. Res. Lett.*, **40** (11), 2742–2747, doi: 10.1002/grl.50532.

703

704

705

Torrence, C., and G. P. Compo, 1998: A practical guide to wavelet analysis. *Bull. Amer. Meteor. Soc.*, **79** (1), 61–78, doi:10.1175/1520-0477(1998)079<0061:APGTWA>2.0.CO;2.

706

707

Treguier, A.-M., and Coauthors, 2006: Internal and forced variability along a section between greenland and portugal in the clipper atlantic model. *Ocean Dynam.*, **56** (5-6), 568–580, doi: 10.1007/s10236-006-0069-y.

708

709

710

Treguier, A. M., and Coauthors, 2014: Meridional transport of salt in the global ocean from an eddy-resolving model. *Ocean Sci.*, **10** (2), 243–255, doi:10.5194/os-10-243-2014.

711

712

Zhang, R., 2010: Latitudinal dependence of atlantic meridional overturning circulation (amoc) variations. *Geophys. Res. Lett.*, **37** (16), L16703, doi:10.1029/2010GL044474.

713
<table><caption>TABLE 1. Model setups and simulations</caption><tr><td>Horizontal Resolution<sup>a</sup></td><td>2°</td><td>1/4°</td><td>1/12°</td></tr><tr><td>Meridional grid refinement</td><td>1/3° (equator)</td><td colspan="2">None</td></tr><tr><td>Nemo version</td><td colspan="2">2.3</td><td>3.4</td></tr><tr><td>Horizontal viscosity operator</td><td>Laplacian</td><td colspan="2">Biharmonic</td></tr><tr><td>Horizontal viscosity (equator)<sup>b</sup></td><td>Specific<sup>c</sup></td><td>1.5 × 10<sup>11</sup> m<sup>4</sup> s<sup>-1</sup></td><td>1.25 × 10<sup>10</sup> m<sup>4</sup> s<sup>-1</sup></td></tr><tr><td>Isopycnal diffusion (equator)<sup>d</sup></td><td>Specific<sup>e</sup></td><td>300 mm s<sup>-1</sup></td><td>100 mm s<sup>-1</sup></td></tr><tr><td>Lateral boundary conditions</td><td>No-slip</td><td colspan="2">Free-slip</td></tr><tr><td>Original forcing dataset</td><td>DFS4<sup>f</sup></td><td colspan="2">DFS4.4</td></tr><tr><td>SSS restoring timescale (days/10 m)</td><td>60</td><td colspan="2">60<sup>g</sup></td></tr><tr><td>AABW restoring</td><td>None</td><td colspan="2">2-year timescale<sup>h</sup></td></tr><tr><td>F runs: simulation name</td><td>ORCA246-G83b</td><td>ORCA025-B83</td><td>ORCA12-MJM88</td></tr><tr><td>F runs: integration period</td><td colspan="3">01/01/1958 - 12/31/2007 (50 years)</td></tr><tr><td>S runs: simulation name</td><td>ORCA246-MJM01</td><td>ORCA025-MJM01</td><td>ORCA12-GJM02</td></tr><tr><td>S runs: integration length</td><td>327 years</td><td>327 years</td><td>85 years</td></tr></table>

<footnote><sup>a</sup>At the equator.</footnote>

<footnote><sup>b</sup>Scales as the square of gridstep.</footnote>

<footnote><sup>c</sup>Cravatte et al. (2007)</footnote>

<footnote><sup>d</sup>Laplacian; scales as the gridstep.</footnote>

<footnote><sup>e</sup>Cravatte et al. (2007)</footnote>

<footnote><sup>f</sup>Brodeau et al. (2010)</footnote>

<footnote><sup>g</sup>Model SSS is smoothed prior to subtraction from Levitus SSS; restoring tends to zero toward continents over a 150km typical scale; corrective freshwater flux limited to 4mm/day.</footnote>

<footnote><sup>h</sup>Within the AABW density range, below 1000m and south of 30°S. See Dufour et al. (2012)'s appendix.</footnote>
<figure><img src="image_68.png" /><figcaption>FIG. 1. Time-averaged $AMOC_{\sigma_2}$ streamfunctions at 1/4° between 1976 and 2007 in the F run (left) and between years 19-50 in the S run (middle); zero contours are shown in white. The right panel shows the time-averaged transport of the $AMOC_{\sigma_2}$ upper cell northward limb in both runs (i.e. maximum of AMOC in density coordinates, refer to Section 2b); a 9-point boxcar filter has been used to smooth the data.</figcaption></figure>
<figure><img src="image_69.png" /><figcaption>FIG. 2. Example of $AMOC_{\sigma_2}$ raw timeseries (Sv, grey/black line) in the 1/4° S run over 327 years at the RAPID latitude (26.5°N). Red and blue lines represent the non-linear trends fitted to, then removed from, the 32-year and 309-year $AMOC_{\sigma_2}$ raw timeseries, respectively. The first 18 years, shown in grey, are discarded from all simulations.</figcaption></figure>
<figure><img src="image_70.png" /><figcaption>FIG. 3. $AMOC_{\sigma_2}$ interannual std (Sv) computed from monthly (plain line) and 5-daily (dashed line) model fields over the last 27 years of the 1/4° S run. A 9-point moving average has been used to smooth the data.</figcaption></figure>
<figure><img src="image_71.png" /><figcaption>FIG. 4. Total (a: years 1976-2007) and intrinsic (b: years 19-50) $AMOC_{\sigma_2}$ interannual std's at 1/4°, and ratio R of their variances (c). The white contours (interval = 2.5 Sv) represent the time-averaged $AMOC_{\sigma_2}$ streamfunctions from Fig. 1; zero contours are shown in black.</figcaption></figure>
<figure><img src="image_72.png" /><figcaption>FIG. 5. Total (a) and intrinsic (b) interannual standard deviations of the $AMOC_{\sigma_2}$ upper-cell northward limb transports at 1/4°, and ratio R of their variances (c). The 9 successive 32-year $s_I$ estimates of intrinsic variability are shown in grey in (b). The plain and dashed black lines in (b) and (c) represent the averages and envelopes of the 9 $s_I$ estimates presented in Section 2f. (d) Correlation between intrinsic interannual $AMOC_{\sigma_2}$ and MHT variabilities at 1/4° over the 9 successive 32-year periods. Plain and dashed lines represent the mean and std of the correlation coefficient over the 9 periods, respectively. Thick plain lines show where the correlations remain significative throughout the 9 periods. (e) same as (c) but for the interannual MHT at 1/4° (envelope of M defined in Section 2f). A 9-point moving average has been used to smooth these data.</figcaption></figure>
<figure><img src="image_73.png" /><figcaption>FIG. 6. Intrinsic-to-Total LF variance ratio R (numbers in %) from $AMOC_{\sigma_2}$ timeseries computed across four sections from our 1/12° simulations. The shaded background shows topography.</figcaption></figure>

<footnote>787</footnote>

<footnote>788</footnote>
<figure><img src="image_74.png" /><figcaption>FIG. 7. Same as Fig. 5 but for $AMOC_z$ at 1/4°, 2° and 1/12°. Two intrinsic std's and $R$ (b and c) at 1/12° resolution are computed from the two 32-year segments available in the S run. 9 and 25-point moving averages were used to smooth the 1/4° and 1/12° results, respectively.</figcaption></figure>
<figure><img src="image_75.png" /><figcaption>FIG. 8. Intrinsic AMOC$_{σ_2}$ anomaly timeseries and Wavelet power spectra from the last 309 years of the 1/4° S run at 29.9°S (a), 26.6°N (b), 54.3°N (c) and 64°N (d). Shown are unfiltered (blue) and band-passed filtered yearly timeseries in the 15-30 (red) and the 50-90 (black) year bands, for comparison with Fig. 10. Wavelet spectra are normalized by the variance of timeseries; white contours indicate 95% confidence levels.</figcaption></figure>
<figure><img src="image_76.png" /><figcaption>FIG. 9. Left panels: frequency-latitude Fourier power spectra of the intrinsic AMOC$_{\sigma_2}$ variability during the nine successive 32-year periods of the 1/4° S run. Spectra are normalized by the AMOC$_{\sigma_2}$ interannual variance at each latitude. Black and white contour lines indicate 95% and 90% confidence levels, respectively. Bottom right panel: probability that these intrinsic variability peaks are significant at 90% over the 9 periods (color, refer to Section 4b); regions of the frequency-latitude plane where total variability spectral peaks in the fully-forced 1/4° F run are significant at 90% (black contours, refer to Section 4c). Dashed white/black line shows the cutoff period of the LOESS detrending operator.</figcaption></figure>
<figure><img src="image_77.png" /><figcaption>FIG. 10. Time-latitude plots of intrinsic AMOC$_{\sigma_2}$ anomalies (Sv) during the 1/4° S run, bandpass-filtered in the 15-30 year band (top) and in the 50-90 year band (bottom); initial and final parts of the timeseries are ignored due to the side effects of the Lanczos filters. Slanted lines and associated numbers indicate ranges of meridional phasespeeds (in cm/s) which dominate within selected latitude ranges; pluses and minuses correspond to northward and southward propagation, respectively.</figcaption></figure>
# Ensembles of eddying ocean simulations for climate

Penduff, T.<fnref n="1" />, B. Barnier<fnref n="1" />, L. Terray², L. Bessières², G. Sérazin<fnref n="1" />, S. Gregorio<fnref n="1" />, J.-M. Brankart<fnref n="1" />, M.-P. Moine², J.-M. Molines<fnref n="1" />, and P. Brasseur<fnref n="1" />

1 CNRS-LGGE, Grenoble, France

2 CERFACS/CNRS, Toulouse, France

Corresponding author:
Thierry.Penduff@legi.grenoble-inp.fr

## 1. Context

The recent IPCC's Fifth Assessment Report (AR5) provides an extensive description of the past evolution, present state, and projected future of the climate system. The modelling contribution to the assessment is based in particular on a large multimodel ensemble of historical and scenario simulations performed within the Coupled Model Intercomparison Project version 5 (CMIP5). The fully-coupled numerical models used for these climate simulations were implemented with a typical resolution of 1-2° in both the atmosphere and ocean components. At these spatial resolutions, atmospheric models spontaneously generate synoptic weather chaotic variability and associated turbulent fluxes that strongly influence larger time and space scales. Simulated ocean dynamics, however, remain laminar in most of these climate simulations: the oceanic "weather" associated with mesoscale eddies is unresolved, and only few of its impacts on larger scales are parameterised, in a rather crude way. Eddying ocean/sea-ice models, implemented at ¼° (and eventually finer) resolutions, will be routinely used in future generations of coupled climate and Earth system models. Indeed, the explicit representation of oceanic eddies improves the consistency of simulated dynamics and yields more realistic results in simulations driven by realistic atmospheric forcing functions (e.g. Penduff et al., 2010). This paper focuses on the nature of interannual-to-decadal variability in eddying ocean/sea-ice simulations with prescribed atmospheric forcing, which greatly differs from that in the laminar regime.

## 2. Intrinsic low-frequency variability in the eddying ocean

Since the pioneering simulations by Holland and Haidvogel (1981) and Cox (1987), idealised studies have shown that the ocean circulation driven by constant atmospheric forcings can spontaneously generate low-frequency (1-10 year) intrinsic variability when resolution (i.e. mesoscale eddy activity) increases. This generic behaviour was proven robust to the addition of seasonal forcing (Sushama et al., 2007; Shimokawa et al., 2010). Two main nonlinear paradigms

have been proposed to explain this phenomenon. Spall (1996), Dewar (2003), Berloff et al. (2007), and Arbic et al. (2012) have shown that intrinsic low-frequency variability may be directly driven by mesoscale turbulence, i.e. through nonlinear interactions between the ambient potential vorticity (PV) field, eddies and the low-frequency variability, or through the temporal inverse cascade of kinetic energy driven by eddy-eddy interactions. Other authors, such as those reviewed in Dijkstra and Ghil (2005), used Dynamical System Theory approaches to show that mesoscale turbulence may "just" favour transitions of the oceanic circulation between multiple large-scale equilibria set by the mean PV field. It is very likely that these two paradigms are mutually compatible and provide complementary views on this phenomenon.

The contribution of low-frequency intrinsic variability in the real ocean cannot be estimated from observations alone since it is fully entangled with atmospherically-forced as well as coupled variability. This contribution is generally estimated from pairs of multi-decadal Ocean General Circulation Model (OGCM) simulations: one "fully-forced" hindcast, driven by the full range of atmospheric timescales over several decades, designed to mimic oceanic observations; and one climatological simulation, driven by a perpetual annual atmospheric cycle, designed to isolate the intrinsic low-frequency variability. Following this approach and Taguchi et al. (2007)'s regional study, Penduff et al. (2011) showed at global scale that when mesoscale eddies are represented (at ¼° resolution), a substantial fraction of the interannual variance of Sea-Level Anomaly (SLA) found in a fully-forced hindcast remains present in the climatological simulation (without interannual forcing). Consistently with idealised studies, most of this intrinsic interannual SLA variance was found in mesoscale-active mid-latitude regions, i.e. in western boundary currents systems, and in the Antarctic Circumpolar Current (ACC).

In the framework of the CHAOCEAN project<fnref n="1" />, these analyses of intrinsic low-frequency SLA variability are being extended to higher resolution simulations (1/12°) and to other climate-relevant variables (sea-surface temperature SST. Atlantic Meridional Overturning Circulation AMOC), with a focus on large spatial scales. Figure 1a shows that at scales larger than about 1000km, a substantial fraction of the SST interannual variance is spontaneously generated by the ocean in the very same regions where air-sea heat fluxes are the largest on Earth. This fraction is about 25-40% in the Gulf Stream and Kuroshio, and reaches 65-80% in the Agulhas Retroreflection region. South of Africa, the intrinsic, large-scale interannual SST standard deviation reaches about 0.5°C. Figure 1b presents the same fraction for interannual AMOC variance: in the South Atlantic up to 40-60% of the fully-forced AMOC variability is actually intrinsic in our eddying simulations; interestingly, the 1/4° and 1/12° results do not differ much, suggesting that so-called "eddying-permitting" resolution captures the essential dynamics at work, at least in mid-latitudes. OGCM simulations are thus consistent with idealised studies: the nature of interannual variability in the eddying regime (1/4° or 1/12°) is different from that in laminar ocean models used in CMIP exercises, where the impact of intrinsic processes on SLA, SST or AMOC variances are barely noticeable (e.g. red line in Figure 1b).

<footnote n="1">1. Four-year research project funded by the Ocean Surface Topography Science Team: http://alturl.com/x8qs8</footnote>
<figure><img src="image_78.png" /></figure>

Figure 1: Panels (a) and (b) show ratios (R in %) between intrinsic and total interannual variances, deduced from seasonally- and fully-forced global ocean/sea-ice Drakkar simulations, respectively. (a) Ratio R is shown for large-scale sea-surface temperature (SST) from global 1/12° simulations: spatial scales smaller than 12° have been filtered out before computing variances. (b) Ratio R for Atlantic Meridional Overturning Circulation (AMOC), shown as an ensemble mean and spread for ¼° simulations (black lines), and as two independent estimates for 1/12° simulations (magenta lines); 2° simulations (red) are almost devoid of intrinsic variability. (c) Normalised wavelet power spectrum of the AMOC intrinsic variability timeseries at 26°N throughout the 327-year seasonally-forced, global ¼° Drakkar OGCM simulation. To appear in Sérazin et al. and Grégorio et al.

### 3. Stochastic character of the low-frequency variability

The idealised studies mentioned above also highlight the chaotic character of the low-frequency, pure intrinsic variability (i.e. simulated with seasonal or constant forcing). In other words, sensitivity to initial conditions and tendency for intermittent, random behaviour in nonlinear regimes are not restricted to mesoscale motions, but also concern low-frequency modulations of the circulation. Figure 1c shows that this property is also found in a realistic context: this wavelet power spectrum, computed from a 327-year 1/4° OGCM seasonally-forced simulation, reveals the intermittent and broadband character of the interannual AMOC pure intrinsic variability at 26°N.

An important question then comes to mind: does the pure (seasonally-forced) intrinsic low-frequency variability remain chaotic when a low-frequency atmospheric forcing is prescribed at the ocean surface? Recent investigations provide preliminary, although contrasting, answers. On the one hand, Pierini (2014) compared idealised shallow-water simulations of the North Pacific driven by constant and interannually-varying winds, and concluded that the low-frequency variability of the Kuroshio is chaotic in the former

case but may become deterministic in the latter. On the other hand, the structure and magnitude of the interannual AMOC intrinsic variance generated by our global 1/4° OGCM under seasonal forcing does not differ much (in structure and intensity) from its counterpart estimated by Hirschi et al. (2013) under interannual forcing with the same model. In other words, interannual forcing might excite pure intrinsic variability modes (as also shown by O'Kane et al., 2013), but the interannually-forced low-frequency oceanic variability is very likely to remain stochastic in the eddying regime.

There are thus various indications that in forced eddying OGCM simulations, the large-scale, low-frequency, mid-latitude ocean variability is less deterministic (i.e. more sensitive to initial conditions) than in the laminar regime, and that eddy-related non-linearities give a stochastic flavour to low-frequency variability. This concerns in particular the large-scale, low-frequency SST and AMOC variability that can have direct and indirect impacts on the atmosphere and climate. This stochastic character advocates for a probabilistic study of the low-frequency oceanic variability in the eddying regime.
4. A global ensemble of multi-decadal oceanic hindcasts

The main goal of the OCCIPUT project<fnref n="2" /> is to provide the first probabilistic description of the global oceanic evolution over the last decades, i.e. based on probability density functions (PDFs). Our plan is to force 50 eddying (1/4°) global ocean/ sea-ice simulations with realistic atmospheric forcings over the period 1958-2013, given perturbed initial conditions. Our integration strategy, including the initialization, spin-up and ensemble generation, is being currently designed using regional reduced-sized ensembles, in close collaboration with the DRAKKAR Group<fnref n="3" />. During the beginning of the integration, we plan to activate within each member a stochastic parameterisation (of e.g. subgrid-scale density gradients, Brankart, 2013) to generate an ensemble spread. Ongoing tests confirm a priori expectations from instability theories: the ensemble spread, i.e. intrinsic variability, emerges at the scale of mesoscale eddies and grows fastest in frontal regions (Figures 2a, 2b, 2c). As expected from inverse cascade studies (Arbic et al., 2014), the space and time scales of intrinsic variability progressively increase: Figure 2d exhibits an inter-member difference that has reached Rossby wave scales (i.e. 1000km and 1-year) after 3 years away from the main fronts, i.e. at 26°N. These tests also show that switching off stochastic perturbations after a few months does not hamper the subsequent growth of the ensemble spread: mesoscale nonlinearities

ultimately control the spin-up of intrinsic variability, as well as (presumably) its cascade toward larger scales and its magnitude. In other words, we plan to stochastically perturb the OCCIPUT ensemble members for 1 year or so, then let the spread grow, reach interannual and regional scales and its non-linearly saturated amplitude. This strategy should provide us with a data set of a few decades over which we plan to describe and study the spatial structure and variability of climate-relevant ensemble statistics in key regions: e.g. ensemble PDF of oceanic heat content, ensemble average (deterministic part) of Kuroshio transport, ensemble spread (stochastic part) of global SST, ensemble covariances (turbulent fluxes), etc.

5. Implementation, model outputs, expected outcomes

The OCCIPUT members are currently adapting the ORCA025 global 1/4° ocean/sea-ice model configuration (Barnier et al., 2006) for ensemble integrations and implementing the system on the CURIE supercomputer<fnref n="4" />. The 50 members will be integrated simultaneously as one executable over a few thousands of processors, so that specific ensemble (inter-member) statistics can be computed online, possibly at every timestep. Such statistics may include e.g. timestep estimates of (exact) turbulent heat fluxes across key sections, deciles of global SST distribution at hourly resolution, etc. The full 3-dimensional 50-member state

<figure><img src="image_79.png" /></figure>

Figure 2: Subsurface temperature (T at 855m depth, °C) in a 3-member, fully-forced North Atlantic NEMO ¼° ensemble simulation with stochastic parameterisation (Brankart et al 2013); only scales larger than 4° are considered. Left panels show 3-year temperature timeseries from the 3 members at locations 1 (a) and 2 (b) indicated on the right. Right panels show the distribution of the inter-member T difference after 3 years as a function of longitude and latitude (c), and as a function of longitude and time along 26°N (d).

<footnote n="2">2. Four-year research project funded by ANR: http://alturl.com/ivfkg</footnote>

<footnote n="3">3. who are developing the ORCA025, global ¼° ocean/sea-ice configuration: http://www.drakkar-ocean.eu/</footnote>

<footnote n="4">4. http://www-hpc.cea.fr/en/complexe/tgcc-curie.htm. An estimated amount of 15 million core hours shall be requested on this machine through a PRACE proposal in fall 2014; we hope to produce the ensemble integration in 2015-2016.</footnote>
vector will be archived as successive monthly averages, and
selected 2-dimensional fields every 5 days; these data will
allow offline computations of additional ensemble statistics
(PDFs, variances, covariances) at 5-daily, monthly, or
interannual timescales of e.g. sea-level anomalies, bottom
pressure, components of the AMOC at 26°N or other
observed latitudes, etc.

The system will also generate ensemble synthetic observations
during the integration, i.e. 50 simultaneous collocated counter-
parts of the ENACT-ENSEMBLE temperature/salinity profile
database (1958-present) or of along-track altimeter timeseries
(e.g. Topex/POSEIDON, Jason, etc). Such data will allow a
thorough, objective observation-based assessment of the
probabilistic simulation using metrics defined in the context
of weather ensemble prediction (e.g. rank histograms, Brier
scores; see Toth et al., 2003). Such synthetic data may also
provide relevant information about the representativeness of
actual observations.

6. Conclusion

Academic models have illustrated the chaotic behaviour of
the ocean circulation at high Reynolds number, not only in
terms of mesoscale turbulence but also in larger (double-
gyre or ACC-like) current systems, up to decadal timescales.
Unlike laminar ocean models used in most current climate
projections, eddying OGCMs also spontaneously generate a
substantial interannual-to-decadal intrinsic variability under
repeated seasonal forcing, with a stochastic character, a
marked signature on AMOC, and on SST in regions where
air-sea fluxes have the largest amplitude. Whether and how
this ocean-driven low-frequency intrinsic variability may
ultimately impact the climate predictability at various time
and spatial scales is an important but unsettled question.

Before answering this question, it is useful to better assess
the stochastic character of the low-frequency ocean
variability, with a focus on climate-relevant indexes. The
OCCIPUT project aims at performing in 2015 a 50-member
ensemble of 1/4° global ocean/sea-ice hindcasts driven
by the same 1958-present atmospheric forcing. The initial
ensemble spread is expected to grow, cascade toward long
space and time scales, and saturate in amplitude. We expect
this eddying ensemble to provide a probabilistic description
of the ocean state and evolution over the last decades, and a
measure of the actual constraint exerted by the atmosphere
on the past interannual-to-decadal ocean variability. New
results are also expected from the computation of ensemble
instead of temporal statistics (variances, covariances)
for the estimation of time-dependant turbulent fluxes,
and from the use of ensemble synthetic observations
to probabilistically assess ocean simulations or the
representativeness of observational datasets.

Acknowledgments

The OCCIPUT project is supported by ANR (ANR-13-
BS06-0007-01), and CHAOCEAN is supported by the CNES/
NASA Ocean Surface Topography Science Team (OST/ST).
The authors also acknowledge support from CNES. The work
presented here used HPC resources from GENCI (TGCC/
CINES/IDRIS) on grants x2013010727 and x2014010727.

References

Arbic, B.K., R.B. Scott, G.R. Flierl, A.J. Morten, J.G. Richman, and J.F. Shriver, 2012: Nonlinear cascades of surface oceanic geostrophic kinetic energy in the frequency domain. J. Phys. Oceanogr., 42, doi:10.1175/JPO-D-11-0151.1, 1577-1600.

Arbic, B.K., M. Müller, J.G. Richman, J.F. Shriver, A.J. Morten, and

R.B. Scott, G. Serazin, and T. Penduff, 2014: Geostrophic turbulence in the frequency-wavenumber domain: Eddy-driven low-frequency variability, J. Phys. Oceanogr, in press.

Barnier, B., G. Madec, T. Penduff, J.-M. Molines, A.-M. Treguier, J. Le Sommer, A. Beckmann, A. Biastoch, C. Böning, J. Dengg, C. Derval, E. Durand, S. Gulev, E. Remy, C. Talandier, S. Theetten, M. Maltrud, J. McClean, B. De Cuevas, 2006: Impact of partial steps and momentum advection schemes in a global circulation model at eddy permitting resolution, Ocean Dynamics, 56(5-6), 543-567.

Berloff, P., A. Mc C. Hogg, W. Dewar, 2007: The Turbulent Oscillator: A Mechanism of Low-Frequency Variability of the Wind-Driven Ocean Gyres. J. Phys. Oceanogr., 37, 2363–2386.

Brankart J.-M., 2013: Impact of uncertainties in the horizontal density gradient upon low resolution global ocean modelling. Ocean Modelling, 66, 64-76.

Cox, M.D., 1987: An eddy-resolving numerical model of the ventilated thermocline: time dependence. J. Phys. Oceanogr., 17, 1044–1056.

Dewar, W.K., 2003: Nonlinear mid-latitude ocean adjustment, J. Phys. Oceanogr., 33, 1057-1082.

Dijkstra, H.A. and M. Ghil, 2005: Low-frequency variability of the large-scale ocean circulation: A dynamical systems approach. Reviews of Geophysics, 43, doi: 10.1029/2002RG000122. issn: 8755-1209.

Grégorio, S., T. Penduff, G. Sérazin, J. Le Sommer, J.-M. Molines, B. Barnier, and J. Hirschi: Intrinsic variability of the Atlantic Meridional Overturning Circulation at interannual-to-multidecadal timescales. In preparation.

Hirschi, J.J.-M., A.T. Blaker, B. Sinha, A. Coward, B. de Cuevas, S. Alderson, and G. Madec, 2013: Chaotic variability of the meridional overturning circulation on subannual to interannual timescales. Ocean Science Discussions, 9(5), 3191-3238. 10.5194/osd-9-3191-2012

Holland, W.R., and D.B. Haidvogel, 1981: On the Vacillation of an
Unstable Baroclinic Wave Field in an Eddy-Resolving Model of the
Oceanic General Circulation. J. Phys. Oceanogr., 11, 557–568.

O'Kane, T.J., R.J. Matear, M.A. Chamberlain, J.S. Risbey, B.M. Sloyan, and I. Horenko, 2013: Decadal variability in an OGCM Southern Ocean: Intrinsic modes, forced modes and metastable states. Ocean Modelling, 69, pp 1-21.

Penduff, T., M. Juza, L. Brodeau, G.C. Smith, B. Barnier, J.-M. Molines, A.-M. Treguier and G. Madec, 2010: Impact of global ocean model resolution on sea-level variability with emphasis on interannual time scales. Ocean Sci., 6, 269–284.

Penduff, T., M. Juza, B. Barnier, J. Zika, W.K.Dewar, A.-M. Treguier, J.-M. Molines, and N. Audiffren, 2011: Sea-level expression of intrinsic and forced ocean variabilities at interannual time scales. J. Climate, 24, 5652–5670.

Pierini S., 2014: Kuroshio Extension bimodality and the North Pacific Oscillation: a case of intrinsic variability paced by external forcing. J. Climate, 27, 448-454

Sérazin, G., T. Penduff, L. Terray, S. Grégorio, B. Barnier, and J.-M. Molines : Spatial scales of the low-frequency intrinsic sea-level variability: a global model study. In preparation

Shimokawa, S., and T. Matsuura, 2010: Chaotic Behaviors in the Response of a Quasigeostrophic Oceanic Double Gyre to Seasonal External Forcing. J. Phys. Oceanogr., 40, 1458-1472.

Spall, M.A., 1996: Dynamics of the Gulf Stream/Deep Western Boundary Current Crossover. Part II: Low-frequency internal oscillations. J. Phys. Oceanogr., 26(10), 2169-2182.

Sushama, L., M. Ghil, and K. Ide, 2007: Spatio-temporal variability in a mid-latitude ocean basin subject to periodic wind forcing, Atmosphere-Ocean, 45(4), 227-250, DOI: 10.3137/ao.450404

Taguchi, B., S.P. Xie, N. Schneider, M. Nonaka, H. Sasaki, and Y. Sasai, 2007: Decadal Variability of the Kuroshio Extension: Observations and an Eddy-Resolving Model Hindcast. J. Climate, 20, 2357–2377.

Toth Z., O. Talagrand, G. Candille, and Y. Zhu. 2003: « Probability and ensemble forecasts », in Forecast Verification: A Practitioner's Guide in Atmospheric Science, Jolliffe I., Stephenson D.B. (eds), Wiley: UK. pp 137–163.