Here, the effects of surface waves on submesoscale instabilities are studied through analytical and linear analyses as well as nonlinear large-eddy simulations of the wave-averaged Boussinesq equations. The wave averaging yields a surface-intensified current (Stokes drift) that advects momentum, adds to the total Coriolis force, and induces a Stokes shear force. The Stokes–Coriolis force alters the geostrophically balanced flow by reducing the burden on the Eulerian–Coriolis force to prop up the front, thereby potentially inciting an anti-Stokes Eulerian shear, while maintaining the Lagrangian (Eulerian plus Stokes) shear. Since the Lagrangian shear is maintained, the Charney–Stern–Pedlosky criteria for quasigeostrophic (QG) baroclinic instability are unchanged with the appropriate Lagrangian interpretation of the shear and QG potential vorticity. While the Stokes drift does not directly affect vorticity, the anti-Stokes Eulerian shear contributes to the Ertel potential vorticity (PV). When the Stokes shear and geostrophic shear are aligned (antialigned), the PV is more (less) cyclonic. If the Stokes-modified PV is anticyclonic, the flow is unstable to symmetric instabilities (SI). Stokes drift also weakly impacts SI through the Stokes shear force. When the Stokes and Eulerian shears are the same (opposite) sign, the Stokes shear force does positive (negative) work on the flow associated with SI. Stokes drift also allows SI to extract more potential energy from the front, providing an indirect mechanism for Stokes-induced restratification.
Submesoscale fronts in the ocean mixed layer are strong horizontal buoyancy gradients (by ≡ M2) that rival the weak vertical stratification (bz ≡ N2). Geostrophic balance implies a geostrophic Richardson number Rig ≡ (N2f2)/M4 ~O(1) (Tandon and Garrett 1994) that depends on the vertical stratification and geostrophic shear M2/f (where f is the Coriolis frequency). It is the O(1) Richardson number and Rossby number (Ro ≡ U/fL; U is the flow velocity, and L is a characteristic length scale) that distinguish the submesoscale from larger-scale motions.
The balance struck between turbulent mixing and restratification processes, such as geostrophic adjustment, restratifying instabilities, and penetrating solar heating, is very important for organisms that live in the mixed layer. Frontal instabilities that drive restratification are able to compete with destabilizing surface forcing in winter and enhance phytoplankton blooms (e.g., Taylor and Ferrari 2011; Mahadevan et al. 2012).
In addition to Ri ~ O(1), submesoscale fronts also have Ro ~ O(1). A large Ro implies that the ageostrophic velocity is significant in the momentum balance. Capet et al. (2008) show that the ageostrophic component of the velocity enables a forward cascade of energy, toward the dissipation scale, rather than geostrophic turbulence, which has an inverse energy cascade. Molemaker et al. (2010) show that submesoscale frontogenesis and instability leads to a forward energy cascade.
These submesoscale flows are typically restricted to the mixed layer of the ocean because strong forcing from wind and strain by mesoscale features creates fast flows over short length scales [where Ro ~ O(1)]. Convection and wind also make the near-surface stratification very weak (Ri ≲ 1). Since submesoscale flows occur at the upper boundary layer, they coexist with wind and wave forcing.
Despite having a partially geostrophically balanced state, these fronts will have shear and are still unstable. The seminal works of Stone (1966, 1970, 1971) showed that geostrophically balanced fronts with constant vertical and horizontal stratification, and with Rig ~ O(1), have three possibly unstable modes: geostrophic instabilities (GI1), symmetric instabilities (SI), and Kelvin–Helmholtz instabilities (KHI). Criteria for all of the instabilities discussed by Stone (1971) have been established for fronts in the absence of wave forcing. Charney and Stern (1962) and Pedlosky (1964) developed criteria for quasigeostrophic (QG) baroclinic instability that depend on the QG potential vorticity (QGPV) and the shear. Although these criteria were derived for QG baroclinic instability (e.g., Eady 1949), the instabilities (GI) that grow into finite-amplitude mixed layer eddies (MLEs) are very similar in growth rate, wavenumber, and vertical structure (Boccaletti et al. 2007). Stone (1966, 1970, 1971) showed that SI only occur when Rig < 1; however, Rig is based strictly on the geostrophic shear. Hoskins (1974) proved a more general result that SI occur only if the Ertel potential vorticity (PV; the absolute vorticity dotted into the buoyancy gradient) is the opposite sign of f. None of these criteria, however, account for the effects of surface gravity waves on the fronts or the instabilities themselves.
The leading-order influence of surface gravity waves (waves hereinafter) on the flow can be addressed by considering the wave-averaged Boussinesq (WAB) equations (Craik and Leibovich 1976; Huang 1979; McWilliams et al. 1997, 2004). The WAB equations incorporate the Stokes drift (a wave-averaged, surface-intensified current) into the momentum equation through the Stokes–Coriolis (Huang 1979) and Stokes shear forces (the combination of the Stokes vortex force and the Stokes-induced perturbation of the Bernoulli effect; N. Suzuki and B. Fox-Kemper 2015, unpublished manuscript). Stokes drift also aids in the advection of momentum, but the Stokes drift itself is not advected. The Stokes shear force is the primary driver of Langmuir cell instabilities (LC; e.g., Craik 1977); however, here modifications to other baroclinic instabilities by Stokes drift are of the greatest interest.
None of the instabilities mentioned above, or the finite-amplitude features that they grow into, are resolved in global climate models (GCMs). With the finest-scale GCM simulations at ≈0.1° (≈10 km) resolution, even large MLEs O(10) km are not resolved. However, the mean effect of MLEs (finite-amplitude GI) has been parameterized for use in GCMs (Fox-Kemper et al. 2011). Since GI extract potential energy from the front (e.g., Haine and Marshall 1998), the net effect of the MLE parameterization is to increase the vertical stratification in the mixed layer. To date, however, no such parameterization has been made for the mean effect of SI. A few parameterizations for LC have been proposed (McWilliams and Sullivan 2000; Smyth et al. 2002; Van Roekel et al. 2012) and implemented in GCMs (Fan and Griffies 2014; Li et al. 2015), all forced by Stokes drift with the net effect of reducing the vertical stratification in the upper ocean. These parameterizations are independent from each other in that the fronts are not accounted for in the LC parameterization, and Stokes drift is not accounted for in the MLE parameterization. This independence implies that the net effects are additive, and any changes to the MLEs due to Stokes drift or changes to LC due to fronts are unaccounted for.
McWilliams and Fox-Kemper (2013) showed that when an initially geostrophically balanced front is encountered by waves (represented as Stokes drift), the Stokes–Coriolis and Stokes shear forces disrupt the geostrophic balance, forcing the front to adjust. For a front with Ro = 0, the adjustment results in an Eulerian, anti-Stokes flow that directly opposes the Stokes drift to balance the Stokes–Coriolis force. For a front with finite Ro (which will not be considered by the linear solutions to follow), the Stokes shear force also alters the hydrostatic balance of the front.
Ursell and Deacon (1950) showed that a canceling Eulerian mean flow arises (anti-Stokes flow) in the presence of the Stokes–Coriolis force. Broström et al. (2014) confirms this anti-Stokes flow and shows that the Stokes–Coriolis force does not inject energy into the mixed layer when considering the vertically integrated momentum and energy equations. McWilliams et al. (2012, 2014) confirmed in nonlinear large-eddy simulations (LES) with Stokes drift and wind forcing that anti-Stokes Eulerian flow is maintained by the leading-order geostrophic and Ekman balance. In simulations of submesoscale fronts with Stokes drift and wind forcing, Hamlington et al. (2014) showed evidence of the anti-Stokes flow from altered PV compared to a simulation without Stokes drift. Anti-Stokes flow has been observed and shown to be important for flows across the inner continental shelf (Lentz et al. 2008; Lentz and Fewings 2012). Fronts in the presence of Stokes drift with constant shear are unstable to a hybrid SI/LC mode (Li et al. 2012). Notably, unlike SI in the absence of Stokes drift, this mode extracts potential energy from the front, thereby restratifying rather than mixing it.
Since the stability criteria for SI and GI depend on the strength of the shear, these criteria are expected to change as the Stokes drift itself vertically varies and also induces the anti-Stokes Eulerian shear. The implications of these effects are of primary interest here, and care is required as Stokes drift does not directly contribute to the vorticity in the wave-averaged equations. Nakamura (1988) showed that although GI and QG baroclinic instability are similar, GI has a reduced wavenumber and growth rate compared to QG baroclinic instability because small Ri fronts have stronger shear.
Surface wind stress induces shear and, in the presence of a front, a frictional PV flux out of (into) the ocean when the winds blow downfront (upfront) (Thomas 2005; downfront implies the wind and geostrophic shear are aligned), which can then induce (suppress) SI. Thomas and Taylor (2010) refer to this driving of SI by wind as forced symmetric instability (FSI) and show that it scales with the Ekman buoyancy flux (EBF; the dot product of the Ekman transport with the horizontal buoyancy gradient). Regions of higher than expected turbulent kinetic energy were observed in the Kuroshio, where the EBF (but not necessarily the wind stress) was large (D’Asaro et al. 2011), suggesting that FSIs were present. Here, a related effect will be studied under downfront (upfront) Stokes drift conditions.
Although the introduction of Stokes drift within a front allows for a new class of instability (LC), the focus of this work is on modifications to baroclinic instabilities (GI and SI) by Stokes drift. With the increased shear due to Stokes drift, how do the criteria for GI and SI change? Which shear—Eulerian, Stokes, or Lagrangian—is relevant for these criteria? How does the Stokes drift influence the PV and thereby trigger or suppress SI? Is the Stokes–Coriolis or Stokes shear force more important for changes to baroclinic instabilities? These questions will be answered in the next several sections.
First, the problem of interacting waves and fronts is described in detail (section 2). Section 3 provides analytically derived stability criteria for SI and GI in idealized settings. The analytic stability criteria are confirmed in more complicated settings with numerical linear stability analysis in section 4. The linear analyses are complemented with nonlinear LES in section 5. The conclusions are discussed in section 6.
2. Problem description
where u is the three-dimensional, Eulerian velocity, uS is the Stokes drift, uL ≡ u + uS is the Lagrangian velocity, f is the Coriolis parameter, p is the pressure divided by a reference density, b is the buoyancy, ν is the kinematic viscosity, and κ is the thermal diffusivity. The term uL,j refers to the jth component of the Lagrangian velocity, and, as usual with Einstein notation, the repeated j index implies a summation over the three spatial components. The Stokes shear force uL,j∇uS,j form of the wave-averaged equations has been chosen here rather than the Stokes vortex force form. These forms are mathematically equivalent. This choice is made because when the Stokes drift is horizontally invariant (which will be assumed) the Stokes shear force only appears in the vertical momentum equation as . Then the effects of Stokes drift are the Stokes–Coriolis force and a vertical force that varies horizontally with the Eulerian mean flow. It is this horizontally varying vertical force that induces the upwelling and downwelling branches of LC. N. Suzuki and B. Fox-Kemper (2015, unpublished manuscript) examine the Stokes shear force in detail and show that this form is equivalent to other forms of the wave-averaged equations.
In the open ocean, Stokes drift has a relatively persistent character, which justifies the assumption of time and horizontally invariant wave forcing. Wave buoy data at Ocean Station Papa (50°N, 145°W; sampled every 30 min for 2 yr; Thomson et al. 2013) reveals an autocorrelation of the Stokes drift on relatively long time scales compared to the growth rates of instabilities considered here, which are O(1) day or less. The surface Stokes drift autocorrelation zero crossing time varies from 9 to 15 days for each vector component, while the integral time scale ranges from 287 to 290 days.2 Similar estimates result from the WAVEWATCH model (Tolman 2009; sampled every 3 h for 1 yr; over the entire ocean, 95% of the zero crossing times are 1.75 days or greater, with an integral time scale of 40 days or greater) are comparable to those of the wind stresses used to drive the wave model (Large and Yeager 2009).
The Stokes drift throughout this work is also assumed to be monochromatic, resulting in a simple exponential profile:
where HS ≡ (2kw)−1 is the e-folding depth of the Stokes drift relating to the wavenumber kw of the waves. A linear Stokes drift profile is considered in section 3c for analytic convenience.
The vorticity equation can be used in place of the momentum equation
where is the (Eulerian) vorticity. Rescaling the momentum equation as in McWilliams (1985), and the Stokes shear term as in McWilliams and Fox-Kemper (2013), gives the dimensionless momentum equation:
where Ro ≡ UL/(fH), Eu ≡ max(1, Ro−1) is the Euler number [p/(ρU2)], , λ ≡ H/HS, and H is the mixed layer depth (see Table 2 below). Note that when Ro < 1, then μλRo is the parameter that compares the Stokes shear force to the buoyancy. This is the small parameter exploited asymptotically by McWilliams and Fox-Kemper (2013).
Scaling the equations
To linearize the equations, several assumptions are made to achieve a steady mean state around which to perturb. In dimensional form, we have the following:
multiple scales of horizontal variation [i.e., , l ≪ L, and similarly for ]. There is only one scale for vertical variation and it is different from the horizontal scales .
The assumption of multiple horizontal scales allows for fast and slow time scales [i.e., ]
Each variable is decomposed into a mean and perturbation (e.g., u = U + u′). The averaging operator is (an average over the small horizontal scales and the fast time scales). For example, , where , and .
The scales above are all chosen based on realistic quantities for the ocean mixed layer. Scales for vertical velocity, buoyancy, and pressure are derived therefrom (see Table 1). The nondimensional numbers, which determine the regime and stability of the problem, are given in Table 2. The complete, rescaled, dimensionless equations are given in the appendix, without approximation.
The complete equation set collapses to several previous linear stability problems that are relevant for ocean dynamics including QG baroclinic instability (Eady 1949), ageostrophic, nonhydrostatic baroclinic instability (Stone 1971), LC instabilities (Leibovich and Paolucci 1981), LC instabilities in an Ekman layer (Gnanadesikan and Weller 1995), and downwave invariant mixed LC/SI instabilities (Li et al. 2012). The complete validation of the present model with respect to the studies above is described in detail in Haney (2015). The combination of the nondimensional parameters required to collapse to each of these regimes is summarized in Table 3.
Although Eqs. (A1)–(A4) can reproduce all the above instabilities as well as those in the parameter space between these disparate regimes, the focus of this work will be the modification to balanced and unbalanced baroclinic instabilities (GI and SI) by the Stokes–Coriolis and Stokes shear forces. As such, inviscid flow will be assumed in the analytic and linear calculations to follow, and a discussion of instabilities in the presence of viscosity will be deferred to subsequent papers.
3. Analytic stability criteria
In a few special cases of fronts in the Lagrangian thermal wind balance (Ek = 0), analytic criteria for GI, SI, and KHI can be derived. A review of what is meant by Lagrangian thermal wind balance is therefore helpful.
McWilliams and Fox-Kemper (2013) showed that a front or filament undergoes an adjustment to accommodate the Stokes–Coriolis and Stokes shear forces. The Stokes shear force alters the vertical momentum equation [Eq. (6)] and scales as μλ (relative to the advection term), whereas the vertical pressure gradient and buoyancy forces scale as Ro−1 when Ro ≪ 1, and therefore Eu ~ Ro−1 [see Eq. (A2)]. Therefore, when Ro−1 ≫ μλ [i.e., μλRo ≡ ε ≪ 1 in McWilliams and Fox-Kemper (2013)], the Stokes shear force can be neglected. Then the only adjustment required in the horizontal momentum equation is to balance the Lagrangian–Coriolis force with the horizontal pressure gradient. Anti-Stokes Eulerian flow (McWilliams and Fox-Kemper 2013) arises because in the absence of any pressure gradient (in small Ro, inviscid flow), the Stokes and Eulerian–Coriolis forces must exactly balance:
With a front, and with the assumption that the Stokes shear force is weak (μλ ≪ Ro−1), then the vertical momentum balance is hydrostatic; one can write the thermal wind balance as
where ∇H = (∂X, ∂Y, 0). The last term on the right-hand side is the anti-Stokes Eulerian shear. In fronts with Ro ~ 1 (μλ ~ Ro−1; which will not be considered here), there is an additional perturbation to the buoyancy and velocity beyond anti-Stokes flow (McWilliams and Fox-Kemper 2013). Note that the Ro ≪ 1 assumption allows for linearized equations, but, as is common in asymptotic theories, the hope is that the linear theory is robust and applies even when this strict assumption is not met. Notably, the conclusions of the linear theory are tested and many are confirmed, with nonlinear simulations in section 5 in which the simulated fronts have Ro ~ 1.
Therefore, in a balanced front with Stokes drift, the Lagrangian shear can be determined by the buoyancy gradient; however, the Eulerian shear will be dramatically different due to the anti-Stokes Eulerian flow. It will be shown in the next few sections that GI, SI, and KHI have stability criteria that depend on the Lagrangian flow, Eulerian flow, or both. Note that the frontal adjustment to the presence of waves cannot alter the PV in the flow (except possibly by advection); therefore, the anti-Stokes contribution to the PV only maintains, fluid parcel by fluid parcel, the structure of PV that was preexisting in the front before the arrival of waves. Of course, diabatic and wind effects may also affect the PV (Thomas 2005), and the winds will lead to waves. For brevity, this paper neglects direct study of viscid, diabatic dynamics. Even so, the connection between winds and waves is not local in space or time, and the wave group velocity is significantly different from the Stokes drift. Thus, waves may arrive and deliver Stokes forces far from any winds or PV anomalies associated with their creation. Last, these changes to the Eulerian flow only represent the changes in the mean state, while the Stokes drift also directly affects instabilities through the Stokes shear force (see sections 4–5).
a. Kelvin–Helmholtz instability
Holm (1996) showed that a sufficient criteria for KHI (Ri < ¼) holds true for the Lagrangian mean Richardson number and that a necessary criteria for KHI is that there must be an inflection point in the Eulerian velocity. Vanneste (1993) showed that since the Eady background state used in Stone (1971) has no inflection point, it is stable to KHI. The mean flows studied here do not have inflection points, and therefore these criteria will not be tested with linear stability. Nevertheless, it is important to keep in mind that in any realistic mixed layer, with a frictional boundary, an inflection point in the Eulerian flow is likely to occur, and therefore flows with a Lagrangian Ri < ¼ will be unstable to KHI.
b. Geostrophic instability
Beginning from the WAB Eqs. (1)–(3), and following a standard derivation of the QG equations [e.g., Pedlosky 1982; or equivalently, assuming the QG limit in Table 3 for Eqs. (A1)–(A4)], the Stokes-modified QG equations are
where , ∇h = (∂x, ∂y, 0), q is the QGPV, β is the change in Coriolis frequency with latitude, and ΨL is a Lagrangian streamfunction that satisfies
Note that Eq. (12) involves a mixture of the Eulerian and Lagrangian streamfunctions. The Eulerian streamfunction appears in the relative vorticity. The Stokes drift only influences the Lagrangian streamfunction that appears in the vortex stretching of the buoyancy surfaces in a (Lagrangian) thermal wind balance. Although , VS = −V may be nonzero, so long as they are horizontally invariant. The perturbation QGPV evolution is given by
where and . The QG perturbations may vary in any direction, and normal-mode form is assumed for the alongfront (x) direction:
The QGPV and buoyancy conservation (on the boundaries) then become
The energy equation is then formed by dotting into Eq. (19). Now, considering the imaginary part of the volume-integrated energy equation,
where ci is the imaginary part of the wave speed. Recall from Eq. (18) that for GI to occur ci must be nonzero. If ci ≠ 0, then at least one of the following must be true [with QL as given by Eq. (12)]:
changes sign in the interior of the domain;
is the opposite sign as at z = 0;
is the same sign as at z = −H;
has the same sign at z = −H and z = 0.
Therefore, the criteria for GI are unchanged so long as the QGPV and shear are appropriately interpreted in their Lagrangian forms as in Eqs.(12), (14), and (15). This implies that once the front reaches Lagrangian thermal wind balance, the stability of the flow to GI is unchanged, and therefore the Stokes drift appears to do very little to alter GI. However, these are only stability criteria that are founded on the Stokes-altered mean state, and they do not give information about what happens to GI under the influence of the perturbation Stokes shear force. This will be explored in section 4.
c. Symmetric instability
An analytic criterion for SI can be derived by following the method of Hoskins (1974), which reveals that the Ertel PV must take the opposite sign of f for SI to occur. Assuming the mixed LC/SI limit in Table 3 for Eqs. (A1)–(A4) and downfront invariance and taking the curl of the velocity yields
Here, we will assume that the buoyancy gradients in the front are constant in the horizontal (M2) and the vertical (N2). This therefore implies that our mean flow, which is in thermal wind balance, has constant geostrophic shear M2/f. Eliminating ω′, u′, and b′, and assuming a normal-mode solution of the form , yields
where τ ≡ tanϕ. Up to this point, any Stokes drift profile may be considered. To obtain stability criteria, the equation is simplified by assuming a Stokes drift profile with constant shear and thus neglecting the last term on the right-hand side. For the exponential Stokes drift profiles considered in the following sections, neglecting this term is unjustified; nevertheless, the assumption is made in order to obtain an analytic solution. For unstable modes, σ2 < 0, and since σ2 is given by a quadratic equation for τ, it is sufficient to show that there exist two distinct real roots of Eq. (26) such that −1 < τ < 1. To obtain real roots, the following [a positive discriminant of Eq. (26)] must be true:
Thus, Hoskins’ original negative (or more precisely, anticyclonic) PV criteria is obtained in Eq. (28), where the anti-Stokes flow contributes to the PV (Q). Note that the PV is multiplied by f to accommodate for either hemisphere:
Since the anti-Stokes Eulerian flow has decoupled the total Eulerian shear from the buoyancy gradient (which, in the absence of ageostrophic Eulerian shear, is related through thermal wind Uz = −BY/f), a criterion on Richardson number for the onset of SI is inadequate.3 A downfront Stokes shear (cosθ > 0) is stabilizing, while an upfront Stokes shear is destabilizing to SI. Recall that downfront winds extract PV from the mixed layer (Thomas 2005), reducing the stability to SI. Therefore, in a front in Lagrangian thermal wind balance, when the winds and waves are aligned, as is often the case (Webb and Fox-Kemper 2015), Stokes shear results in PV of the opposite sign to the PV that is being injected by winds. Although the Stokes drift does not create PV like the wind can, the presence of Stokes drift affects PV through effects on the near-surface velocity and mixing. Such a situation arises profoundly in Hamlington et al. (2014). In that case, two simulations were compared; one with wind only and one with aligned winds and waves (Stokes drift) pointed partially downfront. Both simulations included a front and the same surface buoyancy flux (cooling). In the case with Stokes drift, the PV in the front is significantly more positive, as expected under Lagrangian thermal wind balance. Furthermore, the mixed layer as measured by a sharp change in PV was much deeper than the temperature mixed layer in the wind-only case, suggesting that SI were stronger in the no Stokes case but were still less effective than Langmuir turbulence at mixing away negative PV near the surface.
d. Parcel theory for symmetric instability
To supplement the analytic proof of the SI criterion in section 3c, one can consider the change in potential and kinetic energy when two parcels are switched along a particular path. Haine and Marshall (1998, their appendix B) used parcel switching arguments to show that for an inviscid (no wind) background flow in thermal wind balance (no Stokes), the criterion for SI to exist (fQ < 0) holds true. Following their methods identically, the geostrophic shear can be replaced throughout with an arbitrary Eulerian shear that is a sum of geostrophic and ageostrophic Eulerian components. The Eulerian shear is chosen because the curl of the Stokes drift is absent from the wave-averaged vorticity [Eq. (5)]. Physically, the Stokes drift cannot introduce vorticity into the flow because it is a combination of irrotational wave quantities. Vorticity is only introduced by the imposed mean flow, which then interacts with the waves such that the wave-averaged effect is an interaction between the Stokes drift and the imposed mean vorticity. A consequence of the absence of Stokes vorticity is that the momentum concentration is Eulerian, as in Eq. (1), and that apart from advection, the Stokes drift only acts on the vertical momentum (through the Stokes shear force). Many mathematically equivalent forms of the WAB equations and the momentum concentration are presented in N. Suzuki and B. Fox-Kemper (2015, unpublished manuscript) and Holm (1996). The form chosen here is consistent with absolute Eulerian momentum for a fluid parcel (m = uE − fy). However, because of the anti-Stokes Eulerian components of the background flow, the available kinetic energy and the shapes of constant absolute Eulerian momentum surfaces are changed by the presence of Stokes drift.
A schematic of parcel switching in a front in Lagrangian thermal wind balance with downfront and upfront Stokes drift is show in Figs. 1a and 1b, respectively. Below the influence of the Stokes drift (below ≈0.4), switching parcels along lines of constant absolute Eulerian momentum or constant buoyancy is unstable since the momentum perturbation is amplified. When switching along surfaces of constant absolute Eulerian momentum, a positive vertical velocity perturbation yields a positive buoyancy force, thereby amplifying the positive vertical velocity. In the case with downfront Stokes drift, the momentum surfaces are altered near the surface so that switching along constant buoyancy surfaces is stabilizing. Similarly, switching parcels along constant momentum surfaces in this case results in a stabilizing buoyancy force. The upfront Stokes drift case remains unstable, but none of the cases in the schematic show how unstable these scenarios are. Nevertheless, the instability clearly depends on gradients of absolute Eulerian momentum and buoyancy.
The derivation of a change in kinetic (ΔK) and potential energy (ΔP) of the mean flow follows Haine and Marshall (1998) to find
where s is the slope of the path along which the parcels exchange, and −M2/N2 is the isopycnal slope. The individual sources of energy are the kinetic energy required for parcel switching (Kp; because of the cross-front motion of the parcels), the available kinetic energy from the geostrophic shear (GS) and the ageostrophic Eulerian shear (AS), and the potential energy that needs to be overcome. This can be thought of as mixing (MIX), or restratification, and potential energy loss from the front, if this term is negative. Minimizing Eq. (29) with respect to s gives the energy available for the parcel exchange:
When this quantity is negative, there is available energy for parcel exchange. Relative to the geostrophic shear-only case, the ageostrophic Eulerian shear changes GS by changing the parcel path [Eq. (30)], it contributes available shear (AS), and it forces some mixing or restratification since the parcel path is no longer along isopycnals. Note that, excluding the last term (ES + MIX) and assuming that the only ageostrophic shear is the anti-Stokes Eulerian shear , the expression above is identically fQ [as in Eq. (28)]. Since the last term always contributes to a loss of energy from the mean flow to turbulent motions, the anticyclonic PV criterion derived in the previous section is always stricter. This dissonance will be explored numerically in the next section. Last, this derivation applies in cases forced by Ekman shear (i.e., FSI; Thomas and Taylor 2010; Thomas et al. 2013).
4. Numerical linear stability analysis
The analytic stability criteria (section 3) were derived in the special cases of constant downfront Stokes shear (for SI) and QG perturbations (for GI). However, instabilities in the mixed layer are not constrained this way. Furthermore, while stability criteria predict the onset of SI and GI, they do not address any changes in the structure of these modes due to the presence of Stokes drift, nor do they address the effect of the Stokes shear force on these instabilities. Both of these issues are explored in this section by numerically solving a linear subset of Eqs. (A1)–(A4).
a. Assumptions and setup
To understand the effects of Stokes drift on SI and GI, the viscid effects will be neglected in order to highlight how the familiar, classic versions of these instabilities are altered at leading order by Stokes drift. Beginning from Eqs. (A1) to (A4), assume
Note that here. Assuming Roδ, Roα2 ≪ 1 linearizes and simplifies the equations. Averaging over the small (x, y, t) scales and retaining O(1) and larger terms:
One steady-state solution to the above equations is the Eady-like background state:
Lagrangian thermal wind balance holds, and although VL = 0, V and VS may be nonzero (if the cross-front Eulerian flow is identically an anti-Stokes flow). Since Roμλ ≪ 1, the Stokes shear force μλRoUL ⋅ US can be neglected. However, it is a leading-order term in the perturbation equations (38)–(42). This mean state is depicted in Fig. 2.
The perturbations are given by
The perturbations rely on only five dimensionless parameters: the Rossby number of the perturbations Ro′ ≡ Ro/δ, the aspect ratio of the perturbations α′ ≡ α/δ, the strength μλ and depth scale λ of the Stokes shear μλeλz, and the Lagrangian Richardson number of the mean flow Ri. Note that , putting the perturbed motions in the submesoscale regime (Ro ~ 1, Ri ~ 1). Applying the typical scales in Table 1 yields Ro ≈ 0.1 and Ro′ ≈ 10.
The perturbations are then assumed to have normal-mode form in the x and y directions such that (and similarly for υ′, w′, b′, p′). The problem is posed as a generalized eigenvalue problem with the growth rates σ as eigenvalues and the (complex) vertical structures , , , , and as eigenvectors. The background and perturbation variables are discretized into 50 Chebyshev spectral modes and solved on a Gauss–Lobatto grid using the tau method (see Boyd 2001). Further details and the MATLAB code are given in Haney (2015).
The linear setup is chosen to model only the mixed layer with a rigid lower boundary. No viscosity is used, so no stress conditions are needed. Boccaletti et al. (2007) show that introducing a moving mixed layer base only weakly affects the linear instabilities within the mixed layer, which is assumed to be true here. The nonlinear simulations in section 5 below do include a pycnocline below the mixed layer.
Other useful tools in distinguishing unstable modes are the energy sources of each mode, given by
where is the Lagrangian total derivative; is the kinetic energy of the perturbed flow; PW is the pressure work; ESP is the Eulerian shear production; SSP is the Stokes shear production; BP is buoyancy production; and D is energy dissipation. The energy sources for each mode will help identify the mode. For example, it is well established that, in the absence of Stokes drift, GI extract potential energy from the front (positive BP) thereby reducing the horizontal buoyancy gradient and restratifying, whereas SI extract kinetic energy from the geostrophic shear (positive ESP), thereby reducing it (Haine and Marshall 1998). The energy sources for each mode will also serve as a point of comparison between the linear stability analysis and the nonlinear large-eddy simulations.
c. Stokes-modified baroclinic instabilities
As alluded to in section 3b, Stokes drift has only a weak effect on GI. Since the mean Lagrangian shear in the cases studied here has the same sign at the top and bottom (see section 3b), all cases are unstable to GI. The growth rate is only modestly changed, and a slight shift to higher wavenumbers is evident (Fig. 3). Since GI is driven by the potential energy throughout the mixed layer (e.g., Eady 1949; Haine and Marshall 1998; Boccaletti et al. 2007), it is unsurprising that changes in shear over a small portion of the domain yield only subtle changes in stability. Nakamura (1988) showed that since GI is a result of interacting edge waves for a given wavelength, strong shear reduces the penetration depth of the edge waves. Therefore, longer edge waves are required for them to interact and produce GI. This appears to be the case for Stokes-modified GI, and the relevant shear is the Eulerian shear. With downfront Stokes drift, the anti-Stokes Eulerian shear reduces the total Eulerian shear near the surface, which according to Nakamura (1988) would increase the penetration depth of the edge waves, resulting in smaller GI with larger growth rates. This is consistent with the fact that GI are smaller scale and grow faster with downfront Stokes drift and are larger scale and grow slower with upfront Stokes drift (Fig. 3). However, this effect is weak. The work done on the GI by the Stokes shear force is generally small because GI are driven primarily by the potential energy (BP) rather than the kinetic energy in the shear.
Numerical solutions with exponential Stokes drift confirm that the PV criterion for SI is valid for Stokes drift profiles that decay exponentially with depth as well as the constant shear profiles (not shown) considered in section 3c. In the cases where Stokes drift decays exponentially with depth, the upper part of the domain has PV fixed by the anti-Stokes Eulerian flow, and the lower part has PV that is roughly determined by the geostrophic shear (Fig. 4). These cases also show that neither nor would consistently predict SI. Furthermore, Fig. 4 shows that SI exists in, and only in, the parts of the domain with anticyclonic PV, that is, opposite sign to f.
Furthermore, the anticyclonic PV criterion extends beyond the strictly downfront or upfront Stokes drift cases presented in section 3c. Numerical solutions with cross-front Stokes drift show that the anticyclonic PV criterion holds even with Stokes drift incident at oblique angles. Figure 4 shows cases with Stokes drift rotated π/4 radians to the left of parallel and to the right of antiparallel, with the geostrophic flow to give partial cross-front Stokes drift. Therefore, the criterion for SI with arbitrary Stokes front alignment is as in Eq. (28) with no modification by VS (the cross-front Stokes drift component).
Figure 4 also shows that the alignment of the Stokes-altered SI is across isopycnals as suggested by the altered shape of the momentum surfaces in the parcel switching analysis. The momentum surfaces are less (more) steep where the anti-Stokes Eulerian shear is stronger (weaker or negative) than the geostrophic shear. Therefore, although the instability criteria based on the available energy is necessary but not sufficient, when SI do exist, they align along the paths that tend to maximize their energy.
Haine and Marshall (1998) note that the stability of such a layer depends on the buoyancy gradient within a momentum surface, or the momentum gradient within a buoyancy surface, each giving the anticyclonic PV criterion for instability. Although these are two ways to measure the same mean flow quantity (PV), the alignment of the ensuing instabilities reveals whether they derive energy from the mean buoyancy gradient (potential energy) or from the mean momentum gradient (kinetic energy). In the classical case (waveless case; not shown), SI align mostly with buoyancy surfaces and therefore get their energy from the mean shear. The Stokes-modified SI in Fig. 4 are aligned partially across isopycnals, implying that they get a larger fraction (than classical SI) of their energy from buoyancy production (confirmed in Fig. 5). This is because the stronger vertical shear near the surface adds enough momentum to reduce the penalty for horizontal parcel movement (because of the turning of momentum by the Coriolis force). Therefore, when the momentum surfaces are closer to horizontal, more potential energy may be extracted from the front.
Figure 5 shows vertical profiles of the energy production terms that dominate the SI in Fig. 4. In both cases, BP and ESP are important. Therefore, the Stokes-modified SI get some of their energy from the kinetic energy in the front and some from the potential energy. Therefore, this is a mechanism by which Stokes drift may indirectly drive restratification. This mechanism for restratification is distinct from that discussed in Li et al. (2012) since theirs depends on the nonlinear interaction between multiple tilted LC. Here, the anti-Stokes Eulerian shear maintains negative PV to switch on SI and changes the available kinetic energy allowing SI to do more BP.
The altered PV and available energy due to the anti-Stokes Eulerian flow demonstrate the effects of the mean flow on SI; however, the perturbation Stokes shear force also influences SI. While the Stokes-altered mean flow affects SI through the vertical advection of the mean horizontal momentum (through the term), the Stokes shear force alters the vertical momentum of the perturbed flow. When SI are present, ESP (which comes exclusively from the horizontal momentum equation) is positive . SSP (; which comes exclusively from the vertical momentum equation) is identically the work done by the Stokes shear force. When SI dominates, SSP is only positive when the Stokes and Eulerian shears have the same sign. Therefore, in the case with upfront (downfront) Stokes drift in Fig. 4, the Stokes shear force does negative (positive) work on SI.
The near-surface, Stokes-modified PV layer demonstrates the role of Stokes drift in the stability of balanced and unbalanced, baroclinic motions. Furthermore, this suggests that computing the PV solely based on the geostrophic flow, or the presence of SI based solely or partially on Ri (e.g., Li et al. 2012; Thomas et al. 2013; Hamlington et al. 2014), will result in significant errors. If the Stokes drift and geostrophic flow are downfront (upfront), the Stokes-altered PV suppresses (enhances) SI near the surface. Thus, downfront Stokes drift has the opposite effect on the PV as downfront winds. Additionally, if SI are induced by the Stokes drift, the shape of constant momentum surfaces favors increased BP by SI. This provides a mechanism by which Stokes drift may restratify (rather than mix through LC) the mixed layer.
5. Large-eddy simulations
The previous sections necessarily omit nonlinear effects by assuming small-amplitude perturbations, but what happens when these perturbations are allowed to grow and interact with each other? To complement the analytic and numerical linear stability results in sections 3 and 4, a suite of LES have been run. A limited analysis of these runs is provided here to illustrate the robustness of the preceding theory and linear stability results. These nonlinear simulations confirm many results of the previous sections: Stokes drift sets the PV near the surface through the anti-Stokes Eulerian flow, SI are present in and only in regions where the PV is anticyclonic, and, in regions where SI are dominant, the Stokes shear force does work against the SI.
Three simulations are performed with different initial mean states:
a control case with a negative PV front without Stokes drift;
a case with negative PV below the influence of the Stokes drift; and
a case with positive PV below the influence of the Stokes drift.
The latter two cases were chosen to be similar to the mean states in Fig. 4. The initial PV and buoyancy of each case are shown in Fig. 6, and the LES parameters for each case are given in Table 4. The linear analyses of sections 3 and 4 can be used to predict the presence of SI and the work done by the perturbation Stokes shear force on SI in the LES. These predictions are given in Table 5.
a. Model configuration
The National Center for Atmospheric Research (NCAR) LES model (Moeng 1984; McWilliams et al. 1997; Sullivan and Patton 2011) is used to solve the WAB Eqs. (1)–(3). Horizontal derivatives are computed pseudospectrally, and a second-order finite difference scheme is used for vertical derivatives. This model is described in great detail in McWilliams et al. (1997) and has been used to simulate a realistic mixed layer with fronts and Stokes drift (Hamlington et al. 2014). The side boundaries of the domain are periodic. The top boundary condition is no flux and no stress (i.e., no wind). No wind is applied to attempt to isolate the effects of Stokes drift rather than wind-driven flow. The LES is set up to closely emulate the background state of the linear stability (as in Fig. 4, but without cross-front Stokes drift; VS = 0); however, the periodic boundary condition requires that the two fronts be simulated side by side. This allows for two different alignments of Stokes drift (parallel and antiparallel) with the geostrophic shear in each simulation. A monochromatic Stokes drift is chosen to emulate remotely generated swell with a period of approximately 9 s and an e-folding depth of 10 m. The domain is intentionally restricted in the alongfront direction to suppress GI. Since the influence of Stokes drift on SI is the most robust result of the analytic and linear analyses, this will be the focus of the LES in this work. The LES parameters for each simulation are given in Table 4.
The buoyancy profile is specified as two piecewise continuous fronts with smoothing at the transitions. The fronts are given constant stratification in the horizontal and vertical in order to tightly control Ri (outside of the transition regions). The initial geostrophic velocity is then given by the thermal wind relation . The anti-Stokes Eulerian flow is added to the geostrophic flow to determine the Eulerian flow that the model evolves. The horizontal velocity (both u and v) is then given a small-amplitude, random (white noise) perturbation that is independent of depth. In addition to the two front buoyancy structure, a strongly stratified ( = 10−4 s−2) pycnocline is used at ≈50 m deep to create a more realistic mixed layer than in the linear stability analysis (which had a simple rigid bottom boundary).
In the two cases with Stokes drift, the near-surface PV is negative where the Stokes drift is upfront and positive where it is downfront (Figs. 6b,c). In the case with higher vertical stratification (Ri = 2), the PV is positive at depth, and in the case with lower vertical stratification (Ri = 0.5), the PV is negative at depth (as in the control case; Fig. 6). The two separate frontal regions denoted F1 and F2, which will be referred to throughout the rest of this section, are indicated in Figs. 6, 7, 9, and 12. In both of the cases with Stokes drift, although no wind stress is applied, an Ekman layer develops because the anti-Stokes shear is nonzero at the surface. This sets up a Stokes–Ekman flow (see Gnanadesikan and Weller 1995; Polton et al. 2005; McWilliams et al. 2014) that will be discussed further in a subsequent paper on instabilities that depend on viscosity.
b. Symmetric instabilities
SI are most easily seen in the perturbation cross-front velocity υ′ that alternates in direction. The waveless control case (Ri = 0.5, μ = 0) initially has negative PV throughout both fronts and exhibits SI as expected (Figs. 6a, 7, and 8), which reach finite amplitude in about a day. In this case, since the only mean Eulerian shear is geostrophic, Ri < 1 is equivalent to PV < 0. The PV in this case is initially negative (≈−3 × 10−11 s−3 between the surface and 47 m deep) and is increased by an order of magnitude (to ≈−2 × 10−12 s−3) within a few days of the onset of SI (Fig. 8b). The restoration of PV toward zero is coincident with an increase in volume-integrated vertical stratification (Fig. 8a), implying that SI get some of their energy from BP. The ESP (Figs. 8c,d) shows a strong increase at approximately the same time that the SI appear in the cross-front velocity, and when N2 begins to increase and the PV begins to increase. These features in cross-front velocity, PV, and ESP will be used as a measure of whether SI are present in the cases with Stokes drift.
As predicted in Table 5, SI are present in and only in regions with negative PV. The Ri = 0.5, μ = 2 case (Figs. 6b, 9, 10, 11) exhibits alternating cross-front velocity aligned with isopycnals at depth in both F1 and F2 (Fig. 9). Near the surface, these features only appear in F2 where the Stokes drift contributes to reduce the PV. There are also alternating cross-front velocity features near the surface in F1; however, the vertical velocity at 5 m deep (not shown) indicates that these are a mixed gravitational and LC instabilities. SI affect the mean state in F2 by increasing the stratification and mean PV between 20 and 47 m deep from ≈−6 × 10−11 s−3 initially to ≈−8 × 10−12 s−3 (Figs. 11a,b). In F1, where the SI are weaker, the mean PV below the Stokes layer increases from ≈−4 × 10−12 s−3 initially to ≈2 × 10−11 s−3 (Figs. 10a,b). In F2, at every depth where SI are dominant (below ~10 m), SSP (Figs. 11e,f) is negative as predicted in Table 5 because the Stokes shear opposes the Eulerian shear in this front. ESP (Figs. 11c,d) is positive in this region, consistent with the energetics of SI found in the control case and the linear stability analysis. Within the Stokes layer, no SI are present in F1 because the PV is positive (Figs. 6, 9). The positive SSP and ESP (Figs. 10c–f) near the surface in F1 are due to LC, which display strong vertical velocity perturbations at 5 m (not shown).
SI are absent everywhere below the Stokes drift layer in the Ri = 2, μ = 1 case (Figs. 6c, 12, 13, and 14; as predicted in Table 5) since the PV is positive at depth in both fronts (Fig. 6c). In the downfront Stokes layer (F1), the strongly positive PV is maintained by the anti-Stokes flow, and as expected, no SI are present. The shear production is dominated by positive SSP (Figs. 13e,f), indicating that LC are the dominant instability in this front.
In the upfront Stokes case (F2 of the Ri = 2, μ = 1 case), negative PV is maintained by the anti-Stokes flow, and weak SI are present (Fig. 12). The energetics of F2 (Fig. 14) are not as starkly in agreement with the predictions in Table 5 as in the Ri = 0.5, μ = 2 case (Fig. 11). The ESP (Figs. 14c,d) oscillates but shows no strong increases coincident with SI features in the cross-front velocity. However, the SSP (Figs. 14e,f) is notably negative after about 6 days, indicating that SI, not LC, are dominant. Last, the mean PV between the surface and 20 m increases from ≈−5 × 10−11 s−3 initially to ≈−2 × 10−11 s−3 after 13 days (Fig. 14b) due to the weak SI.
Despite the fact that there are two fronts adjacent to each other, it has been assumed in the analysis that the effects of the fronts on each other are weak. Some evidence for this can be seen in the weak cross-front velocities at depth in Figs. 7, 9, and 12. The strong cross-front velocities near the surface in these cases are due to LC. Although the LES featured here showcase more physics than were discussed, such as LC and Stokes–Ekman layers (see McWilliams et al. 2014), a discussion of these physics will appear in a subsequent paper. Haney (2015) presents a preliminary discussion of these physics.
The presence of Stokes drift alters the stability of the mixed layer through the Stokes–Coriolis force. After a few inertial periods, the Stokes–Coriolis force will be in balance with the Eulerian–Coriolis force, the pressure gradient force, and the shear stress. To leading order in Roμλ, this alters the strength of the Eulerian–Coriolis force and the Eulerian shear while maintaining the same Lagrangian–Coriolis force and shear. This can result in a front with dramatically different Eulerian shear than Lagrangian shear because of the anti-Stokes shear. Therefore, the instabilities that depend most on the Eulerian shear, such as SI, are more affected than those that depend mostly on the Lagrangian shear, such as KHI (Holm 1996; section 3a) and GI (section 3b).
The Charney–Stern–Pedlosky criteria for QG baroclinic instabilities are proven to be amenable to modification and to depend on the properly reinterpreted QGPV and shear. Thus, after the Eulerian flow adjusts to the arrival of waves, the Lagrangian thermal wind, and therefore stability of the mixed layer to GI, is as it was before the waves arrived. If the flow is unstable to GI, an increase (decrease) in Eulerian shear due to anti-Stokes flow appears to reduce (increase) the wavenumber and growth rate of the GI, consistent with the theory established by Nakamura (1988). Also, if the layer is unstable to GI, the work done on the GI by the Stokes shear force is fairly weak because GI extract energy primarily from the potential energy, rather than the shear.
The Hoskins (1974) criterion for SI is proven analytically for a front with Stokes drift that has constant shear. In the presence of anti-Stokes, or any ageostrophic Eulerian shear, the criterion Ri < 1 is not equivalent to the necessary criterion for SI (fQ < 0). This is because the presence of ageostrophic Eulerian shear decouples the relationship between the shear and the horizontal buoyancy gradient, which combine to form the horizontal component of the PV. When the Stokes drift is downfront (upfront), the anti-Stokes Eulerian flow maintains positive (negative) PV, resulting in less (more) favorable conditions for SI. This orientation is opposite that of the frictional injection of PV by winds (Thomas 2005). Therefore, in a common scenario in which the winds and waves are aligned, the PV due to the anti-Stokes flow and frictional injection of PV by winds will be competing effects, and in scenarios with (swell) waves propagating into the wind, these effects will conspire to maintain and inject PV of the same sign. While this rule of thumb aids in remembering the sign of the effects, it is important to remember that the wind injection of PV differs from the indirect way Stokes forces affect PV by changing frontal balances and mixing efficiency.
Parcel switching analysis shows that the anti-Stokes shear changes the path of SI to be closer to surfaces of constant momentum rather than surfaces of constant buoyancy as in inviscid, no Stokes, hydrostatic cases of SI. This results in more cross-isopycnal flow and therefore more restratification (compared to the no Stokes case) done by SI. Although this linear, inviscid result is not confirmed in the nonlinear, viscid LES, this mechanism, as well as the activation of SI by maintaining low PV (which is confirmed by the LES), is a way in which Stokes drift may indirectly cause restratification rather than mixing as is commonly thought of LC.
The analytic and parcel switching results are confirmed with inviscid linear stability analysis and again with viscid nonlinear LES. The linear stability analysis shows that the Stokes shear force does work against SI when the Stokes and Eulerian shears have opposite sign, and this result is confirmed by the LES. The robustness of the Stokes-modified PV, and induced SI, suggests that this is the dominant effect of Stokes drift on baroclinic instabilities in the mixed layer.
Last, there are several unaddressed physical processes that occur in both the linear stability analysis and the LES. These include LC under the influence of a front, the Stokes-modified Ekman layer, and how it modifies SI and GI. These topics will be addressed in a forthcoming paper that focuses on the instability mechanisms that depend on viscosity in the mixed layer.
The authors thank Prof. Gregory Chini for insightful discussions of interactions between symmetric instability and Langmuir circulations, Prof. Eric D’Asaro for suggesting that the submesoscale instabilities of the boundary layer including Stokes forces be explored, and Dr. Peter Sullivan for help with the NCAR LES model. BFK, KJ, AW, and SH were partially supported by NASA NNX09AF38G and NSF OCE 0934737 and 1258907. This work utilized the Janus supercomputer, which is supported by the National Science Foundation (Award Number CNS-0821794) and the University of Colorado Boulder. The Janus supercomputer is a joint effort of the University of Colorado Boulder, the University of Colorado Denver, and the National Center for Atmospheric Research.
The nondimensional, wave-averaged, Boussinesq equations with the choice of scalings in Table 2 are given without approximation below. Approximations to these equations are made in section 4 to solve for the unstable modes of a front with Stokes drift. Other approximations are made by Haney (2015) to reproduce the flow regimes in Table 3:
where, recall from Table 1, , , , and ; , , ; and , with denoting dimensional variables.
This article is included in the Ocean Turbulence Special Collection.
In this work, geostrophic instabilities will refer to instabilities roughly the size of the deformation radius, which need not be strictly geostrophic and are primarily driven by the potential energy in a front. Geostrophic instabilities are often referred to as baroclinic instabilities; however, because symmetric instabilities are also formally baroclinic, the term geostrophic is preferred here because of their similarities with quasigeostrophic baroclinic instabilities. Furthermore, throughout this work baroclinic instabilities will refer to both symmetric and geostrophic instabilities.
The subsurface Stokes drift autocorrelation times were found to vary with depths examined (down to 9 m depth). In general, the zero crossing and zero sum times increased and decreased with depth, respectively, with a 70-day minimum for the latter.
One can write down the equivalent Richardson number criterion , but its usefulness is limited because it is distinct from any other Richardson number since it has a mix of two shears. For example, RiPV is not expected to be relevant for Kelvin–Helmholtz instability. As such, and to avoid confusion over what shear is required in the Richardson number, the PV criterion is highly preferred.