Abstract

The restratification of the oceanic surface mixed layer that results from lateral gradients in the surface density field is studied. The lateral gradients are shown to be unstable to ageostrophic baroclinic instabilities and slump from the horizontal to the vertical. These instabilities, which are referred to as mixed layer instabilities (MLIs), differ from instabilities in the ocean interior because of the weak surface stratification. Spatial scales are O(1–10) km, and growth time scales are on the order of a day. Linear stability analysis and fully nonlinear simulations are used to study MLIs and their impact on mixed layer restratification. The main result is that MLIs are a leading-order process in the ML heat budget acting to constantly restratify the surface ocean. Climate and regional ocean models do not resolve the scales associated with MLIs and are likely to underestimate the rate of ML restratification and consequently suffer from a bias in sea surface temperatures and ML depths. In a forthcoming paper, the authors discuss a parameterization scheme to include the effect of MLIs in ocean models.

1. Introduction

The oceanic mixed layer (ML) is most often described in terms of small-scale vertical processes that reduce the vertical gradients of tracers and momentum and large-scale horizontal motions that stir and mix the lateral gradients. In this description, the small-scale vertical processes are characterized by isotropic eddies with scales on the order of the ML depth (e.g., 100 m) such as Langmuir cells driven by winds and convection cells driven by surface buoyancy fluxes. The horizontal stirring, instead, is thought to be dominated by large-scale ocean currents and mesoscale eddies with scales of O(10–100) km. Dynamics at scales between 100 m and 10 km, the submesoscales, are assumed to be subdominant. Recent work by Spall (1997), Haine and Marshall (1998), Ferrari and Rudnick (2000), and Thomas (2005) has challenged this view: motions in this intermediate range of scales are often leading order in the ML temperature, salinity, and momentum budgets and should be included in theories and models of the upper ocean. The present work focuses on the importance of submesoscale dynamics for ML restratification. We show that baroclinic instabilities develop at the submesoscale and act to continuously restratify the surface ML. Present models of the upper ocean ignore these dynamics and are likely to underestimate the rate of ML restratification.

As an illustration of how these lateral instabilities arise consider the following scenario. A winter storm hits the open ocean, mixing the top 100 m of the water column over a patch of a few hundred kilometers squared. Once subsided, the storm leaves behind a homogenized layer in which horizontal variations of salinity and temperature have survived (Price 1981; Ferrari and Rudnick 2000), yet vertical variations have been virtually erased. An adjustment process begins that restratifies the surface layer by slumping the nearly vertical isopycnals. Restratification is the process whereby warm fluid overrides cold fluid. This trivially happens when the fluid is heated from above (as on summer days), but here restratification occurs by dynamical adjustment. This slumping, initially a simple gravitational overturning, is subsequently modified by rotation leading to a geostrophic adjustment (Rossby 1937, 1938). In this paper we show that that the geostrophically adjusted state is further unstable to submesoscale baroclinic instabilities that continue restratification. In fact, the bulk of the restratification happens after the baroclinic instabilities set in. In due course another storm hits, the stratification is once again erased, and the cycle resumes. The restratification ensuing from the initial gravitational slumping is relatively well understood (e.g., Ou 1984; Tandon and Garrett 1995), but instabilities that arise along the adjusting front are usually neglected. Restratification driven by submesoscale baroclinic instabilities plays an important role in the heat and salt budgets of the ML.

a. Deep and shallow baroclinic instabilities in the ocean

The classical problems of baroclinic instability considered by Charney (1947) and Eady (1949) are successful in explaining the preferred length scale and growth rate of mesoscale disturbances observed in the ocean. These disturbances span a large fraction of the ocean depth and represent the subsidence of the thermocline stratification under the action of gravity and rotation. However, a second class of small-scale baroclinic disturbances appears in the presence of reduced stratification at the boundaries (Blumen 1979). These small-scale instabilities represent the slumping of density fronts within the surface and bottom ML. Nakamura (1988) discusses the properties of submesoscale disturbances for the atmosphere. In this paper we consider the oceanographic case and test the hypothesis that submesoscale baroclinic instabilities (MLIs) play a role in driving ML restratification. The fundamental characteristics of MLIs can be illustrated from ocean observations of shear and stratification. Consider the potential density SeaSoar1 section from the subtropical North Pacific Ocean region shown in Fig. 1. The ML, visible as a layer of weak stratification in the upper 100 m overlying the more stratified interior, is not horizontally homogeneous: there are numerous lateral density gradients down to scales of a few kilometers. Simultaneous ADCP measurements confirm that these lateral gradients are in thermal wind balance. Variations in ML depth are visible as isopycnal domings with scales on the order of the first deformation radius, like that visible between 127° and 125°W, and are likely to be associated with mesoscale eddy heaving.

Fig. 1.

Potential density along a straight section between (32.5°N, 122°W) and (35°N, 132°W), i.e., between the California Current and the middle of the subtropical Pacific gyre, as measured by a sawtooth SeaSoar tow. Longitude is shown in the lower x axis and the corresponding along-track distance in km is shown at the top of the figure. Data are averaged in bins of 3 km in the horizontal by 8 m in the vertical and contoured every 0.2 kg m−3. An ML of weak stratification is evident in the upper 100 m. The ML base is marked by a region of enhanced stratification above the permanent thermocline. The ML is also characterized by lateral density gradients. The data were collected as part of an upper-ocean study of the North Pacific (Ferrari and Rudnick 2000).

Fig. 1.

Potential density along a straight section between (32.5°N, 122°W) and (35°N, 132°W), i.e., between the California Current and the middle of the subtropical Pacific gyre, as measured by a sawtooth SeaSoar tow. Longitude is shown in the lower x axis and the corresponding along-track distance in km is shown at the top of the figure. Data are averaged in bins of 3 km in the horizontal by 8 m in the vertical and contoured every 0.2 kg m−3. An ML of weak stratification is evident in the upper 100 m. The ML base is marked by a region of enhanced stratification above the permanent thermocline. The ML is also characterized by lateral density gradients. The data were collected as part of an upper-ocean study of the North Pacific (Ferrari and Rudnick 2000).

A quasigeostrophic (QG) stability analysis (Gill et al. 1974; Blumen 1979) is performed on a background state computed from a combination of the SeaSoar section in Fig. 1 and Levitus climatology (Fig. 2). Details of the calculations are given in appendix A and the results are shown in Fig. 3. Stratification and shear are weak in the ML, increase suddenly across the ML base, and decrease exponentially through the thermocline. The most unstable eigenmodes correspond to alongfront perturbations and their growth rates are shown as a function of the alongfront wavenumber in Fig. 3. There are two separate classes of instability, as follows.

  1. Deep mesoscale instabilities: These modes have scales greater than 20 km and growth time scales near a month. Their vertical structure, shown in the upper-left panel of Fig. 3, is deep and penetrates to the ocean bottom. These modes are low vertical mode baroclinic instabilities and are the source of the oceanic mesoscale eddy field. These instabilities are known to exert a strong influence on the rate of lateral spreading of tracers in the upper ocean.

  2. Shallow MLIs: The second class of instability is composed of smaller modes, 200 m–20 km, with growth time scales near a day. These modes are trapped in the surface ML (upper right panel of Fig. 3) and energize by slumping ML density fronts; they are the QG approximation to MLIs. They have small vertical scales, implying large shears, and are potentially important in driving ML restratification. The goal of this paper is to assess whether this restratification process is fast enough to compete with the other processes that keep the ML well mixed.

Fig. 2.

(left) Buoyancy frequency N = gρz/ρ0 and (right) vertical shear Uz = y/0 averaged in the horizontal for the SeaSoar section shown in Fig. 1 between 133° and 130°W. The vertical gradients are computed across 8 m, while the horizontal gradients are computed across 10 km. The profiles are extended to the ocean bottom by matching the SeaSoar estimates in the upper 320 m with estimates based on Levitus climatology for the rest of the water column. Details of the calculation are given in appendix A.

Fig. 2.

(left) Buoyancy frequency N = gρz/ρ0 and (right) vertical shear Uz = y/0 averaged in the horizontal for the SeaSoar section shown in Fig. 1 between 133° and 130°W. The vertical gradients are computed across 8 m, while the horizontal gradients are computed across 10 km. The profiles are extended to the ocean bottom by matching the SeaSoar estimates in the upper 320 m with estimates based on Levitus climatology for the rest of the water column. Details of the calculation are given in appendix A.

Fig. 3.

(bottom) Stability analysis of the mean shear and buoyancy shown in Fig. 2. The instability is dominated by two distinct modes: an interior instability with wavelength close to the first internal deformation radius of O(50 km) and an MLI peaking at wavelengths close to the ML deformation radius of O(1) km. (top left) The interior instability has a vertical structure spanning the whole thermocline depth and represents the baroclinic instability that generates the mesoscale eddy field (Eady 1949). (top right) The MLI is confined to the ML and represents restratification due to ageostrophic instability within the ML (Stone 1971).

Fig. 3.

(bottom) Stability analysis of the mean shear and buoyancy shown in Fig. 2. The instability is dominated by two distinct modes: an interior instability with wavelength close to the first internal deformation radius of O(50 km) and an MLI peaking at wavelengths close to the ML deformation radius of O(1) km. (top left) The interior instability has a vertical structure spanning the whole thermocline depth and represents the baroclinic instability that generates the mesoscale eddy field (Eady 1949). (top right) The MLI is confined to the ML and represents restratification due to ageostrophic instability within the ML (Stone 1971).

The differences in spatial and temporal scales between MLIs and interior instabilities can be understood as follows. Baroclinic instability can be characterized by the alongfront velocity scale U and the length scale L associated with a disturbance developing along the front. The scales for the developing instability in terms of the Rossby number Ro = U/fL are

 
formula

For baroclinic instabilities (QG or not) the disturbances grow at a scale near the local deformation radius (Stone 1966):

 
formula

where H is the vertical scale of the mode considered and N and f are the buoyancy and inertial frequencies. Therefore, for these unstable modes,

 
formula

where the bulk Richardson number is defined as Ri = N2H2/U2. Equation (3) determines the Rossby number (and therefore the scale) of the most unstable modes for a given Ri. In the strongly stratified ocean interior, Ri ≫ 1 leads to Ro ≪ 1. This is the familiar QG limit appropriate for deep slow mesoscale eddies with Tf−1 and LUf−1. In the weakly stratified ML, Ri = O(1) selects Ro = O(1) and small (LUf−1) and fast (Tf−1) ageostrophic MLIs. For these modes the QG approximation is inappropriate. Ageostrophic effects increase the growth rates and the spatial scales of the instabilities compared to QG approximation to MLIs, with the fastest-growing-mode scales enlarging as L/LQG1 + Ri−1 (Stone 1966, 1971).

b. Shallow baroclinic instabilities and ML restratification

The literature on shallow baroclinic instabilities in the ML is vast, but largely ignores their effect on ML restratification (Samelson 1993; Barth 1994; Young and Chen 1995; Beron-Vera and Ripa 1997; Spall 1997). Notable exceptions are the works of Jones and Marshall (1997) and Haine and Marshall (1998) that focus on the role of baroclinic instability in restratifying the water column after deep convection at high latitudes. However, the high-latitudes convective chimneys they consider span a large fraction of the water column and there is no clear separation between submesoscale ML dynamics and interior mesoscale dynamics.

The literature on ML restratification, on the other hand, largely ignores the role of baroclinic instabilities. Ou (1984) and Tandon and Garrett (1994, 1995) study the slumping via geostrophic adjustment of lateral density fronts following an intense vertical mixing event. They show that the adjustment leads to an inertially oscillating slanted profile with alongfront velocities in thermal wind balance. Tandon and Garrett (1994) assume that the oscillations continue until another mixing event removes the vertical stratification, leading to a new adjustment process. Young (1994) develops a ML model to capture the essence of this process. In this model fronts can undergo geostrophic adjustment, but mixing is fast enough to prevent further development of alongfront baroclinic instabilities.

Whether or not Young’s approximation is appropriate depends largely on the time scales involved. Consider an ML hit every τ time units by a sudden storm that erases all vertical stratification. The adjustment described by Tandon and Garrett occurs over an inertial time scale f−1 and results in alongfront balanced velocities. Can secondary instabilities develop before the next storm hits? Young (1994) considers only mesoscale instabilities, so that Tτf−1. However, MLIs develop quickly between mixing events, so that τTf−1. The adjustment described by Tandon and Garrett occurs, but the adjusted profile is itself unstable to fast baroclinic instabilities that further restratify the ML before another mixing event occurs. In this paper we show that the bulk of the restratification is associated with MLIs and not with the geostrophic adjustment.

The structure of the paper is as follows. Section 2 introduces the basic ML configuration used in this study. Section 3 illustrates the linear stability analysis of MLIs, while section 4 describes the nonlinear evolution of MLIs and their effect on ML restratification. Section 5 presents the conclusions.

2. Mixed layer instabilities: Linear equations and scalings

The QG analysis presented in the introduction shows that the ML is host to MLIs, smaller and faster than interior mesoscale instabilities. We now extend the analysis to account for the ageostrophic effects that develop as a result of the weak stratification in the ML. The scale separation (Fig. 3) allows an ansatz of fixed large- and mesoscale fields during the submesoscale evolution—effectively a Wentzel–Kramers–Brillouin (WKB) approximation. In other words we consider how the mesoscale may modulate the MLI, but not vice versa. For linear MLI the approximation holds, but full coupling across scales is likely during the nonlinear stages (Spall 1997), and it is a problem for a future study.

The study of baroclinic instability at low Ri is not novel. However, the ocean ML requires a combination of features: a moving interface, possibly with a mean tilt, separating the ML from the stratified interior. A reasonable starting point for this study is the nonhydrostatic Eady-like problem at low Ri formulated by Stone (1971). The dynamics is linearized around the basic state shown in Fig. 4. In the ML a linear vertical shear is in thermal wind balance with a constant horizontal buoyancy gradient. The vertical stratification is constant. This mean state represents conditions in the ML immediately after the gravitational adjustment described by Tandon and Garrett (section 4 shows that excited inertial oscillations are inconsequential). Appropriately, the mean state potential vorticity is 0—a condition believed to occur in the ML after strong mixing (Haine and Marshall 1998; Thomas 2005). The main novelty of the analysis is the addition of a moving interface at z = −H(y) that separates the ML from the stratified interior through a buoyancy jump ΔB and an associated velocity jump ΔU (via the Margules relation). This boundary condition represents the sudden jump in stratification and shear typically found at the ML base and visible in Fig. 2. Equations and boundary conditions are discussed in detail in appendixes B and C.

Fig. 4.

Schematic of the mean state. The ML is represented as a fluid with a linear stratification in y and z. A jet with constant vertical shear U(z) is in thermal wind balance with the horizontal stratification By. The ML base is represented by a density jump between the weakly stratified boundary layer and the stratified ocean interior. The ML base has a mean tilt Hy.

Fig. 4.

Schematic of the mean state. The ML is represented as a fluid with a linear stratification in y and z. A jet with constant vertical shear U(z) is in thermal wind balance with the horizontal stratification By. The ML base is represented by a density jump between the weakly stratified boundary layer and the stratified ocean interior. The ML base has a mean tilt Hy.

Assuming that the buoyancy jump at the ML base is larger than the buoyancy variations within the ML, the nondimensional basic state for the velocity and buoyancy in the ML is

 
formula

Without loss of generality we choose a reference system moving at a constant velocity ΔU, thereby confining the jump in velocity to the ocean interior and keeping only the sheared component in the ML equations. The nondimensionalization follows Stone (1971) and is summarized in Table 1. The external parameters are the inertial time scale f−1, the ML velocity scale U, the mean ML depth H0, and the ML buoyancy frequency N2. The associated nondimensional parameters are the bulk Richardson number, Ri = N2H2/U2, and the Prandtl ratio, δ = f/N. In this nondimensionalization, the vertical stratification is O(1), while the horizontal stratification is of O(Ri−1).

Table 1.

Dimensional scaling for the variables and parameters of the frontal stability analysis. Asterisks represent dimensional variables, f is the inertial frequency, U is the characteristic along-frontal velocity, H0 is the mean ML depth, and N2 is the buoyancy frequency in the ML.

Dimensional scaling for the variables and parameters of the frontal stability analysis. Asterisks represent dimensional variables, f is the inertial frequency, U is the characteristic along-frontal velocity, H0 is the mean ML depth, and N 2 is the buoyancy frequency in the ML.
Dimensional scaling for the variables and parameters of the frontal stability analysis. Asterisks represent dimensional variables, f is the inertial frequency, U is the characteristic along-frontal velocity, H0 is the mean ML depth, and N 2 is the buoyancy frequency in the ML.

Assuming perturbations of the form eikx, the linearized equations for this basic state are

 
formula
 
formula
 
formula
 
formula
 
formula

Deviations from hydrostatic equilibrium are proportional to δ2. In the limit Ri → ∞ and δ ≪ 1, the system (5)(9) recovers the QG Eady (1949) problem. However, in the ML Ri = O(1) and ageostrophic effects modify the instability (Nakamura 1988). Notice that, for a typical ML value of N = 10−3 s−1, δ = O(10−1). Only during times of vigorous mixing, in response to synoptic storms or cold air outbreaks, does stratification become so weak that δO(1). This study focuses on periods when mixing is not actively deepening the ML, and thus we assume that δ is small. Bulk models that assume that the ML is perfectly mixed (δ → ∞) are appropriate to study ML deepening, but not ML restratification.

At the ocean surface, we use a rigid-lid boundary condition (i.e., w = 0 at z = 0). At the ML base, at z = −H(y)/H0, we impose the kinematic boundary condition for a material interface and the dynamic boundary condition of vanishing total pressure gradient to ensure that the ocean interior is nonaccelerated. The complete boundary conditions are derived in appendix C. Under the same approximations used for the dynamics in the ML, the boundary conditions reduce to

 
formula

where η is the nondimensional displacement of the lower boundary from its mean depth.

3. Mixed layer instabilities: Linear stability analysis

Stone (1970) studied the linear stability analysis for the set of (5)(9) with rigid lids (i.e., w = 0 at z = 0 and z = −1) instead of (10). In the absence of a free interface, one can substitute solutions of the form eily+iσt in (5)(9) and solve the eigenproblem for σ. Stone found that there are four different types of instability in the ageostrophic Eady problem, as follows.

  1. Baroclinic mode: The most unstable wavenumber has a large meridional wavenumber and a vanishing zonal wavenumber (k → 0, l = 0).

  2. Convective mode: When Ri ≤ 0 the system is convectively unstable for all wavenumbers. These modes bring the vertical stratification rapidly back to neutral stability (i.e., Ri ≥ 0; see Haine and Marshall 1998).

  3. Symmetric mode: This mode arises when large along-isopycnal angular momentum gradients trigger centrifugal instabilities. This happens when the potential vorticity is negative or 0 < Ri ≤ 1. The most unstable wavenumber has a large meridional wavenumber and a vanishing zonal wavenumber (k → 0, l → ∞). The growth rate increases for increasing l and exhibits ultraviolet catastrophe. The symmetric mode dominates the instability and grows until Ri approaches 1, an effective but minimal restratification (e.g., Thorpe 1998). Numerical simulations confirm that once Ri exceeds 1, the symmetric mode dies, and the baroclinic mode takes over (Haine and Marshall 1998).

  4. Inertial critical layer mode: The existence of a critical layer mode in the short-wave region (k → ∞, lk) is a characteristic of the ageostrophic instability problem (Nakamura 1988). Unlike the geostrophic mode, which is a resonance between two boundary waves, this mode is produced by resonance between a boundary wave and an inertia–gravity wave trapped at an inertia critical level (Jones 1967). We find that similar modes arise from resonance between a boundary wave and a gravity wave along the ML base. However, these modes appear for ML-base slopes steeper than allowed by the approximations used to derive the lower boundary condition in (10) and steeper than expected for mesoscale heaving.

Stone’s analysis suggests that restratification proceeds in two distinct stages. For Ri < 1, convective and symmetric instabilities develop rapidly and stabilize the basic state by bringing Ri to unity. For Ri > 1, the baroclinic mode dominates the instability and drives restratification by releasing the potential energy in the horizontal stratification. If we allowed for curvature in the mean shear, Kelvin–Helmholtz instabilities could also develop for 0 < Ri ≤ 0.25 and modify the transition to Ri = 1. Regardless of the processes that bring Ri to 1, the bulk of restratification is clearly associated with the instabilities that develop for Ri ≥ 1. This conclusion is further supported by the fully nonlinear simulations presented in the next sections. Hence, we focus the linear stability analysis on the baroclinic instabilities that develop at Ri ≥ 1.

The full ageostrophic analysis confirms that the ML baroclinic mode resembles the shallow baroclinic instability in the introduction. The growth rate as a function of zonal wavenumber for Ri = 2 and l = 0 is shown in Fig. 5. Increasing l reduces growth rates as in the QG problem. This mode draws its energy from the potential energy of the basic state and resembles the conventional Eady mode. Stone (1966) derived an analytical approximation for the growth rate as a function of k and the most unstable wavenumber kmax by expanding the eigenvalue problem for small k:

 
formula

with

 
formula

Stone’s solution is a good approximation of the full instability problem (Fig. 5). Ageostrophic effects shift the instability to longer spatial and time scales compared to the QG problem. This effect has a simple explanation (Nakamura 1988). The baroclinic mode requires the coupling of edge waves counterpropagating at the upper and lower boundaries. Larger shears (i.e., small Ri) make the waves shallower and weakens the coupling. In order for the waves to maintain a sufficient depth to interact, a shift to longer scales is necessary because the penetration depth of the waves is proportional to their length scale. Thus, strong shears tend to stabilize the short baroclinic waves through ageostrophic effects. Herein we will refer to ML ageostrophic baroclinic instabilities as MLIs.

Fig. 5.

Growth rate of the baroclinic instability at an ML front between rigid plates with Ri = 2. The three curves correspond to estimates based on the approximate analytical expression of Stone (1966) given in (11), the QG approximation of Eady (1949), and a numerical integration of the eigenvalue problem.

Fig. 5.

Growth rate of the baroclinic instability at an ML front between rigid plates with Ri = 2. The three curves correspond to estimates based on the approximate analytical expression of Stone (1966) given in (11), the QG approximation of Eady (1949), and a numerical integration of the eigenvalue problem.

a. The baroclinic mode and the ML base

The ML base is a deformable interface rather than a solid wall and can affect the MLI by changing the propagation speed of the edge waves at the ML base or by introducing a large potential vorticity gradient associated with the ML base tilt. Analysis of the Ri ≫ 1 limit suggests that these effects substantially modify MLIs (Young 1994; Ripa 1995, 2001). On the other hand, we find that in the ocean ML, where Ri is typically small, the moving interface stabilizes large-scale MLIs but does not affect significantly the most unstable wavenumbers.

Plugging solutions of the form eily+iσt into (5)(10) reduces the problem to an eigenvalue problem in z,

 
formula

with boundary conditions,

 
formula

Solutions to the eigenvalue problem in (12) with homogeneous boundary conditions appear elsewhere (e.g., Moore and Peltier 1987; Nakamura 1988; Fukamachi et al. 1995). The boundary conditions in (13) are less tractable, so (5)(9) are integrated forward in time to find dominant eigensolutions. We solve the eigenvalue problem for two illustrative ML fronts with Ri = 2 and 1/2, Prandtl ratio δ = 0.1. The buoyancy jump at the ML base is set to ΔB = 10. We consider ML-base tilts Hy ranging from −0.05 to 0.05 corresponding to 50 m over 10 km, a generous upper bound for doming due to mesoscale eddies, which are typically an order of magnitude smaller. Also, for |Hy| > 0.1 the bulk Richardson number associated with the density jump at the ML base is less than 1 and shear instabilities are likely to develop at the ML base.

Figure 6 shows the growth rate as a function of zonal wavenumber k and the ML-base tilt Hy for Ri = 2; we set l = 0 for the moment. The fastest-growing mode is around k = 1 and is strongly ageostrophic for all bottom slopes Hy. For U ∼ 0.1 m s−1 and f ∼ 10−4 s−1 this mode has a wavelength of about 6 km and a growth time scale of 0.8 days. There is a sharp high wavenumber cutoff for all Hy. As in the Eady problem, the instability is cured once the wavenumber becomes large enough that the vertical structure of the edge waves can no longer interact. The slight changes in cutoff scale are due to modifications of the phase speed of the edge waves propagating along the tilted boundary (Blumen 1979). There is no ultraviolet catastrophe (i.e., unbounded growth rates with increasing wavenumbers), which emerges in the linear stability analysis of bulk ML models (Ripa 1995; Young 1994); these models cannot reproduce the exponential trapping of the boundary waves that cure the instability at large k (Ripa 2001). The main effect of a tilt in the ML base is to suppresses instabilities at small wavenumbers. The structure of the instability resembles an Eady mode (Fig. 7) with edge waves traveling along the boundaries. Notice that the edge wave at the bottom boundary has contributions from both the buoyancy anomaly along the free interface and the interface displacement.

Fig. 6.

(left) Growth rate for Ri = 2, δ = 0.1, and ΔB = 10, as a function of the ML-base tilt and of alongfront wavenumber k. (right) Growth rate for Ri = ½, δ = 0.1, and ΔB = 10, as a function of the ML-base tilt and of alongfront wavenumber k. The growth rate of the instability vanishes at the k = 2 boundary. In the nondimensionalization used in the paper k = Ro. The contour interval is 0.01.

Fig. 6.

(left) Growth rate for Ri = 2, δ = 0.1, and ΔB = 10, as a function of the ML-base tilt and of alongfront wavenumber k. (right) Growth rate for Ri = ½, δ = 0.1, and ΔB = 10, as a function of the ML-base tilt and of alongfront wavenumber k. The growth rate of the instability vanishes at the k = 2 boundary. In the nondimensionalization used in the paper k = Ro. The contour interval is 0.01.

Fig. 7.

Structure of the most unstable mode for Ri = 2, δ = 0.1, ΔB = 10, k = 0.9, and Hy = 0.04. (top to bottom) Buoyancy b, velocity in cross section of the front (u, w), cross-front velocity υ, and normalized interface displacement η.

Fig. 7.

Structure of the most unstable mode for Ri = 2, δ = 0.1, ΔB = 10, k = 0.9, and Hy = 0.04. (top to bottom) Buoyancy b, velocity in cross section of the front (u, w), cross-front velocity υ, and normalized interface displacement η.

The results for Ri = ½ are qualitatively similar to those for Ri = 2, but the maximum growth occurs at a slightly higher wavenumber and is faster by about 20% (Fig. 6). No Kelvin–Helmholtz instabilities develop, because the shear in the ML is constant (Vanneste 1993). Changes in the density jump at the ML base ΔB have minor effects: for 1 ≤ ΔB ≤ 20, the scale and growth rate of the most unstable wavenumber change by less than 5%. The only appreciable effect of the density jump is to suppress growth at large scales.

Additional calculations with cross-front structure (l ≠ 0) were performed (not shown). For Ri > 1, l ≠ 0 stabilizes the flow and decreases the growth rate, as in QG baroclinic instability (e.g., Pedlosky 1987). Hence, the fastest growth occurs at l = 0. For Ri < 1, the problem is modified by the appearance of symmetric instabilities, which require l ≠ 0. Nonetheless, Fig. 6 captures the fastest-growing baroclinic modes even for Ri = ½.

Figure 6 shows that the ML base tilt tends to suppress long-wave MLIs. This can be explained by inspection of the eigenvalue problem. Let us expand (12)(13) in powers of k2 for k ≪ 1 (equivalent to a small Ro expansion):

 
formula

The moving boundary modifies the instability problem for waves of length comparable to the deformation radius associated with the density jump at the ML base k−1B = RiΔB and of steepness comparable to Hy. For longer and steeper waves the bottom boundary becomes effectively rigid and flat. To focus on such waves we assume that k2kB2Hy in the expansion of the bottom boundary condition. To leading order in small k we obtain

 
formula
 
formula

Instabilities arise for Δ ≤ 0. For a flat and rigid bottom, Hy = kB = 0, this expression reproduces (11). The presence of an interface modifies the instability in two ways. First, the moving interface slightly reduces the growth rates as some energy goes into interfacial waves, but this effect is weak. Second, the bottom slope acts as an effective β and suppresses the instability at low wavenumbers. The low-wavenumber cutoff for Hy ≠ 0 predicted by (16) is evident in Fig. 6. The flow is more stable when Hy < 0 than when Hy > 0, because it corresponds to isopycnal tilts opposing the slope of the ML base. The main point here is that tilts in the ML base and the associated shear can explain the separation between mesoscale and ML instabilities shown in Fig. 3.

Last, we want to comment on an issue that has plagued the literature on the dynamics of ML instabilities. The instability of bulk and subinertial ML models depends on the existence of a tilt in the base of the ML (Young and Chen 1995; Ripa 2001). This instability arises for large Ri and small Ro and persists for arbitrarily large wavenumbers leading to an ultraviolet catastrophe. The growth rates shown in Fig. 6 do not display an ultraviolet catastrophe. As the tilt is increased the growth rate is modified only quantitatively, increasing for positive tilt and decreasing for negative tilt. Furthermore, the sensitivity to tilt decreases as Ri decreases. In conclusion, there are indeed submesoscale instabilities in the ML, but these are the result of ageostrophic baroclinic instability, not the result of QG instabilities interacting with the ML base.

b. The baroclinic mode and frontal width, horizontal strain, and vertical mixing

The preceding analysis has infinite frontal extent, no horizontal strain, and no vertical mixing. The results of introducing these effects may be anticipated from applying results in the literature to MLIs. A finite frontal width Lf introduces a lower bound on l (lπ/Lf) and hence reduces the growth rate from the l = 0 maximum. This effect is better quantified in terms of a Burger number Bu,

 
formula

The larger the Bu, the more stable the front. Eldevik and Dysthe (2002) solve numerically the ML linear instability problem for several combinations of Bu and Ri for a flat base and ΔB = O(10), and find stabilization as Bu goes from 0 to 1: the most unstable wavenumber increases by about 10% while the growth rate drops by nearly a factor of 2. The ML fronts with Bu ≥ 1 have horizontal scales comparable with those of the three-dimensional turbulent eddies that mix the ML and cannot be studied in isolation from ML turbulence. Here we focus on the Bu < 1, which seems to be the most relevant for restratification.

Convergences in the mesoscale horizontal strain field can also modify MLIs. Bishop (1993) shows that linear baroclinic instabilities are stabilized by a two-dimensional strain field (−αx, αy), as their wavenumber grows as eαt and shortens the penetration depth of the boundary waves. Characteristic values of α for mesoscale flows range between 10−7 and 10−6 s−1 (Spall 1997; Ledwell et al. 1993) and are much smaller than the growth rates of MLIs. Thus, we expect MLIs to reach finite amplitude before appreciable modification of the instability. The situation is quite different for convergences driven by winds and frontogenetic processes. These can be as small as the growth rates of MLIs and hence we expect MLIs to be modulated by the surface strain field in the open ocean. We return to this issue in the conclusions.

Fukamachi et al. (1995) study the overall effect of mixing of buoyancy and momentum on MLIs. They represent mixing through enhanced eddy diffusivities and viscosities and find that the growth rates of the baroclinic mode change by less than 10% for vertical viscosities/diffusivities as large as 1.0 m2 s−1. In section 5 we confirm these results with simulations with a more realistic representation of turbulent mixing.

4. Mixed layer instabilities: Nonlinear simulation

The linear stability analysis shows that MLIs are narrow bands with scales of O(1) km and growth rates of O(1) day. This supports our hypothesis that MLIs can develop between mixing events and drive ML restratification. Strong mixing events occur at most daily (convection due to the diurnal cycle) or more rarely (wind- and buoyancy-driven mixing associated with synoptic atmospheric events like storms). Observations support our analysis; rapidly-growing, small-scale, wavelike patterns along ML fronts are commonplace. Flament et al. (1985) show a remarkable example in the Advanced Very High Resolution Radiometer (AVHHR) images of an upwelling front off the California coast. Wave disturbances on the order of a few kilometers are visible along the flanks of the 100-km-wide upwelling filament. In addition, these features have an e-folding time of about 1 day. Munk et al. (2000) show that eddies of a few kilometers in lateral extent populate ocean surface pictures from the Apollo Mission 30.

To illustrate the finite-amplitude development of MLIs and quantify their role in ML restratification, we simulate the adjustment of an ML front in a reentrant channel on the f plane. We use the Massachusetts Institute of Technology general circulation model (MITgcm; Marshall et al. 1997a) in a channel configuration of 100 km × 200 km × 300 m. The ML is 200 m deep and sits over a stratified interior as shown in Fig. 8a. The initial buoyancy profile represents a ML front after the passage of a storm as described in the introduction. It consists of a horizontal temperature gradient of 2.5°C across 50 km. The parameter values are based on observed fronts in the upper ocean (Price 1981; Ferrari and Rudnick 2000). We use a vertical viscosity of 5 × 10−3 m2 s−1 to account for wind-driven turbulence in the ML. With this value the Ekman layer is explicitly resolved by the model mesh grid of 5 m (he2ν/f ≈ 12 m). The front starts from rest and is allowed to adjust under the sole effect of gravity. This is the reference simulation. Different configurations are used to study the dependence of the frontal adjustment on the strength and width of the front. Details about initial conditions and parameters used are given in appendix D.

Fig. 8.

Development of mixed layer baroclinic instabilities along a temperature front undergoing geostrophic adjustment. (a) The initial configuration consists of a lateral temperature front in a well-mixed surface layer on top of stable density stratification. (b) After 10 days the front has tilted from the vertical to the horizontal and wavelike disturbances appear along the front. The tilt of the wave disturbances in the along-channel direction is such as to release the potential energy stored in the horizontal stratification much like in the Eady problem. (c) By day 12 the disturbances are fully nonlinear and start growing in scale as a result of an inverse cascade of energy. (d) At day 17 the disturbances have wrapped up into eddies and frontogenesis develops along the rim of the eddies. The color bar is in degrees Celsius, and the contour interval is 0.25°C.

Fig. 8.

Development of mixed layer baroclinic instabilities along a temperature front undergoing geostrophic adjustment. (a) The initial configuration consists of a lateral temperature front in a well-mixed surface layer on top of stable density stratification. (b) After 10 days the front has tilted from the vertical to the horizontal and wavelike disturbances appear along the front. The tilt of the wave disturbances in the along-channel direction is such as to release the potential energy stored in the horizontal stratification much like in the Eady problem. (c) By day 12 the disturbances are fully nonlinear and start growing in scale as a result of an inverse cascade of energy. (d) At day 17 the disturbances have wrapped up into eddies and frontogenesis develops along the rim of the eddies. The color bar is in degrees Celsius, and the contour interval is 0.25°C.

a. Finite-amplitude development of MLI

Figure 8 show a series of 3D snapshots of the front during the adjustment. Figure 9 shows the corresponding along-channel mean stratification. Within the first day the initially vertical isopycnals start oscillating around the geostrophically adjusted state with N2By2/f2 and UzBy/f, corresponding to Ri = N2/Uz ≈ 1, as predicted by Tandon and Garrett.2 There are only minor changes in stratification during this phase. However, the gravitational adjustment is an important step toward restratification, because it sets the stage for the rapid growth of MLIs that only develop at Ri ≥ 1. After a few days, MLIs are visible as wavelike disturbances along the front (Fig. 8b). The bulk of restratification begins after day 10 when MLIs reach finite amplitude (Fig. 9). The linear phase lasts 10 days, because we seed the initial conditions with infinitesimal perturbations. In the real ocean, fronts would always include large perturbations; hence, the MLIs would reach finite amplitude immediately after the geostrophic adjustment and commence strong restratification.

Fig. 9.

Increase in domain-averaged buoyancy frequency N2 as a result of slumping of the ML front shown in Fig. 8. The initial vertical stratification is 0. The insets show snapshots of the various stages of the along-channel average of buoyancy. The initial slumping oscillates on the inertial period (h 0–24). It is followed by a restratification due primarily to the growth of baroclinic MLIs (days 2–10) and then by the eddies resulting from the nonlinear interaction of the MLIs (day 10 onward). MLI perturbations are infinitesimal until day 10 and thus N2 is seen to simply oscillate around the geostrophic adjusted state with N2by2/f2. Only once MLIs reach finite amplitude does the increase in N2 become significant.

Fig. 9.

Increase in domain-averaged buoyancy frequency N2 as a result of slumping of the ML front shown in Fig. 8. The initial vertical stratification is 0. The insets show snapshots of the various stages of the along-channel average of buoyancy. The initial slumping oscillates on the inertial period (h 0–24). It is followed by a restratification due primarily to the growth of baroclinic MLIs (days 2–10) and then by the eddies resulting from the nonlinear interaction of the MLIs (day 10 onward). MLI perturbations are infinitesimal until day 10 and thus N2 is seen to simply oscillate around the geostrophic adjusted state with N2by2/f2. Only once MLIs reach finite amplitude does the increase in N2 become significant.

MLIs grow following the linear predictions for the first 10 days. Figure 10 shows the growth of the basin-averaged perturbation kinetic energy. Perturbations are defined as departures from a zonal mean along the channel. The e-folding growth rate is σ ≈ 1.1 day−1 within 10% of the predictions of linear stability analysis in (11) with Ri = 1.2. The dominant wavelength of the alongfront wave disturbance shows up as the scale where the spectrum of perturbation kinetic energy peaks (Fig. 11a). The scale is approximately 4 km, very close to the prediction of linear stability analysis. The agreement between linear theory and simulation is very good despite the presence of large inertial oscillations (Fig. 9). This is consistent with previous work suggesting that the balanced motions associated with baroclinic instability are weakly influenced by inertial and superinertial waves (e.g., Reznick et al. 2001).

Fig. 10.

Evolution of the domain-integrated perturbation kinetic energy E(t), where perturbations are departures from a zonal average. Results are shown for a reference simulation with no diurnal buoyancy forcing (thick black line), a simulation with diurnal buoyancy forcing and a convective adjustment scheme (dashed line), and a simulation with diurnal buoyancy forcing and a KPP scheme (dot–dashed line). The evolution of the perturbation kinetic energy is compared to the prediction of linear stability analysis. The curve labeled Stone Max is E(t) = e2σmaxtE(t = 3 days), where σmax = σ(kmax) is the growth rate of the most unstable mode for mixed layer baroclinic instabilities in (11) with Ri = 1.2. Stone Int is E(t) = ∫e2σ(k)tE(k, t = 3 days) dk and takes into account that different wavenumbers have different growth rates. The two estimates converge and agree with the nonlinear simulation between days 5 and 10 when the most unstable wavenumber dominates the perturbation kinetic energy (see Fig. 8).

Fig. 10.

Evolution of the domain-integrated perturbation kinetic energy E(t), where perturbations are departures from a zonal average. Results are shown for a reference simulation with no diurnal buoyancy forcing (thick black line), a simulation with diurnal buoyancy forcing and a convective adjustment scheme (dashed line), and a simulation with diurnal buoyancy forcing and a KPP scheme (dot–dashed line). The evolution of the perturbation kinetic energy is compared to the prediction of linear stability analysis. The curve labeled Stone Max is E(t) = e2σmaxtE(t = 3 days), where σmax = σ(kmax) is the growth rate of the most unstable mode for mixed layer baroclinic instabilities in (11) with Ri = 1.2. Stone Int is E(t) = ∫e2σ(k)tE(k, t = 3 days) dk and takes into account that different wavenumbers have different growth rates. The two estimates converge and agree with the nonlinear simulation between days 5 and 10 when the most unstable wavenumber dominates the perturbation kinetic energy (see Fig. 8).

Fig. 11.

Spectra of the perturbation kinetic energy as a function of zonal wavenumber. (left) Spectra for days 3, 6, and 12 for the reference simulation. Corresponding spectra for simulations with diurnal buoyancy forcing and (middle) a convective adjustment scheme and (right) a KPP scheme. The solid lines mark the most unstable wavenumber predicted by inviscid linear theory.

Fig. 11.

Spectra of the perturbation kinetic energy as a function of zonal wavenumber. (left) Spectra for days 3, 6, and 12 for the reference simulation. Corresponding spectra for simulations with diurnal buoyancy forcing and (middle) a convective adjustment scheme and (right) a KPP scheme. The solid lines mark the most unstable wavenumber predicted by inviscid linear theory.

To assess the role of convective mixing on MLIs, we repeat the frontal adjustment experiment adding a surface buoyancy flux (details in appendix D). We impose a cycle of 200 W m−2 nighttime surface cooling compensated by daytime penetrating short-wave warming, so that the daily averaged heat flux is 0. We run two simulations, one in which convectively forced turbulence is parameterized with the convective adjustment scheme, the other in which it is represented with a K-profile parameterization (KPP; Large et al. 1994). Turbulent convection develops every night and reaches down to the ML base in both simulations. During the day a warm layer forms at the surface and convection stops. The ML depth remains essentially unchanged over the course of the simulations.

The development of MLIs is only weakly affected by the enhanced nighttime turbulence. In the simulation using convective adjustment the growth rate and dominant wavelength of the MLIs remain virtually unchanged (Figs. 10 and 11b). The simulation using KPP is slightly different. First, KPP computes diffusivities based on a local Ri criterion and generates horizontal patchiness in the buoyancy and momentum fields. This patchiness produces a large initial perturbation kinetic energy and the most unstable MLIs emerge already at day 4 (Fig. 10). Second, KPP mixes both momentum and buoyancy, while convective adjustment acts only on buoyancy. As a result the growth rate and dominant wavelength of MLIs are somewhat reduced (Figs. 10 and 11c). The bottom line is that the temporal and spatial scales of MLIs remain within a factor of 2 of the inviscid theory.

We further tested the linear theory by running simulations for a wide range of frontal strengths between 5 × 10−3 and 5 × 10−2 °C km−1 and widths between 4 and 20 km. In all cases, MLIs develop with scales in agreement with linear theory. MLIs are suppressed only in simulations that include diurnal forcing and start from a front in geostrophic balance with Ri ≥ O(100). For such fronts, the growth rate of the baroclinic mode is much longer than a day, that is, much longer than the time scale of vertical mixing of O(1) day.

The agreement between theory and simulations during the linear stages of the adjustment confirms that the frontal collapse is driven by MLIs. However, the bulk of the restratification begins only after day 10, when the instability reaches finite amplitude and causes a fivefold increase of stratification in 6 days (Fig. 9). Once the instability becomes fully nonlinear the horizontal scales of the most energetic eddies are larger than that of the most unstable mode as a result of a turbulent inverse energy cascade. This can be seen in the last two snapshots of the slumping front shown in Fig. 8. The spectra in Fig. 11 confirm that during the first 10 days the increase of perturbation energy is associated with the growth of the most unstable wavenumber, while afterward perturbation energy is transferred to larger scales. The inverse cascade proceeds also in the vertical with eddies reaching scales comparable with the ML depth. To gain insight in the finite-amplitude stages of MLIs, it is very useful to consider in detail the potential vorticity and potential energy budgets of the slumping front.

b. The potential vorticity budget

The full Ertel potential vorticity (PV) is given by P = ωa · ∇b, where ωa = f + ∇ × u is the absolute vorticity. Changes in PV result from convergences/divergences of the PV flux J:

 
formula

The PV flux has an advective component uP and a nonadvective component that arises from diabatic heating ℬ and from frictional forces F. Advective processes, like MLIs, can achieve restratification by rearranging PV in such a way that fbz increases in the ML. We want to investigate how this happens. We consider the reference simulation with no diabatic forcing first.

Integrating (18) over the full horizontal domain and between two vertical levels zt and zb, and for times between 0 and t, we obtain (Thomas 2005)

 
formula

where the overbar denotes an horizontal integral and Jz is the vertical component of J. The term on the lhs of the equation will be referred to as the PV gain. The PV fluxes through the upper and vertical boundaries Jz dominate over the fluxes through the lateral boundaries,3 which are therefore neglected. In Fig. 12 we show the PV gain for three different layers: an upper layer that includes the surface Ekman layer (0–50 m), a middle layer that spans the core of the ML (50–150 m), and a lower layer that spans from the ML base to the stratified interior (150-m bottom). During the first 10 days, when MLIs are linear, PV increases in the upper and lower layers, but not in the middle layer. In this phase MLIs are weak and the advection of PV by the mean jet is 0, because the geostrophic velocity is directed horizontally along PV surfaces. Inertial oscillations are quite linear and do not transport PV integrated in time over a period. Nonadvective PV fluxes are large at the boundaries, as a result of the mechanical stresses and heat fluxes driven by the boundary conditions, but weak in the interior. Hence, PV changes only in the upper and lower layers and the stratification in the core of the ML remains essentially constant (Fig. 9).

Fig. 12.

The PV gain during the frontal collapse for the reference simulation shown in Fig. 8. The PV gain ΔPV is defined as the increase in PV in a specific fluid volume with respect to the initial value. PV gain for the full domain (continuous thick line). PV gain for an upper layer between 0 and 50 m (continuous thin line). PV gain for a lower layer between 150 m and the bottom of the channel (dot–dashed line). PV gain in a middle layer between 50 and 150 m spanning the bulk of the mixed layer (dashed thick line). The gain of fbz, i.e., the gain of PV without the vorticity and baroclinic terms (dashed thin line).

Fig. 12.

The PV gain during the frontal collapse for the reference simulation shown in Fig. 8. The PV gain ΔPV is defined as the increase in PV in a specific fluid volume with respect to the initial value. PV gain for the full domain (continuous thick line). PV gain for an upper layer between 0 and 50 m (continuous thin line). PV gain for a lower layer between 150 m and the bottom of the channel (dot–dashed line). PV gain in a middle layer between 50 and 150 m spanning the bulk of the mixed layer (dashed thick line). The gain of fbz, i.e., the gain of PV without the vorticity and baroclinic terms (dashed thin line).

The PV starts growing in the core of the ML after day 10, once MLIs reach finite amplitude. Fully developed eddies generate a strong circulation in the ML that advects high PV from the stratified interior and the surface Ekman layer (Fig. 12). The injection of high PV filaments along buoyancy surfaces by eddy stirring is evident in snapshots of the PV field (Fig. 13). This injection of high-PV fluid results in ML restratification. In Fig. 12 we compare the middle-layer PV gain with the increase in fbz. Apart from the inertial oscillations for the first 10 days, which appear in fbz but not in the full PV, the two curves are very similar, suggesting that interior PV is a good proxy for stratification. This is not to say that finite Ro and Ri effects are weak: PV is dominated by fbz only in an averaged sense. Locally, at fronts, relative vorticity is important. To summarize, MLIs drive restratification by injecting into the ML highly stratified waters from the interior and the surface Ekman layer.

Fig. 13.

Three snapshots of the PV distribution right after the baroclinic instability reaches finite amplitude for the reference simulation shown in Fig. 8: (left) day 10, (middle) day 10.5, and (right) day 11. The contours are logarithmically spaced to allow for the widely different PV values in the stratified interior and in the surface ML. The dark black lines show two zonal-mean isopycnals that roughly bracket the width of the front.

Fig. 13.

Three snapshots of the PV distribution right after the baroclinic instability reaches finite amplitude for the reference simulation shown in Fig. 8: (left) day 10, (middle) day 10.5, and (right) day 11. The contours are logarithmically spaced to allow for the widely different PV values in the stratified interior and in the surface ML. The dark black lines show two zonal-mean isopycnals that roughly bracket the width of the front.

Figure 12 shows that PV starts growing in the upper and lower layers from day 1. Advective processes can only redistribute PV, hence diabatic or frictional mechanisms must be at work. In the reference simulation there is no surface heating, but we apply a no-stress boundary condition. In order for the shear to vanish at the boundaries, an Ekman shear ∂zue must develop to balance the geostrophic shear ∂zug (Thomas and Rhines 2002),

 
formula

where D is the full-channel depth. The Ekman transport is always directed at 90° to the right of the boundary stress on the fluid, in this case, νzug (Pedlosky 1987). Using thermal wind balance, this implies that the Ekman transport is directed down the surface buoyancy gradient and acts to restratify the fluid. One can indeed see the retratification within the surface and bottom Ekman layers in Fig. 8b. We verified that this process accounts for the PV changes in the upper and lower layers. In the nonlinear stages the frictional PV gain increases, because sharp buoyancy fronts develop at the boundaries and drive larger Ekman transports. The development of frontogenesis at the boundaries is a well-known feature of the nonlinear Eady problem (Hoskins and Bretherton 1972; Garner et al. 1992) and can be seen in Fig. 8d. The PV gain is weaker in the lower than in the upper layer, because PV is destroyed by the no-buoyancy-flux boundary condition at the bottom (Fig. 12).

We repeated the frontal collapse simulation with a reduced viscosity of ν = 5 × 10−5 m2 s−1 and the same vertical resolution. In this case the frictional injection is negligible during the linear evolution and weak afterward. This is possibly an artifact of not resolving the Ekman boundary layer depth. Despite the elimination of the Ekman PV source, the PV injection into the ML remains similar to that shown for the reference simulation. This suggests that the injection of stratified thermocline waters dominates over the injection from the Ekman layer for the range of viscosities considered.

Simulations with diurnal heating are more complicated to analyze because diabatic destruction of PV during convection and diabatic generation of PV during penetrative heating profoundly affect the PV budget. Nevertheless one can still see filaments of high PV being advected into the ML and driving an increase of PV over a few days. Also the eddy kinetic energy and structure of MLIs is very similar to that in simulations with no diurnal forcing. It appears as if the eddy-driven circulation that transports PV is independent of frictional and diabatic mechanisms, while the PV distribution on which this circulation acts is strongly affected by diabatic forcing. An understanding of MLI-driven restratification therefore hinges on determining what controls the eddy-driven circulation. We pursue this idea in the next section.

c. The eddy-driven circulation and the potential energy budget

Held and Schneider (1999) and Plumb and Ferrari (2005) show that baroclinic eddies in a channel drive an overturning circulation in the (y, z) plane that transports tracers, including PV, given by

 
formula

The overbar is a zonal average and primes are departure from that average. The strength of this circulation is proportional to the eddy release of mean potential energy (PE), . Green (1970) and Stone (1972) show that baroclinic eddies tend to release PE from the mean state (i.e., on average > 0). Hence, the overturning circulation typically has the same sign as by and acts to slump the front from the horizontal to the vertical as surmised by Gent and McWilliams (1990). We wish to test whether the release of PE by MLIs is consistent with these arguments.

In Fig. 14 we show the PE extraction from two simulations initialized with fronts of 5 × 10−3 and 5 × 10−2°C km−1 across 20 km. These values span the range of density gradients found in the open ocean (e.g., Ferrari and Rudnick 2000). The PE extraction is always positive and corresponds to surface heat fluxes between 10 and 1000 W m−2. For reference a front of 1 × 10−3°C km−1 generates a vertical heat flux of 100 W m−2 (not shown). The PE release is nearly constant once the eddies reach finite amplitude (i.e., the strength of the eddies achieves a statistically steady state that lasts until the front hits the lateral boundaries). Addition of a diurnal cycle introduces diurnal fluctuations in the release of PE, but the mean trend matches remarkably well that shown in Fig. 14. In addition, the vertical and horizontal structure of the time-averaged remains the same, suggesting that the eddy release of PE is largely independent of diurnal heating and cooling.

Fig. 14.

Vertical eddy fluxes of buoyancy (, left axis) and surface heat flux necessary to change the equivalent amount of buoyancy (cpρ/gαT, right axis) for the strong wide front (dark shading) and the weak wide front (light shading) vs time. See Tables D1 and D2 for details about the two simulations. The range shown is the basinwide mean to the basinwide maximum, so that the mean, median, etc., for only the active front region lies somewhere in the shaded region.

Fig. 14.

Vertical eddy fluxes of buoyancy (, left axis) and surface heat flux necessary to change the equivalent amount of buoyancy (cpρ/gαT, right axis) for the strong wide front (dark shading) and the weak wide front (light shading) vs time. See Tables D1 and D2 for details about the two simulations. The range shown is the basinwide mean to the basinwide maximum, so that the mean, median, etc., for only the active front region lies somewhere in the shaded region.

Table D1. Parameters for the MITgcm simulations.

Table D1. Parameters for the MITgcm simulations.
Table D1. Parameters for the MITgcm simulations.

Table D2. Derived quantities for the MITgcm simulations.

Table D2. Derived quantities for the MITgcm simulations.
Table D2. Derived quantities for the MITgcm simulations.

Before we proceed, it is useful to compare the PE extraction by MLIs with the other terms in the climatological PE budget of the ML: midlatitude surface buoyancy fluxes are O(10–100) W m−2, entrainment O(10) W m−2, and energy dissipation rates O(10–100) W m−2 (Csanady 2000). Heaving by mesoscale eddies generates vertical fluxes in the range of O(10–100) W m−2 (Cronin and Watts 1996). Thus, the release of PE by MLIs is of the same order of magnitude as other ML processes for fronts stronger than about 1 × 10−3°C km−1, a typical ML value. This confirms the results of our numerical simulations that release of PE and ensuing restratification by MLIs is strong enough to compete with the turbulent processes that mix the ML.

We now focus on deriving scalings for the PE release by MLIs. Theories for baroclinic adjustment suggest that baroclinic eddies achieve maximum release of PE by fluxing buoyancy at one-half of the angle of the mean isopycnals (e.g., Pedlosky 1987):

 
formula

MLIs are close to achieving maximum PE extraction as shown in Fig. 15. The eddy flux slope oscillates around half the isopycnal slope for most of the simulation. Oscillations about the mean value are due to inertial oscillations and are inconsequential for the baroclinic adjustment. Only for a couple of days, after the eddies reach finite amplitude, is the ratio larger than one-half and the extraction of PE is somewhat reduced (when the ratio is 1, eddy fluxes are parallel to mean isopycnals and cannot release PE).

Fig. 15.

The median (dark line) and the range between the upper and lower quartile over all grid points of the flux ratio over the isopycnal slope, −bz/by, for the reference simulation without diurnal forcing shown in Fig. 8.

Fig. 15.

The median (dark line) and the range between the upper and lower quartile over all grid points of the flux ratio over the isopycnal slope, −bz/by, for the reference simulation without diurnal forcing shown in Fig. 8.

We can now go back to relating the PE release to the overturning streamfunction. In all the simulations we run, the evolution of ψ tracks that of , because the horizontal buoyancy gradient does not change substantially during restratification—the front slumps from the horizontal to the vertical without becoming wider, hence bz increases but by remains essentially constant. Figure 14 shows that reaches a statistically steady state once the MLIs reach finite amplitude. This is the case in all the simulations we run. Hence, we expect ψ to reach a steady state. We can also infer that ψ is not affected by diurnal forcing, because vertical mixing due to surface heating/cooling does not affect either or by. This confirms our hypothesis that the eddy-driven circulation is independent of diabatic and frictional processes. However, the rate of restratification depends on external forcing, because frictional and diabatic processes modify the background PV state on which ψ acts.

Last, we can use relationship (22) together with (21) to relate ψ to the horizontal eddy fluxes:

 
formula

Green (1970), Gent and McWilliams (1990), and Haine and Marshall (1998) show that eddies generated through baroclinic instability release PE by fluxing buoyancy against the horizontal buoyancy gradient,

 
formula

where K is a positive effective diffusivity, possibly a function of the mean properties of the flow. So we can express the eddy-induced circulation as

 
formula

Applying this kind of closure in the ML has been largely unsuccessful (e.g., Large et al. 1994). Using a canonical value of K = 1000 m2 s−1 for the diffusivity, the closure scheme in (25) produces vertical buoyancy fluxes so large that they instantly restratify the ML everywhere in the ocean. The problem arises from inappropriate choices for the values of isopycnal slope and effective diffusivity. Horizontal mesoscale eddy fluxes are large with diffusivities of O(1000 m2 s−1), but drive overturning circulations proportional to the small isopycnal slopes in the thermocline (Ferrari and McWilliams 2006, manuscript submitted to J. Climate). The ML instabilities have weaker horizontal fluxes, but generate overturning circulations proportional to the steep isopycnal slopes in the ML. Applying (24) to output from our numerical simulations, we find K = O(1–100) m2 s−1 for the range of ML fronts considered. With such values, (25) drives reasonable restratification rates. In summary, MLIs dominate the vertical buoyancy fluxes and the associated release of PE in the ML. Mesoscale eddies drive larger lateral fluxes, but with weak vertical shears and therefore do not contribute much to PE release in the ML. This analysis indicates that ML restratification is primarily a submesoscale phenomenon. In a forthcoming paper we derive a parameterization for restratification driven by MLIs based on these ideas.

5. Conclusions

Mixed layers are traditionally understood in terms of one-dimensional balances: external forcings, whether diurnal or seasonal buoyancy fluxes, prevalent or synoptic winds, drive turbulent motions, which stir the vertical column and maintain a well-mixed surface layer. Horizontal inhomogeneities are traditionally attributed either to surface forcing or to stirring by the large-scale and mesoscale flow fields. Submesoscale dynamics at scales intermediate between ML turbulence and mesoscale eddies are believed to be subdominant. In this paper, we have shown that during ML restratification this view is not correct, because the bulk of PE release is associated with instabilities at the submesoscale.

Observations of submesoscale eddies with scales larger than a few kilometers and yet much smaller than the mesoscale eddy field are ubiquitous in the upper ocean (Pollard and Regier 1992; Munk et al. 2000; Thomas 2005). The main result of this paper has been to show the connection between these eddies and ML restratification. The explanation goes as follows. Lateral buoyancy gradients are generated through surface fluxes or mesoscale straining. Once formed, these gradients are dynamically active and start to slump under the action of gravity. The slumping occurs down the buoyancy gradient with dense waters flowing under light waters. Rotation strongly constrains the efficiency of this process, because thermal wind is established after a pendulum day and further slumping is arrested. However, ageostrophic baroclinic instabilities allow the restratification to continue by generating wavelike disturbances along the front that upset the large-scale thermal wind balance. Because of the weak vertical stratification in the ML, the scale of the instabilities is O(1) km and the growth rate O(1) day, values much smaller than for mesoscale instabilities. The ML baroclinic instabilities rapidly reach finite amplitude and tilt isopycnals from the vertical toward the horizontal (i.e., they restratify the ML). This process is similar to the mesoscale baroclinic adjustment described by Green (1970) and Stone (1972) for the ocean interior, but it proceeds much faster.

The implications of the fast time scale of baroclinic adjustment in the ML are profound. Previous work on ML restratification ignored the role of baroclinic instabilities by arguing that mesoscale restratification acts on time scales too slow to compete with vertical mixing (Young 1994; Tandon and Garrett 1994). MLIs however develop at the submesoscale and are fast enough to restratify between mixing events. Nonlinear numerical simulations show that finite-amplitude MLIs inject high PV thermocline waters into the ML and drive substantial restratification despite the action of vertical mixing. MLI-driven restratification is shown to be a leading-order term in the ML potential energy budget and plays an important role in determining the depth, temperature, and salinity of the ML.

We discussed in detail the properties of MLIs. While the dynamics of baroclinic instability in weakly stratified environments has been studied before, many aspects germane to ML dynamics had not been considered. In particular, we showed that tilts in the ML base suppress instabilities at large wavenumbers, but have otherwise minor effects on the most unstable baroclinic modes that drive the bulk of restratification. Previous studies seemed to suggest that ML-base tilts substantially modify the ML baroclinic adjustment, but were based on QG analysis inappropriate for submesoscale instabilities (Young and Chen 1995; Beron-Vera and Ripa 1997). Despite having limited effect on the properties of MLIs, the presence of a free interface at the ML base has important dynamical implications. Submesoscale instabilities project on the interface and drive large vertical velocities with associated entrainment/detrainment of water masses and nutrients at the ML base. Furthermore, this entrainment/detrainment has impacts on biology, because typical vertical velocities are an order of magnitude larger than those associated with the mesoscale eddies at the depth of the chlorophyll maximum.

A summary of the ML buoyancy budget takes the form

 
formula

where the overbar indicates an horizontal coarse-grain average over the width of the front. Mean and mesoscale stirring dominate the horizontal flux divergence, but MLIs dominate the vertical flux divergence, which leads to restratification. In steady state both terms are balanced by vertical turbulent mixing ∂z. To assess the overall importance of MLIs for the global ocean, we need to estimate how common are fronts strong enough to produce energetic MLIs. A partial answer comes for the ubiquitous presence of MLIs in observations. Munk et al. (2000) and Rudnick (2001) find that eddies in the upper ocean at scales between 1 and 20 km are predominantly cyclonic. Eldevik and Dysthe (2002) show that MLIs drives surface frontogenesis and produce narrow frontal zones of strong cyclonic shear and convergence. Figure 8 confirms that MLIs wind up to produce cyclonic spiral eddies with a corresponding stretching of the frontal zone. In Fig. 16 we compare PDFs of vorticity on scales of 3 km from the reference numerical simulation used in the paper and from ADCP velocity measurements measured in the ML section shown in Fig. 1. The two figures are remarkably similar with thick tails for vorticity larger than f/2. In both cases the skewness disappears at larger scales and below the ML. The ubiquitous skewness in vorticity at the submesoscale supports the idea that MLIs are a common feature in the surface ML. Notice however that frontogenesis and cyclonic shears can also be generated through wind-driven flows and mesoscale strain (Hoskins and Bretherton 1972). Further analysis is needed to univocally identify MLIs in observations.

Fig. 16.

PDF of relative vorticity divided by Coriolis parameter (ζ/f ). (a) Results from the reference numerical simulation of a slumping horizontal density front shown in Fig. 8. The PDF is estimated using ML velocity measurements at day 12. A positive skewness appears as soon as the baroclinic instability enters in the nonlinear stage, and it continues to grow. (b) Results from ADCP measurements for the North Pacific section in Fig. 1. The ADCP data are averaged in bins of 3 km in the horizontal and 8 m in the vertical. The across-track velocity is then differentiated along track to form an estimate of one term in the vertical vorticity. The term involving the across-track difference of along-track velocity is unknown. We consider only measurements in the ML, defined as the layer with potential density within 0.1 kg m−3 from the surface value. The PDF of ζ/f is calculated in bins of 0.02 s−1.

Fig. 16.

PDF of relative vorticity divided by Coriolis parameter (ζ/f ). (a) Results from the reference numerical simulation of a slumping horizontal density front shown in Fig. 8. The PDF is estimated using ML velocity measurements at day 12. A positive skewness appears as soon as the baroclinic instability enters in the nonlinear stage, and it continues to grow. (b) Results from ADCP measurements for the North Pacific section in Fig. 1. The ADCP data are averaged in bins of 3 km in the horizontal and 8 m in the vertical. The across-track velocity is then differentiated along track to form an estimate of one term in the vertical vorticity. The term involving the across-track difference of along-track velocity is unknown. We consider only measurements in the ML, defined as the layer with potential density within 0.1 kg m−3 from the surface value. The PDF of ζ/f is calculated in bins of 0.02 s−1.

Last, in this paper we focused on the role of submesoscale baroclinic instabilities in driving ML restratification. Nurser and Zhang (2000) and, more recently, Lapeyre et al. (2006) suggest that mesoscale instabilities can also act to restratify the upper ocean. However, these studies use coarse-resolution numerical simulations that do not allow for the development of MLIs. We speculate that mesoscale eddies dominate in regions of strong convergence, where MLIs are suppressed (Spall 1997), while MLIs compete and often dominate over mesoscale restratification elsewhere. However, a full answer as to how mesoscale eddies and MLIs interact to drive ML restratification will require numerical simulations that resolve both meso- and submesoscale motions and it is left for a future study.

Acknowledgments

This research was supported by NSF Grants OCE02-41528 (G. Boccaletti and R. Ferrari) and OCE03-36755 (R. Ferrari and B. Fox-Kemper). The authors thank S. Smith for providing the code to compute the quasigeostrophic stability analysis presented in the introduction and L. Thomas for helpful discussions on the PV budget of the ML.

REFERENCES

REFERENCES
Barth
,
J. A.
,
1994
:
Short-wavelength instabilities on coastal jets and fronts.
J. Geophys. Res.
,
99
,
16095
16115
.
Beron-Vera
,
F. J.
, and
P.
Ripa
,
1997
:
Free boundary effects in baroclinic instability.
J. Fluid Mech.
,
352
,
245
264
.
Bishop
,
C. H.
,
1993
:
On the behaviour of baroclinic waves undergoing horizontal deformation. II: Error-bound amplification and Rossby wave diagnostic.
Quart. J. Roy. Meteor. Soc.
,
119
,
241
267
.
Blumen
,
N.
,
1979
:
On short-wave baroclinic instability.
J. Atmos. Sci.
,
36
,
1925
1933
.
Charney
,
J. P.
,
1947
:
The dynamics of long waves in a baroclinic westerly current.
J. Meteor.
,
4
,
135
163
.
Cronin
,
M.
, and
D. R.
Watts
,
1996
:
Eddy-mean flow interactions in the Gulf Stream at 68°W. Part I: Eddy energetics.
J. Phys. Oceanogr.
,
26
,
2107
2131
.
Csanady
,
G. T.
,
2000
:
Air–Sea Interactions.
Cambridge University Press, 239 pp
.
Eady
,
E. T.
,
1949
:
Long waves and cyclone waves.
Tellus
,
1
,
33
52
.
Eldevik
,
T.
, and
K. B.
Dysthe
,
2002
:
Spiral eddies.
J. Phys. Oceanogr.
,
32
,
851
869
.
Ferrari
,
R.
, and
D.
Rudnick
,
2000
:
The thermohaline structure of the upper ocean.
J. Geophys. Res.
,
105
,
16857
16883
.
Flament
,
P.
,
L.
Armi
, and
L.
Washburn
,
1985
:
The evolving structure of an upwelling filament.
J. Geophys. Res.
,
90
,
11765
11778
.
Fukamachi
,
Y.
,
J. P.
McCreary
, and
J. A.
Proehl
,
1995
:
Instability of density fronts in layer and continuously stratified models.
J. Geophys. Res.
,
100
,
2559
2577
.
Garner
,
S. T.
,
N.
Nakamura
, and
I. M.
Held
,
1992
:
Nonlinear equilibration of two-dimensional Eady waves: A new perspective.
J. Atmos. Sci.
,
49
,
1984
1996
.
Gent
,
P. R.
, and
J. C.
McWilliams
,
1990
:
Isopycnal mixing in ocean circulation models.
J. Phys. Oceanogr.
,
20
,
150
155
.
Gill
,
A. E.
,
J. S. A.
Green
, and
A. J.
Simmons
,
1974
:
Energy partition in the large-scale ocean circulation and the production of mid-ocean eddies.
Deep-Sea Res.
,
21
,
499
528
.
Green
,
J.
,
1970
:
Transfer properties of the large-scale eddies and the general circulation of the atmosphere.
Quart. J. Roy. Meteor. Soc.
,
96
,
157
185
.
Griffies
,
S.
, and
R.
Hallberg
,
2000
:
Biharmonic friction with a Smagorinsky-like viscosity for use in large-scale eddy-permitting ocean models.
Mon. Wea. Rev.
,
128
,
2935
2946
.
Haine
,
T. W. N.
, and
J.
Marshall
,
1998
:
Gravitational, symmetric, and baroclinic instability of the ocean mixed layer.
J. Phys. Oceanogr.
,
28
,
634
658
.
Held
,
I. M.
, and
T.
Schneider
,
1999
:
The surface branch of the zonally averaged mass transport circulation in the troposphere.
J. Atmos. Sci.
,
56
,
1688
1697
.
Hoskins
,
B. J.
, and
F. P.
Bretherton
,
1972
:
Atmospheric frontogenesis models: Mathematical formulation and solution.
J. Atmos. Sci.
,
29
,
11
37
.
Jerlov
,
N. G.
,
1976
:
Marine Optics.
Elsevier, 231 pp
.
Jones
,
H.
, and
J.
Marshall
,
1997
:
Restratification after deep convection.
J. Phys. Oceanogr.
,
27
,
2276
2287
.
Jones
,
W. L.
,
1967
:
Propagation of internal gravity waves in fluids with shear flow and rotation.
J. Fluid Mech.
,
30
,
439
448
.
Lapeyre
,
G.
,
P.
Klein
, and
B. L.
Hua
,
2006
:
Oceanic restratification forced by surface frontogenesis.
J. Phys. Oceanogr.
,
36
,
1577
1590
.
Large
,
W. G.
,
J. C.
McWilliams
, and
S. C.
Doney
,
1994
:
Oceanic vertical mixing: A review and a model with a nonlocal boundary layer parameterization.
Rev. Geophys.
,
32
,
363
403
.
Ledwell
,
J. R.
,
A. J.
Watson
, and
C. S.
Law
,
1993
:
Evidence for slow mixing across the pycnocline from an open ocean tracer release experiment.
Nature
,
364
,
701
703
.
Marshall
,
J.
,
A.
Adcroft
,
C.
Hill
,
L.
Perelman
, and
C.
Heisey
,
1997a
:
A finite-volume, incompressible Navier–Stokes model for studies of the ocean on parallel computers.
J. Geophys. Res.
,
102
,
5753
5766
.
Marshall
,
J.
,
C.
Hill
,
L.
Perelman
, and
A.
Adcroft
,
1997b
:
Hydrostatic, quasi-hydrostatic, and nonhydrostatic ocean modeling.
J. Geophys. Res.
,
102
,
5733
5752
.
Moore
,
G. W. K.
, and
W. R.
Peltier
,
1987
:
Cyclogenesis in frontal zones.
J. Atmos. Sci.
,
44
,
384
409
.
Munk
,
W.
,
L.
Armi
,
K.
Fischer
, and
Z.
Zachariasen
,
2000
:
Spirals on the sea.
Proc. Roy. Soc. London
,
456A
,
1217
1280
.
Nakamura
,
N.
,
1988
:
Scale selection of baroclinic instability—Effects of stratification and nongeostrophy.
J. Atmos. Sci.
,
45
,
3253
3268
.
Nurser
,
A. J. G.
, and
J. W.
Zhang
,
2000
:
Eddy-induced mixed-layer shallowing and mixed-layer/thermocline exchange.
J. Geophys. Res.
,
105
,
21851
21868
.
Ou
,
H.
,
1984
:
Geostrophic adjustment: A mechanism for frontogenesis?
J. Phys. Oceanogr.
,
14
,
994
1000
.
Pedlosky
,
J.
,
1987
:
Geophysical Fluid Dynamics.
2d ed. Springer, 710 pp
.
Plumb
,
R. A.
, and
R.
Ferrari
,
2005
:
Transformed Eulerian-mean theory. Part I: Nonquasigeostrophic theory for eddies on a zonal-mean flow.
J. Phys. Oceanogr.
,
35
,
165
174
.
Pollard
,
R.
, and
L. A.
Regier
,
1992
:
Vorticity and vertical circulation at an ocean front.
J. Phys. Oceanogr.
,
22
,
609
625
.
Price
,
J. F.
,
1981
:
Upper ocean response to a hurricane.
J. Phys. Oceanogr.
,
11
,
153
175
.
Price
,
J. F.
,
R. A.
Weller
, and
R.
Pinkel
,
1986
:
Diurnal cycling: Observations and models of the upper ocean response to diurnal heating, cooling, and wind mixing.
J. Geophys. Res.
,
91
,
8411
8427
.
Reznick
,
G. M.
,
V.
Zeitlin
, and
M. B.
Jelloul
,
2001
:
Nonlinear theory of geostrophic adjustment. Part I: Rotating shallow-water model.
J. Fluid Mech.
,
445
,
93
120
.
Ripa
,
P.
,
1995
:
On improving a one-layer model with thermodynamics.
J. Fluid Mech.
,
303
,
169
201
.
Ripa
,
P.
,
2001
:
Waves and resonance in free-boundary baroclinic instability.
J. Fluid Mech.
,
428
,
387
408
.
Rossby
,
C. G.
,
1937
:
On the mutual adjustment of pressure and velocity distributions in certain simple current systems. I.
J. Mar. Res.
,
1
,
15
28
.
Rossby
,
C. G.
,
1938
:
On the mutual adjustment of pressure and velocity distributions in certain simple current systems. II.
J. Mar. Res.
,
2
,
239
263
.
Rudnick
,
D.
,
2001
:
On the skewness of vorticity in the upper ocean.
Geophys. Res. Lett.
,
28
,
2045
2048
.
Samelson
,
R. M.
,
1993
:
Linear instability of a mixed-layer front.
J. Geophys. Res.
,
98
,
10195
10204
.
Smagorinsky
,
J.
,
1963
:
General circulation experiments with the primitive equations: I. The basic experiment.
Mon. Wea. Rev.
,
91
,
99
164
.
Spall
,
M.
,
1997
:
Baroclinic jets in confluent flow.
J. Phys. Oceanogr.
,
27
,
381
402
.
Stone
,
P. H.
,
1966
:
On non-geostrophic baroclinic stability.
J. Atmos. Sci.
,
23
,
390
400
.
Stone
,
P. H.
,
1970
:
On non-geostrophic baroclinic stability: Part II.
J. Atmos. Sci.
,
27
,
721
726
.
Stone
,
P. H.
,
1971
:
Baroclinic stability under non-hydrostatic conditions.
J. Fluid Mech.
,
46
,
659
671
.
Stone
,
P. H.
,
1972
:
A simplified radiative-dynamical model for the static stability of rotating atmospheres.
J. Atmos. Sci.
,
29
,
405
418
.
Tandon
,
A.
, and
C.
Garrett
,
1994
:
Mixed layer restratification due to a horizontal density gradient.
J. Phys. Oceanogr.
,
24
,
1419
1424
.
Tandon
,
A.
, and
C.
Garrett
,
1995
:
Geostrophic adjustment and restratification of a mixed layer with horizontal gradients above a stratified layer.
J. Phys. Oceanogr.
,
25
,
2229
2241
.
Thomas
,
L. N.
,
2005
:
Destruction of potential vorticity by winds.
J. Phys. Oceanogr.
,
35
,
2457
2466
.
Thomas
,
L. N.
, and
P. B.
Rhines
,
2002
:
Nonlinear stratified spin-up.
J. Fluid Mech.
,
473
,
211
244
.
Thorpe
,
S. A.
,
1998
:
Nonlinear reflection of internal waves at a density discontinuity at the base of the mixed layer.
J. Phys. Oceanogr.
,
28
,
1853
1860
.
Thorpe
,
S. A.
,
2005
:
The Turbulent Ocean.
Cambridge University Press, 458 pp
.
Vanneste
,
J.
,
1993
:
The Kelvin–Helmholtz instability in a non-geostrophic baroclinic unstable flow.
Math. Comput. Model.
,
17
,
149
154
.
Young
,
W. R.
,
1994
:
The subinertial mixed layer approximation.
J. Phys. Oceanogr.
,
24
,
1812
1826
.
Young
,
W. R.
, and
L.
Chen
,
1995
:
Baroclinic instability and thermohaline gradient alignment in the mixed layer.
J. Phys. Oceanogr.
,
25
,
3172
3185
.

APPENDIX A

Analysis of a SeaSoar Section

Figure 1 shows a section of potential density taken in the subtropical North Pacific region by towing a CTD with a SeaSoar along a sawtooth pattern in the upper 320 m (Ferrari and Rudnick 2000). The ML is visible as a layer of weak stratification in the upper 100 m overlying the more stratified interior. We wish to study whether this stratification can lead to baroclinic instabilities. We use a high-resolution profile, because we are especially interested in assessing the role of the reduced ML stratification in the stability of the ocean water column. Climatologies are generally too coarse to compute stratification and shears within the ML.

The computation is performed for a 3° latitude band between 133° and 130°W longitude (Fig. 1), where both the vertical stratification and the balanced shear are approximately independent of latitude and can be represented by profiles dependent on depth only (Fig. 2). The vertical buoyancy frequency N2 = −z/ρ0 (g is the acceleration of gravity and ρ is potential density) is computed across 8-m vertical bins and averaged over the 3° latitude range. The vertical shear is computed from thermal wind Uz = y/0; U is the across-SeaSoar track velocity, f is the inertial frequency, ρy is the along-track density gradient computed on scales of 10 km to include both large- and mesoscale gradients.

The SeaSoar profile is confined to the upper portion of the water column. Both N2 and Uz profiles are extended to the ocean bottom by matching the SeaSoar estimates in the upper 320 m with estimates based on Levitus climatology in a patch of 3° × 3° centered around the region considered. The vertical stratification is obtained by taking first-order differences of density at the levels provided in the Levitus climatology. Computation of Uz is a bit more laborious. Levitus has a horizontal resolution of 1° and, hence, does not include mesoscale shears. Ferrari and Rudnick (2000) in their analysis of the SeaSoar profile find that the horizontal spectrum of the density rolls off as k−2 (k is the horizontal wavenumber) for scales between 100 km and 10 m both in the ML and in the interior. For such a spectrum, density gradients are proportional to the square root of the scale at which they are computed (i.e., density gradients across 10 km are 101/2 larger then gradients across 100 km). Hence, we estimated density gradients at 10 km simply multiplying by 101/2 the gradients computed from the Levitus climatology. Using thermal wind, we finally have Uz. The profiles are shown in Fig. 2. We feel that the matching procedure is successful for the present qualitative analysis, because there are no major discontinuities in the profiles at 320 m.

APPENDIX B

Construction of the Basic State

We assume that the ocean interior is motionless. In the ML we impose a basic velocity of the form:

 
formula

Here ΔU is the velocity jump across the ML base at z = −H(y). The second term is a constant shear between the surface and the ML base. The basic buoyancy is given by

 
formula

corresponding to a ΔB density jump across the ML base, a horizontal buoyancy variation B(y), and a constant vertical stratification N2. Integrating the hydrostatic balance from the ML base where perturbation pressure must vanish, we obtain

 
formula

We assume that the basic state is in geostrophic balance:

 
formula

We can make a few simplifying assumptions for typical ML fronts. First, we assume that the density jump across the ML base is larger than horizontal and vertical density variations within the ML (i.e., ΔBB and ΔBN2H0). This is equivalent to assuming that the deformation radius in the ML, RML = NH/f, is smaller than the deformation radius associated with the density jump at the ML base, RΔB = ΔBH/f. Second, fluctuations in ML depth are meant to represent the effect of mesoscale heaving and are therefore smaller than the mean ML depth across the O(1–10) km scales of ML instabilities (i.e., |HH0| ≪ H0; mesoscale isopycnal slopes are typically in the range of 10−4–10−3 and the implied changes in ML depth are on the order of 1 m). With these approximations, the basic state reduces to

 
formula

and we are free to set By and Hy to constants (i.e., we assume that the slope of isopycnals in the ML and the slope of the ML base are constant).

The stability of Kelvin–Helmholtz waves at the base of the ML is a topic of considerable complexity (e.g., Thorpe 2005, chapter 3). Observational and theoretical studies suggest that the ML base is Kelvin–Helmholtz stable if the bulk Richardson number associated with the density and velocity jumps at the ML base is larger than 1 (e.g., Price et al. 1986):

 
formula

We will use this criterion as a guiding line to select stable basic states. In order for Ribase to be greater than 1, we will require that

 
formula

This is easily satisfied in the ML, because the Prandtl ratio f/N is typically of order 10−1, while mesoscale isopycnal slopes are two to three orders of magnitude smaller. In the present study, we will use (B7) to constrain the maximum allowable ML base tilt. Notice that Ribase is the bulk Richardson number associated with the discontinuities at the ML base, while the Ri used in the rest of the paper refers to the shear and stratification within the ML.

APPENDIX C

The Reduced-Gravity Bottom Boundary Condition

The reduced-gravity assumption requires that the total pressure gradient vanishes at the ML base. Linearizing around z = − H(y), we obtain

 
formula

where η is the displacement of the ML base. Using the expression for 𝒫 in (B3), we obtain,

 
formula

The kinematic boundary condition is

 
formula

where wen is the entrainment velocity at the ML base.

Additional simplifications provide consistency with the approximations made in the ML interior. We ignore buoyancy fluctuations within the ML relative to the buoyancy jump at the ML base [i.e., B(y) ≪ ΔB and N2H ≪ ΔB]. Entrainment is also ignored, because we consider a ML after strong mixing has taken place and entrainment is weak (wen = 0). Under these assumptions the boundary conditions reduce to

 
formula

where η is the nondimensional displacement of the lower boundary from its mean depth H(y)/H0. The y dependence in the level at which the lower boundary conditions is applied complicates the analysis. However, displacements of the ML base in this analysis represent heaving by the mesoscale eddy field. Over the scales typical of MLIs (i.e., a few kilometers), mesoscale heaving is of O(1) m or less. Hence, for all but the shallowest MLs, [H(y) − H0]/H0 = O(ɛ) for small ɛ. Expanding the boundary condition (C4) in ɛ, and choosing for consistency ΔB = O−1) and Hy = O(ɛ), the condition reduces to (10).

APPENDIX D

Details of the MITgcm Calculation

The numerical model is the hydrostatic MITgcm (Marshall et al. 1997a, b) with a horizontal resolution of 250 m and a vertical resolution of 5 m. The Coriolis parameter is f = 7.29 × 10−5 s−1. The channel is 300 m deep, and the mixed layer is 200 m deep. The initial velocity is resting, and the initial buoyancy is similar to that in Ou (1984) and Tandon and Garrett (1994):

 
formula

Table D1 shows how the resolution and initial stratification differ between the simulations.

The interior stratification Ni2 is scaled by the Tandon and Garrett (1994) stratification after geostrophic adjustment, By2/f2, so that the stratification ratio between interior and mixed layer is constant. Experiments holding the interior stratification fixed are very similar, so long as the interior stratification is much larger than the mixed layer stratification.

Table D2 uses the results of Tandon and Garrett (1994) and Stone (1970) to produce derived quantities including the maximum mixed layer stratification after geostrophic adjustment (N2) and the corresponding deformation radius (Ld), and Stone’s estimate of the wavelength of the fastest-growing wave (Ls). The deformation radius and the Stone scale are both resolved with better than 3 and 16 grid points, respectively.

The simulations above are typically unforced, with gravitationally unstable density profiles convectively adjusted by means of an increased vertical diffusivity. Some simulations are repeated using a buoyancy forcing at the surface to simulate nighttime convection. In the convectively forced simulations, either convective adjustment or the KPP (Large et al. 1994) are used. The heat flux is given by

 
formula

The nightime cooling strength Q0 is −200 W m−2, the length of daytime heating (th) is 8 h, and the daytime heating parameter (Qd) is set so that no net heat is input at the end of each day. This results in a maximum heating of max(Q) ≈ 718 W m−2. The shortwave heating is allowed to penetrate as Jerlov (1976) water type IA, corresponding to very little turbidity (Jerlov 1976), while the cooling is applied only at the surface.

A number of details of the calculation follow to facilitate reproducibility. A Smagorinsky (1963) harmonic nonlinear horizontal viscosity is used with nondimensional constant 4.0 [implemented on the C-grid as per Griffies and Hallberg (2000)]. Top, side, and bottom boundary conditions are slip. A vertical viscosity of 5 × 10−3 m2 s−1 is used. A third-order, flux-limiting advection of temperature is used, which does not require any explicit diffusion of buoyancy, so none is used (although there is implicit diffusion). All walls are insulating. The time step is 300 s, and an implicit linear free-surface scheme filters faster waves. The initial buoyancy is seeded with a small amount of noise O[2 × 10−8 (m s−2)], which is white in x and y and constant in z.

Footnotes

Corresponding author address: Raffaele Ferrari, Department of Earth, Atmospheric and Planetary Sciences, Massachusetts Institute of Technology, 54-1420, 77 Massachusetts Ave., Cambridge, MA 02139. Email: rferrari@mit.edu

1

SeaSoar is a towed undulating vehicle used to deploy oceanographic monitoring equipment.

2

The calculations of Tandon and Garrett (1994) suggest that during the initial gravitational slumping, before the appearance of MLIs, the front oscillates around a mean state with Ri = (1 − |Byy|max H/2f2)−1. In our simulation H = 200 m and Byy ≈ 7 × 10−12 m−1 s−2, therefore we expect Ri ≈ 1.2. Numerically we find that Ri oscillates between 0.4 and 2.2 around a mean of 1.2.

3

Lateral fluxes become significant only when the slumping front has advanced to reach the meridional walls. We only present results for times before this happens, because lateral walls are a numerical artifact and do not exist in the open ocean.