## Abstract

Baroclinic mixed-layer instabilities have recently been recognized as an important source of submesoscale energy in deep winter mixed layers. While the focus has so far been on the balanced dynamics of these instabilities, they occur in and depend on an environment shaped by atmospherically forced small-scale turbulence. In this study, idealized numerical simulations are presented that allow the development of both baroclinic instability and convective small-scale turbulence, with simple control over the relative strength. If the convection is only weakly forced, baroclinic instability restratifies the layer and shuts off convection, as expected. With increased forcing, however, it is found that baroclinic instabilities are remarkably resilient to the presence of convection. Even if the instability is too weak to restratify the layer and shut off convection, the instability still grows in the convecting environment and generates baroclinic eddies and fronts. This suggests that despite the vigorous atmospherically forced small-scale turbulence in winter mixed layers, baroclinic instabilities can persistently grow, generate balanced submesoscale turbulence, and modify the bulk properties of the upper ocean.

## 1. Introduction

Atmospheric cooling and surface winds frequently mix the surface layer of the ocean. The resulting mixed layer mediates the transfer of heat and momentum between the atmosphere and ocean and thereby affects both the atmospheric climate and the oceanic general circulation. The evolution of the ocean mixed layer has traditionally been understood column by column; atmospheric cooling and wind forcing leads to mixing and deepening of the mixed layer into the thermocline below. It is becoming increasingly clear, however, that lateral exchanges contribute crucially to the dynamical balances of the mixed layer.

Baroclinic instability in the mixed layer, one such agent of lateral exchange, can achieve large vertical buoyancy fluxes by laterally sliding light over dense water, tending to restratify the mixed layer (e.g., Spall 1995; Haine and Marshall 1998; Boccaletti et al. 2007; Fox-Kemper et al. 2008). This restratification modifies the surface properties and thereby feeds back on the surface fluxes of heat and momentum. Baroclinic eddies can also achieve exchanges between the mixed layer and the thermocline below—a process that may be important in subducting heat and atmospheric constituents like carbon into the thermocline as well as bringing nutrients up into the mixed layer, where they can be used for photosynthesis (e.g., Thomas et al. 2008).

The balanced dynamics of baroclinic mixed-layer instabilities suggest they should undergo a seasonal cycle, following the seasonal evolution of the mixed layer itself. The instabilities are expected to be stronger in deep winter mixed layers, in which isopycnals are steep and large amounts of potential energy are available for conversion to kinetic energy. The instability scale, which also scales with the mixed-layer depth, is of order 1–10 km in winter, suggesting these submesoscales should be energized by wintertime baroclinic mixed-layer instabilities. In summer, on the other hand, when mixed layers are shallow, baroclinic instabilities are expected to be much weaker and to occur at much smaller scales.

Observations and models indeed show an energization of submesoscale turbulence (scales of order 1–10 km) in deep winter mixed layers (Mensa et al. 2013; Shcherbina et al. 2013; Sasaki et al. 2014; Callies et al. 2015). The observed characteristics are consistent with a simple model of the balanced dynamics of baroclinic mixed-layer instabilities and the geostrophic turbulence generated by them (Callies et al. 2016). In summer, in contrast, flows 1–10 km in scale carry much less energy and appear to be dominated by thermocline dynamics. It is unclear whether mixed-layer instabilities are damped out in summer because mixing time scales are short in shallow mixed layers or whether they still occur at small scales but are not strong enough to feed a turbulent energy transfer to larger scales and significantly energize the observed 1–10-km scales.

We here seek to go beyond simple arguments based on balanced dynamics because the mixed layer is a place of vigorous atmospherically forced small-scale turbulence. Can baroclinic instabilities grow while the mixed layer is being actively mixed? Or is the instability arrested by mixing and can only develop between mixing events? Is the instability shut off in shallow summer mixed layers?

The goal of this study is thus to delineate the conditions under which baroclinic instability can grow in a turbulent mixed layer. This can help us understand when mixed-layer instabilities are able to energize submesoscales and thereby affect the bulk properties of the upper ocean. We employ an idealized setup that allows the development of both baroclinic instabilities and forced convection in a layer of Boussinesq fluid. The strength of the baroclinic instability is controlled by a prescribed lateral buoyancy gradient, and the strength of convection is controlled by an imposed buoyancy flux at the top and bottom boundaries.

In the limit of weak buoyancy forcing at the top and bottom, it is expected that baroclinic instabilities restratify the mixed layer and thus shut off convection (Mahadevan et al. 2010). This is confirmed in our simple setup. But we find that, conversely, baroclinic growth is not arrested in the limit of strong convection. When baroclinic fluxes are too weak to oppose the destabilization of the forcing, too weak to restratify the layer, baroclinic eddies and fronts still emerge and coexist with the convection that actively mixes the fluid. The baroclinic eddies and fronts organize the convection and modify the bulk properties of the layer.

## 2. Approach

To study the competition between baroclinic instabilities and convection, we depart from previous studies of baroclinic mixed-layer instabilities that considered the transient spindown of a mixed-layer front. In such simulations, the front slumps under the effect of baroclinic eddies and the mixed layer quickly restratifies (e.g., Haine and Marshall 1998; Boccaletti et al. 2007). There is a succession from upright or slantwise convection to baroclinic eddies as the stratification and hence the Richardson number increase, a succession that is expected from linear stability analysis (Stone 1966). Convection and baroclinic eddies do not coexist in these simulations because the restratification by baroclinic eddies shuts off convection.

Mahadevan et al. (2010) attempted to achieve equilibrium between the destabilizing effect of surface forcing and the restratifying effect of baroclinic eddies by blowing a wind down a mixed-layer front. The Ekman transport of such a down-front wind pushes dense over light water, destabilizing the surface layers (Thomas and Lee 2005). The destabilizing effect of this forcing evolves with the flow, however, because whether the Ekman transport is destabilizing depends on the orientation of the front, which under the influence of baroclinic instability starts to meander. The simulations described in Mahadevan et al. (2010) were thus necessarily transient; equilibrium was not achieved, because the increasing misalignment of front and wind stress reduced the effective buoyancy forcing. In the setup described below, we instead use direct buoyancy forcing, which allows us to control the forcing strength and achieve a statistical equilibrium. We can thus study the competition between convection and baroclinic instability in a more controlled environment and reveal that baroclinic instability develops even in the presence of strong convection.

In Callies et al. (2016), the evolution of baroclinic instabilities in the ocean mixed layer was studied using a quasigeostrophic model. Under quasigeostrophic scaling, the restratification by baroclinic eddies is neglected and the mixed-layer stratification remains fixed. This implicitly assumes an equilibrium between the restratification by baroclinic eddies and the destabilizing effects of surface forcing. It also assumes that the balanced flow is affected by the presence of convection only through the mean stratification. One of the goals here is to understand whether such a decoupled evolution of the balanced dynamics in the mixed layer is plausible.

The setup we study here provides a simple example of competing baroclinic instabilities and convective flows in a statistical equilibrium. We impose a fixed meridional buoyancy gradient −*f*Λ in thermal-wind balance with a zonal shear Λ (e.g., Taylor and Ferrari 2010), which provides baroclinicity for the instability ( *f* is the constant Coriolis parameter). In addition, we apply a destabilizing buoyancy flux *F* at the surface and bottom of a domain of depth *H*. For simplicity, there is no representation of the thermocline, so any evolution of the mixed-layer depth is neglected. The restratification by baroclinic eddies for our mixed layer with a rigid bottom, however, is expected to be similar to the case with a thermocline, because instead of folding thermocline water with high potential vorticity (PV) into the mixed layer, such high-PV water is generated at the bottom boundary (Garner et al. 1992).

The artificial buoyancy forcing at the bottom is necessary to allow the system to reach equilibrium. In the real ocean, surface-forced convection can deepen the mixed layer into the thermocline, and equilibrium is not reached. Such a transient evolution means that the system steps through the parameter space explored below. It seems worthwhile to study the equilibrium first, to isolate the parameter dependence, before the more realistic but complicated transient problem is addressed. We posit that the interplay between convection and baroclinic instability in the transient case is qualitatively similar to that exhibited in the equilibrium case studied here.

In this setup, we study the flow under different forcing conditions. The key parameter determining the strength of the forcing is *ε* = *F*/*f*Λ^{2}*H*^{2}. This is the ratio between the destabilizing buoyancy flux and the dimensional part of the Fox-Kemper et al. (2008) scaling for the baroclinic buoyancy flux. We do not include the efficiency factor of Fox-Kemper et al. (2008) because that appears to increase over time in our simulations, and it ultimately depends on the domain size. We vary the parameter *ε* over three orders of magnitude. Over this range, we cover the transition from a weakly forced regime, in which baroclinic instability restratifies the layer and shuts off convection, to a very strongly forced regime, in which convection is essentially upright and the effects of baroclinic eddies negligible. Of particular interest is the intermediate range, where convection persists, but baroclinic instabilities still develop and significantly modify the bulk properties of the fluid.

## 3. Setup

We study the nonhydrostatic Boussinesq equations describing a rotating fluid between two horizontal solid plates at *z* = −*H* and *z* = 0 (Fig. 1). We solve for the perturbations to an imposed background flow with zonal shear Λ and meridional buoyancy gradient −*f*Λ (Taylor and Ferrari 2010). The perturbations satisfy periodic boundary conditions in the horizontal coordinates *x* and *y* in a domain of zonal extent *L*_{x} and meridional extent *L*_{y}. The dynamics are given by the perturbation equations:

where **u** = (*u*, *υ*, *w*) is the velocity vector, *p* is the density-normalized pressure, *b* is buoyancy, *ν* is viscosity, *κ* is diffusivity, and *t* is time. At the boundaries at *z* = −*H* and *z* = 0, we prescribe the buoyancy flux *F*,

and impose no-normal-flow and free-slip conditions: *w* = 0, *u*_{z} = 0, and *υ*_{z} = 0. Note that the free-slip condition on *u* implies small surface and bottom stresses ±*ν*Λ that maintain the background flow against dissipation.

This setup constitutes an energetically consistent system, in which the only energy source is the potential energy supplied by the buoyancy forcing (plus small sources caused by diffusion and viscosity), and the only sink is viscous dissipation. The potential energy equation is

Angle brackets denote a vertical integral over the depth of the domain, the overbar denotes a horizontal average over the domain, and primes denote deviations from this horizontal average. Potential energy is supplied by the buoyancy flux at the boundaries (first term on the right). It can be converted to mean kinetic energy by gravitational slumping (second term on the right) and to eddy kinetic energy by buoyancy production (third term on the right). Downward diffusion of buoyancy also increases the potential energy (fourth term on the right), but this term is generally small for the values of *κ* we use.

To derive the mean kinetic energy budget, we first form the momentum equations of the horizontal-mean velocities. These are free to evolve as they are forced by momentum flux convergences:

The horizontal-mean vertical velocity vanishes, as required by continuity: . We include the imposed background flow in the horizontal-mean flow to obtain a closed energy budget. The mean kinetic energy budget then reads

Gravitational slumping appears as a source of mean kinetic energy (first term on the right). Shear production (second and third terms on the right), when positive, converts mean to eddy kinetic energy. Dissipation of mean kinetic energy (fourth term on the right) is generally small, as is the input of kinetic energy by the surface and bottom stress, which is implied by prescribing a background shear that does not vanish at the boundaries (fifth term on the right). Finally, the eddy kinetic energy budget reads

Buoyancy production (first term on the right) and shear production (second and third terms on the right) appear with opposite signs here. The dissipation term (fourth term on the right) generally dominates the dissipation of energy.

It is worth pointing out a fundamental difference between the potential energy budget of this Boussinesq system and that of an analogous quasigeostrophic system (e.g., Bretherton and Karweit 1975; Salmon 1980; Haidvogel and Held 1980; Larichev and Held 1995), which was used in Callies et al. (2016) to study baroclinic mixed-layer instabilities. In a doubly periodic quasigeostrophic model, there is an infinite reservoir of potential energy available for release by baroclinic instabilities. This is because the effect of vertical buoyancy fluxes on the mean stratification is neglected in quasigeostrophic scaling, such that both the horizontal buoyancy gradient and the stratification are held fixed and baroclinic instabilities can continually grow. In the Boussinesq system, on the other hand, potential energy must be supplied by buoyancy forcing at the boundaries. This is required to supply baroclinic eddies with potential energy available for release. Without this supply, the eddies would restratify the layer until the deformation radius is larger than the finite domain size. At that point, while there remains a horizontal buoyancy gradient, baroclinic eddies are unable to extract any more potential energy and start decaying. In the Boussinesq system, the explicitly captured, atmospherically forced, small-scale turbulence is therefore required to destroy stratification and thereby sustain baroclinic instabilities in the mixed layer.

It should be noted that our setup is different from that studied by Molemaker et al. (2010). They attempted a direct comparison between Boussinesq and quasigeostrophic dynamics, which led them to remove the horizontal-mean momentum flux divergences in (1)–(3), to introduce restoring to a given buoyancy profile and to formulate an evolution equation for the horizontal buoyancy gradient based on energetic arguments. The setup described above is better suited for our purposes, because it has a straightforward energy budget and no interior restoring.

where **x**, **y**, and **z** are the Cartesian unit vectors. In our interpretation of simulation results, we will make use of two fundamental properties of PV dynamics. First, baroclinic instabilities increase the bulk PV by injecting positive PV at the boundaries and folding it into the interior of the fluid (Nakamura and Held 1989). Analogously, baroclinic instability in the real ocean lifts high-PV water from the thermocline into the mixed layer (Garner et al. 1992). Second, convection occurs when *q* < 0 and rapidly restores PV to the marginally stable state *q* = 0 (Hoskins 1974; Emanuel 1994; Haine and Marshall 1998). In the presence of baroclinic shear, this convection is slantwise and the result of symmetric rather than gravitational instabilities. Slantwise convection can produce positive buoyancy stratification (*b*_{z} > 0) by homogenizing buoyancy along tilted absolute momentum surfaces, but it cannot increase the bulk PV above zero. The telltale signs of such slantwise convective adjustment, regions of positive stratification and near-zero PV, have long been known to occur in midlatitude atmospheric fronts (Emanuel 1988). In the ocean, an adjustment from initially strongly negative PV to near-zero PV has recently been observed in the Kuroshio (D’Asaro et al. 2011).

We here want to emphasize the difference between convecting and restratifying states, rather than between upright and slantwise convection. We diagnose restratification using PV: when baroclinic instabilities are strong enough, they increase the PV to *q* > 0, which shuts off convection, including the slantwise kind. Convecting states, on the other hand, have a PV close to zero. We will see, however, that a *q* = 0 state does not mean that no baroclinic eddies are present. It only means that these eddies are too weak to overcome the destruction of PV by the buoyancy forcing.

The system of (1)–(5) is solved numerically using the MITgcm (Marshall et al. 1997), modified to include the additional terms arising from the interaction of the perturbations with the background flow. All simulations have a depth *H* = 100 m, a grid spacing Δ*x* = Δ*y* = Δ*z* = 2 m, and *ν* = *κ* = 10^{−3} m^{2} s^{−1}. We choose harmonic dissipation/diffusion over a Smagorinsky (1963) closure, because it is easier to close the energy budget, and it remains unclear whether a Smagorinsky scheme is appropriate for flows with strong baroclinic fronts. We do not expect the subgrid closure to affect the results discussed in the following.

To reveal the parameter dependence of our setup, we nondimensionalize the system (1)–(5) using the following scales:

This follows Stone (1971), except for the scale of buoyancy and pressure, for which Stone had the imposed stratification as an available scale. The nondimensional system reads

We abbreviate the diffusion operator by *A* = *δ*^{2}(*A*_{xx} + *A*_{yy}) + *A*_{zz}. The boundary conditions are −*σ*^{−1}*αb*_{z }*= ε*, *w* = 0, *u*_{z} = 0, and *υ*_{z} = 0 at *z* = −1 and *z* = 0. The nondimensional parameters of the problem are the ratio of the forcing to the scaling for the baroclinic flux proposed by Fox-Kemper et al. (2008), *ε* = *F*/*f*Λ^{2}*H*^{2}, and Stone’s (1971) nonhydrostatic parameter, *δ* = *f*/Λ. The strength of viscosity is measured by the Ekman number *α* = *ν*/*fH*^{2}, and the Prandtl number is *σ* = 1. Dissipation and diffusion are posited to be unimportant for the bulk behavior as long as they are sufficiently small. The remaining two nondimensional parameters measure the domain size: *λ*_{x} = *fL*_{x}/Λ*H* and *λ*_{y} = *fL*_{y}/Λ*H*. The nondimensional potential energy budget reads

the mean kinetic energy budget is

and the eddy kinetic energy budget is

The nondimensional form of PV is

Note that a state with no stratification and no perturbation to the background flow has a background PV of *q* = −1.

The scales chosen above are appropriate when rotation and baroclinicity play an important role in the dynamics. In the limit of strong forcing (*ε* → ∞), however, the dynamics revert to the classic problem of nonrotating upright convection between two plates (e.g., Emanuel 1994). In this limit, more appropriate scales can be formed from the forcing *F* and depth *H* (e.g., Marshall and Schott 1999):

This nondimensionalization will prove useful when considering strongly forced simulations.

## 4. Experiments

We focus on a series of experiments with *δ* = 1 for two reasons. First, linear stability analysis suggests that with this parameter, the baroclinic instability is already close to the hydrostatic limit *δ* → 0 (Fig. 2; Stone 1971) that is of particular interest for strong ocean fronts, where shears can greatly exceed *f* (e.g., D’Asaro et al. 2011). Second, this parameter choice is practical because the most unstable baroclinic mode then has an aspect ratio close to unity (Stone 1966), which allows us to conveniently resolve the baroclinic mode as well as upright convection. The increasing aspect ratio of the baroclinic instability as *δ* → 0 makes the concurrent simulation of baroclinic and convective motions increasingly challenging, because there is an increasing scale separation between the two types of motion.

We recapitulate the linear properties of the system by performing a linear stability analysis for the baroclinic axis (no meridional variations) following Stone (1971). We consider a mean state with stratification *b*_{z} = 1, such that *q* = 0. This is taken to be representative of the weakly stratified state expected in the limit of strong convection. The growth rates (Fig. 2) show a low-wavenumber geostrophic mode and a high-wavenumber ageostrophic mode (Stone 1970; Nakamura 1988; Molemaker et al. 2005). We focus on the geostrophic mode, because it has a much larger growth rate and represents the balanced dynamics whose interaction with convection we seek to investigate. This mode is most unstable at a zonal wavenumber *k* = 1.13 and has a short-wave cutoff at *k* = 1.68. Because of the finite domain size, we only resolve a discrete set of wavenumbers (Fig. 2). The domain size *λ*_{x} = 10 is chosen such that roughly two wavelengths of the most unstable mode fit into the domain. For comparison, we also perform simulations in narrow domains of zonal extent *λ*_{x} = 2, in which no baroclinic instability can develop (cf. Taylor and Ferrari 2010). Symmetric instabilities, which occur for *q* < 0, only require across-shear variations, so they are well represented in both domains.

In a set of fully nonlinear simulations, we vary the forcing parameter *ε* over three orders of magnitude, covering very weak to very strong forcing conditions. To be able to distinguish a priori between weak and strong forcing conditions, we begin with an unforced simulation (*ε* = 0). This allows us to measure the (dimensionless) buoyancy flux generated by baroclinic instability in the absence of convection. It successfully distinguishes between weak forcing *ε* < *B*, for which restratification occurs, and strong forcing *ε* > *B*, for which convection actively mixes the layer (cf. Mahadevan et al. 2010). To isolate the effect of baroclinic instability, the series of forced simulations are performed in pairs: one with a wide and one with a narrow domain. Baroclinic instability has a short-wave cutoff and cannot develop in the narrow domain, while the instability does fit into the wide domain. All simulation parameters are listed in Table 1.

### a. Unforced simulation

We determine the flux *B* by measuring it in an unforced simulation instead of using the Fox-Kemper et al. (2008) scaling. We resort to this empirical exercise because the scaling does not capture an increase in buoyancy flux as the eddies become increasingly energetic past the initial phase of their development (see Callies and Ferrari 2017, manuscript submitted to *J. Phys. Oceanogr.*). The energy level of the eddies is ultimately set by the finite size of the domain, which translates into a dependence of *B* on the domain size. It is less clear what limits how energetic mixed-layer eddies can become in the real ocean, where scales up to the mesoscale can be energized. We circumvent this question by determining *B* empirically for the domain size used.

We initialize the simulation with a uniform stratification *b*_{z} = 1 and add small random perturbations to *u*, *υ*, and *b* drawn from normal distributions with standard deviations 10^{−7} m s^{−1} and 10^{−7} m s^{−2}. The *b*_{z} = 1 state is marginally stable to symmetric instability, so the development is dominated by the baroclinic mode (Fig. 2). The evolution, once the system has reached finite amplitude and nonlinear terms have become important, is illustrated with surface buoyancy snapshots in Fig. 3. In accordance with the linear stability analysis, the mode that reaches finite amplitude first has a zonal wavelength of 5 (Fig. 3a). This baroclinic wave perturbation then grows in scale, spanning the whole domain size at wavelength 10 (Figs. 3b,c), presumably both because of the linear growth of the wavelength-10 mode and nonlinear energy transfer from the wavelength-5 mode. The domain-filling mode reaches finite amplitude, breaks, and rolls up into a domain-filling vortex (Figs. 3c,d). These baroclinic perturbations drive a vertical buoyancy flux that restratifies the system. The restratification shifts the short-wave cutoff of the instability to larger scales and eventually renders the configuration stable to all modes that fit into the domain. Linear theory predicts a critical value of stratification past which all unstable modes have scales too large to fit into the domain. The shutdown of the vertical buoyancy flux occurs exactly when the critical stratification is reached, confirming that linear theory correctly predicts whether potential energy can be released (cf. Stone 1972; Rhines 1977; Salmon 1978; Fox-Kemper et al. 2008). After the restratification shutdown, the domain-filling vortex persists with roughly constant kinetic energy. We take the maximum *B* = 0.08 as representative of the buoyancy flux that can be achieved by baroclinic instability in the given domain.

This evolution can also be traced in the energy budget shown in Fig. 4. We show all terms in (21)–(23) but here focus on the buoyancy production that dominates the generation of eddy kinetic energy. This term has three distinct peaks at *t* = 44, 62, and 76. The first peak in buoyancy production is associated with the wavelength-5 mode (Fig. 3a), and the third peak is associated with the wavelength-10 mode (Fig. 3c). The second peak has both modes contributing and occurs during the transition from the smaller to the larger mode (Fig. 3b). The third peak is largest and reaches *B* = 0.08. After this peak, the buoyancy production drops, as the restratification has rendered the configuration stable. Subsequently, the kinetic energy tendency is small; the little dissipation the domain-filling vortex experiences is partially offset by a small amount of buoyancy production. Furthermore, the energy budget exhibits the signature of oscillations with a period somewhat longer than inertial that have been triggered by the instability. These oscillations may be baroclinically modified inertial oscillations (Whitt and Thomas 2013) but are not the focus of this paper.

### b. Weak forcing

We now turn to forced simulations and start with a weak forcing *ε* = 0.01, which is much smaller than the *B* = 0.08 estimated from the unforced simulation. The simulation is initialized as before with the *q* = 0 state with *b*_{z} = 1 and no flow anomalies. Convection quickly sets in, but the baroclinic instability develops and restratifies the system. The flow reaches a statistical equilibrium around *t* = 100, characterized by a domain-filling baroclinic vortex (Fig. 5a) and strong stratification greatly exceeding the initial *b*_{z} = 1 (Fig. 6a). The equilibrium is achieved by a balance between the destabilization by the buoyancy forcing and the restratification by the baroclinic mode. There is evidence of convection only in the center of the cyclonic vortex, where the stratification is weakest. The convective plumes reaching the surface can be discerned in the surface buoyancy field, which is otherwise smooth (Fig. 5a).

In the transient development, the peak buoyancy production of about 0.04 is larger than the imposed forcing (Fig. 7). This large buoyancy production draws on the potential energy of the initial state but cannot be sustained by the forcing. The system restratifies, which stabilizes the baroclinic modes and leads to a decay in buoyancy production. The baroclinic instability is not shut off entirely like in the unforced simulation, however, because the destabilizing forcing keeps the system slightly unstable to baroclinic instability. In statistical equilibrium, the buoyancy production thus settles to a value of about 0.01, corresponding to the forcing *ε*. The energy cycle of this equilibrated state is straightforward: potential energy is created by the forcing, potential energy is converted to eddy kinetic energy by buoyancy production, and eddy kinetic energy is dissipated by viscosity. The mean kinetic energy plays no significant role in the energetics.

To further characterize the statistical equilibrium, we show the vertical profiles of the horizontal-mean PV averaged in time over the equilibrated state (Fig. 8a). PV has increased from the initial *q* = 0 and is significantly positive over the bulk of the domain. It is negative only near the top and bottom boundaries, where the boundary conditions force it to be with our choice of *α* = 4.64 × 10^{−4} and *σ* = 1 (Table 1).

This equilibrium state is very different from that achieved in a narrow domain (cf. Taylor and Ferrari 2010). In a domain of width *λ*_{x} = 2, no baroclinic instability can develop, because no unstable mode fits into the domain (Fig. 2). The imposed buoyancy flux is carried by slantwise convection rather than by a baroclinic vortex. The system adjusts to a state that is marginally stable to symmetric instability, that is, one with (Fig. 8a). The stratification is determined by the imposed geostrophic shear through the requirement that and settles at about . This is much weaker a stratification than achieved by baroclinic restratification in the wide domain (Fig. 6a).

### c. Moderate forcing

Next, we consider a moderately forced case, in which the forcing *ε* = 0.1 is slightly larger than the *B* = 0.08 estimated from the unforced simulation. Again, convection sets in very quickly. But a baroclinic mode develops and appears after *t* = 50. The system settles into a statistical equilibrium that exhibits a coexistence of a domain-filling baroclinic vortex and smaller-scale convection (Fig. 5b).

The interplay between the baroclinic and convective flows leads to an interesting state exhibiting significant buoyancy stratification but a PV close to zero (Figs. 6b, 8b). PV is restored to zero by convection, but the convective flows are modified by fronts generated by the finite-amplitude baroclinic instability. At these fronts, the convection is slantwise rather than upright, giving rise to the buoyancy stratification.

The buoyancy stratification is stronger than the expected from slantwise convection operating on the background shear (Fig. 6b). In fact, purely convective flows for this forcing—as simulated in a narrow domain—lead to an even smaller stratification, because momentum fluxes eliminate part of the geostrophic shear. This confirms that the baroclinic vortex significantly modifies the convection.

The horizontal-mean PV is close to zero throughout the bulk of the domain, as expected for slantwise convection (Fig. 8b). There is significant cancellation between a positive and the remaining terms in (24), showing the influence of the baroclinic vortex. This cancellation renders the total PV profile very similar to that of the narrow domain, in which the stratification is much weaker. This similarity in suggests that the flows in both the wide and narrow domains are equally convecting, restoring to a convectively neutral state . The baroclinic eddies do not manage to restratify in the sense of increasing the horizontal-mean PV above zero in the interior of the domain and thereby shutting off (upright or slantwise) convection. But the baroclinic eddies do provide additional shear that modifies the convection and allows some buoyancy stratification.

Another fingerprint of slantwise convection can be found in the energy budget (Fig. 9). The bulk of eddy kinetic energy is generated by buoyancy production, but there is a significant contribution from shear production. This is expected for slantwise convection, which can draw on the kinetic energy of the mean flow (Haine and Marshall 1998; Thomas et al. 2013). The mean kinetic energy is supplied by gravitational slumping that attempts to restore geostrophic balance against the destabilizing effect of the momentum fluxes.

### d. Strong forcing

As the forcing is increased to *ε* = 1, much larger than *B* = 0.08, convection is expected to dominate the flow. We find, however, that the strong convective flows are still modified detectably by baroclinic eddies (Fig. 5c). Baroclinic instabilities develop on top of the strong convection, increase the shear at fronts, and render the convection more slantwise.

The buoyancy stratification clearly shows this effect of the baroclinic instability (Fig. 6c). The stratification that results from this interplay of the baroclinic flow with convection is increased above the weak stratification that results from pure convection in the narrow domain. As expected from *ε* ≫ *B*, the baroclinic instability is not able to achieve true restratification, in the sense of increasing the horizontal-mean PV above zero in the interior of the domain (Fig. 8c). In fact, the PV profiles for the wide and narrow domains are indistinguishable; both are restored to zero by strong convection.

The eddy kinetic energy production is now strongly dominated by buoyancy production, but a significant contribution from shear production is still detectable (Fig. 10). Slantwise convection, while weak, thus still manifests itself in the energy budget.

### e. Very strong forcing

At a very strong forcing of *ε* = 10, the flow is finally in the limit of no rotation and no shear and dominated by strong upright convection (Fig. 5d). The effect of the baroclinic instability is negligible, and the characteristics of the flow in the wide domain are indistinguishable from that in the narrow domain.

The buoyancy profile is now unstratified over the bulk of the domain (Fig. 6d). The profiles are indistinguishable between the wide and narrow domains. The profiles are also very similar to those found in the narrow domain with *ε* = 1 (but not the wide one) if the comparison in this strongly convective limit is done in the alternative nondimensionalization (25)–(26) shown on the upper buoyancy axis (Fig. 6c). Similarly, the PV profiles are indistinguishable between the narrow and wide domains and also match that of the narrow domain with *ε* = 1 in the alternative nondimensionalization (Figs. 8c,d).

The energy budget is now a straightforward conversion of potential energy into eddy kinetic energy by buoyancy production, balanced by dissipation of eddy kinetic energy, as expected for pure upright convection (Fig. 11). There is no significant contribution from shear production.

## 5. Discussion

Our simulations suggest that baroclinic instabilities are remarkably resilient to the presence of convection. When the forcing is weak, the instability’s buoyancy flux is larger than the imposed flux and restratification occurs. Positive PV is injected at the boundaries, and the bulk PV of the fluid is increased above zero, shutting off convection. If the forcing is stronger, exceeding the buoyancy flux that can be generated by baroclinic instability, convection persists and keeps the bulk PV near zero. Baroclinic growth still occurs, however, and modifies the convection by forming baroclinic eddies and fronts, along which baroclinic shear is intensified and convection more slantwise. For very strong forcing, this effect is small and convection becomes indistinguishable from the limit of upright convection. This limit appears to be approached gradually, with no qualitative transition at which the baroclinic instability is shut off.

These results were obtained for a background flow with moderate shear, *δ* = *f*/Λ = 1. We expect, however, that the results carry over to the hydrostatic limit *δ* → 0. As shown by Stone (1971) and in Fig. 2, the linear dynamics are similar for *δ* = 1 and *δ* = 0. In the presence of convection, we expect the baroclinic instability to persist in the *δ* → 0 regime just as well as in the *δ* = 1 case analyzed here. The increasing scale separation between the instability and upright convection as *δ* → 0, which makes the concurrent simulation in this limit so challenging, suggests that strong energy transfer between the processes becomes more difficult. Future work should test this speculation.

One change in the dynamics that is expected as *δ* is decreased is that convection will have a tendency to be more slantwise than in the *δ* = 1 case considered here. The parameter *εδ* is the cube of the ratio of the convective velocity scale (*FH*)^{1/3} to the change in geostrophic flow over the boundary layer Λ*H*. Taylor and Ferrari (2010) and Thomas et al. (2013) argued that this ratio determines whether the geostrophic shear plays a role in the dynamics and renders convection slantwise (*εδ* ≪ 1) or whether convection is upright and overcomes the geostrophic shear (*εδ* ≫ 1). As *δ* is decreased below 1, an increasingly large *ε* is needed to overcome the effect of geostrophic shear to render convection upright. The more realistic *δ* ≪ 1 regime is thus expected to exhibit convection that is more slantwise.

Nonhydrostatic effects become important in the linear dynamics for *δ* ≫ 1 (Stone 1971). Horizontal buoyancy gradients are weak in this regime, and baroclinic instabilities can only play a role when the forcing is weak, because *ε* = *δ*^{2}*F*/*f *^{3}*H*^{2}. It is left for future inquiry to explore whether this regime with weak buoyancy gradients and weak forcing has oceanic relevance or whether baroclinic instabilities are unimportant where buoyancy gradients are weak, in which case the dynamics are dominated by upright or rotating convection (e.g., Julien et al. 1996; Marshall and Schott 1999).

To parameterize the restratification caused by baroclinic mixed-layer instabilities in coarse-resolution ocean models, it is necessary to understand under what atmospheric forcing conditions such restratification occurs. Our analysis shows that the magnitude of the baroclinic flux compared to the forcing distinguishes between convecting and restratifying conditions. As shown in Callies and Ferrari (2017, manuscript submitted to *J. Phys. Oceanogr.*), however, the baroclinic flux depends on the eddy scale, for which we have no good prediction. Understanding this dependence is crucial for the parameterization effort, because it impacts both whether or not restratification occurs and how strong it is when it does occurs.

The finding that baroclinic instabilities can develop even in the presence of strong convection can help us understand the energization of submesoscale turbulence in winter. The persistence of baroclinic instabilities explains how they can develop and produce energetic submesoscale flows even when the mixed layer is actively convecting. The coexistence of baroclinic instability and convection in moderate to strong forcing conditions is consistent with the presence of energetic balanced flows at scales of order 1–10 km and the persistently weak stratification in deep winter mixed layers. Strong atmospheric forcing produces copious amounts of potential energy available for release by baroclinic instabilities and for turbulent redistribution across scales.

The persistence of baroclinic instabilities in the presence of convection also lends credence to the simple model of wintertime submesoscale turbulence developed by Callies et al. (2016). It is at least plausible that the balanced flow evolves with little interference from convection and other small-scale turbulence, in a background provided by a quasi equilibrium between the restratification by baroclinic instabilities and the destabilizing effect of atmospheric forcing. Qualitative agreement between the simple model and submesoscale observations is encouraging. But future work should address whether there are important interactions between the balanced dynamics and small-scale turbulence that are missing from a quasigeostrophic description of submesoscale turbulence in the mixed layer.

The situation in summer is less clear. The fact that baroclinic instabilities can grow in the presence of convection means that they are unlikely to be damped out in shallow summer mixed layers, even if mixing time scales are short. The observed lack of energization of order 1–10-km flows in summer must then be explained by a lack of energy transfer to these scales from the small instability scale. A plausible scenario is that the small amount of potential energy available for release in shallow summer mixed layers is insufficient to significantly energize order of 1–10-km flows. Baroclinic instabilities may successfully grow and restratify the mixed layer in summer, but they quickly exhaust the energy fueling their growth. Unlike in winter, the summertime energy input from atmospheric forcing is weak, leaving the energy throughput too feeble for the instability to transfer significant energy into 1–10-km flows. This scenario should be made quantitative, and its consistency should be tested in a setup that allows for the evolution of the mixed-layer depth.

We finally note that the setup studied here does not allow an assessment of how much energy from submesoscale mixed-layer instabilities is transferred to successively larger scales versus how much is lost to dissipation at small scales (cf. Capet et al. 2008; Molemaker et al. 2010; Sasaki et al. 2014). While the baroclinic eddies have a clear tendency to energize successively larger scales in the transient phase of our simulations, there is no large-scale energy sink, so there cannot be a sustained inverse cascade. The equilibration must instead occur through small-scale dissipation of all of the kinetic energy generated baroclinically. One would have to include motion all the way up to the mesoscale to allow a realistic energy transfer to larger scales. One may speculate, however, that the resilience of baroclinic eddies to the presence of convection might mean that the inverse cascade is more robust than previously thought.

## Acknowledgments

Jean-Michel Campin is thanked for assistance in implementing the MITgcm modifications necessary to prescribe a background flow. Support from the National Science Foundation under Grant OCE-1233832 is gratefully acknowledged.

## REFERENCES

*Numerical Models of Ocean Circulation*, R. O. Reid, A. R. Robinson, and K. Bryan, Eds., National Academy of Sciences, 237–249.

*Atmospheric Convection.*Oxford University Press, 580 pp.

*Marine Modeling*, E. Goldberg et al., Eds., Vol. 6,

*The Sea–Ideas and Observations on Progress in the Study of the Seas*, John Wiley and Sons, 189–318.

*Ocean Modeling in an Eddying Regime*,

*Geophys. Monogr.*, Vol. 177, Amer. Geophys. Union, 17–38, https://doi.org/10.1029/177GM04.

## Footnotes

This article is included in the LatMix: Studies of Submesoscale Stirring and Mixing Special Collection.

This article has a companion article which can be found at http://journals.ametsoc.org/doi/abs/10.1175/JPO-D-17-0029.1

© 2018 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).