Olivier Arzel
email: olivier.arzel@univ-brest.fr

Thierry Huck

Contributions of Atmospheric Stochastic Forcing and Intrinsic Ocean 1 Modes to North Atlantic Ocean Interdecadal Variability 2

Atmospheric stochastic forcing associated with the North Atlantic Oscillation (NAO) and intrinsic ocean modes associated with the large-scale baroclinic instability of the North Atlantic Current (NAC) are recognized as two strong paradigms for the existence of the Atlantic Multidecadal Oscillation (AMO). The degree to which each of these factors contribute to the lowfrequency variability of the North Atlantic is the central question in this paper.

This issue is addressed here using an ocean general circulation model run under a wide range of background conditions extending from a super-critical regime where the oceanic variability spontaneously develops in the absence of any atmospheric noise forcing to a damped regime where the variability requires some noise to appear. The answer to the question is captured by a single dimensionless number Γ measuring the ratio between the oceanic and atmospheric contributions, as inferred from the buoyancy variance budget of the western subpolar region. Using this diagnostic, about two-third of the sea surface temperature (SST) variance in the damped regime is shown to originate from atmospheric stochastic forcing whereas heat content is dominated by internal ocean dynamics. Stochastic wind-stress forcing is shown to substantially increase the role played by damped ocean modes in the variability. The thermal structure of the variability is shown to differ fundamentally between the super-critical and damped regimes, with abrupt modifications around the transition between the two regimes. Ocean circulation changes are further shown to be unimportant for setting the pattern of SST variability in the damped regime but are fundamental for a preferred timescale to emerge.

Introduction

The Atlantic Multidecadal Oscillation (AMO) is a major mode of climate variability explaining nearly 40% of the spatially integrated annual mean sea surface temperature (SST) variance over the North Atlantic [START_REF] Delworth | Decadal to centennial variability of the Atlantic from observations and models[END_REF]). The AMO does not only modulate the climate of the surrounding continents on decadal to multidecadal timescales [START_REF] Zhang | Impact of Atlantic Multidecadal Oscillation on India/Sahel rainfall and Atlantic hurricanes[END_REF][START_REF] Knight | Climate impacts of the Atlantic multidecadal oscillation[END_REF]) but also directly impacts marine ecosystems [START_REF] Edwards | Marine ecosystem response to the Atlantic Multidecadal Oscillation[END_REF] and Arctic sea ice [START_REF] Mahajan | Impact of the Atlantic Meridional Overturning Circulation (AMOC) on Arctic surface air temperature and sea ice variability[END_REF][START_REF] Zhang | Mechanisms for low-frequency variability of summer Arctic sea ice extent[END_REF]. Several physical mechanisms have been put forth to explain the origin of these low-frequency variations, but the diversity of those mechanisms does not allow to provide a clear and robust picture as to which of the ocean or the atmosphere primarily drives the AMO and how it works. More specifically, modelling studies led to the emergence of (at least) two paradigms for the AMO: the first one is related to the integration of the atmospheric white noise by the ocean; the second one has dynamical origins and is related to intrinsic unstable interdecadal ocean modes. The two phenomena probably play a role in the low-frequency variability of the North Atlantic climate as suggested by a number of studies [START_REF] Delworth | Interdecadal variations of the thermohaline circulation in a coupled ocean-atmosphere model[END_REF][START_REF] Delworth | Observed and simulated multidecadal variability in the northern hemisphere[END_REF][START_REF] Dong | Mechanism of interdecadal thermohaline circulation variability in a coupled ocean-atmosphere GCM[END_REF][START_REF] Gastineau | North Atlantic Ocean internal decadal variability: role of the mean state and ocean-atmosphere coupling[END_REF]), but their respective contributions in establishing the pattern and amplitude, and even in determining the very existence, of the AMO remains elusive and model dependent.

The simplest paradigm to explain low-frequency climate variability originates from the seminal work of [START_REF] Hasselmann | Stochastic climate models. Part I: theory[END_REF] who showed that the integration of atmospheric white noise by the ocean along with its large heat capacity gives rise to a reddenned spectrum. This purely thermodynamic response has been invoked by [START_REF] Clement | The Atlantic Multidecadal Oscillation without a role for ocean circulation[END_REF] who questioned the role of ocean circulation changes in the AMO by comparing results from fully coupled models and atmospheric general circulation models coupled to slab-ocean models that do not permit circulation changes.

4

Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

The pattern of SST variability is remarkably similar between the two families of models, leading the authors to conclude that ocean circulation changes are not essential in determining both the pattern and existence of the AMO. Their analysis supports the null hypothesis that the ocean merely integrates the white noise atmospheric forcing of the North Atlantic Oscillation (NAO) to produce a red noise response. Similar conclusions were reached by [START_REF] Schneider | Weather noise forcing of surface climate variability[END_REF] who showed that the null hypothesis is appropriate over much of the World Ocean in the diagnosis of the SST variability in a coupled climate model. The lack of a distinct multidecadal spectral peak in models (at least in the multi-model mean) is in constrast with a number of observations including instrumental measurements [START_REF] Tung | Using data to attribute episodes of warming and cooling in instrumental records[END_REF], tree ring records [START_REF] Delworth | Observed and simulated multidecadal variability in the northern hemisphere[END_REF][START_REF] Gray | A tree-ring based reconsruction of the Atlantic Multidecadal Oscillation since 1567 AD[END_REF]), ice-core records [START_REF] Chylek | Ice-core data evidence for a prominent near 20 year time-scale of the Atlantic multidecadal oscillation[END_REF], and multi-proxy based reconstructions [START_REF] Knudsen | Tracking the Atlantic multidecadal oscillation through the last 8,000 years[END_REF], that show enhanced variability in the 20-80 years range in the Atlantic sector. [START_REF] Dommenget | Analysis of observed and simulated SST spectra in the midlatitudes[END_REF] compared the statistics of large-scale SST variability in the mid-latitudes of the Northern Hemisphere between different coupled models, slab ocean models and observations. In contrast to [START_REF] Clement | The Atlantic Multidecadal Oscillation without a role for ocean circulation[END_REF] these authors concluded that the SST variability in the midlatitudes is significantly different from a red noise response and that processes in the ocean are responsible for these differences. [START_REF] Saravanan | Advective ocean-atmosphere interaction: An analytical stochastic model with implications for decadal variability[END_REF] modified the Hasselmann's model to include steady mean oceanic advection and a spatially variable noise forcing. In contrast to [START_REF] Hasselmann | Stochastic climate models. Part I: theory[END_REF], a preferred timescale is selected by the circulation as long as advective effects dominate thermal damping effects associated with air-sea heat exchanges, leading to a phenomenon called spatial resonance.

The second paradigm relies on the large-scale baroclinic instability of the North Atlantic Current and subsequent westward propagation of unstable planetary waves leading to interdecadal (20-30 yr) oscillations of the Atlantic Meridional Overturning Circulation (AMOC) (Colin de Verdière [START_REF] De Verdière | Baroclinic instability: An oceanic wavemaker for interdecadal variability[END_REF][START_REF] Te Raa | Instability of the thermohaline circulation on interdecadal timescales[END_REF]. Instability occurs at high Peclet numbers through a 5 Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

Hopf bifurcation. The growth rates at bifurcation O(1) year -1 are on the order of effective damping of SST anomalies which explains why this mode could be damped in some coupled models while active in others. [START_REF] Arzel | The internal generation of the Atlantic Ocean interdecadal variability[END_REF] recently studied the bifurcation structure and pattern of this intrinsic mode in the realistic configuration of an ocean general circulation model under prescribed surface fluxes to show that the features previously identified in idealized contexts are robust in a more realistic setting (geometry and physics). In particular, the SST variance now peaks in the western subpolar gyre of the North Atlantic, a feature that is also clearly apparent in observations [START_REF] Deser | Sea surface temperature variability: patterns and mechanisms[END_REF]. The variability disappears for eddy-induced diffusivities O(500 -1000) m 2 s -1 (Huck and Vallis 2001;[START_REF] Arzel | The internal generation of the Atlantic Ocean interdecadal variability[END_REF]) that are in the range of those derived from observations [START_REF] Liu | Adjoint-based estimation of eddy-induced tracer mixing parameters in the global ocean[END_REF][START_REF] Abernathey | Global surface eddy diffusivities derived from satellite altimetry[END_REF] in the subpolar area of the North Atlantic, casting therefore some doubts on the relevance of this second paradigm. These critical values are also in the range of those usually employed in current climate models (see Table 1 in [START_REF] Kuhlbrodt | The influence of eddy parameterizations on the transport of the Antarctic Circumpolar Current in coupled climate models[END_REF] suggesting that stochastic forcing, presumably by the atmosphere, may be needed to sustain this mode in coupled models. Frankcombe et al. (2009) precisely focused on this point to show that atmospheric stochastic forcing leads to oceanic variability in the regime where the intrinsic ocean mode is damped. The effect is strong provided that the noise forcing has a spatial structure (e.g. NAO) and some temporal coherence. What fraction of the variability is driven by this internal mode of variability and the NAO forcing remains however to be determined, in particular in the regime where the internal ocean mode is damped. While some studies show a central role of the NAO forcing in the very existence of North Atlantic climate variability [START_REF] Delworth | Multidecadal thermohaline circulation variability excited by atmospheric surface flux forcing[END_REF][START_REF] Eden | A damped decadal oscillation in the North Atlantic climate System[END_REF][START_REF] Chen | Mechanisms of internally generated decadal-tomultidecadal variability of SST in the Atlantic Ocean in a coupled GCM[END_REF], others point instead to internal ocean dynamics with the noise forcing acting as an amplifier of the variability obtained under climatological surface fluxes [START_REF] Zhu | Interdecadal variability of the meridional overturning circulation as an ocean internal mode[END_REF][START_REF] Gastineau | North Atlantic Ocean internal decadal variability: role of the mean state and ocean-atmosphere coupling[END_REF].

6

Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

The aim of this paper is to investigate in a systematic manner the role played by intrinsic ocean modes in the variability of the Atlantic circulation of an ocean general circulation model subject to atmospheric stochastic forcing. A dynamical system approach is used whereby the characteristics and origins of the variability are systematically assessed against background oceanic conditions.

Different background states are achieved by using different magnitudes of eddy-induced diffusivity, one of the most critical parameter at the relatively low resolution used here. This approach allows us to contrast different oscillatory regimes that have been previously identified in the literature, namely that driven by deterministic dynamics (self-sustained ocean mode) and that excited by atmospheric weather noise (damped ocean mode). Special emphasis will be placed upon the nature and origins of SST variability which is the relevant field in the context of air-sea interactions. The paper seeks to address the following questions: What are the respective contributions of the NAO-like atmospheric stochastic forcing and large-scale baroclinic instability mechanism to the simulated North Atlantic SST and circulation variability? A central aspect is to determine how these contributions depend on background oceanic conditions. Does the spatial pattern of the variability, in particular in terms of horizontal propagation and vertical structure of temperature anomalies, obtained in the regime where the internal ocean mode is active, differ from that obtained in the damped regime? Are oceanic circulation changes fundamental to explain the properties (pattern, amplitude and dominant timescale) of the low-frequency variability?

This paper is organized as follows. Section 2 describes the model and experimental design. The main characteristics of the variability along with its sensitivity to background oceanic conditions are presented in section 3. In section 4, the mechanisms responsible for the maintenance of the variability against all sources of thermal damping are identified and the associated energy sources are quantified. The role of ocean circulation changes is then investigated in section 5. Key findings are summarized and discussed in section 6. Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

Model and experiments

a. The ocean model

The model used for this study is the MITgcm [START_REF] Marshall | A finite-volume, incompressible Navier Stokes model for studies of the ocean on parallel computers[END_REF]) in a configuration identical with that used by [START_REF] Arzel | The internal generation of the Atlantic Ocean interdecadal variability[END_REF]. The only difference lies in the surface heat and momentum fluxes which now include a stochastic part. The ocean model is run at 1 • horizontal resolution and extends from 80 • S to 80 • N. There are 44 levels in the vertical with grid spacing increasing from 10 m at the surface to 250 m at the bottom. Static instability is removed by enhanced mixing (100 m 2 s -1 ). The vertical diffusivity increases downward following [START_REF] Bryan | A water mass model of the world ocean[END_REF] with upper and bottom values of 0.5×10 -4 m 2 s -1 and 1.3×10 -4 m 2 s -1 , respectively. These values are in line with those inferred from large-scale inversion experiments [START_REF] Lumpkin | Global ocean meridional overturning[END_REF], direct measurements [START_REF] Waterhouse | Global patterns of diapycnal mixing from measurements of the turbulent dissipation rate[END_REF]) and more recent robust diagnostic calculations (Arzel and Colin de Verdière 2016). We do not use any mixed layer turbulence parameterization. We use a spatially uniform horizontal Laplacian viscosity ν h of 5×10 4 m 2 s -1 . The Gent-McWilliams (GM, Gent and McWilliams 1990) parameterization of mesoscale eddies is implemented along with the rotated eddy diffusion tensor for isopycnal mixing [START_REF] Redi | Oceanic isopycnal mixing by coordinate rotation[END_REF]. A parameter sensitivity analysis in terms of the eddy-induced turbulent diffusivity K is carried out (Table 1). The isopycnal mixing coefficient is set to 1000 m 2 s -1 in all experiments. The equation of state is that proposed by [START_REF] Jackett | Minimal adjustment of hydrographic profiles to achieve static stability[END_REF], which computes the in-situ density from potential temperature, practical salinity and Boussinesq hydrostatic pressure. Ocean bathymetry is taken from the historical ETOPO1 dataset [START_REF] Amante | ETOPO1: 1 Arc-Minute Global Relief Model: Procedures, Data Sources and Analysis[END_REF]

b. Experimental design

We use flux boundary conditions at the surface for both temperature and salinity, similar to [START_REF] Arzel | The internal generation of the Atlantic Ocean interdecadal variability[END_REF]. The absence of feedback between sea surface salinity (SSS) and freshwater flux justifies the use of a flux formulation for salinity. The use of a flux formulation for temperature resides on the well-established result that on timescales much longer than the atmospheric response time, typically 10 days, atmospheric thermal damping of SST anomalies is relatively weak. [START_REF] Vallis | Stochastic Physics and Climate Modelling, chap. Mechanisms of Climate Variability from Years to Decades[END_REF] estimates this damping timescale to be 4.4 years, which is on the same order as a typical e-folding time of perturbations found in models forced by prescribed surface fluxes (Huck et al. 2001;[START_REF] Arzel | The internal generation of the Atlantic Ocean interdecadal variability[END_REF]. [START_REF] Arzel | The internal generation of the Atlantic Ocean interdecadal variability[END_REF] showed that the addition of a surface restoring flux with a damping timescale α -1 of one year has little influence on the characteristics of the interdecadal variability obtained under deterministic conditions (zero stochastic forcing).

The main effect of thermal damping is to completely damp out the variability near bifurcation, consistent with baroclinic growth rates µ ∼ α, providing a zero net growth of perturbations there.

Away from bifurcation (i.e. towards higher Peclet numbers) µ α in agreement with the stronger circulation leading to a relatively minor impact of surface damping on the variability. On the basis of these results we have chosen to use prescribed surface heat and freshwater fluxes in all numerical experiments. Following [START_REF] Bryan | Parameter sensitivity of primitive equation ocean general circulation models[END_REF], surface buoyancy fluxes are diagnosed from a model integration under restoring boundary conditions, rather than prescribed from observations. The procedure to compute those fluxes is detailed in [START_REF] Arzel | The internal generation of the Atlantic Ocean interdecadal variability[END_REF] but is given here for completeness. For each value of K, the model is first brought to equilibrium through relaxation of the SST and SSS fields toward the World Ocean Atlas climatology (Locarnini et al. 2010;[START_REF] Antonov | World Ocean Atlas[END_REF]

Q(x, y,t) = Q T (x, y,t) + Q NAO (x, y)ζ (t), τ(x, y,t) = τ obs (x, y,t) + τ NAO (x, y)ζ (t),
where Q and τ are the total surface heat and momentum fluxes. The patterns Q NAO (x, y) and τ NAO (x, y) have been obtained by regressing the corresponding annual mean anomalies from [START_REF] Large | The global climatology of an interannually varying air-sea flux data set[END_REF] onto the winter mean (DJFM) NAO index [START_REF] Hurrel | Trends in the North Atlantic Oscillation -Regional temperatures and precipitation[END_REF] multiplied by one standard deviation of the NAO index (Fig. 1). The stochastic forcing is only applied to the North Atlantic. The random discrete timeseries ζ (t) has been built from a first order auto-regressive process with a decorrelation timescale of 10 days. This timescale corresponds to estimates inferred by [START_REF] Feldstein | The timescale, power spectra, and climate noise properties of teleconnection patterns[END_REF] using daily means from the NCEP-NCAR reanalysis. The noise forcing has a sampling frequency (8100 s) corresponding to the time step of the model, yielding a timeseries of 5,760,000 points for a 1,500 years integration under flux boundary conditions (experiments named FTFS, prescribed flux for temperature T and salinity S). The variance of ζ (t) has been adjusted so that the timeseries built from the monthly means of ζ (t) has a variance equal to 1, similar to [START_REF] Herbaut | Response of a simplified oceanic general circulation model to idealised NAO-like stochastic forcing[END_REF]. The sensitivity of the model to the eddy-induced diffusiv- m 2 s -1 . In all experiments, the eddy-diffusivity K is held constant in space. Those values of K span the observed range of eddy diffusivities but do not attempt to capture the strong spatial variations [START_REF] Abernathey | Global surface eddy diffusivities derived from satellite altimetry[END_REF]. All experiments with both stochastic heat and momentum fluxes are repeated with a stochastic heat flux component only (Table 1). The aim of those experiments is to determine the additive effect of stochastic wind-stress forcing on both the characteristics and energy sources of the variability. Additional experiments are designed to determine the precise role of circulation changes in North Atlantic SST variability (details given in section 5).

Results

a. AMOC variability

In all stochastically-forced experiments, a pronounced decadal to multidecadal variability of the Atlantic ocean circulation develops. This can be seen in the timeseries and power spectrum of the AMOC index in Fig. 2 for four different values of K (500, 800, 1200 and 1600 m 2 s -1 ). The AMOC index used here is defined as the maximum value of the annual mean meridional overturning streamfunction below 1000 m and north of 30 • N in the North Atlantic. A clear distinction can be made between the AMOC variability obtained with K < 600 m 2 s -1 from that obtained with K > 600 m 2 s -1 . As shown earlier by [START_REF] Arzel | The internal generation of the Atlantic Ocean interdecadal variability[END_REF], the critical value K = K c = 600 m 2 s -1 corresponds to the existence of a threshold separating a super-critical regime (K < K c ) where the variability spontaneously emerges under deterministic conditions from a damped regime (K > K c ) where the oceanic variability does not emerge in the absence of any noise forcing. In the supercritical regime (K < K c ) the oscillations in the AMOC are large and appear quite regular, thereby producing a distinct spectral peak. In the damped regime (K > K c ), the oscillations have much weaker amplitude and appear less regular with a much broader spectrum. Because density anomalies are dominated by temperature changes rather than salinity changes (not shown), we restrict the description that follows in terms of temperature only. Figure 3 shows the standard deviations of the annual mean SST field as computed from 1000 years of model output from the FTFS experiments for the same four values of K as above. In all cases SST changes are maximum in the western subpolar gyre. Similarly to the AMOC variability (Fig. 2), a clear distinction can however be made between the patterns obtained with K < K c from those obtained with K > K c . For K < K c , SST changes are large in the mid-latitudes, typically between 40 • N and 60 • N, and much weaker elsewhere. There is a significant drop in the amplitude of SST changes around K = K c , in particular in the western subpolar region where the internal ocean mode has its largest fingerprint [START_REF] Arzel | The internal generation of the Atlantic Ocean interdecadal variability[END_REF]. The amplitude of SST changes in the subtropics is in contrast nearly insensitive to K, suggesting that the variability is mostly constrained by the NAO forcing there. As a result, SST changes appear more uniform across the basin in the damped regime with similar amplitudes between the subpolar and subtropical regions.

c. Propagation of SST anomalies

The time evolution of temperature anomalies in relation with the AMOC shows some striking differences between the damped and super-critical regimes. Figure 4 shows composites of SST anomalies as obtained for the same four values of K used above and for periods when the AMOC is maximum (AMOC index larger than the mean plus one standard deviation), when the AMOC anomaly is close to zero and decreasing, when the AMOC is minimum (AMOC index lower than the mean minus one standard deviation), and when the AMOC anomaly is close to zero and increasing. Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

For K < K c the large-scale propagation signals obtained under stochastic surface boundary conditions are almost undistinguishable from those obtained in the deterministic case. With the same ocean model and configuration, [START_REF] Arzel | The internal generation of the Atlantic Ocean interdecadal variability[END_REF] shows that the interdecadal variability in the deterministic case is driven by a large-scale baroclinic instability of the North Atlantic Current (NAC). The strong resemblance between the deterministic and stochastic cases suggests that the same mechanism operate when a noise forcing is present. The effect of the noise forcing on the characteristics of the variability thus appears very limited in the super-critical regime, supporting the idea that the variability mostly originates from internal ocean processes rather than from the NAO forcing in this regime. The time evolution of SST anomalies can be described as follows.

When the AMOC is at its maximum, a prominent SST dipole centered around the mean path of the NAC is present, with a warm anomaly in the east and a cold anomaly in the west, resulting in a stronger than usual NAC via thermal wind balance. As the AMOC decreases, the cold anomaly propagates westward until reaching the western boundary while the warm anomaly splits into two distinct parts on either side of the NAC, one propagating southeastward and the other westward.

By the meantime, a cold anomaly has emerged along the NAC consistent with a reduced poleward heat transport during that period. As time proceeds, this cold anomaly grows up while the western warm anomaly barely evolves. At the AMOC minimum, the situation is exactly opposed to that obtained at the AMOC maximum with an eastern cold anomaly and a western warm anomaly. The subsequent evolution of SST anomalies is similar to that obtained during the decaying phase of the AMOC, but with opposite signs. Central to the existence of the oscillation is the reversal in the sign of the anomalous zonal pressure (temperature) gradient across the NAC. The overall sequence of events is typical of the variability found in idealized models forced by constant surface buoyancy fluxes where westward propagating unstable baroclinic planetary waves grow upon the In the damped regime (K > K c ), the effect of the NAO forcing becomes clearly apparent with the SST anomalies now circulating in a large portion of the North Atlantic from the western subtropical gyre to the mid-latitudes. At midlatitudes, the resemblance with the time evolution of SST anomalies obtained for K < K c is striking. Whether this implies that the SST variability draws its energy from the large-scale baroclinic instability mechanism, as obtained in the super-critical regime, remains to be determined however, and this will be the subject of section 4. The oscillation cycle in this regime is similar to that described previously for K < K c but now large-scale SST signals originating from the subtropics come into play. Subtropical SST anomalies are advected northeastward along the NAC from the Gulf Stream region to the eastern part of the basin at mid-latitudes from where subsequent westward propagation occurs. This pattern of variability is similar to that reported by [START_REF] Eden | North Atlantic interdecadal variability: oceanic response to the North Atlantic Oscillation (1865-1997[END_REF] and [START_REF] Eden | A damped decadal oscillation in the North Atlantic climate System[END_REF] in ocean only simulations either forced by realistic monthly mean surface fluxes associated with the NAO or coupled to a simple stochastic atmosphere model. The fact that similar patterns are obtained is consistent with the result that the internal oceanic variability is damped in [START_REF] Eden | A damped decadal oscillation in the North Atlantic climate System[END_REF]. It should be stressed that subtropical SST anomalies are also present for K < K c under stochastic surface boundary conditions but their amplitude is much smaller than those present at mid-latitudes so that their overall contribution to the North Atlantic SST variability is negligible.

d. Vertical structure of temperature anomalies

To provide further insight into the pattern of the variability, we examine here the vertical structure of temperature anomalies in the subpolar region. In all cases, temperature variability in the western subpolar region (30-60 • W, 40-60 • N) is surface intensified (Fig. 5a) and decreases with depth. In the damped regime (K > K c ), there is a sharp decrease of the temperature variability in the first 100 m and a weaker decrease below as revealed by the vertical derivatives of the stan-14 Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

dard deviations in the inset of Fig. 5a. This relatively strong surface attenuation of temperature changes is consistent with the theoretical vertical scale K v /πν 0 inferred from the heat diffusion

equation ∂ t T = K v ∂ 2 z T
where K v is the vertical mixing coefficient in the upper ocean (assumed uniform and equal to 5×10 -5 m 2 s -1 ) and ν 0 is the characteristic oscillation frequency of the surface temperature. With periods ranging from 25 to 50 years in the damped regime (see below), the vertical attenuation scale varies from 110 to 160 m, in rough agreement with model results.

In the super-critical regime the variability is strongly attenuated around 500 m depth (see inset in

e. Oscillation period

The oscillation period is deduced from the frequency with maximum power in the multi-taper spectrum of North Atlantic average kinetic energy density timeseries. A robust feature across all experiments is a consistent increase of the period with K (Fig. 6). The period typically increases from about 10-20 years in the super-critical regime to about 50 years for the most diffusive case.

Those values are in the range of those inferred from a variety of direct observations and paleoreconstructions [START_REF] Gray | A tree-ring based reconsruction of the Atlantic Multidecadal Oscillation since 1567 AD[END_REF][START_REF] Chylek | Ice-core data evidence for a prominent near 20 year time-scale of the Atlantic multidecadal oscillation[END_REF][START_REF] Knudsen | Tracking the Atlantic multidecadal oscillation through the last 8,000 years[END_REF][START_REF] Tung | Using data to attribute episodes of warming and cooling in instrumental records[END_REF].

The increase of the period with K appears consistent with the decrease of the (westward) phase speed of long baroclinic Rossby waves, given by c = β R 2 d , where R d is in the internal Rossby radius. As K increases, the circulation weakens and so does the northward ocean heat transport.

These changes induce a cooling of the North Atlantic, in particular at mid-latitudes and in the upper 1000 m (typically a 3 • C cooling when K varies from 200 to 1800 m 2 s -1 ). Salinity changes are much weaker in term of their impact on the potential density field. The temperature changes in turn induce a decrease of the stratification below 300 m and an increase above. Note that these changes in the stratification cannot be explained by the direct effect of the GM scheme because of its quasi-adiabatic character. The Rossby radius averaged between 40 with K as the NAC veers more eastward. The oscillation period is therefore found to better scale with 2L/c rather than L/c (as would be the case for a mode 1 propagation). Using length scales L increasing from 2000 km for K = 200 m 2 s -1 to 3000 km for K = 1800 m 2 s -1 , we obtain periods between 21 and 57 years, consistent with those diagnosed from the numerical model. It should be stressed that the rough agreement between the diagnosed oscillation periods in the numerical model and the theoretical values inferred from the phase speeds of long baroclinic Rossby waves does not rule out the possibility that other effects, such as the mean flow and horizontal density gradients, play a role. Determining the contribution of each of these factors to the period is not addressed here.

f. Bifurcation diagrams

Figure 7a shows that the mean AMOC strength is strongly impacted by K with a decrease of about 60% when the diffusivity increases from 200 to 1800 m 2 s -1 . This sensitivity was rationalized by [START_REF] Marshall | The dependence of the ocean's MOC on mesoscale eddy diffusivities: A model study[END_REF] using scaling laws built upon the strong interplay between the AMOC changes, Southern Ocean upwelling and strength of the abyssal cell emanating from Antarctica. In addition to the weakening of the circulation, the North Atlantic Current tends to veer more eastward as K increases (Fig. 3). Because stronger vertical shears lead to larger growth rates of (large-scale) baroclinic instability, and because the stabilizing influence of β (the meridional gradient of planetary vorticity) is maximum in the zonal direction [START_REF] Pedlosky | Geophysical Fluid Dynamics[END_REF], the simulated changes in the circulation when K increases lead to a damping of the internal ocean mode. A critical threshold is indeed confirmed and clearly visible at K = K c = 600 m 2 s -1 for all quantities under deterministic conditions (Fig. 7b-c-d). This threshold has the nature of a super-critical Hopf bifurcation, where the amplitude of oscillations in the vicinity of the bifurcation increases with the square root of the distance from the bifurcation with the Peclet number as the control parameter Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

(Colin de Verdière and Huck 1999; [START_REF] Arzel | The internal generation of the Atlantic Ocean interdecadal variability[END_REF]. For K > K c , no variability emerges under deterministic conditions since baroclinic growth rates are too weak to overcome the large damping rates associated with eddy mixing rates: the internal ocean mode is damped in this regime. For a given value of K, the annual mean AMOC strength in the RTRS experiments (where noise forcing is absent) is very close to that obtained in the FTFS runs. This shows that rectification of the long-term annual mean flow strength by stochastic forcing does not occur in our model, or at least is of minor importance.

The amplitude of the variability in the FTFS experiments is measured in terms of changes in North Atlantic kinetic energy density, AMOC strength, and western subpolar SST (the region where SST changes are maximum, Fig. 3) and is illustrated in Fig. 7b-c-d. In the super-critical regime, the effect of the noise forcing on the variability is relatively weak away from the bifurcation and strong near the bifurcation. In the damped regime, a low-frequency variability emerges unlike the deterministic case with peak-to-peak AMOC variations of 1-3 Sv depending on K and noise forcing characteristics. Interestingly, the amplitude of SST changes in the western subpolar gyre in this regime are relatively insensitive to eddy mixing rates unlike the amplitude of circulation changes that consistently decrease with increasing diffusivities. This behaviour is at odds with the common understanding that the amplitude of SST changes are positively correlated with the amplitude of circulation changes, as is the case in the super-critical regime for instance. This suggests that subpolar SST changes become decoupled from the circulation anomalies in the damped regime, an hypothesis that will be further explored in section 5. Those bifurcation diagrams further show that the AMOC, SST and extra-tropical (north of 20 • N) kinetic energy variability in the damped regime (K > 600 m 2 s -1 ) are mainly driven by noise heat flux forcing (see also the AMOC timeseries in Fig. 2) in agreement with [START_REF] Delworth | Multidecadal thermohaline circulation variability excited by atmospheric surface flux forcing[END_REF] with the stochastic wind component having a small amplifying effect. 18 Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

Energy sources of the variability

To which physical process does the SST variability mostly owe its existence? Is it primarily related to atmospheric stochastic forcing or large-scale oceanic baroclinic instability or a combination of both? How do the energy sources associated with each of those two processes depend on the background state? The analysis of the SST patterns in the previous section provides a possible answer to these questions and suggests that the physical mechanism driving the SST variability is not the same across experiments: the NAO forcing is the leading process in the damped regime

(K > K c ) whereas intrinsic ocean dynamics is dominant for K < K c .

a. Method

In order to provide a quantitative estimate of the contribution of each of these two processes (i.e. atmospheric versus oceanic energy source) in the variability we refer to the buoyancy variance budget which has proven to be a powerful tool to infer the origins of the variability. Such an approach has been previously and successfully applied to the interdecadal climate variability problem in either oceanic (Colin de Verdière and Huck 1999; [START_REF] Arzel | The different nature of interdecadal variability of the thermohaline circulation under mixed and flux boundary conditions[END_REF][START_REF] Arzel | The internal generation of the Atlantic Ocean interdecadal variability[END_REF] or coupled models [START_REF] Arzel | On the origin of interdecadal oscillations in a coupled ocean-atmosphere model[END_REF][START_REF] Arzel | Abrupt millennial variability and interdecadal-interstadial oscillations in a global coupled model: sensitivity to the background climate state[END_REF][START_REF] Buckley | On the relationship between decadal buoyancy anomalies and variability of the Atlantic meridional overturning circulation[END_REF][START_REF] Jamet | Oceanic control of multidecadal variability in an idealized coupled GCM[END_REF][START_REF] Gastineau | North Atlantic Ocean internal decadal variability: role of the mean state and ocean-atmosphere coupling[END_REF] with complexities ranging from idealized to fully coupled and realistic. We consider the linearized buoyancy variance equation

1 2 ∂ b 2 ∂t = -u h b .∇ h b -w b ∂ z b - 1 2 u.∇b 2 + b Q b + b D b , (1)
where the overbar denotes a time-mean average and the prime the perturbation. cludes the eddy transport velocity associated with the GM scheme. The GM term destroys buoyancy variance and is therefore not relevant to the growth of perturbations. The GM term is included in the dissipation term D b in the form of a skew flux. The objective here is to focus on the energy sources of the variability, that is the positive terms on the rhs of (1). Focusing on the buoyancy b rather than the temperature variance equation is suitable since temperature variations dominate the buoyancy changes for all experiments (not shown). The first term on the rhs of ( 1) is a source of buoyancy variance when transient buoyancy fluxes u h b are oriented down the mean buoyancy gradient, where u h is the horizontal Eulerian velocity. This configuration is typical of baroclinic instability for which potential energy is extracted from the mean stratification. This term has been pinpointed as the primary source of the variability in many ocean-only and coupled models (see references above). Associated with baroclinic instability is a conversion of potential to kinetic energy of perturbations through the positive exchange term w b . Under such unstable conditions, the second term in (1) is always negative (provided ∂ z b > 0 in stably stratified waters) and is therefore a sink of buoyancy variance. The third term represents the spatial redistribution of buoyancy variance by the three dimensional background flow u. It plays an important role at the regional scale by decreasing or increasing the variance, but cannot be at the very origin of the variability at a global scale since its global average is identically zero. The fourth term is a source of buoyancy variance when the surface buoyancy anomalies and the surface buoyancy flux anomalies We next take the spatial average (denoted by angle brackets below) of ( 1) over the western subpolar area (40-60 • N, 30-70 • W) where maximum SST changes consistently occur in all experiments. We define the quantities

Q b = g 0 α T Q /C o
S A =< b Q b >, S O = -< u h b .∇ h b > and R O = -< 1 2 u.∇b 2 >.
To objectively determine which of the ocean or the atmosphere explains the most the growth of interdecadal oscillations against dissipation, we concentrate in what follows on the ratio Γ = S O /S A between the oceanic and atmospheric energy sources: the origin of the growth of perturbations in the region of interest will be ascribed to internal ocean dynamics when Γ 1 and to the NAO forcing when Γ 1, while Γ = O(1) corresponds to cases where the ocean and atmosphere play equal roles in the growth of perturbations. The input of buoyancy variance by the mean currents is evaluated against the atmospheric energy source by computing the ratio Λ = R O /S A .

b. Growth of sea surface temperature variance

Focusing first on the origin of the growth of SST variance, we see that Γ 1 in the supercritical regime whereas Γ 1 in the damped regime (Fig. 8a). This shows that the NAO forcing is the leading process for generating surface buoyancy (temperature) variance in the damped regime whereas internal ocean processes associated with large-scale baroclinic instability is the leading one in the super-critical regime. The decrease of Γ with K can only be explained by a reduction in S O (the internal generation of buoyancy variance in the ocean) since S A is nearly insensitive to the eddy diffusivity K (Figs. 8b and9). The result that the covariance term S A =< b Q b > be nearly independent on K could be unexpected since the amplitude of SST changes in the damped regime is significantly less than in the super-critical regime (Fig. 7d). However the correlation (strictly equal to the normalized covariance) between SST anomalies and surface heat flux anomalies over the region of interest is considerably larger in the damped (r = 0.5 -0.6) than in super-critical regime (r = 0.2 -0.3, not shown). One explanation for such a behaviour is to note that the kinetic 21 Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

energy density variability is significantly lower in the damped regime (Fig. 7b). Low anomalous oceanic advection tend to keep the noise-forced SST anomalies in the forcing region, a process that favors relatively high correlations between the forcing and the SST field. This increase of the correlation with K compensates for the decrease in the SST variance leading to an almost unchanged covariance term S A across the range of values of K explored here. Near the bifurcation, 1) indicating that the oceanic and atmospheric energy sources contribute almost equally to the growth of SST variance in the western subpolar region. To sum up, this analysis reveals that although the internal ocean mode is clearly excited by the noise forcing in the damped regime (S O > 0), its role in the existence of SST variability in the northern North Atlantic in this regime is much weaker than that associated with the NAO forcing. At this point, it is important to recall that this analysis says nothing about the role of the ocean in setting the oscillation period (a question that will be tackled in the next section) but instead provides firm answers about the physical processes sustaining the interdecadal oscillations against all sources of thermal damping.

Γ = O(
The same conclusions hold when focusing at a specific location in the western subpolar area. damped regime, an effect that will be confirmed below and shown to be caused by the presence of a stochastic wind forced component. Similarly to the statistics averaged over the subpolar box, the atmospheric production term in the damped regime is larger than its oceanic counterpart at virtually all locations in the western subpolar region.

We finally mention that S O is similar between the deterministic and stochastic cases in the supercritical regime (as indicated by stars in Fig. 8a) demonstrating the limited ability of the noise forcing to increase the creation of buoyancy variance by internal ocean dynamics in this regime.

Advection of buoyancy variance by the mean circulation tends to extract surface buoyancy variance from the western supolar gyre in the super-critical regime but deposits surface buoyancy variance in the damped regime (Fig. 8a). The amplitude of the terms R O and S O is similar in the damped regime (at least for some K values), but their combined effect still remains much smaller than the energy input associated with the NAO forcing.

c. Growth of upper ocean heat content (UOHC) variance

When the buoyancy variance budget is carried out over the upper 1000 m rather than over the forcing layer (10 m thick), the surface forcing contribution is reduced by a factor 100. The oceanic energy source (Fig. 8e) is also reduced but much less, typically by a factor ∼ 4 in the supercritical regime and up to a factor ∼ 40 (∼ 7) in the damped regime when stochastic wind-stress forcing is present (absent). Consistently larger values of Γ are therefore obtained in this case whereas the sensitivity to K remains unchanged (Fig. 8d). The key here is that Γ becomes now larger than one over a large portion of the damped regime in contrast to what has been obtained previously for the surface buoyancy variance budget (where Γ 1). When averaged over the The analysis therefore demonstrates that fundamentally different mechanisms govern the SST and UOHC variability in the damped regime: the NAO forcing is the leading process for maintaining SST variability whereas UOHC variability is mostly sustained by internal ocean dynamics.

d. Effect of a stochastic wind component

The presence of stochastic momentum fluxes does not alter the above conclusions but has nevertheless a substantial effect on the internal generation of buoyancy variance in the ocean. Its effect is strong in the damped regime and in particular near the surface and negligible in the super-critical regime. Figure 8b shows that the presence of a stochastic wind component increases the oceanic term S O at the surface by a factor ranging from O(1) at bifurcation to about 20 for the most diffusive case compared to experiments using only stochastic surface heat fluxes. The oceanic term S O results from the interaction of transient buoyancy fluxes and time mean horizontal buoyancy gradients. These latter are very similar between the cases with and without stochastic wind forcing (not shown). As a result, the much stronger value of S O obtained when stochastic winds are present can only be caused by the much larger transient buoyancy fluxes. This feature is illustrated in Fig. 10 where the meridional contribution v b (the largest contribution to the total buoyancy fluxes) is shown at the surface for K = 1000 m 2 s -1 for cases with and without stochastic wind-stress forcing. Averaging over a greater depth considerably reduces the differences between the two cases (not shown), thereby highlighting the central role of anomalous Ekman velocities in increasing the internal generation of buoyancy variance at the surface.

With stochastic surface heat fluxes only, the NAO forcing explains about 90% (computed as the ratio S A /(S O + S A )) of the total production of surface buoyancy variance by the ocean and atmosphere in the damped regime (Fig. 8c). This ratio falls to about 65% in the presence of a stochastic wind component. Therefore it is estimated that about one third of the creation of Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

SST variance by the ocean-atmosphere system can be directly ascribed to the noise excitation of damped ocean modes. Focusing on the creation of UOHC variance, the relative contribution of the atmosphere (Fig. 8f) is close to zero in the super-critical regime and increases with K in the damped regime. For K > 1300 m 2 s -1 the dominant energy source for the UOHC variance switches from the atmosphere for noise heat flux forcing only to the ocean when both stochastic heat and wind-stress forcing are applied. For smaller values of K the ocean provides the largest contribution in all cases.

The role of circulation changes

This section examines the role that ocean circulation changes have in determining the amplitude, pattern and timescale of North Atlantic SST variability in the FTFS experiments.

In the super-critical regime (K < K c ), transient buoyancy fluxes associated with (westward) planetary wave propagation are at the heart of the existence of the variability: these fluxes are the process by which unstable waves extract energy from the mean flow to grow against all dissipative processes (Colin de Verdière and Huck 1999). It is therefore not surprising to see that in this regime circulation changes are central to the variability as we shall confirm below. In the damped regime (K > K c ), the buoyancy variance budget analysis shows that the NAO forcing is essential in maintaining the SST variability against all sources of dissipation. It is thus tempting in this case to expect circulation changes to be of minor importance, at least for determining both the amplitude and pattern of the SST variability.

To determine the role of ocean circulation changes in North Atlantic SST variability, we compare the reference experiments where the circulation is free to evolve to experiments with prescribed oceanic velocities from the climatological seasonal cycle diagnosed from the RTRS runs. In these experiments, the circulation is decoupled from the buoyancy field which is thus passively advected 25 Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

by the seasonally varying prescribed circulation but can still respond to atmospheric stochastic forcing. The noise forcing includes only a heat component (Table 1). Adding a stochastic windstress component in those experiments has no effect on the ocean circulation (which is prescribed by definition) and thereby on the oceanic tracer field. Figure 11 compares the amplitude of SST variations in the subpolar box (40-60 • N, 30-70 • W) between cases with and without circulation changes against the eddy diffusivity K. When the circulation is prescribed, the subpolar SST variance increases with eddy mixing rates: the larger the eddy diffusivity K, the slower the circulation and the larger the SST response consistent with the Hasselmann's theory modified by the addition of steady mean oceanic advection [START_REF] Saravanan | Advective ocean-atmosphere interaction: An analytical stochastic model with implications for decadal variability[END_REF]. In the super-critical regime, circulation changes substantially increase the subpolar SST variance compared to the case with prescribed oceanic currents. In the damped regime, a significant fraction (typically between 70% and 85%) of the subpolar SST variance obtained when both stochastic heat and wind-stress forcing are present is captured by the pure thermodynamic response without circulation changes. As expected, the SST patterns strongly project onto the NAO forcing when the circulation is prescribed (Fig. 12), with the leading EOF explaining about 70% of the spatially integrated annual mean SST variance. The comparison of the SST patterns between the prescribed and free circulation cases further indicates the minor (strong) impact of changes in ocean circulation on the leading pattern (Fig. 12) and amplitude of SST variability (Figs. 3 and right panels in Fig. 12) in the damped (super-critical) regime. We finally note that the pure thermodynamic response obtained with a prescribed circulation shows maximum SST variance in the western subpolar gyre. These SST changes add up to the internally-generated SST changes that also reach their maximum in this area.

We now investigate whether oceanic circulation changes are essential in setting the oscillation period. Figure 13 shows power spectra of western subpolar SST for three different values of K Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

covering both the super-critical and damped regimes. In all cases changes in ocean circulation are essential to produce a preferred timescale in the system. The purely thermodynamic response obtained with a prescribed circulation is consistent with a red noise response which demonstrates that the spatial resonance mechanism put forth by [START_REF] Saravanan | Advective ocean-atmosphere interaction: An analytical stochastic model with implications for decadal variability[END_REF] does not operate in our simulations. The case with K = 1600 m 2 s -1 and prescribed circulation (Fig. 13c) does indicate enhanced power in the 40-50 years range, as in the free circulation case, but is not statistically significant. As discussed in section 3e, Fig. 13 clearly shows that the peak period increases with K in agreement with Fig. 6.

Summary and discussion

Understanding the ocean's response to atmospheric stochastic forcing requires to separate explicitely the thermodynamic contribution from the dynamical one, the latter being associated with either self-sustained ocean modes or a noise excitation of damped ocean modes. This issue becomes fundamental when applied to the North Atlantic interdecadal climate variability problem, and more specifically to the Atlantic Multidecadal Oscillation (AMO) for which the precise roles of the ocean and atmosphere continue to be fiercely debated [START_REF] Clement | The Atlantic Multidecadal Oscillation without a role for ocean circulation[END_REF][START_REF] Zhang | Comment on "The Atlantic multidecadal oscillation without a role for ocean circulation[END_REF][START_REF] Cane | Low-pass filtering, heat flux, and Atlantic multidecadal variability[END_REF][START_REF] Zhang | On the persistence and coherence of subpolar sea surface temperature and salinity anomalies associated with the Atlantic Multidecadal Oscillation[END_REF]. In this paper a method has been proposed to objectively determine the contribution of atmospheric stochastic forcing and internal ocean dynamics to the North Atlantic SST and circulation variability, in the limit of no feedback to the atmosphere. Numerical simulations of an ocean general circulation model have been carried out at a 1 • horizontal resolution under prescribed surface fluxes including a climatological seasonal forcing and NAOrelated stochastic surface fluxes. The analysis was carried out across a range of eddy-induced diffusivities that was chosen to be sufficiently large to explore the physics of two contrasting regimes: a super-critical regime where intrinsic oceanic variability spontaneously develops in the Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

absence of any noise forcing and a damped regime where the oceanic variability requires some atmospheric noise to show up.

A buoyancy variance budget in the western subpolar region is used to objectively determine which of the ocean or atmosphere primarily sustains interdecadal oscillations against dissipation.

Our results demonstrates that the fraction of the variability explained by the ocean and atmosphere is a strong function of background oceanic conditions. In the super-critical regime, intrinsic ocean dynamics is the determining factor for all aspects of the low-frequency variability with stochastic forcing having a relatively weak impact except near bifurcation, in agreement with Frankcombe and Dijkstra (2009). In the damped regime, the analysis provides evidence of a stochastic excitation of the intrinsic ocean mode. Despite this clear stochastic excitation however, the maintenance of the SST variability in this regime is shown to be mostly caused by the NAO forcing. In contrast, upper ocean heat content (0-1000m) variability in the damped regime is mostly sustained by internal ocean dynamics. Caution must therefore be granted when interpreting low-frequency variability in terms of SST alone or upper ocean heat content alone.

Stochastic wind-stress forcing is shown to substantially increase the internal generation of buoyancy variance in the ocean. The effect is strong in the damped regime and near the surface and is shown to be caused by the much stronger transient buoyancy fluxes in relation with anomalous Ekman velocities in the western subpolar area. Without stochastic wind-stress forcing the growth of surface buoyancy variance caused by atmospheric stochastic fluxes is between one and two orders of magnitude larger than its oceanic counterpart. With stochastic wind-stress forcing, the atmospheric energy source is only about twice larger than the oceanic energy source. To put this another way, our results indicate that in the damped regime about 90% (65%) of the entire production of surface buoyancy variance is accomplished by the atmosphere when stochastic wind-forcing is absent (present). Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

The transition from the self-sustained to the damped regime produces changes in the spatial structure of the variability that are consistent with baroclinic instability. In the super-critical regime, the SST signal is strong and intensified at mid-latitudes and features a zonal dipolar structure centered around the mean path of the North Atlantic Current. In the damped regime, the SST pattern has a much broader latitudinal extent and features a basin-scale dipole extending from the western subtropical gyre to the subpolar area, in good agreement with the large-scale spatial pattern of the NAO forcing. Temperature anomalies in the super-critical regime are baroclinic with a clear westward phase shift with depth and are relatively deep. Temperature anomalies in the damped regime do not exhibit such a vertical structure and are concentrated in the thermocline. We further note that the SST variability be primarily stochastically-forced (damped regime) or internally-generated (super-critical regime) does not modify the region of peak SST variance, which is always found in the western part of the subpolar gyre.

Ocean circulation changes are shown to be unimportant for establishing the leading pattern of SST variability in the damped regime but are fundamental to select a preferred timescale in the system. Hence the spatial resonance mechanism [START_REF] Saravanan | Advective ocean-atmosphere interaction: An analytical stochastic model with implications for decadal variability[END_REF] does not occur in our simulations. The amplitude of the variability in the damped regime is to a large extent (from 70% to 85% depending on K and with stochastic wind-stress forcing) imposed by the pure thermodynamic oceanic response to atmospheric stochastic forcing. In the super-critical regime by contrast, ocean circulation changes are central to all aspects of the variability, as expected. [START_REF] Clement | The Atlantic Multidecadal Oscillation without a role for ocean circulation[END_REF] showed that ocean circulation changes are unimportant for establishing the pattern and amplitude of the North Atlantic low-frequency SST variability in fully coupled climate models. [START_REF] Clement | Response to Comment on "The Atlantic Multidecadal Oscillation without a role for ocean circulation[END_REF] and [START_REF] Colfescu | Internal atmospheric noise characteristics in twentieth century coupled atmosphere-ocean model simulations[END_REF] further argue that changes in oceanic heat transport convergence plays a minor role on interdecadal timescales in coupled climate models.

The present results suggest that this behaviour is consistent with damped interdecadal internal 29 Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

ocean modes in fully coupled models with the NAO forcing providing the main energy source for the growth of SST variance.

The above model results reveal that a clear dichotomy exists in the characteristics and leading patterns of the variability between the super-critical and damped regimes, that is remarkably captured by a single dimensionless number Γ measuring the ratio between the oceanic and atmospheric energy sources, as inferred from the buoyancy variance budget of the western subpolar region. The abrupt change in Γ around the stochastic Hopf bifurcation ( 1 in the super-critical regime and 1 in the damped regime) strongly suggests that it is a very useful quantity to objectively separate the two regimes, at least in the limit of no feedback to the atmosphere. In any case, applying this diagnostic to coupled climate model configurations would be certainly very informative about the profound nature of the variability (either sustained by atmospheric noise or driven by deterministic dynamics), as for instance recently done by [START_REF] Gastineau | North Atlantic Ocean internal decadal variability: role of the mean state and ocean-atmosphere coupling[END_REF]. Addressing this issue using observations remains unfortunately very difficult if not impossible because of the too short instrumental record compared to the timescales of the AMO and the too low spatial coverage, in particular at depth.

Nevertheless the comparison of our model results with the statistics of the observed ocean temperature record can give some hints on the relative importance of the two mechanisms, stochastic forcing and internal ocean mode, and eventually tell if the real ocean belongs to either the supercritical or damped regime. First, the stochastic forcing is based on the actual amplitude of the atmospheric NAO forcing, so we can expect the amplitude of the oceanic response to be fairly well constrained. In contrast, the amplitude of the internal ocean mode critically depends on model parameters, here the strength of eddy-diffusivity, that is not sufficiently constrained (and varying spatially) to infer the mode amplitude. The standard deviation of annual-mean SST in observations (detrended to get rid of the warming signal) is globally stronger than in the model for the damped 30 Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

regime. The peak amplitude, around 1 • C, has a similar intensity and location, east of Newfoundland, around 50 • N-45 • W, but the pattern is more widely spread over the whole subpolar gyre (the comparison in the subtropical region is probably not relevant because of the importance of airsea coupling). The vertical structure of the temperature variability in the western subpolar gyre, based on annual anomalies of the World Ocean Atlas [START_REF] Levitus | World ocean heat content and thermosteric sea level change (0-2000 m), 1955-2010[END_REF]) is also suggesting that the variability in the damped regime is too weak by a factor of 4. The characteristic sign change of temperature anomalies on the vertical in the super-critical regime is not seen in observations. However observations extend only to 700 m, a depth close to that where the sign change is found in the model. On the other hand, EOFs of SST anomalies (taken from the HADISST data set, [START_REF] Rayner | Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century[END_REF] show a dipole pattern in the meridional direction more similar to the damped regime, whereas the internal ocean mode shows a dipole pattern in the zonal direction maximum around 50 • N as was shown in Fig. 12. As a whole the comparison with observations is not fully conclusive and does not allow to rule out any of the two candidate mechanisms. Very likely, the actual ocean regime is close to the bifurcation such that the internal ocean mode strengthens the response to stochastic forcing at the surface, and increases the variability in the thermocline. If the ocean mode is super-critical, its amplitude is probably similar to the oceanic response to stochastic forcing at the surface, as found in coupled model simulations by [START_REF] Gastineau | North Atlantic Ocean internal decadal variability: role of the mean state and ocean-atmosphere coupling[END_REF].

Our experimental setup has several simplifying assumptions, the most critical one being the absence of air-sea coupling. First the effective damping rate of SST anomalies by air-sea fluxes is on the order of that associated with the large-scale baroclinic instability mechanism near bifurcation. It is thus expected that the bifurcation structure of interdecadal variability be preserved under coupling with the atmosphere with the transition between the two regimes occuring at slightly higher horizontal Peclet numbers [START_REF] Arzel | The internal generation of the Atlantic Ocean interdecadal variability[END_REF]. [START_REF] Barsugli | The basic effects of atmosphere-ocean thermal coupling on midlatitude variability[END_REF] showed that the effect of (local) ocean-atmosphere coupling is to reduce internal damping of temperature anoma-Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

lies causing greater thermal variance in both the ocean and atmosphere compared to an uncoupled situation. Whether the same amplifying effect applies to the covariance terms of the buoyancy variance equation remains to be determined, since thermal coupling between the ocean and atmosphere does not only affect the variance of each quantity (in particular the oceanic temperature and oceanic currents), but also the correlation between these quantities.

The eddy-induced diffusivity K has been used here to place the ocean state into either the damped or super-critical regimes. The coefficient K was deliberately chosen to be spatially uniform. Observationally based studies show however that eddy mixing rates are highly variable in space [START_REF] Liu | Adjoint-based estimation of eddy-induced tracer mixing parameters in the global ocean[END_REF][START_REF] Abernathey | Global surface eddy diffusivities derived from satellite altimetry[END_REF] with values ranging from O(10 2 ) to O(10 4 ) m 2 s -1 at mid-latitudes, including the Gulf Stream region. Coupled climate models traditionally use the [START_REF] Visbeck | Specification of eddy transfer coefficients in coarseresolution ocean circulation models[END_REF] parameterization or related schemes for representing the spatial heterogeneity of K. Unfortunately it is not possible from these studies to relate the different variability regimes (e.g. [START_REF] Delworth | Multidecadal thermohaline circulation variability excited by atmospheric surface flux forcing[END_REF][START_REF] Gastineau | North Atlantic Ocean internal decadal variability: role of the mean state and ocean-atmosphere coupling[END_REF] to either the magnitude or spatial distribution of K because a myriad of other aspects come obviously into play (mean flow structure, surface forcing, parameterizations, etc). [START_REF] Arzel | The internal generation of the Atlantic Ocean interdecadal variability[END_REF] carried out ocean-only experiments without atmospheric stochastic forcing and with such spatially variable K, and found the circulation to be in the super-critical regime. However, the local values of K in those experiments are not in full agreement with observational estimates. Finally, our model configuration uses a relatively low spatial resolution and does not represent mesoscale eddies. These latter do not only impact the mean current positions but also strongly interact with the larger scales. Oceanic mesoscale turbulence can force strong interannual to decadal fluctuations of the AMOC [START_REF] Le Roux | Intrinsic and atmospherically forced variability of the AMOC: insights from a large-Ensemble ocean hindcast[END_REF] and induce an inverse cascade of kinetic energy toward the larger spatial scales and lower frequencies [START_REF] Sérazin | Inverse cascades of kinetic energy as a source of intrinsic variability: a global ogcm study[END_REF]. [START_REF] Huck | Multidecadal variability of the overturning circulation in presence of eddy turbulence[END_REF] investigated the nature of the multidecadal variability in the presence of eddy turbulence using an 32 Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

idealized ocean model configuration. Mesoscale eddies were shown to strongly rectify the mean circulation but the generic mechanism driving the variability was found to be identical to that obtained at coarse resolution. In view of the buoyancy variance budget investigated in the present study, it remains to determine the impact that the oceanic mesoscale has on the internal generation of temperature variance at large scales and multidecadal periods. Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

LIST OF FIGURES

  interpolated onto the model grid using a simple gaussian filter with a width of 100 km. The model uses a climatological seasonal wind stress (Large and Yeager 2009) averaged over the years 1949-2006. 8 Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

  . The restoring procedure occurs on a monthly timescale in order to mimick seasonal variations of the 9 Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.surface buoyancy flux. The temperature and salinity restoring timescales are fixed to 10 days and 6 months respectively. These experiments, termed RTRS (Restoring T Restoring S), start from the same initial condition corresponding to the end state of a previous 6000 yr long model integration. Each RTRS run is 1200 years long, which is sufficient to reach a new equilibrium. Monthly mean surface heat and freshwater fluxes (Q T and Q S respectively) are diagnosed from the equilibrated states of each RTRS experiment to form a synthetic seasonal cycle. Stochastic surface fluxes are then added to the climatological surface heat flux Q T and observed seasonal surface wind-stress τ obs as follows

  ity K is assessed by performing twelve experiments with values of K ranging from 200 to 1800 10 Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

11

  Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

  mean circulation and stratification (Colin de Verdière and Huck 1999; te Raa and Dijkstra 2002). 13 Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

Fig. 5a )

 5a Fig.5a). There is also a clear secondary maximum at 3000 m depth in the super-critical regime,

  • N and 60 • N and between 60 • W and 30 • W has been obtained by solving the Sturm-Liouville eigenvalue problem for the vertical structure of the vertical velocity. The calculation indicates a decrease from 20 km to 15 km as K increases from 200 to 1800 m 2 s -1 , consistent with the decrease of the stratification below 300 m. Fig. 3 shows that the perturbations propagate westward from the NAC in the form of monopoles (mode 1/2) and that the zonal extent L over which this propagation occurs increases 16 Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

  The third-order term associated with advection of buoyancy variance by the disturbed flow is between one and three orders of magnitude smaller than -u b .∇b for all values of K (not shown) and has been dropped during the linearization procedure. Here the velocity u is the Eulerian velocity and ex-19 Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

(

  with g 0 is the acceleration of gravity at the sea surface, α T is the spatially varying surface thermal expansion coefficient, C o is the specific heat capacity of the first top model layer and Q is the anomalous surface heat flux) are positively correlated. The dissipation term b D b , which contains contributions from eddy-induced, vertical and isopycnal mixing processes, is a sink of buoyancy variance and is always negative. 20 Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

Figure 9

 9 Figure9shows that the oceanic production term -u h b .∇ h b is very localized along the NAC in the super-critical regime with values much larger than b Q b . The oceanic term is at least an order of magnitude larger than its atmospheric counterpart at virtually all locations in the subpolar re-

  upper 1000m, advection by the mean flow always extracts buoyancy variance from the western subpolar gyre (not shown) and thereby acts to reduce the growth of perturbations in this region. 23 Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

Fig. 1 .Fig. 4 .Fig. 6 .

 146 Fig. 1. The anomalies in turbulent (sensible + latent) surface heat flux (a) and wind-stress (b) associated with a positive NAO phase. The patterns are obtained by regressing the annual mean surface flux anomalies from[START_REF] Large | The global climatology of an interannually varying air-sea flux data set[END_REF] onto the normalized station-based winter mean (DJFM) NAO Index[START_REF] Hurrel | Trends in the North Atlantic Oscillation -Regional temperatures and precipitation[END_REF] and multiplying the patterns by one standard deviation of the NAO index. Positive fluxes of the surface heat flux are directed out of the ocean. . . . . . . . . . . . . . . . . . . 46 Fig. 2. (top) AMOC index (Sv, 1 Sv = 10 6 m 3 s -1 ) timeseries for four different values of K covering both the super-critical and damped regimes. Experiments are carried out with stochastic heat flux forcing only (gray) and with the addition of a stochastic wind-stress component (black). (bottom) Estimation of power spectra of the AMOC index timeseries with both stochastic heat and wind-stress forcing applied. The calculation is based on a multi-taper technique with 3 tapers. The smooth solid lines are the power of a red noise spectrum with the same AR1 (first order auto-regressive) coefficient as the data, and the dashed lines are the 95% confidence limits. The analysis is based on 1500 years of annual mean model output. . . . 47 Fig. 3. Standard deviations of annual mean SST anomalies in the stochastic (heat and wind-stress) FTFS experiments for four different values of the eddy-induced diffusivity K. Note the different colorscales between the super-critical (K < K c ) and damped (K > K c ) regimes. Longterm mean ocean currents averaged in the upper 250 m are superimposed. The calculation is based on 1000 years of model output. . . . . . . . . . . . . . . . 48 Fig. 4. Composites of SST anomalies and mean upper ocean (250 m) currents associated with four phases of the AMOC under deterministic conditions for K = 500 m 2 s -1 (top), and under stochastic forcing (heat and wind-stress) for K = 500, 800, 1200 and 1600 m 2 s -1 (from second row downward). The AMOC is maximum in the first column, close to its mean value and decreasing in the second column, minimum in the third column, and close to its mean value and increasing in the fourth column. The corresponding AMOC timeseries have been discussed previously and shown in Fig. 2. Different colorscales have been used between the super-critical and damped regimes because of the much lower amplitude of SST anomalies in the latter. The analysis is based of 1000 years of annual mean model output. . . . . . 49 Fig. 5. Vertical structure of temperature anomalies in the stochastic (heat and wind-stress) FTFS experiments. Standard deviation (top) of horizontally-averaged temperature anomalies in the western subpolar area (30-60 • W, 40-60 • N). First (thickness-weighted) EOF (bottom) of horizontally averaged temperature anomalies over the North Atlantic Current (50-55 • W, 25-35 • W). In average, the first EOF explains about 90% of the total variance in the super-critical regime, and 74% in the damped regime. The light (dark) grey shading indicates the spread across the super-critical (damped) regime (centered over the mean profiles ± one standard deviation). The inset in the top panel shows the vertical derivative of standard deviations of temperature anomalies in the first 1000 m. The calculation is based on 1000 years of annual mean model output. . . . . . . . . . . . . . . . . . . . . . 50 Fig. 6. Dominant timescale of the variability as a function of the eddy diffusivity K in both the deterministic and stochastic cases. The period is computed from a multi-taper spectral analysis of the North Atlantic average kinetic energy density timeseries. . . . . . . . . . 51 Fig. 7. Statistics of key indices as a function of the eddy-induced diffusivity K under deterministic and stochastic boundary conditions and for cases with (open circles) and without stochastic surface wind-stress forcing (open squares). (a) Mean strength of the AMOC (Sv) in the RTRS and stochastic FTFS experiments. The index is computed as the maximum value of

FIG. 1 .

 1 FIG. 1. The anomalies in turbulent (sensible + latent) surface heat flux (a) and wind-stress (b) associated with a positive NAO phase. The patterns are obtained by regressing the annual mean surface flux anomalies (1949-2006) from Large and Yeager (2009) onto the normalized station-based winter mean (DJFM) NAO Index (Hurrel 1995) and multiplying the patterns by one standard deviation of the NAO index. Positive fluxes of the surface heat flux are directed out of the ocean.

  FIG. 2. (top)AMOC index (Sv, 1 Sv = 10 6 m 3 s -1 ) timeseries for four different values of K covering both the super-critical and damped regimes. Experiments are carried out with stochastic heat flux forcing only (gray) and with the addition of a stochastic wind-stress component (black). (bottom) Estimation of power spectra of the AMOC index timeseries with both stochastic heat and wind-stress forcing applied. The calculation is based on a multi-taper technique with 3 tapers. The smooth solid lines are the power of a red noise spectrum with the same AR1 (first order auto-regressive) coefficient as the data, and the dashed lines are the 95% confidence limits. The analysis is based on 1500 years of annual mean model output.

  FIG. 3. Standard deviations of annual mean SST anomalies in the stochastic (heat and wind-stress) FTFS experiments for four different values of the eddy-induced diffusivity K. Note the different colorscales between the super-critical (K < K c ) and damped (K > K c ) regimes. Long-term mean ocean currents averaged in the upper 250 m are superimposed. The calculation is based on 1000 years of model output.

FIG. 6 .

 6 FIG. 6. Dominant timescale of the variability as a function of the eddy diffusivity K in both the deterministic and stochastic cases. The period is computed from a multi-taper spectral analysis of the North Atlantic average kinetic energy density timeseries.

FIG. 7 .

 7 FIG. 7. Statistics of key indices as a function of the eddy-induced diffusivity K under deterministic and stochastic boundary conditions and for cases with (open circles) and without stochastic surface wind-stress forcing (open squares). (a) Mean strength of the AMOC (Sv) in the RTRS and stochastic FTFS experiments. The index is computed as the maximum value of the overturning streamfunction below 1000 m and north of 30 • N in the North Atlantic. (b) Amplitude of North Atlantic kinetic energy density (J m -3 ) averaged in the upper 500 m and north of 20 • N. (c) Amplitude of AMOC variations (Sv). (d) Amplitude of SST changes averaged in western subpolar area (30-70 • W, 40-60 • N). The amplitude of the variability in (b-c-d) has been estimated from a composite analysis of the last 1000 years of each experiment.

FIG. 9 .

 9 FIG. 9. Surface patterns of atmospheric (top) and oceanic (bottom) energy sources (×10 -14 m 2 s -5 ) for four different values of eddy diffusivity K covering both the super-critical and damped regimes, and for cases with both noise surface heat and momentum fluxes applied. The amplitude of b Q b barely varies with K whereas -u h b .∇ h b experiences a strong decrease from the super-critical to the damped regime. The same colorscale is applied for the top panels and lower right two panels. The streamlines indicate the long-term mean upper ocean (250m) currents.

FIG. 10 .

 10 FIG. 10. Meridional transient buoyancy flux v b (10 -6 m 2 s -3 ) at the surface for K = 1000 m 2 s -1 (damped regime). The top panel has stochastic heat flux forcing only whereas the bottom panel also includes a stochastic wind forced component. The streamlines indicate the long-term mean upper ocean (250m) currents.

  FIG. 11. Impact of ocean circulation changes on the amplitude of SST changes in the western subpolar region (30-70 • W, 40-60 • N). The amplitude of the changes is estimated from a composite analysis of the last 1000 years of each experiment. The crossed thick line corresponds to the deterministic case. Black dots (open circles) correspond to the prescribed (free) circulation case with stochastic surface heat flux only. Open squares correspond to the free circulation case with stochastic surface heat and wind-stress forcing.

  FIG.12. Impact of circulation changes on SST patterns for two values of the eddy-induced diffusivity, namely K = 500 m 2 s -1 in the super-critical regime (upper) and K = 1000 m 2 s -1 in the damped regime (bottom). Shown is the leading (area-weighted) EOF of annual mean SST anomalies obtained when the circulation is free to evolve (left) and prescribed to a repeating seasonal cycle diagnosed from the RTRS runs (middle). The right panels show the standard deviations of the SST field for the prescribed circulation case. The analysis is based on 1000 years of model output. The streamlines indicate the long-term mean upper ocean (250m) currents.

1 FIG. 13 .

 113 FIG. 13. Power spectra of the western subpolar SST index, defined as the average of SST in the region (30-70 • W, 40-60 • N). Shown are the results obtained under stochastic heat and momentum fluxes for 3 different values of eddy-induced diffusivities K covering both the super-critical and damped regimes. The blue (red) lines correspond to cases where the circulation is free to evolve (prescribed). Estimation of power spectra is based on a multi-taper technique with 3 tapers. The smooth solid lines are the power of a red noise spectrum with the same AR1 (first order auto-regressive) coefficient as the data, and the dashed lines are the 99% confidence limits. The analysis is based on 1000 years of model output without any temporal filtering.

Accepted for publication in Journal of Climate. DOI 10.1175/JCLI-D-19-0522.1.

Acknowledgments. We thank Edwin Schneider and two anonymous reviewers for valuable comments that helped to improve the manuscript. We also thank A. Colin de Verdière for comments on an early draft of this manuscript. This research was supported by the MesoVarClim project funded through the French CNRS/INSU/LEFE program. The authors acknowledge the Pôle de Calcul et de Données Marines (PCDM at Ifremer, Brest) for providing DATARMOR computational resources for the MITgcm simulations. We thank the MITgcm group for making the model publicly available.

LIST OF TABLES Table 1. List of experiments. The low-frequency variability arising in the stochastic FTFS experiments is assessed against the value of the eddy-induced turbulent diffusivity K (m 2 s -1 ) and in the presence or absence of stochastic wind-stress forcing (noise heat flux forcing is always present). The role of circulation changes in North Atlantic SST variability is studied through the use of experiments with stochastic heat flux forcing only and prescribed oceanic circulation (details of the method given in section 5). The results obtained under stochastic forcing are also compared to those obtained by [START_REF] Arzel | The internal generation of the Atlantic Ocean interdecadal variability[END_REF] 200, 300, 400, 500 600,700,800,1000200, 300, 400, 500 600,700,800, 1200200, 300, 400, 500 600,700,800, , 1400200, 300, 400, 500 600,700,800, , 1600200, 300, 400, 500 600,700,800, , 1800200, 300, 400, 500 600,700,800, 200, 300, 400, 450, 500 525, 550, 575, 600 (Arzel et al., 2018) ) TABLE 1. List of experiments. The low-frequency variability arising in the stochastic FTFS experiments is assessed against the value of the eddy-induced turbulent diffusivity K (m 2 s -1 ) and in the presence or absence of stochastic wind-stress forcing (noise heat flux forcing is always present). The role of circulation changes in North Atlantic SST variability is studied through the use of experiments with stochastic heat flux forcing only and prescribed oceanic circulation (details of the method given in section 5). The results obtained under stochastic forcing are also compared to those obtained by [START_REF] Arzel | The internal generation of the Atlantic Ocean interdecadal variability[END_REF]