This study focuses on the description of an oceanic variant of the Charney baroclinic instability, arising from the joint presence of i) an equatorward buoyancy gradient that extends from the surface into the ocean interior and ii) reduced subsurface stratification, e.g., as produced by wintertime convection or subduction. We analyze forced-dissipative simulations with and without Charney baroclinic instability (C-BCI). In the former C-BCI strengthens near-surface frontal activity with important consequences in terms of turbulent statistics: increased variance of vertical vorticity and velocity, increased vertical turbulent fluxes. Energetic consequences are explored. Despite the atypical enhancement of submesoscale activity in the simulation subjected to C-BCI and contrary to several recent studies the downscale energy flux at submesoscale en route to dissipation remains modest in the flow energetic equilibration. In particular, it is modest vis a vis the global energy input to the system, the eddy kinetic energy input through conversion of available potential energy, and the classical inverse cascade of kinetic energy. Linear stability analysis suggests that the southern flank of the Gulf Stream may be conducive to oceanic Charney baroclinic instability in spring, following mode water formation and upper ocean destratification.

Introduction

There has recently been a great interest in near-surface ocean turbulence properties in the length scale range from a few hundreds of meters up to a few tens of kilometers. At these scales oceanic flows exhibit frontal structures associated with significant deviations from geostrophy. This range of turbulence scales (referred to as sub-mesoscale) is responsible, for transferring energy up to larger scales [START_REF] Sasaki | Impact of oceanic-scale interactions on the seasonal modulation of ocean dynamics by the atmosphere[END_REF], down to small scales (Capet et al. 2008c;[START_REF] Molemaker | Balanced and unbalanced routes to dissipation in an equilibrated eady flow[END_REF], for mediating air-sea interactions (Thomas and Taylor 2010) and perhaps contributing to interior ocean mixing and the maintenance of the ocean thermohaline structure [START_REF] D'asaro | Enhanced turbulence and energy dissipation at ocean fronts[END_REF]. Submesoscale turbulence is also important through its associated vertical fluxes of properties between the near-surface and the interior, for example by contributing to the injection of deep nutrients into the euphotic layer [START_REF] Levy | Impact of sub-mesoscale physics on production and subduction of phytoplankton in an oligotrophic regime[END_REF][START_REF] Levy | Bringing physics to life at the submesoscale[END_REF].

Recent studies suggest that a useful distinction exists between two different forms of upper ocean submesoscale activity resulting from either stirring of surface buoyancy gradients by interior mesoscale eddies [START_REF] Nurser | Eddy-induced mixed layer shallowing and mixed layerthermocline exchange[END_REF]Lapeyre et al. 2006;Thomas and Ferrari 2008) or from the growth of mixed layer eddies and fronts through unstable processes (Fox-Kemper et al. 2008;Fox-Kemper and Ferrari 2008;[START_REF] Mahadevan | Rapid changes in mixed layer stratification driven by submesoscale instabilities and winds[END_REF]Capet et al. 2008a;Mensa et al. 2013;[START_REF] Skyllingstad | Baroclinic frontal instabilities and turbulent mixing in the surface boundary layer. part i: Unforced simulations[END_REF]. The former is a consequence of baroclinic instability (BCI) involving the ocean interior and leading to mesoscale turbulence that stirs the inhomogeneous surface buoyancy field. This yields frontogenesis essentially around mesoscale eddies while other parts of the ocean remain devoid of fronts. The latter arise when laterally heterogeneous boundary layers subjected to sufficient destabilization by the atmosphere undergo mixed layer instability.

Mixed layer eddies exert a strong influence on the evolution of the boundary layer by contributing to its rapid restratification, e.g., following a mixing event. On the other hand they are considered unimportant for the transport (advective or diffusive) of properties between the surface mixed layer and the ocean interior [START_REF] Boccaletti | Mixed layer instabilities and restratification[END_REF]Capet et al. 2008b;[START_REF] Callies | The role of mixed-layer instabilities in submesoscale turbulence[END_REF]. Although typically less numerous the fronts associated with interior mesoscale activity are presumably more important for property exchanges between the surface and interior because their vertical extension is usually much larger than the height of the mixed layer [START_REF] Ascani | Physical and biological controls of nitrate concentrations in the upper subtropical North Pacific Ocean[END_REF].

Schematically, two types of upper ocean frontal intensification may thus be associated with either: interior BCI having a near-surface frontogenetic expression because buoyancy gradients are present at the surface; or mixed layer BCI whose small horizontal length scale can be interpreted as resulting from weak near-surface stratification, hence a deformation radius of a few kilometers [START_REF] Boccaletti | Mixed layer instabilities and restratification[END_REF]. In practice the turbulent dynamics of the upper ocean is more complex though. One obvious reason is that both types of processes described above tend to coexist [START_REF] Callies | The role of mixed-layer instabilities in submesoscale turbulence[END_REF]Capet et al. 2008d). A more original reason exposed in this study is that other baroclinic instability processes can exist that couple the unstable behaviors of the ocean surface and interior. [START_REF] Roullet | Properties of steady geostrophic turbulence with isopycnal outcropping[END_REF] consider a turbulent regime in which the primary instability arises from the joint presence of a surface buoyancy gradient and an interior velocity shear. This is the so-called Charney instability [START_REF] Charney | The dynamics of long waves in a baroclinic westerly current[END_REF]) that is perhaps best understood with the help of the Charney-Stern necessary condition for baroclinic instability [START_REF] Charney | On the stability of internal baroclinic jets in a rotating atmosphere[END_REF]. For a quasi-geostrophic (QG) zonal flow characterized by buoyancy anomaly, background Brunt-Väisälä frequency, and velocity field denoted respectively b(y, z), N 2 and U(y, z), the meridional gradient of QG potential vorticity is:

∂ y Q = β -f ∂ z s b + f 2 N 2 ∂ z U δ 0 - f 2 N 2 ∂ z U δ -H (1)
In (1) f is the Coriolis parameter, β is the planetary vorticity gradient, s b is the slope of buoyancy surfaces (s b = -b y /N 2 ), H is the water depth assumed constant and δ is the Dirac delta function.

Relative vorticity has been neglected for simplicity. Surface and bottom boundary conditions for buoyancy can be reexpressed in terms of infinitely thin sheets of a distinct potential vorticity form [START_REF] Bretherton | Critical layer instability in baroclinic flows[END_REF]Lapeyre and Klein 2006;[START_REF] Roullet | Properties of steady geostrophic turbulence with isopycnal outcropping[END_REF], hence the Dirac delta functions in (1).

A necessary and most frequently sufficient condition for BCI to develop is that ∂ y Q changes sign in the vertical. Phillips BCI arises from a sign change of interior PV first plus second rhs terms).

Eady BCI results from the fact that the third and fourth term (associated with upper and lower interface PV sheets) are of opposite sign [START_REF] Eady | Long waves and cyclone waves[END_REF][START_REF] Molemaker | Balanced and unbalanced routes to dissipation in an equilibrated eady flow[END_REF]. In its original atmospheric version Charney BCI corresponds to the situation where the stretching term is constant and the sign reversal of Q y occurs because βf ∂ z s b and the lower PV sheet contribution are of opposite signs. A variant situation more relevant to the upper ocean will be considered here where the stretching term varies with depth and β is negligible. The oceanic Charney baroclinic instability (hereafter C-BCI) will refer to the unstable mode arising from the surface Dirac contribution and underlying interior PV gradient (≈ -∂ z s b ) having opposite signs.

In the ocean, situations conducive to C-BCI have been associated with negative s b , i.e., with easterly flows as found on the equatorward flank of subtropical gyres (Tulloch et al. 2011;[START_REF] Gill | Energy partition in the large-scale ocean circulation and the production of mid-ocean eddies[END_REF]. There, the stretching term is expected to be negative while surface buoyancy decreasing poleward is associated with a surface PV gradient that is positive. However, surface intensified westerly currents with positive s b can also be subject to C-BCI provided that the isopycnal tilt s b increases toward the surface. These are the conditions analyzed in [START_REF] Roullet | Properties of steady geostrophic turbulence with isopycnal outcropping[END_REF]. Visual inspection of the thermohaline structures shown in Klein et al. (2008a), Spall and Richards (2000), [START_REF] Levy | Impact of sub-mesoscale physics on production and subduction of phytoplankton in an oligotrophic regime[END_REF], [START_REF] Lima | Biological response to frontal dynamics and mesoscale variability in oligotrophic environments: Biological production and community structure[END_REF] and [START_REF] Akitomo | Baroclinic instability and submesoscale eddy formation in weakly stratified oceans under cooling[END_REF] indicates that process studies on submesoscale dynamics make frequent implicit use of flow conditions characterized by ∂ z s b > 0 and hence subjected to C-BCI (again, strictly speaking, the condition is ∂ z s b > β but β will be negligible in practice as we justify below).

In the present study, our main goal is to gain insight into a submesoscale turbulent regime fueled by C-BCI for an eastward flow, and in particular into its energetics. We compare two configurations with the same degree of frontality but that differ in their upper ocean stratification so that one has a mean state conducive to C-BCI while the other does not. A reference is provided by a third setup that has no surface buoyancy gradient and whose turbulence is only due to interior Phillips-like BCI [START_REF] Pedlosky | Geophysical Fluid Dynamics[END_REF]Lapeyre et al. 2006). The three setups are constructed so that their differences are as limited as possible to facilitate their inter comparisons.

The paper is organized as follows. Methods are described in Section 2. Linear analyses of the initial and equilibrated turbulent states are performed in Sec. 3. The simulations for the three regimes and two numerical resolutions (eddy resolving with ∆x = 8 km and submesoscale rich with ∆x = 1 km) are analyzed and compared in terms of turbulent statistics (Sec. 4) and energetics (Sec. 5). The possibility that a finite depth surface mixed layer and its associated mixed layer eddies disrupt C-BCI dynamics is subsequently investigated using one additional numerical simulation (Sec. 6). We then explore the significance of the Charney BCI regime to the real ocean based on in situ observations (Sec. 7). Some concluding remarks and perspectives are offered in the final section.

Methods

a. Setup designs

In this study, the PE simulations for three baroclinic flows (S1,S2,S3) are analyzed and compared at two different numerical resolutions, one typical of eddy resolving models (∆x = 8 km) and one that is submesoscale rich (∆x = 1 km). In all cases the same computational domain is used: a zonal reentrant channel that is 512 km long, 2040 km wide and 4000 m deep with flat bottom.

The Overall, S2 corresponds to a more stratified upper ocean with N 2 ≈ 0.4-0.7 × 10 -4 s -2 in the upper 200 m (compared to N 2 ≈ 0.1-0.2 × 10 -4 s -2 in S1) but note that the meridional buoyancy gradient is the same at all depths in these two setups (Fig. 1c). As we will see below, the differences in upper ocean N 2 lead to distinct instability properties for S1 and S2, with S1 being subjected to C-BCI while S2 is not. Hydrographic observations presented in Sec. 7 suggest that S1 and S2

bear some resemblance with real ocean flow regimes. Meridional buoyancy gradients in S3 are concentrated in the subsurface and vanish at surface. Lapeyre et al. (2006) analyze in details the differences between two simulations that resemble S1 and S3. Our focus here is on the S1-S2

differences, although S3 occasionally provides a useful "no surface front" reference.

b. Numerical framework PE numerical integrations are performed using the Regional Ocean Modeling System (ROMS) configured for a β -plane at 45 o N latitude (f=1.0313×10 -4 s -1 and β = 1.6186×10 -11 m -1 s -1 ).

The version of the code we use is the one referred to as UCLA-ROMS in [START_REF] Shchepetkin | Correction and commentary for "Ocean forecasting in terrain-following coordinates: Formulation and skill assessment of the regional ocean modeling system[END_REF]. For simplicity, only one prognostic state variable (buoyancy-like) is used.

The code uses a third-order upstream biased scheme for advection which provides lateral diffusivity/viscosity (Shchepetkin andMcWilliams 1998, 2005). Baseline simulations use linear bottom friction with a coefficient r d = 1.5 × 10 -3 m s -1 . Because the intensity of bottom drag is a key parameter in the equilibration of baroclinic flows [START_REF] Rivière | Effects of bottom friction on nonlinear equilibration of an oceanic baroclinic jet[END_REF][START_REF] Arbic | Baroclinically unstable geostrophic turbulence in the limits of strong and weak bottom ekman friction: Application to midocean eddies[END_REF] we

will also analyze a sensitivity run to the value of r d . Simulation S1 rd is identical to S1 except that it has r d = 1.5 × 10 -4 m s -1 , on the low end of commonly used linear bottom frictions.

Our main interest is in the statistical properties of their turbulence once a quasi-equilibrated state is reached. Quasi-equilibrium is achieved by restoring zonally averaged and meridionally smoothed buoyancy and zonal velocities toward their respective reference states. Thus, instabilities growing out of the zonal mean flow can evolve without direct interference by the forcing. We choose a restoring time scale equal to 50 days. Meridional smoothing is performed using a Hanning window of width 50 km which ensures that submesoscales are not directly forced [START_REF] Roullet | Properties of steady geostrophic turbulence with isopycnal outcropping[END_REF]. with values equal to respectively 10 -5 and 10 -4 m 2 s -1 . In addition, diffusion is enhanced locally depending on the local Richardson number and occurrence of static instability, as described in [START_REF] Large | Oceanic vertical mixing: A review and a model with a nonlocal boundary layer parameterization[END_REF]. We note this contribution with a Ri subscript. The real ocean is generally overlaid by a well mixed layer in response to atmospheric forcing. In section 6, we test the importance of mixed layer dynamics and its interplay with Charney BCI. To do so, we artificially produce an upper ocean boundary layer in S1 and S2 by prescribing enhanced diffusion κ (b) and viscosity κ (u) above a level z mld = -65 m. The simulations are respectively named S1 ml and S2 ml .

Equilibrated

Overall we thus have:

κ (b,u) = κ (b,u) b + κ (b,u) Ri + κ (b,u) mld with κ (b,u) mld = 10 -2 1 -tanh( z + 5 h trans ) 1 + tanh( z -z mld h trans )
in S1 ml and S2 ml and 0 otherwise. Note that κ (b,u)

mld nearly vanish at the surface as observed in the ocean. For simplicity, the same vertical transition length scales h trans is chosen equal to 10 m at the surface and bottom of the mixed layer.

c. Experimental layout

Analyses are carried out on equilibrated solutions without and with a mixed layer. The simulations are performed as follows. Initial zonal velocities are in geostrophic balance with the zonally invariant initial buoyancy field constructed as described above. Flow destabilization is triggered by a very small random perturbation added to the initial buoyancy field. The solutions are integrated over a four year period. For simplicity a year is assumed to last 360 days. The first three years are discarded as spin up because this is the time period required for eddy kinetic energy (EKE hereafter; kinetic energy will be abbreviated as KE) to reach a plateau. The 36 model outputs stored every 10 days during the fourth year are used for analyses of the quasi-equilibrated solutions. At the end of the fourth year, near-surface vertical mixing is enhanced in S1 and S2 to create a mixed layer and the simulation is pursued during 420 day. Interaction between the Charney baroclinic instability regime and mixed layer dynamics is evaluated over the last 360 days (36 outputs every 10 days) in these simulation referred to as S1 ml and S2 ml .

Simulation S1 rd with reduced bottom drag is restarted from S1 instantaneous state after completion of the 3 year spin up. It is further spun up for 6 months to let turbulence adjust and is analyzed over the subsequent year.

d. Analysis framework

We use and combine spatial and temporal forms of averaging to obtain statistical robustness. • refers to spatial averaging over the entire domain. Zonal invariance is key model symmetry. We will rely on zonal averaging denoted • zo to decompose the flow into mean and turbulent parts:

X(x, y, z) = X zo (y, z) + X (x, y, z). Also, to increase statistical reliability, temporal averaging < • > can be performed over the 36 model outputs that are available for each simulation.

Combining these averaging operators, we define the root mean square (rms) of a variable X X rms =< X 2 > and the standard deviation

X std =< X 2 >
The definition of eddy kinetic energy follows:

EKE = 1/2 (u std ) 2 + (v std ) 2 .
where u and v are the zonal and meridional velocity components. The frontogenetic tendency presented in section 4 is defined as [START_REF] Giordani | Sensitivity of cyclogenesis to sea surface temperature in the northwestern atlantic[END_REF]:

F sh = -(∂ x ρ) 2 ∂ x u + (∂ y ρ) 2 ∂ y v + ∂ x ρ ∂ y ρ (∂ y u + ∂ x v) ,
where potential density ρ is related to buoyancy by b = -g ρ 0 ρ, with ρ 0 a reference density that we choose equal to 1000 kg m -3 .

Our model outputs are 10 days apart which is not long enough to consider them independent realizations. We find that mesoscale energy levels tend to fluctuate on time scales from weeks to a few months (fluctuations on longer scales are also present albeit with reduced amplitude). When computing uncertainty we will assume that our 360 days runs correspond to only 12 independent realizations, hence uncertainty on estimates of < X > will be 1/ √ 12 X std . However we note that our analyses are for equilibrated solutions with very limited drift in energy between start and end states, so that fluxes presented in Sec. 5 closely balance each other.

Linear analysis

Linear stability analysis is performed to characterize the unstable modes of the three setups. The method follows [START_REF] Roullet | Properties of steady geostrophic turbulence with isopycnal outcropping[END_REF] with a linearization of the QG equations about the reference states. Relative vorticity is neglected. The key inputs are the northern and southern density profiles and width of the baroclinic zone. Linear stability of the equilibrated states produced by model integrations and presented in details in the next section are also analyzed. To do so, time-and zonal-mean density profiles 250 km north and south of the jet axis (at y=1000 km) are used.

Results for a shorter across-jet distance are not significantly different. Growth rates are strongly reduced for northern and southern density profiles taken at the extreme ends of the domain but they do not reflect, in this case, the turbulent behavior of the most energetic region near the jet axis.

Results in terms of growth rate and vertical structure of the most unstable modes are shown in Fig. 3. The quantity IQ y = z -H ∂ y Q dz, where ∂ y Q is given in (1), is also shown for the initial/reference and model equilibrated states to connect instability modes with the Charney-Stern theorem. Vertically integrating ∂ y Q allows us to explicitly represent the Dirac contribution associated with the surface PV sheet. Note that sign changes of the meridional PV gradient correspond to extrema for its vertical integral (represented with open circles in Fig. 3). Below 500 m depth all setups have similar IQ y profiles with the extremum around 1000 m being responsible for the low wavenumber unstable mode (e-folding time scale around 20 days, wavelength of the order of 200 km, i.e., ≈ 5-6 deformation radii). In addition, S1 initial state possesses a prominent shorter mode with a most unstable wavelength at about 18km (associated e-folding time scale is 12 days).

Eigenvector's modulus reveal the contrasted nature of the low and high wavenumber modes. The former corresponds to the gravest baroclinic mode, with a unique zero-crossing of the real and complex parts of the eigenvectors about 1000 m depth (not shown). The latter is strongly surface trapped with negligible signature below 400-500 m depth. It results from the presence of an IQ y minimum below the surface and a prominent discontinuity at the surface induced by the Dirac density contribution. This unstable mode is similar in essence to the Charney mode [START_REF] Charney | The dynamics of long waves in a baroclinic westerly current[END_REF].

Note that in the upper 1500 m, IQ y changes with depth due to vortex stretching systematically dominate over the β contribution.

The precise choice of some parameters used in the construction of the initial states can now be more easily justified. In particular, the finite δ b s value chosen for S2 (see appendix and Table 2) corresponds to the minimal amount of surface stratification that needs to be added to S1 buoyancy profiles to make IQ y monotonic above 1000 m depth, thereby suppressing the high wavenumber Charney mode.

Anticipating on Sec. 4 we present the linear instability analysis performed on the mean turbulent states obtained in our numerical simulations, i.e., on states rectified by turbulent fluxes. Results

indicate that flow turbulence affects the low wavenumber interior modes significantly less than the Charney mode of S1 which is strongly modified by turbulent rectification (maximum growth rate reduced by 40 % and significant shift of the corresponding most unstable wavenumber; compare gray and black curves in Fig. 3b). In agreement with the nonlinear results we present below and also the conclusions of [START_REF] Roullet | Properties of steady geostrophic turbulence with isopycnal outcropping[END_REF] we attribute this to the efficiency with which the Charney mode converts available potential energy (APE) into EKE, thereby producing an intense near-surface restratification that reverses the sign of Q y in the upper 70-80 m and regularizes the near-surface behavior of IQ y (compare the near surface structure of IQ y for the basic and model mean states in Fig. 3a). This near-surface sign reversal of Q y means that modes with short enough vertical/horizontal scales are no longer unstable because Q y is one-signed in their depth range of influence. A short wave cut off is indeed present for the S1 nonlinear state (see Fig. 3b) while it is absent for the reference state, as in the original Charney problem [START_REF] Charney | The dynamics of long waves in a baroclinic westerly current[END_REF][START_REF] Pedlosky | Geophysical Fluid Dynamics[END_REF].

Turbulence properties in quasi-equilibrium

The differences between the solutions for the three setups in terms of eddy kinetic energy, vertical velocity rms, vertical eddy flux of buoyancy and frontogenetic tendency are revealed in Figs. 4

and 5 for two resolutions typical of eddy-resolving (∆x = 8 km, 50 vertical levels) and submesoscale rich (∆x = 1 km, 200 vertical levels) simulations.

Mesoscale resolving solutions exhibit limited differences in terms of relative vorticity (Figs. 4 a-c at the surface) and w rms depth distribution (fig. 5d) despite a noticeable signature of upper ocean surface frontogenesis in S1, S2 above ∼ 500 m depth (see Fig. 5h). Effects of frontogenesis on domain-averaged vertical buoyancy flux are also limited although less so than for w rms . Comparison of S1 and S2 reveals that Charney BCI brings only modest differences, most noticeably for < w b > above 200m depth (compare S1 and S2 profiles in Fig. 5f).

Because the fine scales are energized by near-surface frontogenesis present in S1-2 and by high wavenumber Charney instability in S1, the solution most (respectively, least) affected by the res-olution increase is S1 (resp. S3). In S1, a dramatic change with resolution is the emergence of a well marked < w b > maximum at 100 m depth for ∆x = 1 km. At that depth < w b > increases 4-fold compared to ∆x = 8 km. The < w b > near-surface increase with resolution in S2 is more limited (about two-fold). By construction, the basic states of S1 and S2 have identical meridional buoyancy gradients and we have verified that this is also the case for their rectified mean states.

Their eddy kinetic energy is also very similar (Fig. 5a). Therefore near surface frontogenesis resulting from the stirring of mean buoyancy gradients by mesoscale activity has no reason to differ.

The difference in w rms and < w b > thus arises from the expression of the Charney mode in S1

when resolution allows it. Between 0 and 300 m depth, i.e., roughly the ocean layer where the PV meridional gradient is negative in S1 (see Fig. 3a) the Charney mode is thus associated with over 60 % of the vertical buoyancy flux, provided resolution is fine enough to permit the representation of fine scale instabilities (and their nonlinear consequences). As noted in Sec. 3 this flux has important consequences in terms of instability of the mean rectified state.

The visual aspect of the instantaneous relative vorticity and w fields is in agreement with the statistical description. In S2, frontogenesis manifests itself through the intensification of a limited number of fine-scale filamentous structures mainly located in between mesoscale eddies. Thus the typical submesoscale features in S2 are fronts and filaments whose alongfront length scale seems of the order of 100 km or more, i.e., indicative of mesoscale stirring. In contrast, in S1 the central part of the computation domain is more densely populated with submesoscale features: numerous submesoscale eddies with diameters less than 20 km and a marked cyclonic preference [START_REF] Munk | Spirals on the sea[END_REF]; fronts/filaments broken up or undulating with associated length scales of tens rather than hundreds of kilometers. The S1/S2 differences are most evident at the surface but they are still noticeable at 110 m depth, where S1 has its < w b > maximum and the linear Charney mode is still quite influential (Fig. 3b inset). Vertical velocity fields at 110 m also exhibit similar differences, S1 vertical velocities being richer in fine-scale structures.

For S3, increased resolution helps better resolve vorticity filamentation arising from stirring by interior mesoscale activity (compare Figs. 4c andf). However, no surface intensification occurs owing to homogeneous surface buoyancy. Note though the significant enhancement of the deep vertical buoyancy flux maximum when going from 8 to 1 km resolution (compare Fig 5e andf).

This surprising result is further discussed in the final section.

Energetics

In this section we investigate how the solutions energetically equilibrate. Our main focus is the strength of the downscale kinetic energy transfer due to advection. The importance of such downward energy transfers in the overall equilibration of the ocean is debated (Capet et al. 2008c;Klein et al. 2008a;[START_REF] Molemaker | Balanced and unbalanced routes to dissipation in an equilibrated eady flow[END_REF][START_REF] Marchesiello | Submesoscale turbulence in Tropical Instability Waves[END_REF][START_REF] Pouquet | Geophysical turbulence and the duality of the energy flow across scales[END_REF]Barkan et al. 2015). One issue is that the studies with the largest downscale KE transfers are those for idealized flows whose relevance to the ocean are unclear, in particular because they do not have upper ocean intensified stratification. On the other hand, realistic flows tend to have limited or negligible forward cascades of KE. We see S1 as a good candidate to maximize the importance of the downscale KE flux (given the intensity of its submesoscale frontal activity) while remaining close to realistic ocean conditions.

We start by examining the 2D horizontal velocity spectra computed as in Capet et al. (2008c) except that no procedure is applied to periodize the fields because our reentrant channel setup implies periodicity in the x direction by construction while the weakness of the flow near the northern and southern walls makes the issue of periodicity in the y direction irrelevant (tapering in the y direction does not affect our results, not shown). Power spectrum densities for surface horizontal velocities (Fig. 6) are generally consistent with an energization of submesoscale turbulence by both C-BCI, and general stirring of the surface buoyancy field. S1 has the shallowest surface EKE spectrum which roughly follows a -5/3 power law1 (Fig. 6) whereas the spectral slope of S3 is 10 -4 rad m -1 which is roughly consistent with the range of influence of the Charney mode as determined from the linear stability analysis. Nevertheless S1 also exhibits a similar steepening of its KE spectrum as depth increases, even well within the range of influence of the C-BCI mode (e.g., at 110 m). We interpret this tendency as a combination of reduced frontogenetic efficiency with depth (also felt by S2, see Fig. 5h) and reduced energization of the submesoscale flow by C-BCI (see the eigenvector vertical structure shown in Fig. 3b from the linear analysis).

∼ -3 to -3.5,
Underlying these differences in EKE spectrum are major differences in APE to EKE conversion < w b > as already noticed (see Fig. 5c for their overall magnitude). To better rationalize the behavior of EKE as a function of depth we find useful to examine the fraction of APE to EKE conversion C(z) =< w b > that takes place in the submesoscale range. Fig. 7 represents the ratio C sm /C as a function of depth, where

C sm (z) =< k>k sm Re w * b dk > , (2)
the caret is a horizontal 2D Fourier transform, the symbol Re denotes the operator that selects the real part and wavenumber integration is done in k shells (k = |k h | is the horizontal wavenumber).

We choose k sm = 3 × 10 -4 rad m -1 (20 km wavelength) and note that this value corresponds approximately to the maximum unstable wavenumber for the C-BCI mode in S1 (see Fig. 3b).

It is also the lower bound of the wavenumber range where a downscale EKE flux is found in S1 as we will see below. For both S1 and S2 C sm /C is maximum at the surface and decreases rapidly with depth. A factor 2 to 3 reduction is found at 20 m depth in S1 (the reduction is even more abrupt in S2). We interpret this rapid decrease as a consequence of reduced frontogenetic efficiency with depth (also felt by S2, see Fig. 5h) while the depth dependance of the amplitude of the C-BCI linear mode would predict a much gentler reduction with depth (see the eigenvector vertical structure in Fig. 3b).

The near-surface confinement of submesoscale energization in S1 (and even more so in S2) has important implications for the flow dissipation. In the remainder of this section, we present some analyses suggesting that, even in the submesoscale energized S1, downward KE transfers (by advection) are of limited significance in terms of flow equilibration. Spectral EKE transfer are estimated as (Capet et al. 2008c):

Σ Π (k) = -< 0 z=-200m k max k Re u * h •( u h •r)u h dk dz > (3)
where u h is the instantaneous horizontal velocity vector field and temporal averaging over the 36 model outputs is performed. Although the precise value (-200 m) is arbitrary the bottom of the integration depth range is chosen so as to evaluate the strength of the downscale KE flux in a layer that is rich in submesoscale but also sufficiently thick to be of some significance (as opposed to a layer of a few tens of meters, i.e., the typical width of oceanic surface boundary layers).

In agreement with previous studies Σ Π takes both positive and negative values in S1 and S2

(Fig. 8) with a sign change for k 2-3 × 10 -4 rad m -1 which corresponds to wavelengths of the order of 25 km, well resolved at ∆x = 1 km. Positive values are located at high wavenumbers and reflect the EKE leakage from the balanced flow toward smaller scales and dissipation. Negative values at low wavenumbers correspond to the well-known inverse cascade found in rotating stratified flows [START_REF] Charney | Geostrophic turbulence[END_REF]). Our focus is on the forward cascade and some issues complicate the interpretation of this inverse cascade in our simulations (limited statistical robustness and limited domain size in the zonal direction) but we note that the magnitude of the latter is about 40 times larger than that of the former.

For the forward cascade, the largest values for Σ Π are found in S1 (∼ 4 × 10 -2 mW m -2 ) while high wavenumber Σ Π are indistinguishable from 0 in S3 (not shown). In S2 Σ Π is positive over roughly the same range as in S1 but with a peak value 7 times smaller. Useful comparisons can be made between the values of Σ Π in S1 and two other flux terms essential to the energy budget, 1-APE to EKE conversion < w b >, which is the primary source of EKE in the simulations, and 2-the amount of energy F I injected by the restoring terms (see Sec. 2). In practice APE input through buoyancy F b I strongly dominates over KE input through velocity restoring (see table 1).

F b I is computed as in [START_REF] Molemaker | Balanced and unbalanced routes to dissipation in an equilibrated eady flow[END_REF], i.e., we do not make use of the quasi-geostrophic approximation to compute APE.

The main energy route as it emerges from this analysis is as follows (see table 1): the forcing provides large amounts of APE in the upper ocean, most of which (0.72 mW m -2 ) is converted into KE in the submeso and mesoscale range. The forward energy cascade at submesoscale only participates marginally to the redistribution and dissipation of that energy with a flux that is less than 5% of the KE input. A much larger upscale flux takes care of bringing injected energy at horizontal and vertical scales large enough that bottom friction can act on them. Note that the maximum intensity of the upscale flux can be larger than the total conversion term owing to net fluxes of pressure work into the 0-200 m layer 2 . Integrated vertically down to the bottom at z = -h (max) though we have verified thatmin(Σ h (max)

Π

) < w b > (see Table 1). Note also that the intensity of the full depth forward cascade max(Σ h (max)

Π

) is an even smaller fraction of the full depth APE release, of the order of 1 %.

Another interesting remark concerns the respective roles played by the inverse and forward 5b) is 200 × 0.07 = 14 m 3 s -2 . Assuming approximate equipartition of energy we obtain an estimate for mesoscale energy ∼ 28 m 3 s -2 . This leads to a decay time scale for the mesoscale 20 years (in contrast using the maximum intensity of the inverse cascade diagnosed in S1 -1.6 mW m -2 leads to a physically consistent 200 days time scale). This further confirms that the direct KE cascade in S1 is of small magnitude. This is even more true for S2

which has comparable levels of mesoscale energy but a downscale KE flux about 7 times smaller than S1.

To make sure that our conclusion does not overly depend on the relative efficiency of the large scale energy sink we have performed a sensitivity run S1 rd identical to S1 except that the bottom drag coefficient was reduced by a factor 10 (r d = 1.5 × 10 -4 m s -1 ). The simulation is analyzed at equilibrium as detailed in Sec. 2c. It is important to realize that reducing bottom drag has profound consequences on the whole energetics, in ways that may not be consistent with the functioning of the real ocean. More intense eddy activity permitted by the reduction of viscous forces leads to a very large increase in release of APE which, in turn, leads to an increase in APE injection because of the type of full depth restoring we use. In the real ocean no external forcing maintains interior buoyancy gradients. Buoyancy restoring in S1 rd is responsible for an unrealistic 26 mW/m 2 of domain-average APE injection. And domain average EKE levels near the surface approach 2600 cm 2 s -2 (1.5 increase compared to S1), which is on the upper limit of what can be locally found in the real ocean [START_REF] Ducet | A comparison of surface eddy kinetic energy and reynolds stresses in the gulf stream and the kuroshio current systems from merged topex/poseidon and ers-1/2 altimetric data[END_REF]. Even then, the downscale EKE cascade offers a dissipation pathway to just above 6 % of the APE to EKE conversion (see table 1), i.e., a fraction similar to that found for S1. We tend to see this 6 % figure as an upper bound and will further elaborate on this in the conclusion.

Effect of a mixed layer

The simulations analyzed so far have no well-mixed boundary layer. This may be an important limitation. We indeed wonder if a conflict exists between the perturbed motions needed to amplify the Charney mode versus mixed layer instability modes. Essential to our Charney instability is the subsurface layer where Q y is negative (and IQ y decreases, see Fig 3a). If a mixed layer occupies a significant fraction of that layer C-BCI growth and the associated turbulent fluxes may be reduced by the growth of mixed layer eddies.

The numerical methods we use to perform linear instability analysis do not withstand vanishingly small N 2 in the upper ocean. Thus, we explore the interplay between Charney BCI and mixed layer instability through nonlinear simulations. Specifically, we compute solutions for S1

and S2 with a mixed layer artificially created by increasing the vertical diffusion and viscosity coefficients within ∼ 65 m from the surface (see Sec. 2b and c for details). These simulations are called S1 ml and S2 ml . A reduction of < w b > is also noticeable below the mixed layer compared to the simulations without mixed layer but it concerns a limited depth range [100-130] m and is modest in amplitude, especially in S1 ml . In the case of S2 ml interaction between mixed layer and subsurface dynamics leads to negative < w b > in the subsurface. This is because the subsurface destratification tendency that accompanies the mixed layer restratification (Lapeyre et al. 2006) cannot be compensated by other processes (Phillips instability has a deeper center of action). Overall, subsurface vertical buoyancy fluxes remain much stronger in S1 ml than in S2 ml , i.e., the mixed layer eddies do not cancel the C-BCI effects, even for the relatively large mean mixed layer depth we have chosen (65 m).

Charney BCI in the real ocean

We now discuss the plausibility that the C-BCI regime exists in the real ocean (just like Phillipstype instability has been linked with the dynamical regimes in the Southern Ocean, Gulf Stream . . . ). General occurrence of C-BCI in the global ocean is studied in [START_REF] Smith | The geography of linear baroclinic instability in earth's oceans[END_REF] and Tulloch et al.

(2011) but their main focus is on westward flows. Consistent with our numerical experiments, we focus on C-BCI for near-surface intensified eastward flows. For such flows ∂ y b is negative, ∂ z U is positive and so is the surface PV sheet gradient term in (1). Satisfying the Charney-Stern condition for C-BCI implies that the stretching term, -∂ z s b , be negative below the surface, i.e., that the isopycnal slope increases when approaching the surface. One particularly favorable situation is when ∂ z N 2 < 0. This can be readily seen by re-expressing the stretching term as:

-f ∂ z s b = + f ∂ z b y N 2 = - f N 2 ∂ z (-∂ y b) + f 2 N 2 ∂ z U 1 N 2 ∂ z N 2 (4)
Due to the nature of frontogenetic processes, ocean fronts tend to be intensified near the surface so ∂ z (-∂ y b) is generally positive and we have rearranged the first rhs term to reveal its negative sign. ∂ z N 2 < 0 below the surface implies that the second rhs term is also negative. ∂ z N 2 < 0 also causes the 1/N 2 and 1/N 4 factors to increase toward the surface so that stretching is expected to dominate over β in (1). In S1, stratification weakens upward in the upper 300 m (Fig. 1b) and this is what explains the presence of C-BCI (the second term in ( 4) is overwhelmingly responsible for the negative upper ocean Q y , not shown).

In most parts of the ocean, robust ∂ z N 2 < 0 are restricted to a thin "transition layer", typically O(10) m, located between the mixed layer base (above which ∂ z N 2 ∼ 0) and the relatively shallow core of the pycnocline (e.g., see transition layer thickness statistics in Johnston and Rudnick 2009, their Fig. 7a). Because the depth range of influence of C-BCI roughly coincides with the layer where ∂ z N 2 < 0 (and Q y < 0) C-BCI will be unimportant when this layer is very thin.

A place where one expects to find thick layers with ∂ z N 2 < 0 are mode water regions. By definition, these regions are characterized by the presence of pycnostads, i.e., , weak upper ocean N 2 and deep N 2 maxima, hence negative ∂ z N 2 (and Q y ) over the depth range of the pycnostad.

C-BCI can thus be important in such regions. On the other hand, in places and at times where mode waters are actively formed through convection, we expect more energetic instabilities (convective/symmetric/mixed layer instability) to dominate upper ocean dynamics over a depth range comparable to that where ∂ z N 2 is negative. This leads us to hypothesize that C-BCI may impact mode water dynamics when: 1-winter time convection/subduction has strongly destratified the water column and led to ∂ z N 2 < 0 over hundreds of meters; 2-calm weather conditions prevail so the surface mixed layer is relatively shallow, i.e., we are in the type of situation studied numerically in Sec. 6 where C-BCI is important.

We now show partial evidence that the spring time southern flank of the Gulf Stream may offer a particular instance where 1-and 2-are satisfied and C-BCI plays a role in the seasonal restratification of the subsurface ocean. Wintertime destratification is clearly visible in the depth range 0-200 m while the permanent pool of North Atlantic mode water at greater depths resides, down to about 500m depth. This translates in a monotonic increase of the isopycnal slope roughly above that depth as revealed by the IQ y profile in Fig. 10. Above 200 m depth s b and IQ y change rapidly with depth which is due to the second term in (4) (not shown). Combined to the northward surface density gradient this is a priori conducive to C-BCI and the IQ y vertical profile is similar to that of S1 (see Fig. 3a), albeit with reduced excursions away from the planetary IQ y profile.

The maximum growth rate curve differs from that found in Sec. 3 for S1 in that it does not exhibit two well-separated instability modes. Nevertheless, the extended range of instability beyond k= 10 -4 rad m -1 strongly suggest that Charney type BCI instability is regionally important, as the Charney-Stern criterium predicts. The vertical structure of the eigenvector found for k= 5 × 10 -4 rad m -1 (i.e., about 13 km wavelength) is also consistent with Charney BCI upper intensification.

The maximum growth rate, obtained at that wavenumber reaches 0.04 days -1 . It is about half that of the most unstable mode in S1 and also half of the growth rate found for a Charney BCI mode in the Gulf Stream region by [START_REF] Smith | The geography of linear baroclinic instability in earth's oceans[END_REF]. Note that, in contrast to [START_REF] Smith | The geography of linear baroclinic instability in earth's oceans[END_REF], we have restricted ourselves to the zonal mean part of the flow and our growth rate has a regional significance. A more local analysis with y 1 and y 2 chosen closer to each other dot not qualitatively change our conclusions but leads to significantly higher maximum growth rates comparable to or stronger than those in [START_REF] Smith | The geography of linear baroclinic instability in earth's oceans[END_REF].

Stratification evolution between April and July is confined in the upper 200 m where it yields a dramatic reduction of the jump between interior and surface potential vorticity (compare left and right columns in Fig. 10). As a result, the IQ y profile is approximately monotonic in the upper ∼ 200 m and high wavenumber baroclinic instability is shut down in July. In the upper 30-50m we expect air-sea heat exchanges to contribute most of the restratification and PV increase. At greater depths, the role of short-wave radiative fluxes and vertical mixing should be limited and our study suggests that rapid restratification may arise from Charney-type BCI. A dedicated examination of the Oleander data for signs of turbulence intensification by C-BCI could be useful to make progress but C-BCI effects, e.g., on spectral EKE distribution variability, may be subtle and hidden by other processes [START_REF] Callies | km)[END_REF].

In the vicinity of the north Pacific Subtropical Countercurrent (STCC) the buoyancy field also exhibits positive ∂ z s b with an equatorward buoyancy gradient. It is therefore consistent with the presence of a Charney baroclinic instability mode (see Fig. 2a in [START_REF] Qiu | Seasonal mesoscale and submesoscale eddy variability along the North Pacific Subtropical Countercurrent[END_REF]. As in the Gulf Stream region the time period when this mode can be influential is limited to spring after baroclinic modes associated with ∼ 100 m thick wintertime mixed layers have disappeared but subsurface stratification remains weak.

In the Southern Ocean, mode water formation regions and intense wintertime destratification also occur along some parts of the ACC. Alongstream heterogeneities and the relative lack of observations there make it more difficult to investigate the existence of C-BCI. C-BCI may also play a role on mesoscale dynamics, e.g., in the vicinity of upper ocean anticyclone which generally have strong fronts and ∂ z s b < 0 in their periphery.

Conclusions

Numerous flavors of parallel flow instabilities have been found relevant to explain ocean turbulence. The upper ocean is prolific in such instabilities because frontal intensification is particularly efficient and also because atmospheric forcings inject large amounts of available potential energy there. Owing to these specificities, baroclinic instabilities associated with the release of APE in the vicinity of fronts are especially relevant to the upper ocean. This study focuses on the description of an oceanic variant of the Charney baroclinic instability (C-BCI), arising from the joint presence of i) an equatorward buoyancy gradient that extends from the surface into the ocean interior and ii) reduced subsurface stratification, e.g., as produced by wintertime convection or subduction. (2015). In direct numerical simulations for an idealized wind and buoyancy forced periodic flow, they find a vigorous forward EKE cascade associated with frontal instabilities. Dissipation occurs preferentially at small scale in these simulations with a secondary role played by bottom friction, i.e., the dissipation route toward larger scales.

In Capet et al. (2008c) and [START_REF] Marchesiello | Submesoscale turbulence in Tropical Instability Waves[END_REF] downscale KE fluxes rapidly vanish below 50 m depth (e.g., see Fig. 2 bottom panel in Capet et al. 2008c) and their quantitative importance is limited. For example the intensity of the KE forward cascade diagnosed in Capet et al. (2008c) implies a dissipation of the order of 1.5 × 10 -9 W kg -1 in the mixed layer, a value that is typically very small in regard to that induced by air-sea interactions. It is also not very effective in terms of mesoscale energy dissipation because it is so confined into the mixed layer.

The present study allows us to place these contradictory results into perspective. Comparison between S1 and S2 demonstrate the major sensitivity of the downscale energy cascade to the upper ocean stratification. For a given large scale buoyancy gradient, baroclinic flows with reduced stratification undergo larger downscale energy cascades. S1 with intense levels of upper ocean KE, a large scale north south buoyancy gradient squared ∼ 1.5 × 10 -8 s -2 and a stratification characterized by N 2 ∼ 0.1-0.2 s -2 , on the lower end of what can be found in the mid-latitude upper ocean [START_REF] Emery | Geographic and seasonal distributions of brunt-väisälä frequency and rossby radii in the north pacific and north atlantic[END_REF], may offer a useful upper bound estimate for the intensity of that cascade. An important conclusion of this study is that even in the favorable situation S1 the downscale advective flux of KE remains a small fraction (≈ 0.05) of the other fluxes involved in the energetics of the system (energy input through the restoring terms mainly as APE, APE to KE conversion, and KE flux toward large scales). Modest sensitivities to, e.g., large scale dissipation intensity through bottom drag may change this number somewhat but do not cast doubt on our general conclusion that the downscale route is of limited significance. Sensitivity to horizontal resolution below ∆x = 1 km would also need to be explored but this is not expected to dramatically affect the intensity of the forward cascade (Capet et al. 2008c). Irrespective of such a sensitivity, our main conclusion on the dependence of the downscale energy flux intensity upon subsurface stratification should hold.

In the light of these results (and in agreement with the conclusions of [START_REF] Jouanno | Dissipation of the energy imparted by mid-latitude storms in the Southern Ocean[END_REF] we attribute the substantial downscale cascades present in [START_REF] Molemaker | Balanced and unbalanced routes to dissipation in an equilibrated eady flow[END_REF] and Barkan et al. (2015) to unrealistically low N 2 and/or N 2 /M 2 . This is evident in Barkan et al. (2015). Isopycnals in their wind and buoyancy forced flow intersect both the bottom and surface over most of the baroclinic zone, which rarely occurs in the open ocean (see their figure 3; their N 2 is 5-10 times weaker than in the real ocean). Likewise, the basic state in [START_REF] Molemaker | Balanced and unbalanced routes to dissipation in an equilibrated eady flow[END_REF] has a M 2 /N 2 dimensional value which is uniform by definition of the Eady flow and equal to 0.2, i.e., much larger than typically found in the real ocean.

Our general results are also at odd with the global estimates of the downscale energy flux obtained by [START_REF] Brüggemann | Routes to dissipation under different dynamical conditions[END_REF]. Based on a set of numerical forced-dissipative experiments [START_REF] Brüggemann | Routes to dissipation under different dynamical conditions[END_REF] establish a relation between the bulk Richardson number of baroclinic flows and the fraction of the APE release < w b > dissipated at small scale as a consequence of the submesoscale energy cascade. Extrapolations based on an eddy permitting realistic global ocean simulation then yields an estimate of small scale dissipation due to submesoscales in the range 0.2-0.5 TW, i.e., a substantial fraction of the energy input to the general circulation.

We question the validity of the extrapolation performed by [START_REF] Brüggemann | Routes to dissipation under different dynamical conditions[END_REF] because it relies on Eady-type experiments to infer ocean interior downscale energy transfers. By construction, surface and bottom surface dynamics is influential at all depths in an Eady flow which makes it a priori ill-suited to explore interior ocean dynamics. It is well known that near-surface and ocean interior significantly differ in their frontogenetic behavior (e.g., MacVean and Woods 1980; [START_REF] Klein | Three-dimensional stirring of thermohaline fronts[END_REF]. Fig. 7 exemplifies these surface-interior differences and confirms unambiguously that not all APE releases are equally suited to feed a downscale energy cascade at submesoscale.

The reasons why the simulations analyzed in [START_REF] Marino | Resolving the paradox of oceanic large-scale balance and small-scale mixing[END_REF] and [START_REF] Pouquet | Geophysical turbulence and the duality of the energy flow across scales[END_REF] develop a substantial downscale energy cascade in contrast to ours may be due to the specifics of their numerical setting. Indeed, we note that their simulations are for a triperiodic ocean (with no boundaries, hence also no near-surface frontal intensification) devoid of large scale energy sink, subjected to a 3D isotropic stochastic forcing, and energetically unequilibrated.

Although the energetic analysis of S1 indicates that the downscale route to dissipation remains limited for the overall equilibration of realistic baroclinic flow, the role of such energy transfers in the local equilibration and dissipation of ocean fronts is presumably important as found, e.g., by The analysis of multiple < w b > vertical profiles for different baroclinic flows and numerical resolutions (eddy resolving and submesoscale rich) raises several issues regarding the existing parameterizations (GM90 and FK08) and provides some limited guidance on how to improve them.

As previously found by [START_REF] Ramachandran | Enhancement in vertical fluxes at a front by mesoscale-submesoscale coupling[END_REF] turbulent motions in the mixed layer and in the interior layer below can couple to produce eddy fluxes which differ from those that would be found if these layers did not interact. In particular, mixed layer and subsurface flows in weakly stratified environments conspire to boost mixed layer restratification. This further pleads for a relaxation of the assumption made in FK08 that the base of the mixed layer behaves like a lid preventing surface-subsurface exchanges.

Contrary to common wisdom, our analyses also reveal that subsurface vertical eddy buoyancy fluxes have not necessarily converged in eddy resolving simulations and they can further increase when refining horizontal resolution, at least down to ∆x = 1 km. This occurs when the flow is susceptible to fine-scale instabilities such as the Charney baroclinic instability modes we considered. A preliminary exploration suggests that this is not very common but may affect subduction regions and dynamics that are important to the global ocean and climate functioning.

More problematically, our simulations show that deep < w b > associated with interior Phillipstype baroclinic instability are also enhanced when going from ∆x = 8 km to ∆x = 1 km. We have at present limited insight to offer about why that may be, given that the deep APE release takes place at scales that should be well resolved at ∆x = 8 km. We are also hesitant on the significance of this result. Although the tendency to intensifying interior baroclinic instability with increased resolution may be physical the forcing used in our simulations amplifies this tendency in a way that is not. Indeed, in our solutions energy injection naturally increases as eddies drive the ocean mean state closer to rest (full depth energy injection is 30 % higher at 1km than at 8km, e.g., in solution S3 with only Phillips baroclinic instability). The energetics of the real ocean interior is fundamentally different.

Because it results from the coupling between surface dynamics (driven by surface buoyancy gradients) and subsurface dynamics (driven by potential vorticity anomalies present below the surface) the Charney baroclinic instability is a process that we find conceptually interesting when contemplating the parameterization of submesoscale subsurface eddy fluxes. The Charney problem suggests for example that both the background stratification N 2 and its vertical derivative ∂ z N 2 be involved in any parameter sweep exploring the sensitivity of < w b > to the flow characteristics. Unfortunately, no suitable analytical expression exists to relate linear growth rates for Charney BCI to environmental parameters, even when N 2 and the interior PV gradient are con-stant [START_REF] Pedlosky | Geophysical Fluid Dynamics[END_REF]. Ways to circumvent this issue numerically can certainly be found if the Charney problem confirms to be a good prototype for near-surface/subsurface baroclinic coupling and property exchanges.

APPENDIX

Construction of initial/reference buoyancy profiles for S1-3

Three northern and southern buoyancy profiles (respectively b N (z) and b S (z)) are needed to construct the initial/reference states corresponding to S1-3. b N (z) and b S (z) are constructed as the sum of four terms: a small depth-independent background stratification that applies equally to the northern and southern profiles and guarantees static stability; an exponential buoyancy profile which enhances near-surface stratification without affecting North-South gradients either; a distorted hyperbolic tangent buoyancy profile responsible for the interior meridional buoyancy gradient and also for a stratification asymmetry between the lower (less stratified) and upper (more stratified) pycnocline, but with no surface meridional buoyancy difference; a hyperbolic tangent buoyancy profile added to the southern profile only which has its inflection point close to the surface and is the only source of surface North-South buoyancy difference. Precisely we have:

b N,S (z) = b (min) +N 2 b (z+h (max) ) +δ b s 1 e exp z s -z z s + 1 2 δ b (1) N,S   1 + tanh   Z N,S (z) -z (1) N,S δ z (1) N,S     + 1 2 tanh(1) δ b (2) N,S 1 + tanh z (2) -z z (2) (A1)
where,

Z N,S (z) = z (1) N,S + z -z (1) N,S    1 + 0.5   z -z (1)
N,S + |zz

(1)

N,S | 1.3 δ z N,S   2    0.5 , δ b (1) N = δ b (1) S 1 + tanh Z S (0) -z (1) S δ z (1) S / 1 + tanh Z N (0) -z (1) N δ z (1) N
, and δ b

(2) N = 0 z denotes depth, h (max) is the ocean depth (4000 m), b (min) = 2.775 × 10 -1 m s -2 , N 2 b = 9.8 ×10 N (in bold) is calculated so that the first hyperbolic tangent term in (A1) does not contribute to the North-South buoyancy difference at the surface. for the solutions S1 (solid black) and S2 (solid gray). S1 uncertainty is also represented as dotted lines based on the assumption that the 360 days run correspond to 12 independent realizations. In the inset the vertical scale is refined in the wavenumber range [10 -4 -2 10 -3 ] rad m -1 to better appreciate the forward cascade.

  Fig.1).

  solutions for the three setups S1-3 and at two different resolutions are presented in this study: ∆x = 8 km and 40 levels ; ∆x = 1 km and 200 levels. These baseline simulations are performed using a simple treatment of vertical mixing, which is possible because air-sea fluxes are absent. Constant background vertical diffusion for buoyancy κ

  cascades in fluxing the additional KE energy input < w b > present in S1 compared to S2, which is ≈ 0.5 mW m -2 . Most of it is being dealt with by the inverse energy cascade (whose increase from S2 to S1 is slightly more than 0.5 mW m -2 as a plausible consequence of changes in fluxes of pressure work through the z=-200 m surface) and not by the forward cascade whose intensity in S1 is only a small fraction (∼ 10 %) of that amount.All these results are not affected by the fact that we include the quiescent northern and southern edges of the domain in the flux calculations. Estimates over a restricted central area of width 680 km lead to a similar ratio (6 %) between the intensity of the downscale KE transfer and the KE input. (When computing Σ Π over an area restricted in the meridional direction, we make use of the 1D hanning window to prevent non-periodicity from affecting our calculation. The obtained flux is then multiplied by a factor 2 to compensate for variance reduction due to the tapering procedure, as inCapet et al. 2008c).Overall, the downscale KE cascade translated into a local dissipation rate leads to a modest 2 × 10 -10 W kg -1 , averaged over the depth range 0-200 m. To put this number into perspective we can derive a time scale for the decay of mesoscale energy in the depth range 0-200 m through the downscale KE route. Vertically integrated KE in this depth range for the mesoscale resolving solution (∆x = 8 km, see Fig.

Fig. 9

 9 Fig. 9 represents the time-averaged vertical buoyancy flux for S1 ml and S2 ml along with those

  Fig. 10 shows monthly mean density structure across the Gulf Stream derived from streamline averaging of all Argo profiles available over the period 2005-2012 between longitudes 40 o W and 60 o W. April and July density sections are represented as well as the integrated meridional PV gradient between two reference profiles located respectively at y 1 =1000 km and y 2 =1650 km south of the Gulf stream axis. The month of April was chosen because it corresponds to post wintertime convection conditions with reduced subsurface stratification and thinning mixed layers, of the order of a few tens of meters. Therefore, the vertical buoyancy fluxes (and restratifying tendency) in the upper 100-200 m of the water column are not driven by mixed layer eddies anymore.

  Background conditions conducive to C-BCI have been widely used in previous process studies of submesoscale flows. C-BCI strengthens near-surface frontal activity with important consequences in terms of turbulent statistics: increased variance of vertical vorticity and velocity, increased vertical turbulent flux of buoyancy which, in turn, implies enhanced nonlinear transfers of KE both upscale and downscale. Dedicated analyses allow us to offer insight into the energetics that underlies submesoscale intensification of flows conducive to C-BCI.The relevance of C-BCI to the real ocean is an important issue that we have explored in a way complementary to[START_REF] Smith | The geography of linear baroclinic instability in earth's oceans[END_REF] by focusing on eastward flows. Linear stability analysis suggests that regions where seasonal episodes of convection and/or subduction of mode waters creates weakly stratified subsurface conditions will be susceptible to Charney baroclinic instability. A numerical sensitivity experiment indicates that mixed layer eddy dynamics does not fundamentally affect the subsurface expression of Charney instability modes. Submesoscale rich realistic simulations of the Gulf Stream region and in particular its southern flank might be useful to further explore Charney instability in a fully realistic context.Notwithstanding the importance of C-BCI in the real ocean the present investigation offers useful theoretical lessons on two subjects, the significance of the downscale energy cascade and the parameterization of eddy induced fluxes of properties in the ocean.a. Lessons on the intensity of the downscale kinetic energy cascade in the oceanThe possibility that submesoscale allows meso-and large-scale flows to equilibrate through a forward cascade of kinetic energy en route to dissipation has been explored for almost two decades. Theoretical arguments support the breakdown of the QG inverse cascade paradigm[START_REF] Charney | Geostrophic turbulence[END_REF] and the existence of the forward route for primitive equations or non-hydrostatic flows[START_REF] Mcwilliams | From stirring to mixing of momentum: Cascades from balanced flows to dissipation in the oceanic interior[END_REF][START_REF] Molemaker | Baroclinic instability and loss of balance[END_REF][START_REF] Müller | Routes to dissipation in the ocean: The 2D/3D turbulence conundrum[END_REF]) but existing studies disagree on its strength and overall energetic significance. Idealized studies tend to indicate that it is significant.[START_REF] Molemaker | Balanced and unbalanced routes to dissipation in an equilibrated eady flow[END_REF] demonstrate that the QG and non-hydrostatic solutions of an Eady flow strongly differ in their submesoscale regime, with implications for the spin down and forced-dissipative equilibration of the solutions. As the submesoscale is progressively resolved small scale dissipation decreases dramatically in QG whereas it is maintained at finite levels in non-hydrostatic solutions. In equilibrated solutions about 1/2 (resp. 1/3) of the kinetic (resp. mechanical, i.e., kinetic and potential) energy dissipation involves advective fluxes toward small scales in the submesoscale range. Similar findings and conclusions are obtained byBarkan et al.

D'

  Asaro et al. (2011) and[START_REF] Johnston | Elevated mixing at a front[END_REF]. Overall though, our study emphasizes the major role played by the inverse kinetic energy cascade even at submesoscale, thereby confirming results fromKlein et al. (2008b) and more recently[START_REF] Sasaki | Impact of oceanic-scale interactions on the seasonal modulation of ocean dynamics by the atmosphere[END_REF].b. Lessons on limitations of and needs for eddy flux parameterizationsThe vertical structure of the APE release < w b > is a key characteristic of baroclinic instability processes. Models with insufficient resolution to produce baroclinic eddies through baroclinic instability rely on parameterizations to mimic the effects of these eddies in terms of fluxes of properties, including vertical fluxes of buoyancy < w b >. Typically, the[START_REF] Gent | Isopycnal mixing in ocean circulation models[END_REF] parameterization (hereafter GM90) is used to represent the effects of mesoscale eddies in non-eddy resolving simulations, and with reduced intensity in eddy-permitting models (∆x 20 km). For ∆x < 10-20 km the processes accounted for by GM are well resolved and GM should be turned off. Mixed layer turbulent motions have much smaller scales than mesoscale eddies so that their effects become well resolved only for ∆x 1 km. In simulations at coarser resolution the Fox-Kemper et al. (2008) parameterization (FK08) accounts for vertical fluxes by mixed layer eddies under the key assumption that submesoscale frontal circulations are well confined into the mixed layer.

Fig. 1 .Fig. 2 .Fig. 3 .Fig. 5 .Fig. 6 .Fig. 7 .Fig. 8 .Fig. 9 .Fig. 10 .

 1235678910 Fig. 1. Northern and southern vertical profiles of buoyancy (left), Brunt-Väisälä frequency (center) and lateral buoyancy gradient modulus (right) at the center of the domain . Continuous (resp. dashed and dotted) lines correspond to S1 (resp. S2 and S3) profiles. Northern profiles for S1 and S3 are identical. Lateral buoyancy gradients for S1 and S2 are also identical. . . . . 47 Fig. 2. Cross-section of initial/reference buoyancy b re f for S1 (top), S2 (middle) and S3 (right). Contour interval is 2.5 × 10 -3 m s -2 . . . . . . . . . . . . . . . . . 48 Fig. 3. Left (a,c,e): IQ y vertical profiles for the initial (black) and mean equilibrated (gray) state of S1 (a), S2 (c) and S3 (e). S1 and S2 surface values of IQ y accounting for the Dirac contribution are shown with asterisks. Note the discontinuity associated with these contributions. Extrema of IQ y profiles for the initial state are shown with open circles. The contribution of planetary vorticity alone is also shown (dashed line). Right (b,d,f): Growth rate of infinitesimal perturbation as a function of wavenumber for S1 (b), S2 (d) and S3 (f). Growth rate maxima are indicated with a plus sign. Vertical profiles of the corresponding eigenvector's modulus are shown between -2000 m depth and the surface (insets). Note the surface trapping of the shortest mode in the case of S1 (dashed line, black for the reference state, gray for the ROMS mean state). . . . . . . . . . . . . . . . . . . . 49 Fig. 4. Surface (a-c, d-f) and 110 m depth (g-i) vertical vorticity (normalized by the Coriolis frequency f) in S1 (top), S2 (middle) and S3 (bottom). Vertical velocity at 110 m is also shown (right, j-l, in m/day). All panels correspond to simulations with horizontal resolution ∆x = 1 km except surface vorticity panels a-c which are for ∆x = 8 km. All panels represent instantaneous fields at t=180 days. Because no surface intensification is present in S2 vorticity at surface and 110 m are virtually identical, hence the similarity between panels f) and i). . . . 50 Fig. 5. Vertical profiles of EKE (a,b), w rms (c,d), < w b > (e,f) and frontogenetic tendency due to horizontal advection F sh (g,h) for S1 (solid black), S2 (solid gray) and S3 (dashed) for ∆x = 1km (left panels) and 8km (right panels). Estimation uncertainty is indicated with markers (crosses for S1, + for S2 and open circles for S3). It is computed for each domain average variable X as ± 1/ √ n X std where we suppose that we only have 12 independent realizations, hence n=12 (see Section 2.4.) . . . . . . . . . . . . . . . . . . 51 Fig. 6. Time averaged KE spectra for the horizontal velocity u h surface fluctuations as a function of horizontal wavenumber magnitude, k = |k h |, i.e., the 2D spectrum is azimuthally integrated in k shells. Spectra for S1 (black), S2 (dark gray) and S3 (light gray) at the surface (left) and 110m depth (right) are represented. For comparison straight lines indicate -5/3 (dotted) and -3 (dot-dash) spectrum slopes. . . . . . . . . . . . . . . . . . . 52 Fig. 7. Fraction of the conversion from APE to EKE C = < w b > that is achieved at spatial scales smaller than k sm = 3 × 10 -4 rad m -1 , as determined from Fourier spectral analysis (see text and eq (2) for details). Profiles for S1 (resp. S2) are represented with a solid (resp. dashed) line. . . . . . . . . . . . . . . . . . . . . . . . . . 53 Fig. 8. Time averaged KE transfer function Σ Π (k) (in mW m -2 , vertical integration over the upper 200 m) for the solutions S1 (solid black) and S2 (solid gray). S1 uncertainty is also represented as dotted lines based on the assumption that the 360 days run correspond to 12 independent realizations. In the inset the vertical scale is refined in the wavenumber range [10 -4 -2 10 -3 ] rad m -1 to better appreciate the forward cascade. . . . . . . . . . 54

  FIG. 5. Vertical profiles of EKE (a,b), w rms (c,d), < w b > (e,f) and frontogenetic tendency due to horizontal advection F sh (g,h) for S1 (solid black), S2 (solid gray) and S3 (dashed) for ∆x = 1km (left panels) and 8km (right panels). Estimation uncertainty is indicated with markers (crosses for S1, + for S2 and open circles for S3). It is computed for each domain average variable X as ± 1/ √ n X std where we suppose that we only have 12 independent realizations, hence n=12 (see Section 2.4.)

  FIG. 9. Vertical profiles of < w b > for simulations S1 ml (solid black) and S2 ml (solid gray). Profiles for S1 (dashed black) and S2 (dashed gray) are repeated from Fig.5. Estimations of uncertainty are also given as in Fig.5based on the assumption that 12 independent realizations are obtained during our 360 days long simulations.

  as expected for interior QG turbulence. S2 spectral slope is intermediate. Spectra for S3 at the surface and 110 m depth are not significantly different given estimate uncertainty (computed from the spectrum variance between model outputs, not shown) and neither are 110m depth spectra for S2 and S3 in the range k 8 × 10 -5 rad m -1 . As discussed in, e.g.,Klein et al.

	(2008a) S2 energization by surface dynamics is progressively restricted to smaller wavenumbers
	with increasing depth. S1 EKE level at 110 m depth is still significantly higher than in S2 for k

  -8 s -2 and other parameter values depend on jet side and/or configuration case, as listed in table 2. dynamics -i. model formulation and initial experiments. Deep Sea Res., 47 (I), 1261-1301. Thomas, L., and R. Ferrari, 2008: Friction, frontogenesis, and the stratification of the surface mixed layer. J. Phys. Oceanogr., 38, 2501-2518. Thomas, L., and J. Taylor, 2010: Reduction of the usable wind-work on the general circulation by forced symmetric. Geophys. Res. Lett., 37, L18 606. Tulloch, R., J. Marshall, C. Hill, and S. Smith, 2011: Scales, growth rates and spectral fluxes of baroclinic instability in the ocean. J. Phys. Oceanogr., 41, 1057-1076.

	The difference between S1 and S2 arises from the exponential term, absent in S1. S1 and S3 differ
	because the final term in (A1) is absent in S3.

ton

TABLE 1 .

 1 Maximum value of S1 advective KE flux Σ Π reached in the submesoscale range; maximum upscale flux intensity reached in the mesoscale range; total energy input F I through zonal restoring (on velocities and buoyancy); contribution to this energy input due to buoyancy restoring F b I alone; conversion from APE to EKE ( < w b >) and rate of change δ t EKE accounting for EKE drift between the solution end state and initial state. Numbers in the first row (resp. second row) are for the reference S1 experiment (resp. the sensitivity experiment with reduced bottom drag S1 rd . In each box numbers before and after the '/' separator refer to integrated values from surface to 200 m and bottom respectively. All numbers are in mW m -2 .

	max(Σ Π )	-min(Σ Π )	F I	F b I	< w b >	δ t EKE
	0.04 / 0.055	1.7 / 3.6	0.82 / 8.4 0.78 / 8.2 0.72 / 5.9 -0.05 / -0.03
	0.12 / 0.85	4.6 / 12	2.7 / 22	3.1 / 26	2.0 / 25	+0.06 / 1.38

This energy distribution should not be mistakenly interpreted as the consequence of an inertial range for surface KE energy. To the contrary it is an additional KE input into the submesoscale range that is responsible for the difference between S1 and S2 power spectra. Both spectra reflect a complex balance between nonlinear fluxes scales and sinks/source terms that are active over the entire wavenumber range.

The fact that pressure forces are responsible for an upward flux of energy is surprising. We relate this to our mean state restoring procedure which leads to mechanical energy injection directly into the ocean interior. This is in contrast with the energization of the real interior ocean by winds through a downward flux of pressure work[START_REF] Roquet | On the patterns of wind-power input to the ocean circulation[END_REF]).

Growth rate of infinitesimal perturbation as a function of wavenumber for S1 (b), S2 (d) and S3 (f). Growth rate maxima are indicated with a plus sign. Vertical profiles of the corresponding eigenvector's modulus are shown between -2000 m depth and the surface (insets). Note the surface trapping of the shortest mode in the case of S1 (dashed line, black for the reference state, gray for the ROMS mean state).

Acknowledgments. Most numerical simulations were performed through computing allocations PRACE RA1616 and GENCI GEN1140 on Curie. We acknowledge financial support from Agence Nationale pour la Recherche through ANR project SMOC. ARGO data were collected and made freely available by the Coriolis project and programmes that contribute to it.

LIST OF TABLES