The dominant processes governing ocean mixing during an active phase of the Madden–Julian oscillation are identified. Air–sea fluxes and upper-ocean currents and hydrography, measured aboard the R/V Revelle during boreal fall 2011 in the Indian Ocean at 0°, 80.5°E, are integrated by means of a large-eddy simulation (LES) to infer mixing mechanisms and quantify the resulting vertical property fluxes. In the simulation, wind accelerates the mixed layer, and shear mixes the momentum downward, causing the mixed layer base to descend. Turbulent kinetic energy gains due to shear production and Langmuir circulations are opposed by stirring gravity and frictional losses. The strongest stirring of buoyancy follows precipitation events and penetrates to the base of the mixed layer. The focus here is on the initial 24 h of an unusually strong wind burst that began on 24 November 2011. The model shows that Langmuir turbulence influences only the uppermost few meters of the ocean. Below the wave-energized region, shear instability responds to the integrated momentum flux into the mixed layer, lagging the initial onset of the storm. Shear below the mixed layer persists after the storm has weakened and decelerates the surface jet slowly (compared with the acceleration at the peak of the storm). Slow loss of momentum from the mixed layer extends the effect of the surface wind burst by energizing the fluid at the base of the mixed layer, thereby prolonging heat uptake due to the storm. Ocean turbulence and air–sea fluxes contribute to the cooling of the mixed layer approximately in the ratio 1:3, consistent with observations.
The intraseasonal variability of sea surface temperature (SST), winds, and outgoing radiation in the equatorial eastern Indian Ocean is dominated by the 30–90-day period of the Madden Julian oscillation (MJO; Hendon and Glick 1997; Wheeler and Hendon 2004; de Szoeke et al. 2015). The active phase of the MJO features an area of increased westerly wind anomalies (Madden and Julian 1971, 1972) near the surface and strong precipitation at the equator, propagating eastward at ~5 m s−1 (Zhang 2005). These westerly wind anomalies, or wind bursts, last only 1–3 days (Zhang 2013). Despite their short duration, wind bursts can account for the majority of the climatological momentum flux into the ocean. The mean wind forcing is weak ≲0.05 Pa (Wyrtki 1973; Schott and McCreary 2001) and highly variable. During a westerly wind burst (WWB), surface stress is typically ~0.5 Pa and can reach values >0.8 Pa in 1-min averages. In addition to dominating the momentum exchange, the active phase of the MJO changes the nature of the air–sea heat and freshwater exchanges due to increased latent heat flux and precipitation at the surface (de Szoeke et al. 2015).
The intense surface fluxes associated with a westerly wind burst change the physics of heat, salt, and momentum transport within and below the ocean mixed layer. In the MJO suppressed phase, heat exchanges are dominated by a diurnal cycle of strong daytime solar heating ~400 W m−2 and weak ≲100 W m−2 cooling with weak precipitation and winds (Moum et al. 2014; de Szoeke et al. 2015). In addition to direct heat, salt, and momentum fluxes, mixed layer transports are partially mediated by surface waves. Typical waves in the equatorial Indian Ocean are too shallow, 1–2 m high with periods >10 s (Young 1999; Chen et al. 2002; Sterl and Caires 2005), to drive significant motion in the mixed layer. During a westerly wind burst, in addition to strong momentum flux, there is strong precipitation, net surface cooling of 300 W m−2 and wave heights increasing to ≳1.5 m with periods as short as 3 s (Moum et al. 2014). The strong fluxes freshen, cool, and accelerate the surface water, while steep waves, shear, and buoyancy fluxes energize mixing (Fig. 1) and exchange salt, heat, and momentum with deeper water previously inaccessible to air–sea interaction. Observations during a westerly wind burst show the potential influence of a barrier layer on the exchange of heat across the mixed layer base (Chi et al. 2014; Moum et al. 2014).
In this study, our main goals are 1) to identify the most important subsurface mixing mechanisms and 2) to assess the relative importance of subsurface and atmospheric processes in determining the surface cooling that damps the storm. To this end, we explore atmosphere–ocean feedbacks and subsurface mixing processes in the MJO active phase using a large-eddy simulation (LES).
After initial applications to the atmospheric boundary layer (Deardorff 1972), LES was adapted for the ocean in the 1990s (Skyllingstad and Denbo 1995) and has been used in numerous studies since. Skyllingstad et al. (1999, hereinafter S99) explored the limits of the technique by comparing LES-derived turbulence statistics with microstructure measurements in the context of a westerly wind burst. Several modeling studies using idealized (McWilliams et al. 1997; Wang and Müller 2002; Harcourt and D’Asaro 2008; Grant and Belcher 2009; Noh et al. 2011) and empirical (Kukulka et al. 2009, 2010) forcing have contributed to the general understanding of upper-ocean physics. LES of shallow mixed layers driven by strong wind forcing show that turbulence near the surface is driven mainly by shear associated with the Stokes drift of the surface waves (Langmuir turbulence; McWilliams et al. 1997). Langmuir circulation transports momentum and buoyancy more quickly than a simple shear-driven mixed layer (Kukulka et al. 2009, 2010), rapidly mixing away shear near the surface (Noh et al. 2011). Away from the surface, Langmuir turbulence relies on the inertia of vertical motions of Langmuir cells (Grant and Belcher 2009) because the driving Stokes shear decays rapidly with depth, unlike turbulence driven by convective or shear instability, which need not depend on depth (Thorpe 2004). The influence of the Langmuir turbulence is present in much of the mixed layer with the effect strongest in shallow mixed layers with weak underlying stratification (Noh et al. 2011). Langmuir turbulence may deposit momentum at the mixed layer base, setting the stage for enhanced shear production (Kukulka et al. 2010).
Observational estimates of surface heat flux, precipitation, and wind stress can be used to drive ocean LES of a dynamically evolving mixed layer and inform the analysis of the subsurface mixing mechanisms of a specific event. S99 did this for a westerly wind burst observed in the western Pacific using measurements from the Coupled Ocean–Atmosphere Response Experiment (COARE). Here, we use the same approach to examine mixing a wind burst observed in the equatorial Indian Ocean at 0°, 80.5°W during the Dynamics of the Madden–Julian Oscillation (DYNAMO) field campaign (Yoneyama et al. 2013; Moum et al. 2014; de Szoeke et al. 2015).
We will first discuss how the boundary and initial conditions were extracted from the observational record (section 2) and estimate the relative influence of convection and Langmuir turbulence on the mixed layer before summarizing the numerical model (section 3). The simulation is used to partition turbulent energy production among Stokes-, shear-, and buoyancy-driven production (section 4b). Special attention is given to shear instability below the mixed layer (section 4c). Vertical transports of momentum, heat, and salt are compared to the surface forcing (section 4d). Results are compared with the COARE wind burst (section 5), and conclusions are summarized in section 6
2. Initial ocean state and air–sea fluxes during the wind burst
Initial conditions for the LES model are defined as follows: Vertical profiles of horizontal velocity are obtained from the acoustic Doppler current profiler at the nearby equatorial Research Moored Array for African–Asian–Australian Monsoon Analysis and Prediction (RAMA) mooring (McPhaden et al. 2009) at 80.5°E. We assume uniform currents above 12-m depth due to the unreliability of near-surface data. The Chameleon vertical microstructure profiler provides profiles of salinity and temperature (Moum et al. 1995, 2014; Pujiana et al. 2015). Potential density is calculated using the Gibbs Seawater package version 3.02 (McDougall and Barker 2011) and referenced to surface pressure. Depth Conservative Temperature Θ and Absolute Salinity SA, as functions of density at surface pressure ρ(Θ, SA, P = Patm), are calculated by interpolating in density using a cubic polynomial. The observed currents and hydrography reveal a strong tidal signature that is not included in the model. To remove the resulting bias, isopycnal averages of velocity, temperature, salinity, and depth are computed over the 5 days preceding the wind burst.
After the MJO suppressed phase that preceded the wind burst, the upper ocean was stable and the currents were weak (Fig. 2). A layer of stable stratification near z = −20 m was mostly due to a stable salinity gradient of gβ∂zS ~ 4 × 10−4 s−2 embedded in weaker temperature gradient gα∂zT ~ 10−4 s−2. This halocline coincides with the upper edge of a meridional current extending from 20 to 50 m, carrying salty water into the Southern Hemisphere (Fig. 2). Below the strong halocline, the salinity continues to increase, but temperature is the dominant source of density stratification. Salinity reaches a maximum at 50 m, below which the temperature gradient increases to gα∂zT ~ 4 × 10−4 s−2. The near-surface current was dominated by its eastward component and exhibited very little shear except at z = −20 m (Fig. 2b, solid line).
Momentum fluxes are taken from 1-min average wind speed observations measured aboard R/V Revelle. Surface heat and salt fluxes include contributions from the observed net precipitation and evaporation P − E and the surface heat flux J0 calculated using the observed winds, humidity, radiative fluxes, and air–sea temperature differential. Observations are converted to surface fluxes using COARE 3.5 revision 3 (de Szoeke et al. 2015). The processed fluxes are then low-pass filtered at 1 h.
In situ air–sea fluxes (de Szoeke et al. 2015) and subsurface profiles are available from the DYNAMO project from September 2011 to January 2012 while R/V Revelle was on station at 0°, 80.5°E. The strongest measured heat flux at the base of the ocean mixed layer of the DYNAMO record is associated with the wind burst at the end of November 2011 (Chi et al. 2014). We choose to model the 24 November 2011 westerly wind burst because it exemplifies critical atmosphere–ocean feedbacks and was brief enough to be accessible to high-resolution modeling.
The model is initialized on yearday 328 (24 November) at 0000 UTC, shortly before local sunrise. The previous suppressed phase of the MJO was characterized by weak winds (τ ≲ 0.1 N m−2) and strong daytime heating (~−400 W m−2). On 24 November, the zonal wind stress increased rapidly to τ ~ 0.5 N m−2 (Fig. 3a) with peaks as high as 0.83 N m−2 in 1-min averages. This wind event was accompanied by significant precipitation (15 mm h−1) and surface cooling (+400 W m−2) (Fig. 3b). The surface buoyancy flux (Fig. 4b) was dominated by its thermal component (~10−7 m2 s−3) except during rain squalls (~10−6 m2 s−3) when the saline component was as great or greater.
During the wind burst, the spectral distribution and vertical attenuation of the shortwave radiative heat flux J were measured (C. Ohlmann 2011, personal communication) and found to be nearly constant. These were consistent with the Paulson and Simpson (1977) formula:
with coefficient A1 = 0.69 with penetration depths λ1 = 1.1 m and λ2 = 23 m.
Surface wave effects are parameterized as a function of Uwind, with the surface wind at 19.5 m, by assuming an equilibrium sea state (Pierson and Moskowitz 1964; Li and Garrett 1993; Harcourt and D’Asaro 2008). In this approximation, the Stokes drift velocity uS and the e-folding depth are given by
During the wind burst, the surface Stokes drift remained near 0.2 m s−1 (Fig. 4a), while its vertical e-folding scale was 2–4 m. Both the turbulent Langmuir number [where uS is the Stokes drift speed (|us|)] and the Hoenikker number Ho = 4B0ρ0L/uSτ (where B0 is the net surface buoyancy flux) were less than unity during the storm, indicating quantitatively that Langmuir turbulence is likely to have been a factor (Li and Garrett 1995; McWilliams et al. 1997).
3. The LES model
The capabilities and limitations of upper-ocean LES were established by S99, who simulated the ocean response to a westerly wind burst observed in the equatorial Pacific (Smyth et al. 1996a,b) and carried out statistical comparisons between the modeled turbulence and concurrent microstructure observations. Statistically, the turbulent kinetic energy dissipation rate ε was found to agree very well with the microstructure measurements under two conditions: First, a spinup period of a few hours is required to produce realistic dissipation rates. Second, ε can be underestimated in strongly stratified layers where the model grid fails to resolve the Ozmidov scale. The resolution requirement for accurate turbulent fluxes is much less stringent since flux-carrying motions are resolved explicitly.
Our ocean LES model is essentially the same as that used by S99.1 The model equations include the surface waves by the inclusion of Stokes drift , Coriolis effect, and buoyancy b = −g(ρ′/ρ0) due to variations in potential density ρ′ = ρ(T, S, Patm) − ρ0 from temperature T, salinity S, and atmospheric pressure Patm using the equation of state from United Nations Educational, Scientific and Cultural Organization (UNESCO 1981). The model does not include mixing due to surface wave breaking. The governing equations for velocity u, temperature T, and salinity S are
where εijk is the alternating tensor (not to be confused with the turbulent dissipation rate defined below), δij is the identity tensor, p = P/ρ0 is the normalized pressure, Ω is the planetary rotation, and νt is the eddy viscosity from the subgrid-scale model of the respective field. Profiles of the Stokes drift uS and radiative heat flux J are obtained from observations using Eqs. (1)–(3).
The eddy viscosity νt is estimated using an anisotropic Smagorinsky closure (Ducros et al. 1996; Wilcox 2006) detailed in the appendix. An additional hyperviscosity term α∇12ui is included for numerical stability to remove variance at the grid scale, which is not absorbed by the Smagorinsky eddy viscosity. Both the turbulent Prandtl number Pr and Schmidt number Sc are assumed to be 0.6. The LES is conducted in a 256 m × 256 m × 60 m horizontally periodic domain with 0.5-m cubic cells. Equations (4)–(7) are advanced in time as in S99.
4. Simulation results
a. Spatial organization of turbulence
In our LES, subsurface turbulence develops immediately after the sharp increase in wind stress (Fig. 5). A few hours into the wind burst, the vertical velocity field shows considerable complexity, but at least two distinct, coherent flow geometries are evident (Fig. 5b). Near the surface, we see periodic bands of upwelling and downwelling. These are coherent over 50–150 m and are spaced at 5–10-m intervals. A range of orientations is visible, but the longest bands are oriented ~20° from the zonal, as expected for Langmuir cells (Leibovich 1983; McWilliams et al. 1997; Thorpe 2004). At the base of the surface mixed layer (~30-m depth) are upwelling and downwelling bands oriented at ~70° north of zonal with a wavelength of 128 m (half the domain extent). In this section, we will diagnose the driving mechanisms of the modeled turbulence and quantify the resulting vertical fluxes. In the process, we will explain the patterns seen in Fig. 5.
b. Evolution of turbulent kinetic energy
As we are interested in turbulent mixing processes, the simplest metric is the mean turbulent kinetic energy tke (Fig. 5a). To derive the governing equation for , we decompose velocity ui, pressure p, buoyancy b, and diffusivity νt into a horizontal mean and a perturbation . The momentum equation [Eq. (4)] is multiplied by and averaged over x and y, resulting in
is the viscous dissipation rate and
is the subgrid-scale flux.
We distinguish between four generation/dissipation mechanisms and three transport mechanisms of turbulent kinetic energy. Surface waves interact with turbulent eddies through the shear of the Stokes drift. We refer to this process as Stokes production; StP = ∂zus〈u′w′〉x,y. Where Stokes production dominates, the turbulence is referred to as Langmuir turbulence (McWilliams et al. 1997). Similarly, the mean shear can exchange energy with turbulent eddies through the mechanism of shear production; . The conversion of potential energy to turbulent kinetic energy is quantified in the buoyancy production term BP = 〈b′w′〉x,y. Eddy viscosity dissipates turbulent kinetic energy at the rate ε. The model distinguishes between the resolved advection of turbulent kinetic energy 〈w′tke〉x,y, pressure work 〈w′p′〉x,y, and transport at subgrid scales 〈sgs〉x,y.
The sources and sinks are isolated by taking the vertical mean of Eq. (8):
In this vertically averaged form, tke is generated by Stokes and shear production in roughly equal proportion (Fig. 6b, red and dark blue curves), with the former dominating early in the wind burst and the latter dominating later. The vertically averaged buoyancy production (Fig. 6b, light blue curve) was negative, fluctuating between 10% and 20% of the shear production. This is consistent with typical oceanic values of the flux Richardson number (Osborn 1980). The hyperviscosity term is small, indicating that the Smagorinsky subgrid model is effectively absorbing the downscale energy cascade.
As the simulation progresses, turbulence spreads downward; by the end of the 30-h simulation, turbulence has spread throughout the upper 40 m. The evolution of tke features several maxima (Fig. 5a) originating at the surface and extending as deep as ~30 m (Fig. 6a). These coincide with extrema in all of the vertically integrated production terms (Fig. 6b). Close examination of Fig. 5a shows that these peaks correspond to wind maxima. The wind stress changes by ~25% from hour to hour (Figs. 3a, 5a) and these variations influence tke production through multiple routes. In the Li–Garrett parameterizations [Eqs. (2) and (3)], both the speed and the penetration depth of the Stokes drift increase with the wind and therefore so does the Stokes production term (Fig. 6b, dark blue; Fig. 7a). Wind also directly drives the mean shear and hence the shear production (Fig. 6b, red; Fig. 7b). The final tke maximum (day 329, hour 6) corresponded to a shift in the wind direction, with the result that the Stokes production acted to reduce tke in competition with the positive shear production.
Langmuir turbulence is confined the upper 5 m throughout the simulation (Fig. 7a). The fact that Langmuir turbulence appears so rapidly is due in part to the assumption that the waves are always in equilibrium with the wind [Eqs. (2) and (3); Li and Garrett 1993]. In a more realistic model, the wave field, and the attendant Langmuir turbulence, might require more time to become established after the onset of strong winds. Some of the intense tke generated in the Langmuir turbulence is advected downward and deposited at 3–12-m depth (Fig. 8a) where some is lost, via negative shear production, to the acceleration of mean flow (Fig. 8b). The reversal of StP near the end of the simulation (Fig. 6b) is also evident here.
The increased surface stress accelerated the near-surface current so that the shear at its base (Fig. 9a) descended rapidly to about 20-m depth, generating a layer of positive SP over the same layer (Fig. 7b). The strongest shear coincided with a sharp pycnocline (maximum of N2; Fig. 9b), which was due to a warm fresh surface layer (Fig. 2a and accompanying discussion). Turbulence mixed the stable stratification, doing work against gravity and generating negative BP in the same layer (Fig. 7c). Dissipation ε is of comparable magnitude to generation throughout the simulation (Fig. 7d) and suggests turbulence rapidly adjusts to the relatively slowly changing driving fluxes and mean flow.
Turbulent kinetic energy is transported vertically via resolved advection, pressure work, and subgrid-scale processes. Resolved advection of kinetic energy 〈w′tke〉 (Fig. 8a) serves to move tke downward from the strong production region at the surface and away from regions of shear production late in the simulation. The pressure work 〈w′p′〉 (Fig. 8b) is limited to the region of strong Stokes drift influence. It has variable direction. It is consistently upward below regions of strong Stokes production and becomes positive late in the simulation as the Stokes production changes sign. The subgrid-scale diffusive flux is small and generally diffuses tke downward sgs (Fig. 8c).
c. Marginal shear instability below the mixed layer
Below the Stokes penetration depth, turbulence is governed largely by a competition between shear and buoyancy. These influences are quantified here using the squared shear
and the squared buoyancy frequency
The ability of turbulence to overcome gravity is consistent with the fact that the gradient Richardson number Ri = N2/S2 was small in this layer. The transfer of kinetic energy from the mean shear to small eddies require that Ri be smaller than a critical value that is approximately ¼ (Miles 1961; Howard 1961; Thorpe and Liu 2009). That condition is satisfied in a layer surrounding 20-m depth (Fig. 9c). At day 328.67, the time of the snapshot shown in Fig. 5, the shear magnitude S was at maximum at 28-m depth. At this depth, the direction of the shear vector was ~20°. The large upwelling and downwelling bands evident in Fig. 5 are oriented perpendicular to the shear direction, just as one would expect for shear instability (Smyth et al. 2013). The mechanics of shear and buoyancy production are discussed further in section 4c.
Marginal instability is a diagnostic property of sheared, stratified turbulence and is readily identified by the statistics of the gradient Richardson number. Specifically, Ri fluctuates about a critical value (often approximated as ¼) due to competing effects of forcing and dissipation (e.g., Smyth and Moum 2013). Based on Ri, the evolving upper ocean can be segregated into three distinct regimes (Fig. 9c):
During the strong solar heating and weak winds prior to the onset of the westerly wind burst, the full water column is stable (Ri > ¼; dark red regions of Fig. 9c).
At night, between precipitation events, the model develops convectively unstable stratification (Ri < 0; blue regions in Fig. 9c) in the upper 5 m. Although negative Ri indicates sheared convection, the dominant mechanism near the surface is Stokes production.
As the surface current accelerates, the shear descends and accumulates at the base of the mixed layer. The increasing shear at the base of the mixed layer is evident by the region of dynamical instability (Ri ≲ ¼; yellow and orange regions of Fig. 9c) that thickens and descends following the onset of the wind burst. This coincides with the largest values of the shear production term in the tke equation (Figs. 7b, 10b). The proximity of the mean Ri to ¼ in that region indicates marginal instability.
We next examine the statistics of Ri conditioned on the shear production rate. Figure 10 shows the fraction of time–depth regions defined by specified ranges of SP and Ri. SP bins are chosen so that each bin contains 2.5% of the time–depth points, while Ri bins are logarithmically spaced. The fraction of time–depth points in each SP/Ri bin is shown in Fig. 10b. The region of significant SP is conveniently defined to include the uppermost 35% of values, whose time–depth boundary is contoured in Fig. 7b. The Ri abundance shows a marked shift in this region (Fig. 10). The cumulative fraction of the upper 35% of SP values (Fig. 10a, solid line) shows the Ri clustered near a central value slightly in excess of ¼. In regions of weak SP, the Ri distribution (Fig. 10a, dashed line) is broader and centered near unity, suggesting that shear instability is rare.
These characteristics of the Ri distribution correspond well with observations of marginal instability associated with deep cycle turbulence in the equatorial Pacific, another example of forced dissipative turbulence driven by the dynamic instability of sheared, stratified flow (Smyth and Moum 2013).
d. Turbulent fluxes of heat, salt, and momentum
We focus on vertical fluxes across three particular surfaces. The first is the air–sea interface z = 0, where fluxes were inferred from the DYNAMO observations. The other two are alternative definitions of the lower boundary of what we casually call the mixed layer in our model. Both are defined using the density difference Δρ between the depth of interest and the surface. The first choice is Δρ = 0.01 kg m−3 (Moum et al. 2014). Above this depth, vertical fluid motion is essentially unhindered by buoyancy, so that the layer’s temperature is nearly the SST. For clarity, we will refer to this layer as the surface layer (SL); it is equivalent the “diurnal mixed layer” employed, for example, by Smyth et al. (1996a,b, 2013). We will also consider the choice Δρ = 0.1 kg m−3. This depth is relatively stable and is in particular resistant to rapid shoaling during rain events. It has been used recently to explore the subsurface heat budget during DYNAMO (Chi et al. 2014). We refer to the corresponding layer Δρ < 0.1 kg m−3 as the mixed layer (ML).
Prior to the wind burst, the turbulent heat flux is negative as sun-warmed surface water is mixed downward (Fig. 11). With the onset of strong winds, the heat flux adjacent to the surface changes from negative to positive, indicating a downward flux of surface water cooled, mainly, by the latent heat release due to wind-enhanced evaporation. After wind onset, a region of strong downward heat flux forms at the surface and descends rapidly to ~15-m depth and then gradually to ~35 m. This coincides with the ML base and also with regime 3, identified in the previous subsection, in which shear instability is active. The resulting turbulence transports cool water from the thermocline upward, exchanging it with water warmed at the surface prior to the wind burst. At the base of the SL, the modeled turbulent flux is approximately ⅓–½ of the surface flux (Fig. 11b, solid curves). This ratio is consistent with the observational estimate by Moum et al. (2014). At the base of the ML, the modeled heat flux is usually slightly larger, with values fluctuating in the range −300 ± 100 W m−2. For comparison, in the observational analyses of Chi et al. (2014), a smoothed estimate of this flux (based on a budget residual) decreased from zero to around −300 W m−2 over several days after the beginning of the storm (Fig. 12).
Like the heat flux, the salt flux shows a maximum near the ML base (see below). The flux is positive (i.e., upward), since salty deep water is mixed with fresh surface water. Near-surface flux events correspond to rain. The SL and ML bases bracket the salt flux maximum so that values across those surfaces are nearly equal. On average the salt flux in that layer is slightly higher than the surface value, suggesting that some salt is being mixed up from below, most likely from the barrier layer around 20-m depth (Fig. 2a, dashed curve), consistent with observations (Moum et al. 2014).
Sporadic rain events associated with the MJO active phase replenish the saline stratification that existed prior to the wind burst. The renewed fresh layer causes the SL base to shoal as a stratified fresh layer at the surface is formed (e.g., asterisk in Figs. 6 and 12). The fresh layer insulates the underlying water from the surface momentum flux, leading to a rapid decrease of turbulence around 20-m depth in Fig. 6 near the initial halocline. Smyth et al. (1997) have described several similar events observed during COARE. The fresh layer persists until near-surface turbulence is able to reentrain the fluid that had been part of the mixed layer before the downpour, leading to a resumption of turbulence around 20-m depth.
The model resolution is chosen so that the resolved turbulent momentum flux dominates the momentum budget (i.e., the subgrid-scale fluxes are relatively small). As the wind event begins, an area of strong momentum flux forms near the surface driven by Stokes production (Fig. 7a) and accelerates the upper 10–20 m (Fig. 13). A large fraction of the momentum is deposited near the ML base (Kukulka et al. 2010). During these events, the resulting shear is the cause of the low Richardson number and high shear production in this region (section 4c) and also of the intense heat and salt fluxes discussed above.
There is little lag between the surface stress and the momentum flux through the mixed layer. Below the ML base, the turbulent momentum flux lags the surface stress, as illustrated by the curved plumes of intense momentum flux in Fig. 13. While the present LES covers only the first 24 h of the wind burst, observations suggest that this deep turbulence remains active long after the wind subsides, potentially impacting the subsequent switch to the suppressed phase of the MJO (Moum et al. 2014).
5. Comparison with the COARE wind burst
Direct measurements of turbulence during WWBs are rare due to the technical difficulty of making small-scale measurements in heavy weather. Such observations have been made once previously, in the western Pacific warm pool as part of COARE. Here, we will compare the two cases. The COARE wind burst lasted about 3 days, and S99 modeled the second day. The observation site was slightly south of the equator (2°S) but was close enough that the inertial period was long compared with the duration of the WWB, so that Coriolis effects were negligible over the modeled interval. During that time, the wind was mostly westerly (veering to northwesterly during the last few hours) with magnitude around 0.15 N m−2. Because wind stress was significant prior to the modeling period, the initial surface mixed layer was about 50 m deep.
In comparison, the DYNAMO wind burst was brief and intense, with wind stress rising from 0 to 0.6 N m−2. Because the wind burst followed a lengthy calm period associated with the MJO suppressed phase (Moum et al. 2014), the upper ocean was well stratified and the surface mixed layer was only a few meters deep.
Both the COARE and DYNAMO simulations lasted ~24 h. In each case, the net surface heat flux changed from a daytime value around −400 W m−2 to a night time value +400 W m−2, turbulent kinetic energy production in the upper few meters was dominated by Stokes production (Fig. 7a, S99), and shear production was strong at the mixed layer base as it descended in late afternoon and early evening before relaxing to a quasi-stationary state. There was an order of magnitude difference in the quasi-equilibrium value of shear production: 10−7 W kg−1 in COARE versus 10−6 W kg−1 in DYNAMO. This difference is consistent with the combination of stronger winds and shallower mixed layer in DYNAMO.
In COARE, buoyancy production was significantly positive down to 30–40-m depth (S99, their Fig. 7c), whereas in DYNAMO the convectively unstable layer was restricted to the upper 2–3 m (Fig. 9). This is consistent with the difference in Hoenikker number: Ho ~ 1 for COARE and Ho ~ 0.1 for DYNAMO. At the mixed layer base, in both cases, buoyancy production was smaller than shear production, typically by a factor consistent with a flux Richardson number in the usual range of 0.2–0.3. Like shear production, the turbulent kinetic energy dissipation rate at the base of the nocturnal mixed layer was an order of magnitude larger in DYNAMO (10−6 W kg−1) than in COARE (10−7 W kg−1). As in COARE, the DYNAMO modeled dissipation rates were consistent with observations, typically to within a factor of 2. A fraction of the turbulence observed in DYNAMO was biogenic and therefore has no counterpart in the LES (Pujiana et al. 2015). At the mixed layer base that fraction was a few tens of percent at most, and it therefore does not alter the approximate agreement noted here.
6. Conclusions and discussion
We have analyzed the upper-ocean response to a westerly wind burst that occurred in the equatorial Indian Ocean in boreal fall 2011 as part of an active MJO phase. During the long period of relative calm and strong solar heating preceding the wind burst, a layer of strong stable stratification was established at the base of a shallow surface mixed layer. In the simulation, the upper ocean reacts to the strong surface forcing with enhanced turbulent kinetic energy production and expansion of the mixed layer. The resulting Langmuir turbulence mixed a layer extending to ≳5 m. Within this layer the vertical kinetic energy is consistent with the surface layer Langmuir number LaSL scaling of Harcourt and D’Asaro (2008) with an observed range of LaSL = 1–2.5.
Subsequently the wind-driven mixed layer deepens by entrainment to about 30 m. Turbulence at and below the mixed layer base is consistent with generation by shear instability. Convective turbulence is not a factor. Instead, the Langmuir- and shear-driven turbulence generates potential energy at a rate consistent with a flux Richardson number of 0.1–0.2 typical of ocean mixing (Osborn 1980; Moum 1996). The near surface is dominated by Stokes production, while shear production dominates turbulent kinetic energy production near the mixed layer base (Fig. 7). Mixed layer depth (Figs. 11, 12, 13) is initially comparable to the penetration depth of the Stokes drift (Fig. 3) due to stratification developed prior to the MJO active phase (Fig. 2). During the wind burst, surface buoyancy and momentum fluxes (Fig. 3) are quickly communicated into the mixed layer by Langmuir turbulence before shear or convective instabilities can grow as expected due to the low Hoenikker and turbulent Langmuir numbers (Figs. 4, 11, 12, 13). As the mixed layer extends beyond the Stokes penetration depth, mixing shifts into a region where shear dominates production of turbulent kinetic energy.
In the simulation, the heat flux across the mixed layer base is ⅓–½ as strong as the sum of latent, radiative, and sensible fluxes at the surface. This is consistent with concurrent microstructure observations (Moum et al. 2014). Interpretation of the microstructure measurements has been complicated by the presence of fish near the ship (Pujiana et al. 2015), which added significantly to the observed turbulence levels. The ichthyogenic component of mixing was strongest near 60-m depth, suggesting that observations at the base of the mixed layer (around 30-m depth) may not have been significantly contaminated. The agreement between the observed estimate of the turbulent heat flux (Moum et al. 2014) and the present results tends to support that view. Both the model and observations (Chi et al. 2014; Moum et al. 2014) show a ~0.5°C drop off in sea surface temperature. It is now clear that this cooling was controlled in substantial measure by subsurface processes. The large downward heat flux occurs despite of a preexisting saline barrier layer at the top of the thermocline (Chi et al. 2014; Moum et al. 2014).
An intriguing possibility is that storm-driven changes to sea surface temperature may be rapid enough to feed back onto the storm itself. It is established that changes in sea surface temperature can modulate MJO properties, but a clear causal relationship has yet to be demonstrated (Lau and Waliser 2012, their chapter 7). Some theories of the MJO attempt to include these as surface fluxes (Emanuel 1987; Neelin et al. 1987) or parameterize them as damping terms (Chang 1977). If the MJO is indeed sensitive to SST changes on this time scale, then MJO models must account for subsurface turbulent processes transporting heat through the mixed layer base.
Relative to the COARE wind burst previously observed in the Pacific (S99), the DYNAMO event involved a brief, intense forcing applied to a previously stable upper ocean. The result was relatively strong turbulence concentrated in a shallow layer in which shear from Stokes drift and the mean flow dominated over the convection that was important in the COARE observation.
The significant contribution of surface waves to the surface turbulent kinetic energy budget highlights the need for accurate characterization of the wave state. Inclusion of wave breaking (e.g., Sullivan et al. 2007) could further improve the fidelity of the model. Direct measurements of the sea surface height were made a part of the DYNAMO campaign (Moum et al. 2014) and may be used in place of Eqs. (2) and (3) for future studies of near-surface mixing.
This manuscript benefited from conversations with many members of the DYNAMO project. In particular, we thank Jim Moum and Aurelie Moulin for help interpreting near-surface observations; Jim Edson, Simon de Szoeke, and June Marion for their help with the observed surface fluxes; and Christopher Zappa for productive discussions about surface waves. We thank Carter Ohlmann for access to his penetrative radiation data prior to publication. The work was supported by National Science Foundation Grant OCE-1129419.
The subgrid-scale fluxes use a Smagorinsky (first order) closure model that gives the eddy velocity,
in terms of the variance of the (high-pass filtered) velocity and the grid spacing Δx (Ducros et al. 1996; Wilcox 2006). In the presence of large-scale gradients, the unfiltered velocity overestimates the local kinetic energy. To filter out large-scale gradients, we use a discrete Laplacian whose vertical second derivative term has been neglected to remove the influence of strong vertical gradients supported by stratification:
This is analogous to the method of Ducros et al. (1996) for bounded flows in which the wall-normal direction is neglected and a two-dimensional Laplacian is used to remove large-scale gradients. The resulting filter transfer function
is of the order of unity at the Nyquist wavenumber of the simulation [the average value of the filter transfer function for a cubic grid; Δx = Δy at is 1.3]. The Laplacian is iterated three times,
to strongly suppress the influence of large scales. In addition to Smagorinsky subgrid-scale viscosity, a high-order hyperviscosity term is included in Eq. (4) with a constant coefficient α = 0.001. This prevents aliasing by damping the smallest scales of motion.
Current affiliation: University of Michigan, Ann Arbor, Michigan.
To estimate the e-folding depth of the Stokes drift, S99 used visual observations of the dominant swell instead of Eq. (3). Had Eq. (3) been used, the e-folding scale would have been smaller and the Langmuir turbulence would have been stronger and concentrated even more tightly at the surface. That would not have affected the conclusions.