Olivier Arzel
email: olivierarzel@univ-brest.fr

Thierry Huck

Antoine Hochet

Alexandre Mussa

Internal ocean dynamics contribution to North Atlantic interdecadal variability strengthened by ocean-atmosphere thermal coupling

Identifying the primary drivers of North Atlantic interdecadal climate variability is crucial for improving climatic prediction over the coming decades. Here the effect of thermal coupling on the leading energy sources of the interdecadal variability of the ocean-atmosphere system is examined by means of a stochastically-forced idealized coupled model. The effect of coupling is quantified from a comparison of the buoyancy variance budget of coupled and uncoupled model configurations. The simplicity of the model allows us to contrast the effect of coupling between a super-critical regime where the deterministic ocean dynamics drive the variability and a damped regime where noise forcing is central to its existence. The results show that changes in surface buoyancy fluxes act as a sink of temperature variance in the super-critical regime, and only become a source in the strongly damped regime. By contrast, internal ocean dynamics associated with the interaction of transient buoyancy fluxes with mean buoyancy gradients always act as a source of interdecadal variability. In addition to the reduced thermal damping effect in coupled integrations, thermal coupling with the atmosphere is shown to significantly increase the role of internal ocean dynamics in the variability, in particular in the regime where interdecadal modes are damped. Only for oceanic background states in the strongly damped regime do changes in surface buoyancy fluxes play a leading role in the upper ocean variability. A stochastically-forced coupled box model is proposed that captures the basic effect of thermal coupling on atmospheric and oceanic energy sources of variability.

2 SIGNIFICANCE STATEMENT: The purpose of this study is to better understand the impact of ocean-atmosphere thermal coupling on the leading energy sources of Atlantic interdecadal variability. Increasing our understanding of the physical mechanisms driving climate variability at interdecadal timescales is important to improve climate prediction. We show that the effect of ocean-atmosphere thermal coupling is to substantially increase the role of internal ocean dynamics in the low-frequency variability of the upper ocean heat content and sea surface temperature.

Atmospheric stochastic forcing only becomes the primary driver of the oceanic temperature variability in the large dissipative limit, when internal ocean modes are strongly damped.

Introduction

Multidecadal variability of North Atlantic sea surface temperatures has been observed in both the instrumental record [START_REF] Delworth | Decadal to centennial variability of the Atlantic from observations and models[END_REF][START_REF] Deser | Sea surface temperature variability: patterns and mechanisms[END_REF][START_REF] Tung | Using data to attribute episodes of warming and cooling in instrumental records[END_REF] and a number of climatic proxy reconstructions [START_REF] Mann | Global-scale temperature patterns and climate forcing over the past six centuries[END_REF][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].

Atlantic Multidecadal Variability (AMV, [START_REF] Sutton | Atlantic Multidecadal Variability and the U.K. ACSIS program[END_REF]) is characterized by alternating basinwide decadal-scale warming and cooling periods with maximum variance at subpolar latitudes. [START_REF] Clement | The Atlantic Multidecadal Oscillation without a role for ocean circulation[END_REF] suggested that the AMV could result from the thermodynamic response to atmospheric stochastic forcing and that ocean circulation variations are unimportant. This mechanism seems however incompatible with the paleorecord that exhibits statistically significant multidecadal periods above a red-noise background (e.g. [START_REF] Gray | A tree-ring based reconsruction of the Atlantic Multidecadal Oscillation since 1567 AD[END_REF][START_REF] Wang | Internal and external forcing of multidecadal Atlantic climate variability over the past 1,200 years[END_REF]. [START_REF] O'reilly | The signature of low-frequency oceanic forcing in the Atlantic Multidecadal Oscillation[END_REF] showed instead that that ocean circulation variations are key to reproduce the sign of the observed correlation between sea surface temperatures (SST) and the AMV. [START_REF] Zhang | A review of the role of the Atlantic Meridional Overturning Circulation in Atlantic Multidecadal Variability and associated climate impacts[END_REF] reviewed the observational and modelling evidence for the origins of the AMV and concluded that internal ocean dynamics is a key driver. The existence of such internal variability was recently questioned by [START_REF] Mann | Multidecadal climate oscillations during the past millennium driven by volcanic forcing[END_REF] who found instead a prominent role of volcanic forcing and anthropogenic aerosols. Which of internal ocean dynamics, atmospheric stochastic forcing or external radiative forcing lies at the very origin of the observed AMV remains therefore much debated. The instrumental record is unfortunately too short compared to the timescale of the variability and too sparse in terms of spatial coverage to obtain a definitive answer to this question and climate models have yet to reach a consensus [START_REF] Zhang | A review of the role of the Atlantic Meridional Overturning Circulation in Atlantic Multidecadal Variability and associated climate impacts[END_REF]. The present study will solely focus on the variability generated internally in the ocean-atmosphere system. Even in this con-text numerical models are inconclusive, in part because they show that the relative contributions of the ocean and atmosphere to the variability are a strong function of poorly constrained critical parameters, such as turbulent eddy diffusivities associated with unresolved scales [START_REF] Arzel | Contributions of atmospheric stochastic forcing and intrinsic ocean modes to North Atlantic ocean interdecadal variability[END_REF]. In addition ocean-atmosphere interactions certainly have an impact on the amplitude of the variability, at least at low-frequencies [START_REF] Barsugli | The basic effects of atmosphere-ocean thermal coupling on midlatitude variability[END_REF], but the net effect of those interactions on the leading energy sources of the variability have yet to be found. This study will thus focus on the effect of thermal coupling on the primary drivers of the interdecadal variability of the extratropical ocean-atmosphere system.

The effect of ocean-atmosphere coupling on the decadal or longer variability is traditionnaly evaluated by comparing fully coupled ocean-atmosphere simulations to ocean-only experiments where the surface heat, momentum and freshwater fluxes are diagnosed from the coupled run. Under such flux boundary conditions atmospheric damping on sea surface temperature (SST) anomalies is reduced to zero and the amplitude of oscillations is maximized. This approach has been used in many modelling studies and across a large diversity of model configurations and complexities. [START_REF] Delworth | Multidecadal thermohaline circulation variability excited by atmospheric surface flux forcing[END_REF] and [START_REF] Gastineau | North Atlantic Ocean internal decadal variability: role of the mean state and ocean-atmosphere coupling[END_REF] performed such experiments with comprehensive coupled general circulation models (GCMs) to show that coupling with the atmosphere was not essential to the variability, and that atmospheric stochastic forcing associated with mid-latitudes disturbances was necessary to sustain interdecadal oscillations against dissipation.

One drawback of such an approach is that the surface fluxes used to drive ocean-only integrations contain information about the response of the atmosphere to SST changes. The comparison of those ocean-only runs with coupled integrations does not therefore allow to isolate the effect of the feedback of the atmosphere on the SST anomalies, but instead to assess the entire effect of surface heat flux damping on SST anomalies. Another approach consists in forcing the ocean with the time history of surface air temperatures (SAT) and winds diagnosed from an atmospheric-only integration itself forced by climatological SSTs from the coupled simulation [START_REF] Barsugli | The basic effects of atmosphere-ocean thermal coupling on midlatitude variability[END_REF]. These "uncoupled" ocean-only experiments therefore lack the feedback of the atmosphere onto the SST anomalies. A direct comparison with the coupled runs thus allows to unambigously isolate the basic impact of this feedback (i.e. thermal coupling) on the variability of the coupled system. The absence of atmospheric feedback in uncoupled experiments makes the phasing of SST and SAT anomalies less coherent, leading to a smaller persistence of SST anomalies in uncoupled than in coupled integrations. The basic effect of thermal coupling is therefore to enhance the temperature variance in both the ocean and atmosphere, and more specifically at low-frequencies [START_REF] Barsugli | The basic effects of atmosphere-ocean thermal coupling on midlatitude variability[END_REF]. [START_REF] Farneti | Mechanisms of interdecadal climate variability and the role of ocean-atmosphere coupling[END_REF] performed uncoupled atmosphere-only experiments in the context of interdecadal climate variability and reached the same conclusions. The spectral peak at interdecadal periods in the SAT was shown to be only present when prescribed, time-varying, SST from the coupled run were used, indicating an influence of either an oceanic interdecadal mode or some form of coupled interaction. [START_REF] Wu | North Atlantic decadal variability: air-sea coupling, oceanic memory, and potential Northern Hemisphere resonance[END_REF] run an uncoupled experiment with their realistic coupled model to show the critical role of coupling in sustaining North Atlantic Ocean decadal variability. [START_REF] Arzel | Contributions of atmospheric stochastic forcing and intrinsic ocean modes to North Atlantic ocean interdecadal variability[END_REF] used a realistic ocean general circulation model forced by combination of steady surface buoyancy and momentum fluxes and North Atlantic Oscillation-related stochastic forcing to show that internal ocean dynamics plays a leading role in the growth of SST variance in the super-critical regime (i.e. the regime where the multidecadal variability develops without any noise forcing, at relatively low diffusivities), but that the atmospheric stochastic forcing represents the major energy source for the variability in the damped regime (i.e. the regime where the variability requires some atmospheric noise to emerge, at relatively high diffusivities). Therefore the leading energy source of the variability is a strong function of critical parameters (such as eddy diffusivities) and background climate conditions, at least in the forced case. How ocean-atmosphere thermal coupling modifies this behaviour remains to be determined. The present study will assess the impact of thermal coupling on the leading energy sources of the variability and is therefore a natural extension of [START_REF] Arzel | Contributions of atmospheric stochastic forcing and intrinsic ocean modes to North Atlantic ocean interdecadal variability[END_REF] who used an ocean-only configuration. The aim of the present study is thus to progress in the understanding of the physical processes driving Atlantic climate variability at interdecadal timescales. The main purpose is to quantify the impact of oceanatmosphere thermal coupling on the drivers of temperature variance at interdecadal timescales. As such the present study is also a natural extension of [START_REF] Barsugli | The basic effects of atmosphere-ocean thermal coupling on midlatitude variability[END_REF] who focused their analysis on the impact of coupling on the temperature variance of the ocean-atmosphere system.

A joint objective is to determine how thermal coupling impact the energy sources of the variability in both the super-critical and damped regimes. Different variability regimes will be obtained by systematically varying the magnitude of the turbulent eddy diffusivity K which has been shown to be one of the most critical parameter to the interdecadal variability of the overturning circula-tion (Huck et al. 1999a). For each value of K the effect of thermal coupling is inferred from a comparison of coupled, uncoupled and forced integrations. Given the relatively large number of millenial-scale experiments required to achieve this task a computationally efficient coupled model is required. Following the views [START_REF] Held | The gap between simulating and understanding in climate modeling[END_REF] we believe that the use of a model hierarchy such as the one proposed here is necessary to help gain understanding of the mechanisms driving climate variability at interdecadal timescales. The model comprises a planetary geostrophic oceanic component coupled to a stochastically-forced energy balance model. The geometry is idealized with a two-hemisphere sector of sphere and a southern periodic channel. The model captures the large-scale features of the Atlantic circulation and for some parameter range exhibit interdecadal oscillations under deterministic forcing conditions. The cause of the variability was shown to be a large-scale baroclinic instability of the mean flow driving interdecadal oscillations with the oscillation period set by the zonal transit time of long baroclinic planetary waves (Colin de Verdière and Huck 1999;[START_REF] Te Raa | Instability of the thermohaline circulation on interdecadal timescales[END_REF]. This mode of variability has been shown to persist in realistic configurations of ocean-only [START_REF] Sévellec | The leading, interdecadal eigenmode of the Atlantic meridional overturning circulation in a realistic ocean model[END_REF][START_REF] Arzel | The internal generation of the Atlantic Ocean interdecadal variability[END_REF]) and coupled models [START_REF] Ortega | Reconciling two alternative mechanisms behind bi-decadal variability in the North Atlantic[END_REF][START_REF] Muir | Evidence for the AMOC interdecadal mode related to westward propagation of temperature anomalies in CMIP5 models[END_REF][START_REF] Gastineau | North Atlantic Ocean internal decadal variability: role of the mean state and ocean-atmosphere coupling[END_REF].

The outline of this paper is as follows. Section 2 provides a description of the idealized coupled ocean-atmosphere model as well as the experimental procedure used to quantify the effect of ocean-atmosphere thermal coupling. The main characteristics of the variability are presented in section 3, and its sensitivity to horizontal turbulent diffusivity is presented in section 4 for the coupled, uncoupled and forced cases. Section 5 quantifies the leading energy sources of the interdecadal variability and their sensitivity to K using an approach based on the buoyancy variance budget. The model hierarchy is finally extended to a coupled stochastically-forced oceanatmosphere dynamical system in section 6 that reproduces with very high fidelity the results obtained with the 3D model. The results are then summarized and discussed in section 7.

Model and experiments

a. The coupled model

The 3D ocean model is based on the planetary geostrophic equations, valid for time scales much larger than the inertial period and spatial scales much larger than the internal Rossby radius of deformation [START_REF] Salmon | A simplified linear ocean circulation theory[END_REF][START_REF] De Verdière | Buoyancy driven planetary flows[END_REF]. The domain is a flat-bottom sector of a and Sarachik (1990) (Fig. 1a).

The atmospheric model is a standard, two-dimensional, dry energy balance model. It closely follows the scheme developed by [START_REF] Fanning | An atmospheric energy-moisture balance model: climatology, interpentadal climate change, and coupling to an ocean general circulation model[END_REF]. The freshwater forcing is imposed and follows the latitudinal profile shown in Fig. 1b. The atmospheric eddy temperature diffusivity K a is uniform and fixed to 1.3×10 6 m 2 s -1 . The atmosphere-ocean heat exchange coefficient λ falls within the range of observed values and is fixed to 30 W m -2 K -1 is all experiments. Both the incoming solar radiation S at the top of the atmosphere and the planetary albedo α p follow annual mean latitudinal profiles given by [START_REF] North | Energy balance climate models[END_REF] and [START_REF] Graves | New parametrizations and sensitivities for simple climate models[END_REF], respectively. In the absence of atmospheric dynamics, stochastic forcing is introduced in the atmospheric temperature equation as

Q N (x, y, t) = σ N Q NAO (x, y)ζ(t), (1)
This term parameterizes the stochastic forcing associated with the divergence of eddy temperature fluxes which are typically enhanced at mid-latitudes. There the observed transient sensible and latent heat fluxes amount to about P = 3.5 PW [START_REF] Peixoto | Physics of climate[END_REF]. Using an eddy length scale L eddy of 1000 km and a zonal length scale L x of 25,000 km for the length of a latitude circle at mid-latitudes, we obtain σ N = P/L x L eddy = O(100) W m -2 for the amplitude of the eddy forcing. We then specify σ N = 100 W m -2 in all stochastically-forced experiments. The pattern Q NAO is similar to that used by [START_REF] Herbaut | Response of a simplified oceanic general circulation model to idealised NAO-like stochastic forcing[END_REF]. It mimicks the observed North Atlantic

Oscillation pattern with a dipolar structure centered at mid-latitudes and with maximum values in the west of the basin. 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. To assess the effect of noise coherence additional experiments using a white noise forcing and a longer temporal coherence (30 days) have been conducted. The noise forcing has a sampling frequency (0.05 days) corresponding to the oceanic timestep of the model. The variance of ζ(t) is 1. Note finally that the stochastic timeseries are stricly identical in all stochastically-forced integrations.

b. Experimental design

Three different model configurations (coupled, uncoupled, forced) are used to assess the impact of ocean-atmosphere thermal coupling on the variability. For each of those three configurations the model is run with and without atmospheric stochastic forcing. The latter case is said to be deterministic. This procedure is repeated for a wide range of oceanic horizontal eddy diffusivities K between 500 and 2600 m 2 s -1 . 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]. This approach allows us to explore the effect of ocean-atmosphere coupling in two different variability regimes, namely a super-critical regime where the variability spontaneously develops without atmospheric stochastic forcing and a damped regime where the variability requires some noise to emerge. Deterministic runs are integrated for 5000 years and are initialized from a resting ocean with uniform temperature (4 • C) and salinity (35 psu). Stochastic runs are integrated for 1000 years and are initialized from the end state of deterministic runs.

(1) Coupled. The model is integrated in its standard coupled mode (denoted by superscript C)

where the atmospheric temperature T a and the surface buoyancy (heat and freshwater) flux B O toward the ocean are given by

C a dT C a /dt = F(T C a , T C o ) + Q N , (2)
B C O = G(T C a , T C o ), (3)
where C a is the atmospheric heat capacity and T o the sea surface temperature. The operators F and G include both turbulent and radiative fluxes,

F(T a , T o ) = λ(T o -T a ) + o σT 4 o -a σT 4 a -p σT 4 a + S a s (1 -α p ) + K a ∇ 2 T a , (4)
B O (T a , T o ) = λ(T a -T o ) -o σT 4 o + a σT 4 a + S (1 -a s )(1 -α p ) ( 5
)
where σ is the Stefan-Boltzmann constant, o = 0.96 is the longwave emissivity of the surface, a = 0.85 is the longwave emissivity of the atmosphere, p = 0.63 is the planetary emissivity to space, a s = 0.3 is the shortwave absorptivity. The terms on the right hand side of (4) are, from first to last, the sensible heat exchange with the surface, the upward longwave emission from the surface, the downward longwave emission from the atmosphere, the outgoing longwave radiation to space, the absorbed solar radiation within the atmosphere and the diffusive heat transport. The last term on the right hand side of (5) is the absorbed solar flux at the surface.

(2) Uncoupled. (a) The atmosphere is first run with the fixed 100 years averaged SST distribution diagnosed from deterministic coupled integrations. This SST field is defined here as the climatological SST field T clim o . (b) The ocean model is then forced by the surface fluxes computed using the time history of atmospheric temperatures from phase (a) and the actual SST field from the uncoupled integration. Those two phases are run synchronously with the same coupling frequency as the coupled model (0.05 days). Hence

C a dT U a /dt = F(T U a , T clim o ) + Q N , (6)
B U O = G(T U a , T U o ), (7)
where the superscript U stands for uncoupled. In those uncoupled runs the atmospheric temperature can only change in response to the stochastic forcing, all other forcings including the SST field being kept constant in time. The atmosphere therefore does not see the SST variations that can occur in response to either atmospheric stochastic forcing (if any, embedded in T U a ) or changes due to the existence of an intrinsic interdecadal ocean mode. In the coupled integrations by contrast the atmosphere can vary in response to both atmospheric stochastic forcing and changing

SSTs. The comparison of the coupled and uncoupled runs therefore allows to unambigously isolate the impact of the atmospheric feedback on SST anomalies.

(3) Forced. The main motivation behind these ocean-only forced experiments is to determine whether the variability obtained in the coupled model can be interpreted as a manifestation of a self-sustained ocean mode or as an excitation of a damped ocean mode. The procedure here is very similar to that used for building uncoupled runs except for the previous phase (b) where the surface heat flux to the ocean B F O is now computed using the climatological SST T clim o rather than the actual SST obtained in those ocean-only integrations, thermal damping on SST anomalies that is at the exact opposite of what is expected when thermal coupling with the atmosphere is suppressed.

B F O = G(T U a , T clim o ), (8)

Interdecadal variability

In what follows we focus on the stochastic integrations and contrast the variability patterns obtained for the coupled, uncoupled and forced systems. The comparison is made for the canonical horizontal eddy diffusivity value K = 1000 m 2 s -1 .

Figure 2 shows the time-mean meridional overturning streamfunction in the stochastic coupled run, as well as the timeseries of the MOC index over a 1000 years period for the coupled, uncoupled and forced systems. The MOC index is defined as the maximum value of the meridional overturning streamfunction north of 30 • N and below 850m. Despite the relative simplicity of the model, in terms of both geometry and physics, the observed North Atlantic meridional overturning circulation is well reproduced with a thermally direct cell intensified in the Northern Hemisphere with a peak value close to that observed (18 Sv, Colin de Verdière et al. 2019), and a weaker (5 Sv) thermally-indirect cell at depth reminiscent of the Antarctic Bottom Cell. Quasi-regular oscillations with interdecadal periods (∼25 years) are found in all cases, with peak-to-peak amplitudes reaching 38% (56%) of the mean transport in the coupled (forced) case, but only 12% in the uncoupled system. Thermal coupling with the atmosphere, as inferred from the comparison between the coupled and uncoupled systems, therefore acts as a strong amplifier of interdecadal oscillations in agreement with [START_REF] Farneti | Mechanisms of interdecadal climate variability and the role of ocean-atmosphere coupling[END_REF] and [START_REF] Wu | North Atlantic decadal variability: air-sea coupling, oceanic memory, and potential Northern Hemisphere resonance[END_REF]. As originally stated by [START_REF] Barsugli | The basic effects of atmosphere-ocean thermal coupling on midlatitude variability[END_REF], coupling with the atmosphere makes SST and SAT anomalies to vary more in phase in the coupled case, thereby reducing the damping on SST anomalies compared to the uncoupled case. In the forced case the damping of SST anomalies is at its minimum (but non-zero because of the presence of mixing processes) leading to much greater variability.

This feature can be seen in the standard deviations of SST anomalies in Fig. 3. For all cases SST variance is maximum in the northwest corner of the basin and is about three times greater in the coupled compared to the uncoupled case. There is virtually no variability in the Southern Hemisphere: the presence of the southern channel induces a relatively weak upper circulation that does not favor the baroclinic instability mechanism [START_REF] Arzel | On the origin of interdecadal oscillations in a coupled ocean-atmosphere model[END_REF]). The spectral peaks (Fig. ). In the uncoupled case, the NAO signature along with its dipolar structure is clearly evident in the first EOF of SST anomalies (Fig. 5). Therefore thermal coupling does not only increase the amplitude of the mode but also significantly alter the leading pattern of the variability, at least for this specific value of K = 1000 m 2 s -1 . The second EOF of the uncoupled case exhibits strong similarities with the leading EOF of the coupled and forced cases, suggesting an implication of internal ocean dynamics in the uncoupled variability.

Sensitivity to K

Figure 6 presents the sensitivity of some key quantities to the horizontal diffusivity K in the coupled, uncoupled and forced systems, and under both deterministic and stochastic forcing conditions. We first see that the strength of the mean overturning increases with K in all cases. This behaviour sharply contrasts with results obtained with models where the eddy-induced diffusivity rather the horizontal diffusivity is varied [START_REF] Marshall | The dependence of the ocean's MOC on mesoscale eddy diffusivities: A model study[END_REF][START_REF] Arzel | Contributions of atmospheric stochastic forcing and intrinsic ocean modes to North Atlantic ocean interdecadal variability[END_REF]. The reason for this discrepancy may be caused by the "Veronis effect" whereby horizontal diffusion induces large diapycnal fluxes once the isopycnals are tilted by coastal upwellings, in particular along the western boundary [START_REF] Veronis | The role of models in tracer studies, chap[END_REF]Huck et al. 1999b). Quite remarkably the strength of the mean overturning for a given value of K keeps nearly the same value whatever the configuration used (coupled, forced, uncoupled) and irrespective of the stochastic forcing being present or not. Hence rectification of the mean flow by stochastic forcing does not occur in our simulations as opposed to [START_REF] Frankcombe | Coherent multidecadal variability in North Atlantic sea level[END_REF]. This implies that the different variability characteristics associated with the coupled, uncoupled and forced systems, in particular when stochastic forcing is present, cannot be ascribed to modifications of the mean circulation.

The transition from a steady circulation to an interdecadal oscillatory behaviour was shown to occur through a Hopf bifurcation as the horizontal diffusivity (Colin de Verdière and Huck 1999) or ocean-atmosphere heat exchange coefficient [START_REF] Chen | Interdecadal variability in a hybrid coupled ocean-atmosphere model[END_REF] is decreased. Here Hopf bifurcations under deterministic conditions are successively found at K = 210, 1140 and 1560 m 2 s -1 for the uncoupled, coupled and forced cases respectively (note that only the last two bifurcations are shown in Fig. 6). This sequence of bifurcations is consistent with the magnitude of atmospheric damping acting on SST anomalies which is at its largest in the uncoupled case and at its lowest in the forced case, as discussed in Section 2b. The suppression of variability for large K is not caused by changes in the circulation, which favors increasingly baroclinically unstable situations as K increases (Fig. 6a), but instead by the direct damping effect of dissipation on perturbations. When stochastic forcing is present interdecadal variability exists over the full range of diffusivities explored here, but only as long as the noise forcing has some temporal coherence, in agrement with [START_REF] Frankcombe | Coherent multidecadal variability in North Atlantic sea level[END_REF] and [START_REF] Arzel | Contributions of atmospheric stochastic forcing and intrinsic ocean modes to North Atlantic ocean interdecadal variability[END_REF]. Additional spectral analyses of the MOC and SST (averaged in the region 70 -50 • W, 46 -74 • N) indices show that the interdecadal variability is always statistically significant in both fields, even for the most diffusive cases. The effect of stochastic forcing on the variability is large in the damped regime (for diffusivity values larger than the critical value at bifurcation) and much weaker in the super-critical regime (for diffusivity values weaker than the critical value at bifurcation). The comparison of coupled and uncoupled configurations of the model clearly shows that thermal coupling between the ocean and atmosphere increases the temperature variance in both fluids, in agreement with [START_REF] Barsugli | The basic effects of atmosphere-ocean thermal coupling on midlatitude variability[END_REF], as well as the amplitude of ocean circulation changes. The novel aspect is that the amplifying effect of thermal coupling on the variability is much stronger in the super-critical than in the damped regime. For SST for instance, thermal coupling typically quadruples the amplitude of the variability in the super-critical regime (from 0.55 to 2.31 • C for K = 500 m 2 s -1 ) but less than doubles it in the damped regime (from 0.28 to 0.45 • C for K = 2600 m 2 s -1 ). Let us finally mention that the amplitude of the changes in surface heat flux (diagnosed from a composite analysis identical to that used for SST and SAT in Fig. 6) in the region of maximum SST variance are relatively independent of K and amount to about 25 W m -2 (not shown), an amplitude comparable to the observed annual mean changes in turbulent fluxes associated with the NAO (e.g. Fig. 1 in Arzel and Huck 2020).

Energy sources

We now turn to the analysis of the effect of thermal coupling on the leading energy sources of the variability. These energy sources are identified as the major terms driving the growth of buoyancy variance in the ocean against all sources of damping. As such the buoyancy variance budget provides a quantitative estimate of the contribution of the ocean and atmosphere to the growth of perturbations. Such an approach has been previously and successfully applied to the interdecadal climate variability problem in either stand-alone ocean models (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]Buckley et al. 2012;[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.

a. Method

The linearized buoyancy variance equation is given by 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 , (9)
where the overbar denotes a time average and the prime the perturbation. Integrated over the domain and weighted appropriately, it provides an equation for the available potential energy in the quasi-geostrophic approximation [START_REF] Huang | Mixing and available potential energy in a Boussinesq ocean[END_REF]. 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. The first term on the right hand side 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). Baroclinic instability is associated with a conversion of potential to kinetic energy of perturbations through the positive exchange term w b . Under such unstable conditions, the second term in (9) 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 the basin 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

Q b = g 0 α T Q /C O (

b. Results

The buoyancy variance budget is applied to the northwestern corner of the basin (70-50 • W, 46-74 • N) which roughly coincides with the region of maximum temperature (buoyancy) changes in all experiments (Fig. 3). The buoyancy changes are largely controlled by temperature variations so that the buoyancy variance budget presented hereafter is to a large extent identical to the temperature variance budget (not shown). We define the quantities

S A =< b Q b >, S O = -< u h b .∇ h b >
where the brackets denote volume averaging. The analysis entirely focuses on the budget within the uppermost 250 m of the ocean model. Because our region of averaging encompasses the bulk of the variability the advection of buoyancy variance into or out of the region by the mean flow is always an order of magnitude smaller than S O (not shown).

Figure 7ab presents the sensitivity of S O and S A to K for the coupled, uncoupled, and forced systems, and under both deterministic and stochastic forcing conditions. We first see that S O is positive for all model configurations and for all values of K. Consistent with the amplitude of the variability, the largest values of S O are obtained for the forced system, and the lowest values for the uncoupled one. The oceanic energy source term in the coupled case is about an order of magnitude larger than that in the uncoupled one. The effect of thermal coupling is therefore to significantly increase the creation of buoyancy variance by internal ocean dynamics for all values of K, with the strongest impact in the super-critical regime (Fig. 7c).

Focusing now on S A in the coupled case, we see that this term is strongly negative in the supercritical regime and even beyond (up to K = 1400 m 2 s -1 , Fig. 7b). In this range of K values, surface heat flux anomalies therefore act to damp the variability. For larger K values, S A becomes positive, which along with S O contributes to the production of buoyancy variance in the northwestern corner of the domain. The change in the sign of S A with K can be rationalized as follows. Density anomalies are largely controlled by temperature changes so that b ∝ T . Using this result, and assuming that the upwelling and downwelling longwave radiative fluxes nearly cancel out, the covariance term between the surface buoyancy flux anomalies and the buoyancy anomalies can be approximated as S A ∝< λ(T T a -T 2 ) >, where λ is the air-sea heat exchange coefficient (recall that both the solar and freshwater forcings are steady). Now in the super-critical regime, SST anomalies are relatively large implying that the covariance term S A is dominated by -λ < T 2 > which is always negative. Larger negative values of S A are obtained in the coupled case in the super-critical regime because of the much greater oceanic temperature variance T 2 compared to the uncoupled case. In the damped regime, SST anomalies have a much weaker magnitude and the leading term in S A is λ < T T a > which is positive since both SST and SAT anomalies tend to vary in phase. In the uncoupled system, the same behaviour occurs but the change in S A occurs for lower values of K, consistent with the lower SST variance compared to the coupled case.

Interestingly the change in the sign of S A in the uncoupled case occurs near the Hopf bifurcation of the coupled system. Whether this occurs by pure coincidence or not remains to be found.

Consistent with the larger SST variance in the coupled system compared to the uncoupled one, thermal coupling acts to increase the sink of buoyancy variance associated with surface fluxes across a large range of K values, up to 1700 m 2 s -1 beyond the bifurcation at 1140 m 2 s -1 , but increase the production of buoyancy variance for the most diffusive cases (Fig. 7c).

In the forced case surface flux anomalies are always a source of buoyancy variance (Fig. 7).

We also note that S A in this case is nearly independent on K despite the significant changes in the amplitude of the variability between the super-critical and damped regimes. An explanation for this behaviour was provided by [START_REF] Arzel | Contributions of atmospheric stochastic forcing and intrinsic ocean modes to North Atlantic ocean interdecadal variability[END_REF]. In the damped regime, the kinetic energy variability is much weaker than in the super-critical regime. 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. The increase of the correlation with K compensates for the decrease in the SST variance leading to an almost unchanged covariance term S A as K increases. Compared to the forced system, the coupled system decreases the production of buoyancy variance associated with both the surface fluxes and internal ocean dynamics.

Which of the ocean or atmosphere represents the major energy source for the variability? How does the contribution of the atmosphere and ocean to the production of temperature variance depend on thermal coupling and oceanic eddy diffusivities? To answer these questions we now tends beyond the super-critical regime to a much larger range of K values when thermal coupling with the atmosphere is enabled. Of course stochastic forcing is required to generate variability in the damped regime, but the growth of upper ocean heat content (250 m) anomalies is however mostly explained by changing ocean currents rather than by the direct thermodynamical response of the mixed layer to the noise forcing. Importantly this behaviour occurs despite the fact that the leading pattern of the variability in the damped regime bears some resemblance with the imposed NAO forcing (Fig. 8). The imprint of the internal ocean mode in this damped regime is clearly seen in the second EOF of SST anomalies (Fig. 8) which appears to be similar to the leading pattern obtained in the super-critical regime (Fig. 5). We further note that the pattern correlation between the coupled and uncoupled systems in the damped regime is high (Fig. 8), a behaviour that sharply contrasts with that obtained in the super-critical regime (Fig. 5). In the forced configuration, internal ocean dynamics is the primary driver of the upper ocean heat content (0-250 m)

variability for all values of K (Fig. 7d).

The results obtained here have been obtained using a very specific timescale of noise coherence of 10 days. Additional experiments not presented here reveal that increasing this timescale to 30 days has the effect of increasing the production of buoyancy variance by both surface fluxes (S A ) and internal ocean dynamics (S O ), the effect being much more pronounced in the damped regime. As a whole the increase in S A is larger than the increase in S O implying that a more coherent noise tends to increase (decrease) the role of the surface fluxes (internal ocean dynamics) in the variability. However the ocean is still the dominant contributor to the growth of buoyancy variance in the northwestern corner of the basin and in the upper 250m for diffusivity values up to 2200 m 2 s -1 (not shown).

Figure 9 illustrates finally how the contribution of the ocean and atmosphere to the temperature variability of the coupled system changes as the averaging depth increases from 50 m (the depth of the forcing layer) to 1150 m (roughly the thermocline depth), the 250 m averaging depth chosen previously being intermediate. Quite clearly the greater the averaging depth the stronger the contribution of internal ocean dynamics to the variability, irrespective of the values of K. For instance changes in surface heat fluxes represent the leading energy source of the SST variability in the strongly damped regime, for K > 1700 m 2 s -1 . Upper ocean heat content variability (0-1150 m) in this range is however mostly constrained by the production of temperature variance associated with internal ocean dynamics. The reason for this behaviour is that the S A term is only nonzero in the forcing layer (upper 50 m) whereas S O monotonically decreases with depth, con-sistent with the result that the temperature anomalies mostly project onto the first baroclinic mode whose vertical extent is O(1000)m. Performing the variance budget over a greater depth has the effect of decreasing S O but at a much lower rate than S A , implying a greater role of internal ocean dynamics in the variability.

A coupled ocean-atmosphere box model

To better understand the way thermal coupling reinforces the internal oceanic contribution to the variability we propose herein a conceptual two degree of freedom dynamical system first introduced by Colin de Verdière and [START_REF] De Verdière | A 2 degree of freedom dynamical system for interdecadal oscillations of the ocean-atmosphere[END_REF] to study interdecadal oscillations of the oceanic overturning under deterministic conditions. The strength of this conceptual model resides in the fact that it captures the behaviour obtained in GCMs, and in particular the Hopf bifurcation at the transition between the steady and oscillatory regimes [START_REF] Arzel | The internal generation of the Atlantic Ocean interdecadal variability[END_REF]). This model is modified here to include atmospheric stochastic forcing (Fig. 10).

a. Model description

A detailed derivation of the model from first principles can be found in the Appendix. We simply provide here a brief description of the model and governing equations before going on to its analysis in the following sections. The model represents advection, dissipation and exchange of heat within and between the ocean and atmosphere in a single-hemisphere configuration (Fig. 10).

The meridional redistribution of heat by advective processes between the tropical and polar boxes in the ocean is assumed to be entirely accomplished by the meridional overturning circulation. The oscillations we are studying here arise under constant surface wind-stress forcing. As such, heat transport by the gyre circulation does not play a fundamental role in the variability and will therefore be discarded, as opposed for instance to [START_REF] Marshall | North atlantic climate variability: phenomena, impacts and mechanisms[END_REF]. The overturning circulation is assumed to vary in quadrature with the anomalous temperature contrast, in agreement with results from previous studies (Huck et al. 1999a;[START_REF] Te Raa | Instability of the thermohaline circulation on interdecadal timescales[END_REF]. This feature is also clearly apparent here with the rate of change in the MOC being highly correlated (r = 0.75) with the anomalous SST contrast between the subtropics (20 • N-48 • N) and the subpolar area (48 • N-74 • N) in all model configurations, irrespective of the oceanic state being in the damped or super-critical regime. The time lag between meridional circulation anomalies and changes in the meridional

Q AO Q AO Q S Q L Q L Q S ψ F F eddies eddies OCEAN ATMOSPHERE equator pole
ψ = κx + µψ -γψ 3 , (10a)
where x is the perturbation oceanic temperature difference between the tropical and polar boxes and ψ the perturbation oceanic meridional overturning circulation. There is no evolution equation of the atmospheric temperature since on the long, decadal, timescales of interest here, the atmosphere is in instantaneous equilibrium with its fluxes. Under this approximation, the atmospheric temperature becomes diagnostic and the atmosphere enters the problem only through the parameters δ, β and N 0 . More specifically, δ sums up the oceanic and atmospheric contributions to oceanic thermal damping (oceanic and atmospheric eddy diffusivity, turbulent air-sea heat exchange, longwave radiative feedback) and β represents the amplitude of the atmospheric stochastic forcing N 0 seen by the ocean. The last cubic term has been introduced to stabilize the system at large amplitudes. We choose γ = 300 in all experiments. Following the previous approach, three different configurations of the box model impacting the coefficients δ and β are considered. The first one is the standard case where ocean-atmosphere thermal coupling is enabled, with coefficients

δ = α(2K O + λ(1 -σ)) and β = ασσ N /x (11)
with α = τ/C O and σ = λ/(λ + B + 2K A ). Here τ = 1 year is the scale for time, C O = 4 × 10 9 J K -1 m -2 is the heat capacity of the thermocline whose depth is fixed to 1000 m, x = 20 • C is the mean meridional temperature contrast between the tropical and polar boxes,

λ = 30 W m -2 K -1 is the air-sea heat exchange coefficent, K O = 1 W m -2 K -1 (corresponding to 1000 m 2 s -1 ), K A = 2.275 W m -2 K -1
is the atmospheric eddy diffusivity (corresponding to 1.3×10 6 m 2 s -1 ), and B = 1.7 W m -2 K -1 is the longwave feedback at the top of the atmophere. The second case is the forced ocean-only case where value of the baroclinic growth rate µ the coupled system is integrated for 10,000 years using a fourth order Runge-Kutta scheme and the statistics are computed over the second half of the integration. The integration is done for the coupled, uncoupled and forced cases under both deterministic and stochastic conditions. Hopf bifurcations occur successively at µ = 0.016, 0.057 and 0.25 for the forced, coupled and uncoupled cases. Over the range of values of µ considered here, the deterministic uncoupled system is always stable.

δ = 2αK O and β = ασσ N /x (
The third case is the uncoupled one where the feedback of the atmosphere on oceanic temperature anomalies is suppressed where

δ = α(2K O + λ) and β = ασσ N /x (13)
b. Bifurcations, oscillation period and power spectra

Using the values of δ for the coupled, uncoupled and forced cases, we see that as µ increases, oscillations first emerge in the forced case (µ C = 0.016), then in the coupled case (µ c = 0.057), and lastly in the uncoupled case (µ C = 0.25). The transition from the non-oscillatory to the oscillatory regime occurs through a genuine Hopf bifurcation (see Appendix for details). This sequence of bifurcations is presented in Fig. 11 for the amplitude of the changes in the meridional temperature contrast, as obtained from a 10,000 years numerical integration of the box model. This sequence of bifurcations is consistent with the results from the 3D model (Fig. 6), demonstrating the relevance of this box model to study the physics of interdecadal oscillations of the overturning circulation and the effect of ocean-atmosphere thermal coupling as we shall now see. Period (year) Fig. 12. Theoretical power spectra (A8) of temperature anomalies for µ = 0 and for the coupled, uncoupled and forced cases (solid lines). When circulation changes are inhibited (ψ = 0) the canonical red noise response is obtained (dashed lines) with a temperature spectra given by |

X(ν)| 2 = β 2 | N(ν)| 2 )/(4π 2 ν 2 + µ 2 c ).
All the spectra are normalized by the peak value of the forced case.

Figure 12 shows the theoretical power spectra (A8) for the meridional temperature contrast at µ = 0 for the coupled, uncoupled and forced cases and for cases where the circulation is free to interact with the temperature field (solid lines) and prescribed to its background value (dashed lines). If the circulation is not allowed to vary, technically by setting ψ = 0 in (10a), the canonical response of [START_REF] Hasselmann | Stochastic climate models. Part I: theory[END_REF] is obtained: a red spectrum at high frequencies that levels out when ω < µ c . Circulation changes are thus essential in producing a preferred timescale in the system, consistent with realistic ocean-only simulations of the MITgcm [START_REF] Arzel | Contributions of atmospheric stochastic forcing and intrinsic ocean modes to North Atlantic ocean interdecadal variability[END_REF]. The effect of thermal coupling is to reduce the damping of oceanic temperature anomalies compared to the uncoupled case. There is therefore more variability in the coupled case compared to the uncoupled case. This result is independent of the circulation being prescribed or not. The forced response exhibits the highest variance since the damping of temperature anomalies is even weaker in this case, with an effective damping timescale τδ -1 of 63 years (compared to 18 and 4 years in the coupled and uncoupled cases respectively). These results are in good agreement with those deduced from the 3D model (Fig. 4).

c. Energy sources

Following previous studies and the present approach, insight into the physical mechanisms driving the variability can be obtained by multiplying the temperature equation by x and averaging over a long time period (denoted by an overbar),

1 2 ẋ 2 = -2ψ x + βx N 0 -δx 2 (14)
Similar to (9), the growth of temperature variance can either originate from internal ocean dynamics (-2ψ x ) or from the direct effect of stochastic forcing on the oceanic temperature (βx N 0 ).

The last term is always negative and represents a sink of temperature variance. The results for the coupled and uncoupled cases are presented in Fig. 13. As can be seen the box model reproduces with very high fidelity the results obtained with the full 3D model (Fig. 7). The major points are that thermal coupling (1) significantly increases the production of temperature variance by internal oceanic processes, (2) significantly increases the destruction of temperature variance by air-sea heat fluxes, but in the most weakly unstable cases where the effect of coupling is to reinforce the production of temperature variance by air-sea heat fluxes, (3) and as a whole significantly reinforces the role of internal ocean dynamics in the variability, in particular in the damped regime.

Summary and discussion

The impact of ocean-atmosphere thermal coupling on the primary drivers of interdecadal variability has been investigated using an idealized coupled model. The main advantage of such a model is its ability to explore a greater parameter regime than that would be possible with a comprehensive coupled GCM. The effect of thermal coupling on the primary drivers of the variability was quantified across a wide range of oceanic eddy diffusivities from a comparison of the buoyancy variance budget of coupled and uncoupled integrations.

The results first confirm that the primary effect of thermal coupling is to reduce the internal damping of temperature anomalies due to surface heat fluxes [START_REF] Barsugli | The basic effects of atmosphere-ocean thermal coupling on midlatitude variability[END_REF]. This reduced damping leads to a greater thermal variance in both the ocean and atmosphere as well as greater ocean circulation changes in coupled integrations compared to uncoupled ones, in agree-29 ment with [START_REF] Wu | North Atlantic decadal variability: air-sea coupling, oceanic memory, and potential Northern Hemisphere resonance[END_REF]. One new aspect is that the amplifying effect of thermal coupling on the variability is much stronger in the super-critical than in the damped regime. The presence of an internal ocean mode in the super-critical regime therefore provides a significant positive feedback on the amplitude of the variability through reduced thermal damping by surface heat fluxes.

Thermal coupling is also shown to significantly alter the spatial pattern of the variability, in particular in the super-critical regime. When the interdecadal internal ocean mode is damped however, the pattern correlation between the coupled and uncoupled leading EOFs is large. The analysis presented by [START_REF] Wu | North Atlantic decadal variability: air-sea coupling, oceanic memory, and potential Northern Hemisphere resonance[END_REF] falls within this latter regime suggesting that the internal ocean mode in their model is damped.

The buoyancy variance budget then reveals that internal ocean dynamics always acts as a source of temperature variance, irrespective of the coupling with the atmosphere being present or not or the oceanic state being in the super-critical or damped regimes. On the contrary changes in surface buoyancy fluxes always act to damp the variability in the super-critical regime. Only for the most diffusive, weakly unstable states, does the atmosphere act as a source of energy for the variability through the stochastic forcing exceeding surface heat flux damping. These behaviours are present in both the coupled and uncoupled configurations of the model, but with different magnitudes. Overall thermal coupling with the atmosphere is shown to significantly reinforces the role of internal ocean dynamics in the interdecadal variability of the coupled oceanatmosphere system, in particular in the damped regime. This result constitutes the major finding of our study. Specifically it is shown that the leading role of internal ocean dynamics in the variability extends beyond the super-critical regime to a much larger range of diffusivity values when thermal coupling with the atmosphere is enabled. Put another way, thermal coupling with the atmosphere significantly widens the range of diffusivity values over which internal ocean dynamics drives the variability. Of course stochastic forcing is required to generate variability in the damped regime.

The growth of upper ocean heat content anomalies in this regime is however mostly constrained by changing ocean currents associated with the noise excitation of the internal ocean mode rather than by the direct thermodynamical response of the mixed layer to the noise forcing. Importantly this behaviour occurs despite the fact that the leading pattern of the variability in the damped regime bears some resemblance with the imposed NAO forcing. The same conclusion applies to the underlying mechanism driving interdecadal sea surface temperature variability. These results are in line with the more general statement that the ocean drives the mid-latitude North Atlantic variability on decadal and longer timescales [START_REF] Bjerknes | Atlantic air-sea interaction[END_REF][START_REF] Gulev | North Atlantic ocean control on surface heat flux on multidecadal timescales[END_REF][START_REF] O'reilly | The signature of low-frequency oceanic forcing in the Atlantic Multidecadal Oscillation[END_REF]. Our results are also consistent with [START_REF] Garuba | On the relative roles of the atmosphere and ocean in the Atlantic multidecadal variability[END_REF] who showed on the basis of coupled and partially-coupled integrations of a realistic coupled climate model that the AMV is mostly driven by ocean circulation variability.

A stochastic coupled ocean-atmosphere model is finally proposed that captures the basic effect of ocean-atmosphere thermal coupling, as obtained in the 3D model. The box model builds upon a stochastically-forced atmospheric component coupled to a dynamical ocean. The presence of the latter significantly modifies the response of the coupled system to stochastic forcing, first by introducing an additional and significant source of temperature variance in the system, and second by selecting an interdecadal timescale as opposed to pure thermodynamical models [START_REF] Hasselmann | Stochastic climate models. Part I: theory[END_REF]; [START_REF] Barsugli | The basic effects of atmosphere-ocean thermal coupling on midlatitude variability[END_REF][START_REF] Clement | Response to Comment on "The Atlantic Multidecadal Oscillation without a role for ocean circulation[END_REF]. The fact that the effect of coupling is captured by a simple two-degree of freedom dynamical system suggests that the proposed mechanism does not depend on the details of the 3D model and is therefore a robust feature of the ocean-atmosphere system.

There are of course several fundamental aspects of the real climate system that are missing from our idealized framework. Perhaps the most critical is the absence of a dynamical atmospheric component. With a full dynamical atmosphere, NAO related stochastic forcing would adjust its amplitude and spatial pattern to the changing ocean surface conditions and feed back onto the SST field such that coupled air-sea modes or resonant behaviours could emerge [START_REF] Weaver | On the variability of the thermohaline circulation in the GFDL coupled model[END_REF][START_REF] Wu | North Atlantic decadal variability: air-sea coupling, oceanic memory, and potential Northern Hemisphere resonance[END_REF][START_REF] Ortega | Reconciling two alternative mechanisms behind bi-decadal variability in the North Atlantic[END_REF]. The absence of a seasonal cycle and the limited poleward extent of the ocean domain does not allow the emergence of sea ice although a simple one-layer thermodynamic sea ice component identical to that used by Colin de Verdière and te Raa (2010) is present. The variability studied here is typically enhanced at mid-latitudes along the North Atlantic Current [START_REF] Ortega | Reconciling two alternative mechanisms behind bi-decadal variability in the North Atlantic[END_REF][START_REF] Arzel | The internal generation of the Atlantic Ocean interdecadal variability[END_REF][START_REF] Gastineau | North Atlantic Ocean internal decadal variability: role of the mean state and ocean-atmosphere coupling[END_REF][START_REF] Arzel | Contributions of atmospheric stochastic forcing and intrinsic ocean modes to North Atlantic ocean interdecadal variability[END_REF][START_REF] Gastineau | North Atlantic Ocean internal decadal variability: role of the mean state and ocean-atmosphere coupling[END_REF]) so that direct sea ice effects are unimportant. It would however be interesting to assess the effect of thermal coupling on the variability of coupled climate models presenting a source of variability related to Arctic-Atlantic interactions [START_REF] Jungclaus | Arctic-North Atlantic interactions and multidecadal variability of the meridional overturning circulation[END_REF][START_REF] Escudier | A 20-year coupled ocean-sea ice-atmosphere variability mode in the north atlantic in an aogcm[END_REF]. The lack of oceanic turbulence is another limitation of our model. With resolved mesoscale eddies the explicit horizontal diffusivity used here becomes irrelevant and another critical parameter, such as vertical diffusivity used by [START_REF] Farneti | Mechanisms of interdecadal climate variability and the role of ocean-atmosphere coupling[END_REF] for instance, would be required to explore the physics of the variability in both the damped and super-critical regimes. The presence of oceanic mesoscale eddies does not modify the generic mechanism of baroclinic instability captured in our simple model [START_REF] Huck | Multidecadal variability of the overturning circulation in presence of eddy turbulence[END_REF], but simply acts as a sink of temperature variance at interdecadal timescales [START_REF] Hochet | Direct temporal cascade of temperature variance in eddy-permitting simulations of multidecadal variability[END_REF][START_REF] Hochet | Energy transfers between multidecadal and turbulent variability[END_REF]. As such the leading role of internal ocean dynamics in the variability advocated in the present study is probably overestimated compared to a situation where mesoscale eddies are resolved. Additional studies based on strongly eddying quasi-geostrophic coupled models point to a strong coupling between the ocean and atmosphere on decadal and longer timescales [START_REF] Kravtsov | A highly nonlinear coupled mode of decadal variability in a mid-latitude ocean-atmosphere model[END_REF][START_REF] Martin | Drivers of atmospheric and oceanic surface temperature variance: a frequency domain approach[END_REF]). These quasi-geostrophic models cannot capture the large-scale baroclinic instability mechanism that is at play in our model, because it develops in regions of steep isopycnal slopes. Future work based on primitive ocean models at eddying resolutions will need to be pursued to assess which mechanism among the few listed here is the most relevant for driving temperature variance from mid to subpolar latitudes as well as the role of thermal coupling which is still debated (e.g. [START_REF] Weaver | On the variability of the thermohaline circulation in the GFDL coupled model[END_REF][START_REF] Timmermann | Northern hemispheric interdecadal variability: a coupled air-sea mode[END_REF][START_REF] Delworth | Multidecadal thermohaline circulation variability excited by atmospheric surface flux forcing[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]. The generic large-scale baroclinic instability mechanism driving internal oceanic variability in the planetary geostrophic ocean component of the present model (Colin de Verdière and Huck 1999;[START_REF] Huck | On the robustness of the interdecadal modes of the thermohaline circulation[END_REF] has been shown to be robust in comprehensive ocean climate models [START_REF] Sévellec | The leading, interdecadal eigenmode of the Atlantic meridional overturning circulation in a realistic ocean model[END_REF][START_REF] Ortega | Reconciling two alternative mechanisms behind bi-decadal variability in the North Atlantic[END_REF][START_REF] Arzel | The internal generation of the Atlantic Ocean interdecadal variability[END_REF][START_REF] Gastineau | North Atlantic Ocean internal decadal variability: role of the mean state and ocean-atmosphere coupling[END_REF]. We thus feel confident that the mechanisms presented here might play a role in more realistic contexts. As a final note we wish to stress the extreme usefulness of the buoyancy variance budget combined with the coupled/uncoupled modelling approach to identify the primary drivers of the variability and unravel the role of thermal coupling. Applying such a budget in coupled and uncoupled configurations of comprehensive GCMs and realistic climate models will undoubtedly improve our understanding of the mechanisms driving climate variability on interdecadal timescales in those models.

a. Model equations

Radiative imbalance at the top of the atmosphere (Fig. 10) arises due to differences between incoming solar flux Q S (assumed constant) and outgoing longwave radiation Q L which is linearized around the atmospheric temperature T a , i.e. Q L = A + BT a . The term F acting on the atmospheric layer represents the dynamical part of the forcing which is taken to be stochastic. The turbulent surface heat flux Q AO (assumed positive downward) between the ocean and atmosphere is equal to λ(T a -T o ) where λ is the air-sea heat exchange coefficient and T o the oceanic temperature. Using λ = 30 W m -2 K -1 , an atmospheric heat capacity C a = ρ a C pa h a = 7 × 10 6 J m -2 K -1 (assuming an atmosperic scale height h a = 7 km, and standard values for air density ρ a and specific heat capacity of air C pa ), we obtain an atmospheric damping timescale C a /λ of O( 1) day. Since we are interested in timescales of the order of 10 years or longer, it is legitimate to assume that the atmosphere is in instantaneous equilibrium with the ocean and the stochastic forcing [START_REF] Saravanan | Advective ocean-atmosphere interaction: An analytical stochastic model with implications for decadal variability[END_REF]. The heat budget of the atmospheric box i then reads

Q i S -Q i L -Q i ao + K A (T j a -T i a ) + F i = 0 (A1)
where j = 3-i and K A parameterizes the turbulent heat exchange by large-scale atmospheric eddies between the tropical (i = 1) and the polar (i = 2) boxes. The heat budget of the oceanic box i is ocean depth), we obtain a damping timescale τδ -1 of oceanic temperature anomalies of about 18 years. This timescale is much longer than the typical e-folding times of O(1) year of unstable planetary waves computed from OGCM studies [START_REF] Arzel | The internal generation of the Atlantic Ocean interdecadal variability[END_REF] suggesting the potential for spontaneous oscillations developping in ocean-only integrations to survive to coupling with the atmosphere. Notice that the magnitude of the noise forcing seen by the ocean is proportional to σ which decreases with decreasing λ values and increasing K A values. To close the system, an equation for the anomaly in the strength of the circulation is needed. With an attempt to encapsulate the ideas developed in the previous sections, at least two features must be represented. The first one is the baroclinic instability mechanism, which is parameterized here through the use of a linear growth rate µ. The second one is the apparent phase lag that exists between the meridional while the amplitude of temperature variations is given by ( | X(ω)| 2 dω) 1/2 . With κ = 0.0315, a period of 25 years is obtained in the coupled, uncoupled and forced cases, similar to the 3D model.

C O ∂T i ∂t = Q i ao + C O a i ψ(T j -T i ) + K O (T j -T i ), (
It is readily seen that the period in the deterministic case 2π[2κ -(µ + µ c ) 2 /4] -1/2 closely follows that obtained under stochastic forcing as long as µ < µ c 1. At bifurcation µ = µ c the period of the variability predicted by the deterministic solution is exactly recovered by (A10) but the linear assumption breaks down since | X(ω 0 )| 2 becomes infinite. The stochastic linear solutions become invalid when approaching the Hopf bifurcation and are mostly useful in the most diffusive, weakly unstable cases where µ µ c .

Fig. 1 .

 1 Fig. 1. (a) Zonal (solid) and meridional (dashed) surface wind-stress. (b) Imposed freshwater forcing (evaporation minus precipitation), (c) spatial pattern of the stochastic forcing (c.i. 10 W m -2 ) entering the atmospheric thermodynamic balance. The amplitude in the subpolar gyre (100 W m -2) is twice that of the subtropical gyre (50 W m -2 ) similar to[START_REF] Herbaut | Response of a simplified oceanic general circulation model to idealised NAO-like stochastic forcing[END_REF]. The sign change in the stochastic forcing roughly occurs at the intergyre boundary at 47 • N.

  Fig. 2. (top) Time-mean meridional overturning circulation (MOC) in the stochastic coupled integration for K = 1000 m 2 s -1 . (bottom) Time-series of the MOC index in the stochastic experiments. The MOC index is defined as the maximum of the meridional overturning streamfunction north of 30 • N and below 850m. The stochastic experiments are run for 1000 years and start from the end state of the corresponding deterministic solutions. These later solutions are obtained from a 5000 years long integration which was necessary to reach statistical equilibrium.

Fig. 3 .

 3 Fig.3. SST standard deviations diagnosed from the annual mean ouput of SST anomalies over a 100 years period for the stochastic case with K = 1000 m 2 s -1 for the coupled (left), uncoupled (middle) and forced (right) configurations.

4

  ) in both the MOC index and SST anomalies (computed as the weighted average in the region of maximum SST variance, 70-50 • E, 46 • N-74 • N) greatly exceed the red-noise background (but barely for the uncoupled case which just passes the 95% confidence level) indicating that interdecadal oscillations cannot be simply interpreted as an integration of the atmospheric noise by the oceanic mixed layer. The first and second leading patterns of variability have been estimated in terms of Empirical Orthogonal Functions (EOFs) and are presented in Fig.5. The patterns in the coupled and forced systems are very similar but with a greater amplitude for the latter. The first EOF features a monopole intensified in the northwestern corner of the basin, whereas the second EOF features a monopole of oppositve sign centered in the subtropics and with a much weaker amplitude. These two EOFS vary in quadrature (not shown) with the first EOF leading the second one by about 6 years, thereby covering a full cycle in about 25 years. A lagged regression analysis between SST anomalies and the MOC index highlights the emergence and growth of perturbations in the eastward extension of the western boundary current in the northwestern corner of the basin (not shown). The perturbations then propagate westward toward the western boundary where they subduct and quickly dissipate (not shown). This spatio-temporal organization of the variability is similar to several previous idealized ocean-only model studies using fixed flux boundary conditions without noise forcing[START_REF] Greatbatch | An interdecadal oscillation in an idealized ocean basin forced by constant heat flux[END_REF][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]

Fig. 4 .

 4 Fig. 4. Power spectra of (top) the MOC index (Sv 2 years) and (bottom) SST anomalies ( • C 2 years) timeseries averaged in the northwest corner (70-50 • W and north of 46 • N) in the stochastic experiments with K = 1000 m 2 s -1 and for the coupled (blue), uncoupled (red) and forced (black) systems. The calculation uses a multitaper technique with 3 tapers and is based on 1000 years of annual mean model output. The smooth dashed lines indicate the red noise 95% confidence level.

Fig. 5 .

 5 Fig. 5. First (top) and second (bottom) EOFs of annual mean SST anomalies along with the explained variance (%) in the stochastic experiments with K = 1000 m 2 s -1 and for the coupled (left), uncoupled (middle) and forced (right) cases.

Fig. 6 .

 6 Fig. 6. Statistics of key indices as a function of the horizontal eddy diffusivity K under deterministic and stochastic forcing conditions for the coupled, uncoupled and forced cases. (a) Time-mean MOC index. The index is defined as the maximum value of the meridional overturning streamfunction north of 30 • N and below 850m. (b) Amplitude of MOC variations (Sv). (c) Amplitude of SST changes ( • C) averaged in the northwestern corner (70 -50 • W, 46 -74 • N), which corresponds to the region of largest temperature changes in all experiments. (d) Same as (c) but for the SAT ( • C). The amplitude of the variability in (b-c-d) has been estimatedfrom a composite analysis of the last 1000 years of each experiment where the max (min) of the timeseries is computed as the time average of the values greater (smaller) than the mean plus (minus) one standard deviation.Note that only the changes greater than zero are shown. The uncoupled deterministic system does not have variability in the range of diffusivity values shown here. When the stochastic forcing is white (orange dotted line) the amplitude of the variability in the damped regime is greatly reduced compared to a situation where a noise forcing has a finite temporal coherence.

  with g 0 the acceleration of gravity at the sea surface, C O the specific heat capacity of the forcing layer and Q the anomalous surface heat flux) are positively correlated. The dissipation term b D b , which contains contributions from horizontal and vertical mixing processes (including convective mixing), is a sink of buoyancy variance at the basin scale and will not be considered. In what follows, we will restrict our attention to the two only terms in (9) that can take positive values at the basin scale, namely the atmospheric energy source associated with surface buoyancy fluxes anomalies and the oceanic energy source associated with the interaction of transient buoyancy fluxes with mean horizontal buoyancy gradients.

Fig. 7 .

 7 Fig. 7. Buoyancy variance budget in the northwestern corner of the domain (70-50 • W, 46-74 • N) and in the uppermost 250m as a function of the horizontal diffusivity K. Top panels show the oceanic S O = -< u h b .∇ h b > (left) and atmospheric S A =< b Q b > (right) energy source terms under both deterministic and stochastic boundary conditions and for the coupled, uncoupled and forced systems. The uncoupled case under deterministic forcing conditions does not exhibit variability in the range of diffusivities considered here which explains why it is absent. The bottom left panel shows the changes (coupled minus uncoupled) in atmospheric and oceanic source terms caused by coupling with the atmosphere. Note that a log vertical scale is used for both positive and negative values in the upper and bottom left panels. The bottom right panel shows the contribution R = S O /(S O + S +A ) (%) of internal ocean dynamics to the production of buoyancy variance in the northwestern corner of the domain for the coupled (R C ), uncoupled (R U ) and forced (R F ) cases, with S + A = S A if S A > 0 and 0 otherwise (see text for details). In all panels the vertical dashed lines indicate the position of the Hopf bifurcation of the coupled (blue, K = 1140 m 2 s -1 ) and forced (black, K = 1560 m 2 s -1 ) cases.

Fig. 8 .

 8 Fig. 8. Leading patterns of SST variability ( • C) and explained variance (%) in the damped regime shown here for the specific diffusivity value K = 1400 m 2 s -1 . Shown are the first two EOFs and the first EOF in the stochastic coupled and uncoupled systems, respectively.

Fig. 9 .

 9 Fig. 9. Contribution R = S O /(S O + S + A ) (%, see text for details) of internal ocean dynamics to the production of buoyancy variance in the northwestern corner of the domain (70-50 • W, 46-74 • N) as a function of the horizontal diffusivity K. The calculation is done for the stochastic coupled (R C ) case only and in the upper 50, 250 and 1150m. The vertical dashed (dotted) line indicates the corresponding position of the Hopf bifurcation at K = 1140 m 2 s -1 whereas the horizontal dashed line corresponds to the pivotal value where S A = S O .

Fig. 10 .

 10 Fig. 10. Geometry of the single-hemisphere coupled ocean-atmosphere box model. Two oceanic boxes are coupled to two atmospheric boxes in the meridional plane. Radiative fluxes at the top of the atmosphere are the solar forcing Q S and outgoing longwave radiation Q L . The ocean and atmosphere exchange heat through turbulent sensible fluxes Q AO . The term F acting on the atmospheric layer represents the dynamical part of the forcing and is taken to be stochastic. Mixing by turbulent motions between the tropical and polar boxes within the atmosphere and ocean is taken into account. Advective heat exchange between the two oceanic boxes is entirely accomplished by the meridional overturning circulation Ψ.

Fig. 11 .

 11 Fig. 11. Standard deviation of the anomalous meridional temperature contrast in the box model. For each

Fig. 13 .

 13 Fig. 13. Temperature variance budget in the coupled-ocean-atmosphere box model, as obtained from a 10,000 years long numerical integration. Shown are the oceanic (S O = -2x Ψ ) and atmospheric energy source (S A = βx N 0 ) terms for both the coupled and uncoupled systems under both deterministic and stochastic forcing conditions as a function of the growth rate µ, the differences in S O and S A between the coupled and uncoupled systems and the ratio R = S O /(S O + S + A ) measuring the fraction of the total production of temperature variance explained by internal ocean dynamics. In the calculation of R only positive values of S A are retained so that S + A = 0 wherever S A < 0.

  A2)where ψ (m 3 s -1 ) is the overturning strength, C O = ρ o C po h is the heat capacity of the ocean where ρ o , C po and h are the density, specific heat of water and ocean depth respectively, and K O parameterizes the turbulent heat exchange by mesoscale eddies between the two oceanic boxes. Following Colin de Verdière and[START_REF] De Verdière | A 2 degree of freedom dynamical system for interdecadal oscillations of the ocean-atmosphere[END_REF] we concentrate on the variability of the temperature differences between the two boxes. Introducing x = T 1 -T 2 for the ocean and y = T 1 A -T 2 A for the atmosphere and assuming equal areas a i (a i = a = 4 × 10 6 km 2 ) for the two boxes allows us to writeS -By -λ(yx) -2K A y + N = 0 (A3) and ẋ = αλ(yx) -2ψx -2K O αx,(A4)where α = τ/C O and ψ is scaled by ha/τ where τ = 1 year. The terms S and N are the differences of solar flux and noise forcing between the tropical and polar boxes respectively. The noise forcing is written as N = σ N N 0 where σ N is the amplitude of the forcing and N 0 obeys to a first order autoregressive process with standard deviation equal to 1 and decorrelation time τ N = 10 days.Linearizing now the thermodynamic equations (A3 and A4) around a time mean state (denoted by an overbar), substituing (A3) into (A4), and nondimensionalizing the oceanic temperature anomalyx by x the mean meridional temperature contrast, we obtainẋ = -2ψ -δx + βN 0 , (A5)where ψ is the anomalous overturning circulation. It has been shown from 3D numerical experiments(Huck et al. 1999a) that temperature anomalies are mostly constrained by anomalous advection rather than mean flow effects. Hence the term -2ψx is neglected. The coefficients δ and β are given byδ = α(2K O + λ(1 -σ)) and β = ασσ N /x (A6) with σ = λ/(λ + B + 2K A ). Using B = 1.7 W m -2 K -1 , K A = 2.275 W m -2 K -1(equivalent to an atmospheric eddy diffusivity K A × a/C a = 1.3 × 10 6 m 2 s -1 ), K O = 1 W m -2 K -1 (equivalent to an oceanic eddy diffusivity of K O × a/C O = 1000 m 2 s -1 , with C O = 4 × 10 9 J m -2 K -1 ,for a 1000m

Acknowledgments. We thank the three anonymous reviewers for their helpful suggestions. We are grateful to Alain Colin de Verdière who provided helpful criticism on an early version of this manuscript. The authors acknowledge the Pôle de Calcul et de Données Marines (PCDM at Ifremer, Brest) for providing DATARMOR computational resources.

Data availability statement. The data that support the findings of this study are stored on IFRE-MER (Brest, France) servers and can be made available upon reasonable request to the corresponding author.

APPENDIX

The coupled ocean-atmosphere box model

temperature gradient and the MOC. These aspects encourage us to follow Colin de Verdière and [START_REF] De Verdière | A 2 degree of freedom dynamical system for interdecadal oscillations of the ocean-atmosphere[END_REF] and use the same dynamics, that is

where the last cubic term has been introduced to stabilize the system at large amplitudes. The sensitivity of model solutions is studied in terms of the baroclinic growth rate µ.

b. Properties of the solutions

In the deterministic case (β = 0) and in the small amplitude limit, it can easily be shown (Colin de Verdière and Huck 2000) that the eigensolutions are unstable for µ > µ c = δ provided that δ < √ 2κ.

At µ = µ c a super-critical Hopf bifurcation occurs. Physically this means that when the growth rate of pertubations becomes larger than all sources of thermal damping, the oscillations can grow and eventually settle into a finite amplitude limit cycle. In the stochastic regime, the statistics of the variability can be deduced from the temperature spectra in the small amplitude limit,

where N0 (ω) is the Fourier transform of the noise forcing and

For timescales much longer than O(10) days, the stochastic forcing is essentially white and N0 (ω) = 1. In this case an analytical expression for the most energetic timescale T 0 = 2πω -1 0 of the variability can be deduced from (A8) and will be valid as long as it is much longer than the atmospheric spin down timescale, )