## Abstract

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 low-frequency 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 supercritical 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-thirds 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 supercritical 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 time scale to emerge.

## 1. 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 (Delworth et al. 2007). The AMO does not only modulate the climate of the surrounding continents on decadal to multidecadal time scales (Zhang and Delworth 2006; Knight et al. 2006) but also directly impacts marine ecosystems (Edwards et al. 2013) and Arctic sea ice (Mahajan et al. 2011; Zhang 2015). 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 us 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, modeling 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 (Delworth et al. 1993; Delworth and Mann 2000; Dong and Sutton 2005; Gastineau et al. 2018), but their respective contributions in establishing the pattern and amplitude, and even in determining the very existence, of the AMO remain elusive and model dependent.

The simplest paradigm to explain low-frequency climate variability originates from the seminal work of Hasselmann (1976), who showed that the integration of atmospheric white noise by the ocean along with its large heat capacity gives rise to a reddened spectrum. This purely thermodynamic response has been invoked by Clement et al. (2015), 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. 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 Schneider and Fan (2007), 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 multimodel mean) is in contrast with a number of observations including instrumental measurements (Tung and Zhou 2013), tree-ring records (Delworth and Mann 2000; Gray et al. 2004), ice-core records (Chylek et al. 2011), and multiproxy based reconstructions (Knudsen et al. 2011), that show enhanced variability in the 20–80-yr range in the Atlantic sector. Dommenget and Latif (2002) compared the statistics of large-scale SST variability in the midlatitudes of the Northern Hemisphere between different coupled models, slab ocean models and observations. In contrast to Clement et al. (2015), 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. Saravanan and McWilliams (1998) modified Hasselmann’s model to include steady mean oceanic advection and a spatially variable noise forcing. In contrast to Hasselmann (1976), a preferred time scale 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 and Huck 1999; te Raa and Dijkstra 2002). Instability occurs at high Peclet numbers through a Hopf bifurcation. The growth rates at bifurcation $O\u2061(1)$ yr^{−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. Arzel et al. (2018) 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 (Deser et al. 2010). The variability disappears for eddy-induced diffusivities $O\u2061(500\u20131000)$ m^{2} s^{−1} (Huck and Vallis 2001; Arzel et al. 2018) that are in the range of those derived from observations (Liu et al. 2012; Abernathey and Marshall 2013) 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 Kuhlbrodt et al. 2012) 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 (Delworth and Greatbatch 2000; Eden and Greatbatch 2003; Chen et al. 2016), others point instead to internal ocean dynamics with the noise forcing acting as an amplifier of the variability obtained under climatological surface fluxes (Zhu and Jungclaus 2008; Gastineau et al. 2018).

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 time scale) 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.

## 2. Model and experiments

### a. The ocean model

The model used for this study is the MITgcm (Marshall et al. 1997) in a configuration identical to that used by Arzel et al. (2018). The only difference lies in the surface heat and momentum fluxes that 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 Bryan and Lewis (1979) with upper and bottom values of 0.5 × 10^{−4} and 1.3 × 10^{−4} m^{2} s^{−1}, respectively. These values are in line with those inferred from large-scale inversion experiments (Lumpkin and Speer 2007), direct measurements (Waterhouse et al. 2014) 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 (Redi 1982). 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 Jackett and McDougall (1995), which computes the in situ density from potential temperature, practical salinity and Boussinesq hydrostatic pressure. Ocean bathymetry is taken from the historical ETOPO1 dataset (Amante and Eakins 2009) 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.

### b. Experimental design

We use flux boundary conditions at the surface for both temperature and salinity, similar to Arzel et al. (2018). 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 time scales much longer than the atmospheric response time, typically 10 days, atmospheric thermal damping of SST anomalies is relatively weak. Vallis (2009) estimates this damping time scale 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; Arzel et al. 2018). Arzel et al. (2018) showed that the addition of a surface restoring flux with a damping time scale *α*^{−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., toward 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 Bryan (1987), 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 Arzel et al. (2018) 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; Antonov et al. 2010). The restoring procedure occurs on a monthly time scale in order to mimic seasonal variations of the surface buoyancy flux. The temperature *T* and salinity *S* restoring time scales are fixed to 10 days and 6 months, respectively. These experiments, termed “restoring *T* restoring *S*” (RTRS), 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:

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 (1949–2006) from Large and Yeager (2009) onto the winter mean [December–March (DJFM)] NAO index (Hurrell 1995) multiplied by one standard deviation of the NAO index (Fig. 1). The stochastic forcing is only applied to the North Atlantic. The random discrete time series *ζ*(*t*) has been built from a first-order autoregressive (AR1) process with a decorrelation time scale of 10 days. This time scale corresponds to estimates inferred by Feldstein (2000) 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 time series of 5 760 000 points for a 1500 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 time series built from the monthly means of *ζ*(*t*) has a variance equal to 1, similar to Herbaut et al. (2002). The sensitivity of the model to the eddy-induced diffusivity *K* is assessed by performing 12 experiments with values of *K* ranging from 200 to 1800 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 (Abernathey and Marshall 2013). 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).

## 3. 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 time series 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 Arzel et al. (2018), the critical value *K* = *K*_{c} = 600 m^{2} s^{−1} corresponds to the existence of a threshold separating a supercritical 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.

### b. Patterns of temperature variability

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 midlatitudes, typically between 40° 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 (Arzel et al. 2018). 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 supercritical 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.

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, Arzel et al. (2018) 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 supercritical 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 mean circulation and stratification (Colin de Verdière and Huck 1999; te Raa and Dijkstra 2002).

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 midlatitudes. 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 supercritical 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 midlatitudes from where subsequent westward propagation occurs. This pattern of variability is similar to that reported by Eden and Jung (2001) and Eden and Greatbatch (2003) 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 Eden and Greatbatch (2003). 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 midlatitudes 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 standard deviations in the inset of Fig. 5a. This relatively strong surface attenuation of temperature changes is consistent with the theoretical vertical scale $K\upsilon /\pi \nu 0$ inferred from the heat diffusion equation $\u2202tT=K\upsilon \u2202z2T$ where *K*_{υ} 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 supercritical regime the variability is strongly attenuated around 500-m depth (see inset in Fig. 5a). There is also a clear secondary maximum at 3000-m depth in the supercritical regime, much less pronounced in the damped regime, that coincides with the depth at which temperature anomalies are exported southward from the convective region along the deep western boundary current. The vertical structure of temperature anomalies is deduced from an empirical orthogonal function (EOF) analysis of horizontally averaged annual-mean temperature anomalies over a region encompassing the mean path of the NAC (50°–55°N, 25°–35°W). In the supercritical regime, the temperature anomalies are strongly phase-shifted on the vertical and change sign around 600-m depth (Fig. 5b). This should not be surprising since the North Atlantic Current is baroclinically unstable in this regime and a (westward) vertical tilt of buoyancy anomalies is exactly what is required for the waves to be unstable. Such a vertical organization of temperature anomalies is not captured by the first EOF in the damped regime. The second EOF (accounting for 21% of the variance, not shown) does, however, capture a clear sign change around a depth of 250 m in very good agreement with the same EOF in the supercritical regime (not shown). This suggests that the interdecadal mode characteristic of the supercritical regime is excited by the noise forcing in the damped regime, in agreement with Frankcombe et al. (2009). A remarkable feature is the radically different flavors of temperature variability between the two regimes (Fig. 5). Clearly, the characteristics of the temperature variability in the western subpolar area vary abruptly around the critical threshold at *K* = *K*_{c}. The same is true in the eastern part of the basin at midlatitudes (not shown).

### e. Oscillation period

The oscillation period is deduced from the frequency with maximum power in the multitaper spectrum of North Atlantic average kinetic energy density time series. 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 supercritical 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 paleo-reconstructions (Gray et al. 2004; Chylek et al. 2011; Knudsen et al. 2011; Tung and Zhou 2013). 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=\beta Rd2$, 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 midlatitudes 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 terms 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° and 60°N and between 60° 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 to 15 km as *K* increases from 200 to 1800 m^{2} s^{−1}, consistent with the decrease of the stratification below 300 m. Figure 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 with *K* as the NAC veers more eastward. The oscillation period is therefore found to better scale with 2*L*/*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 Marshall et al. (2017) 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 (Pedlosky 1987), 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 (Figs. 7b–d). This threshold has the nature of a supercritical 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 (Colin de Verdière and Huck 1999; Arzel et al. 2018). 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 Figs. 7b–d. In the supercritical 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 (1 Sv ≡ 10^{6} m^{3} s^{−1}) 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 behavior is at odds with the common understanding that the amplitudes of SST changes are positively correlated with the amplitude of circulation changes, as is the case in the supercritical regime, for instance. This suggests that subpolar SST changes become decoupled from the circulation anomalies in the damped regime, a hypothesis that will be further explored in section 5. Those bifurcation diagrams further show that the AMOC, SST and extratropical (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 time series in Fig. 2) in agreement with Delworth and Greatbatch (2000) with the stochastic wind component having a small amplifying effect.

## 4. 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

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; Arzel et al. 2006, 2018) or coupled models (Arzel et al. 2007, 2012; Buckley et al. 2012; Jamet et al. 2016; Gastineau et al. 2018) with complexities ranging from idealized to fully coupled and realistic. We consider the linearized buoyancy variance equation

where the overbar denotes a time-mean average and the prime the perturbation. The third-order term associated with advection of buoyancy variance by the disturbed flow is between one and three orders of magnitude smaller than $\u2212u\u2032b\u2032\xaf\u22c5\u2207b\xaf$ for all values of *K* (not shown) and has been dropped during the linearization procedure. Here the velocity **u** is the Eulerian velocity and excludes 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 $Db\u2032$ 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 $uh\u2032b\u2032\xaf$ 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\u2032b\u2032\xaf$. Under such unstable conditions, the second term in (1) is always negative (provided $\u2202zb\xaf>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\xaf$. 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 $Qb\u2032=g0\alpha TQ\u2032/Co$ (where *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\u2032Db\u2032\xaf$, which contains contributions from eddy-induced, vertical, and isopycnal mixing processes, is a sink of buoyancy variance and is always negative.

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 $SA=\u2329b\u2032Qb\u2032\xaf\u232a$, $SO=\u2212\u2329uh\u2032b\u2032\xaf\u22c5\u2207hb\xaf\u232a$ and $RO=\u2212\u2329\u2061(1/2)u\xaf\u22c5\u2207b\u20322\xaf\u232a$. 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 $\Gamma =O\u2061(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 supercritical 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 and 9). The result that the covariance term $SA=\u2329b\u2032Qb\u2032\xaf\u232a$ be nearly independent of *K* could be unexpected since the amplitude of SST changes in the damped regime is significantly less than in the supercritical 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 supercritical regime (*r* = 0.2–0.3, not shown). One explanation for such a behavior is to note that the kinetic energy density variability is significantly lower in the damped regime (Fig. 7b). Low anomalous oceanic advection tends 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, $\Gamma =O\u2061(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.

The same conclusions hold when focusing at a specific location in the western subpolar area. Figure 9 shows that the oceanic production term $\u2212uh\u2032b\u2032\xaf\u22c5\u2207hb\xaf$ is very localized along the NAC in the supercritical regime with values much larger than $b\u2032Qb\u2032\xaf$. The oceanic term is at least an order of magnitude larger than its atmospheric counterpart at virtually all locations in the subpolar region in this regime. The covariance term $b\u2032Qb\u2032\xaf$ has a much broader spatial structure with subpolar and subtropical centers of action corresponding to those of the NAO forcing. Beyond the bifurcation, the oceanic generation of buoyancy variance along the NAC falls drastically and becomes more uniformly distributed across the western subpolar region. As such the patterns of $b\u2032Qb\u2032\xaf$ and $\u2212uh\u2032b\u2032\xaf\u22c5\u2207hb\xaf$ in the western subpolar region in the damped regime bear some resemblance as opposed to what occurs in the supercritical regime. This resemblance suggests that the NAO forcing projects similarly onto the atmospheric and oceanic production terms in the subpolar area in the 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 subpolar gyre in the supercritical 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 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 of 100. The oceanic energy source (Fig. 8e) is also reduced but much less, typically by a factor of ~4 in the supercritical regime and up to a factor of ~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 upper 1000 m, 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. The analysis therefore demonstrates that fundamentally different mechanisms govern the SST and upper-ocean heat content (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 supercritical 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 $\upsilon \u2032b\u2032\xaf$ (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 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 supercritical 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.

## 5. The role of circulation changes

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

In the supercritical 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 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 wind stress 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 (Saravanan and McWilliams 1998). In the supercritical 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 (Fig. 3 and right panels in Fig. 12) in the damped (supercritical) 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* covering both the supercritical and damped regimes. In all cases changes in ocean circulation are essential to produce a preferred time scale 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 Saravanan and McWilliams (1998) 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-yr 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.

## 6. Summary and discussion

Understanding the ocean’s response to atmospheric stochastic forcing requires us to separate explicitly 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 (Clement et al. 2015; Zhang et al. 2016; Cane et al. 2017; Zhang 2017). 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 NAO-related 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 supercritical regime where intrinsic oceanic variability spontaneously develops in the 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 demonstrate that the fraction of the variability explained by the ocean and atmosphere is a strong function of background oceanic conditions. In the supercritical 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–1000 m) 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).

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 supercritical regime, the SST signal is strong and intensified at midlatitudes 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 supercritical 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 (supercritical 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 time scale in the system. Hence the spatial resonance mechanism (Saravanan and McWilliams 1998) 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 supercritical regime by contrast, ocean circulation changes are central to all aspects of the variability, as expected. Clement et al. (2015) 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. Clement et al. (2016) and Colfescu and Schneider (2017) further argue that changes in oceanic heat transport convergence plays a minor role on interdecadal time scales in coupled climate models. The present results suggest that this behavior is consistent with damped interdecadal internal 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 supercritical and damped regimes, which 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 supercritical 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 Gastineau et al. (2018). Addressing this issue using observations remains unfortunately very difficult if not impossible because of the too-short instrumental record compared to the time scales 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 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 air–sea coupling). The vertical structure of the temperature variability in the western subpolar gyre, based on annual anomalies of the *World Ocean Atlas* (Levitus et al. 2012) 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 supercritical 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 dataset; Rayner et al. 2003) 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 us 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 supercritical, its amplitude is probably similar to the oceanic response to stochastic forcing at the surface, as found in coupled model simulations by Gastineau et al. (2018).

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 occurring at slightly higher horizontal Peclet numbers (Arzel et al. 2018). Barsugli and Battisti (1998) showed that the effect of (local) ocean–atmosphere coupling is to reduce internal damping of temperature anomalies 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 supercritical 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 (Liu et al. 2012; Abernathey and Marshall 2013) with values ranging from $O\u2061(102)$ to $O\u2061(104)$ m^{2} s^{−1} at midlatitudes, including the Gulf Stream region. Coupled climate models traditionally use the Visbeck et al. (1997) 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., Delworth and Greatbatch 2000; Gastineau et al. 2018) 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.). Arzel et al. (2018) carried out ocean-only experiments without atmospheric stochastic forcing and with such spatially variable *K*, and found the circulation to be in the supercritical 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 (Le Roux et al. 2018) and induce an inverse cascade of kinetic energy toward the larger spatial scales and lower frequencies (Sérazin et al. 2018). Huck et al. (2015) investigated the nature of the multidecadal variability in the presence of eddy turbulence using an 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.

## 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.

## REFERENCES

*J. Geophys. Res. Oceans*

*Salinity*. Vol. 2,

*World Ocean Atlas 2009*, NOAA Atlas NESDIS 69, 184 pp.

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*Tellus*

*Climate Dyn.*

*J. Climate*

*J. Atmos. Sci.*

*J. Phys. Oceanogr.*

*J. Geophys. Res.*

*J. Climate*

*J. Climate*

*Climate Dyn.*

*Geophys. Res. Lett.*

*Science*

*Science*

*Climate Dyn.*

*J. Phys. Oceanogr.*

*J. Climate*

*Climate Dyn.*

*J. Climate*

*Ocean Circulation: Mechanisms and Impacts—Past and Future Changes of Meridional Overturning*.

*Geophys. Monogr.*, Vol. 173, Amer. Geophys. Union, 131–148.

*Annu. Rev. Mar. Sci.*

*Climate Dyn.*

*J. Climate*

*J. Climate*

*J. Climate*

*PLOS ONE*

*J. Climate*

*Geophys. Res. Lett.*

*J. Phys. Oceanogr.*

*J. Geophys. Res.*

*J. Phys. Oceanogr.*

*Geophys. Res. Lett.*

*Tellus*

*J. Phys. Oceanogr.*

*Tellus*

*J. Climate*

*J. Phys. Oceanogr.*

*Science*

*J. Atmos. Oceanic Technol.*

*Climate Dyn.*

*Geophys. Res. Lett.*

*Nat. Commun.*

*Ocean Modell.*

*Climate Dyn.*

*J. Climate*

*Geophys. Res. Lett.*

*J. Phys. Oceanogr.*

*Temperature*. Vol. 1,

*World Ocean Atlas 2009*, NOAA Atlas NESDIS 68, 184 pp.

*J. Phys. Oceanogr.*

*J. Climate*

*J. Geophys. Res.*

*Ocean Modell.*

*Geophysical Fluid Dynamics*. 2nd ed. Springer-Verlag, 710 pp.

*J. Geophys. Res.*

*J. Phys. Oceanogr.*

*J. Climate*

*J. Atmos. Sci.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*Proc. Natl. Acad. Sci. USA*

*Stochastic Physics and Climate Modelling*, T. Palmer and P. Williams, Eds., Cambridge University Press, 1–33.

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*Proc. Natl. Acad. Sci. USA*

*Geophys. Res. Lett.*

*Geophys. Res. Lett.*

*Science*

*Climate Dyn.*