## Abstract

Understanding the processes that control the evolution of the ocean surface boundary layer (OSBL) is a prerequisite for obtaining accurate simulations of air–sea fluxes of heat and trace gases. Observations of the rate of dissipation of turbulent kinetic energy (*ε*), temperature, salinity, current structure, and wave field over a period of 9.5 days in the northeast Atlantic during the Ocean Surface Mixing, Ocean Submesoscale Interaction Study (OSMOSIS) are presented. The focus of this study is a storm that passed over the observational area during this period. The profiles of *ε* in the OSBL are consistent with profiles from large-eddy simulation (LES) of Langmuir turbulence. In the transition layer (TL), at the base of the OSBL, *ε* was found to vary periodically at the local inertial frequency. A simple bulk model of the OSBL and a parameterization of shear driven turbulence in the TL are developed. The parameterization of *ε* is based on assumptions about the momentum balance of the OSBL and shear across the TL. The predicted rate of deepening, heat budget, and the inertial currents in the OSBL were in good agreement with the observations, as is the agreement between the observed value of *ε* and that predicted using the parameterization. A previous study reported spikes of elevated dissipation related to enhanced wind shear alignment at the base of the OSBL after this storm. The spikes in dissipation are not predicted by this new parameterization, implying that they are not an important source of dissipation during the storm.

## 1. Introduction

The ocean surface boundary layer (OSBL) is a critical interface in the Earth system, through which heat, freshwater, momentum, and trace gases are fluxed between the atmosphere and the ocean (Belcher et al. 2012; Rippeth et al. 2014). Because of its importance, there is strong interest in understanding the processes that determine the characteristics and evolution of the OSBL (e.g., Kilbourne and Girton 2015; Aijaz et al. 2017).

Mixing by turbulence in the OSBL tends to produce a layer with relatively weak vertical gradients in temperature and salinity, which will be referred to as the well-mixed layer (WML). The WML is separated from the deeper ocean by a stratified transition layer (TL). Current shear across the transition layer may be associated with near-inertial waves (NIW) (Plueddemann and Weller 1999), which are a ubiquitous feature of the surface ocean (Pollard 1980; D’Asaro 1985), and are a significant energy source driving turbulent mixing in the ocean (Alford 2003; Watanabe and Hibiya 2002). The generation of inertial motions in the OSBL is highly intermittent, with storms providing an important source of energy (D’Asaro 1985; Large and Crawford 1995).

The shear associated with the NIW is concentrated across the stratified transition layer (Pollard and Millard 1970; D’Asaro 1985), which is often in a state of marginal stability, with a Richardson number *O*(1), so that the shear may result in the generation of turbulence in the TL (Johnston and Rudnick 2009; Rippeth et al. 2005, 2009). Observations reported by Burchard and Rippeth (2009), Lenn et al. (2011), Brannigan et al. (2013), and Lincoln et al. (2016) suggest that the generation of turbulence within the transition layer is a result of surface wind and current shear alignment which produces enhanced shear at the base of the WML. Using a one-dimensional model, Plueddemann and Farrar (2006) showed that the energy input into the NIW is balanced by the downward propagation of NIW energy and the local dissipation due to shear generated turbulence within the TL.

Generation of turbulence by shear across the TL is particularly efficient when the rate at which near-surface winds rotate due to the motion of a storm, matching the inertial period of the wind-driven currents (Large and Crawford 1995; Dohan and Davis 2011; Chen et al. 2015). This resonance condition allows large shears to build at the base of the well-mixed layer, which can lead to the growth of a stratified shear layer below the WML, without significant impact on the thickness of the WML (Dohan and Davis 2011; Johnston et al. 2016; Chen et al. 2016).

Skyllingstad et al. (2000) used large-eddy simulation (LES) to look at shear production in the TL for resonant and nonresonant situations. Compared to resonant conditions, the predicted current shear and turbulence were significantly smaller for nonresonant conditions. Grant and Belcher (2011) have derived a parameterization for the magnitude of the maximum shear production and dissipation at the base of the WML, due to resonant wind forcing.

Here we present measurements of the dissipation rate *ε*, temperature, and salinity of the water column, obtained by an ocean microstructure glider over a period of 9.5 days. The measurements were obtained during the process cruise of the Ocean Surface Mixing, Ocean Submesoscale Interaction Study (OSMOSIS) project in the northeast Atlantic in September 2012.

During the period of the observations a significant storm occurred. The aim of this study is to investigate the processes responsible for the evolution of the OSBL during the storm. To achieve this aim we will combine the profiles of *ε* from the microstructure glider, with supporting data, to test a new parameterization for OSBL mixing.

This paper is arranged as follows: section 2 gives a description of the observational campaign together with the methods used to collect the data. Section 3 provides a description of the main experimental results, section 4 describes the turbulence measurements in the transition layer, and section 5 concludes the paper with a discussion of the key results.

## 2. Observations and modeling

The observations used in this study were collected during a multidisciplinary cruise aboard the Royal Research Ship *Discovery* (Allen et al. 2013), as part of the NERC OSMOSIS project. The cruise took place in the vicinity of the Porcupine Abyssal Plain (PAP) observatory (48.69°N, 16.19°W; Fig. 1), which is to the southwest of the United Kingdom. The site is representative of the open ocean, with a water depth of ~4800 m. The measurements were made from the evening of 17 September to the evening of 27 September 2012 (year days 260–270).

The observations to be discussed were obtained using 1) an ocean microstructure glider (OMG), 2) a TRIAXYS directional wave buoy, 3) shipborne measurements of meteorological data, and 4) water velocity profiles from the vessel mounted ADCP.

### a. Ocean microstructure glider

The OMG was a Teledyne Webb Research Slocum coastal electric glider equipped with an unpumped Sea-Bird CTD sensor and a Rockland Scientific International MicroRider microstructure package. The OMG microstructure package samples shear microstructure from which estimates of the dissipation rate *ε* of turbulent kinetic energy (TKE) were determined (Fer et al. 2014; Palmer et al. 2015). During the OSMOSIS cruise, the glider profiled between the sea surface and ~100-m depth, capturing 1420 profiles over 9.5 days. A profile was obtained approximately every 10 min, with a 20-min gap for data upload every 10 profiles.

Spikes in the raw OMG microstructure shear data were removed by hand, after which the dissipation was calculated in bins of approximately 1-m vertical resolution following the methods outlined in Merckelbach et al. (2010), Fer et al. (2014), and Palmer et al. (2015). Close to the ocean surface the dissipation estimates can be contaminated by glider motions induced by surface waves. To account for this, the near-surface portion of the glider dissipation profiles have been truncated to exclude data from the surface to the deepest of 1) the significant wave height, 2) the point where the glider speed drops below one standard deviation from the median (for this deployment), and 3) the point where the glider pitch changes by more than one standard deviation from the median value (for this deployment). Typically the cutoff depth is about 5 m.

Conductivity, temperature, and depth (CTD) were provided by standard payload sensors (Sea-Bird Electronics) housed in the central section of the OMG and are used to calculate salinity and density. CTD data were collected at 1 Hz during periods when the MicroRider was operative. Errors in salinity and density may occur due to inconsistencies between temperature and conductivity sensors, which are partly attributable to the physical separation of sensors and different response times, both of which can be simply corrected for. The raw temperature data from the unpumped CTD sensor was low-pass filtered (using a third-order Butterworth filter), subsequently corrections to account for thermal inertia in the conductivity cell were made following the methods of Lueck and Picklo (1990) using modified parameters according to Palmer et al. (2015).

Temperature and salinity were then calibrated against CTD profiles obtained from the ship. Between year days 269 and 270.4, the unpumped CTD sensor on the glider failed, and for this period temperatures have been obtained from the MicroRider probe. To ensure the MicroRider temperatures were consistent with the CTD temperatures a regression was made between the two instruments when they showed the same structure, between the deep water (75–100 m) and surface waters (1–20 m), and was used to reconstruct temperature when the CTD failed. It was not possible to reconstruct the missing salinity data.

### b. The TRIAXYS wave buoy

The TRIAXYS directional wave buoy was deployed on 7 September and recovered on 27 September 2012, providing spectral energy data from 0 to 0.64 Hz in frequency bins of 0.01 Hz, with directional dependence resolved to 3° divisions, captured every 20 min. The Stokes drift $Us0$ is the integral of the third moment of the energy spectrum and is estimated following Webb and Fox-Kemper (2011), namely,

where *ω* is the wave frequency and *θ* the directional angle. In practice the wave buoy has a physical cutoff frequency after 0.64 Hz, such that the highest-frequency components of the wave field are not resolved. Thus, for each time stamp and direction, a best-fit tail is extrapolated from the cutoff frequency using a minus fifth-order power law (Phillips 1977). This directional patch for the unresolved tail is added to the resolved spectra and similarly integrated over frequency. Nondirectional Stokes drift data used in this study are taken as the absolute values of the total directional Stokes drift, given at 20-min intervals.

### c. Ship data

Atmospheric data were sampled throughout the cruise using the ship’s continuous recording instrumentation. Wind speed, direction, atmospheric pressure, air temperature, relative humidity, and upwelling and downwelling shortwave irradiances were all measured at a height of 18 m above mean sea surface. In all cases raw data was recorded at 10-s intervals. Quality control, despiking, and smoothing was applied to all data following Inall and Audsley (2012). The *u* and *υ* components of the wind were smoothed, and obvious spikes removed manually. The remaining data was then interpolated onto a regular grid and a 120-s median smoothing window applied.

The surface (air) friction velocity was calculated using

where $u*a$ is the friction velocity on the air side of the air–ocean interface. The drag coefficient *C*_{D} and 10-m wind speed *W*_{10} were obtained iteratively by applying a log-law boundary layer to adjust for measurement height (Beardsley and Pawlowicz 1999).

The surface buoyancy flux *B*_{0} was calculated from the net surface heat flux *H*_{0} as

where *C*_{T} = 1.6 × 10^{−4} K^{−1} is the thermal expansion coefficient, *g* = 9.81 m s^{−2} is gravitational acceleration, *ρ* is the water density, and *C*_{p} = 3993 J kg^{−1} K^{−1} is the specific heat capacity of water.

The net surface heat flux was calculated using

with shortwave radiation (SW) from the total incident radiation (TIR) sensors on board, the longwave (IR) radiation obtained from reanalysis data (National Oceanography Centre/University of Southampton 2008), and the sensible (SH) and latent (LE) heat fluxes obtained using the TOGA COARE 2.0 algorithm (Fairall et al. 1996). To account for shading of the irradiance sensors, values of the TIR were created by taking the maximum value recorded by the port and starboard sensors.

Currents were determined with an RDI Ocean Surveyor 75-kHz vessel-mounted ADCP, configured to sample over 120-s intervals with 96 bins of 8-m length, giving a standard deviation of 1.1 cm s^{−1}. The instrument calibration and calculations of the GPS accuracy are documented in the D381 cruise report (Allen et al. 2013). During some high wind/wave events data dropout was apparent in the ADCP current profile data, probably due to cavitation below the ship’s hull. These data were identified and masked for any 120-s epoch that was more than 35% incomplete between the surface and 200-m depth; this systematically removed most of the bad data, but the data around these incidents should be treated with caution.

### d. Modeling

The dissipation rates in the transition layer obtained from the OMG glider will be compared with a simple parameterization of the dissipation due to shear production. The parameterization is based on results from LES, and the derivation of the parameterization is given in appendix A.

A simple bulk model of the OSBL is used to determine the inertial currents in the OSBL, and the evolution of the mixed layer depth. The model is described in appendix B. The thickness of OSBL is assumed to increase through entrainment, with two parameterizations of entrainment considered. The first assumes that entrainment is driven by a combination of convective and Langmuir turbulence (this will be referred to as the Langmuir model). Parameterization of entrainment due to Langmuir turbulence have been proposed by Grant and Belcher (2009), McWilliams et al. (2014) and Li and Fox-Kemper (2017). The second parameterization assumes that entrainment is due to a combination of convective and conventional shear turbulence (referred to as the shear model). The parameterization for shear turbulence is taken from Grant and Belcher (2009), which is similar to the parameterization due to Li and Fox-Kemper (2017).

It is generally thought that Langmuir turbulence is important in the OSBL (McWilliams et al. 1997; D’Asaro 2014), and so in the main part of the paper the results from the Langmuir model will be shown. Additionally, the results from a shear model will then be considered in section 5.

The model is forced using ERA-Interim data, which includes wave data. The friction velocity, surface Stokes drift and buoyancy fluxes from ERA-Interim are in very good agreement with the estimates obtained from the onboard meteorological measurements and the Stokes drift obtained from the TRIAXYS wave buoy. Differences in the surface fluxes obtained from ERA-Interim and from the ship data contribute to the uncertainties in making comparisons between the model and the observations. In lieu of formal estimates of the surface flux errors, it was decided to use ERA-Interim to force the models, and where necessary, the ship based flux to derive estimates of entrainment fluxes from the observations.

The initial temperature and salinity profiles that are used to initialize the bulk model were obtained from the glider and resolved to a grid of 1 m. These profiles were used to provide the temperature and salinity structure below the OSBL, which was assumed to remain constant through the storm period. The initial depth of the WML was taken to be 35 m, the same as the average mixed layer depth obtained from the glider profiles at the beginning of the storm.

## 3. Surface forcing and the evolution of the OSBL

Figures 2a and 2b show time series of the surface friction velocity, Stokes drift and the surface buoyancy flux obtained from the ship and buoy data over the full 9.5 days of the glider deployment. Time–depth cross sections of temperature and the turbulent dissipation rate are shown in Figs. 2c and 2d, where the black dotted line shows the mixed layer depth determined from the temperature profiles, defined as the level at which the temperature is 0.2°C lower than the temperature at a depth of 10 m (de Boyer Montégut et al. 2004). The black line shows the depth of the base of the transition layer, which is the depth below the MLD of the deepest isopycnal that lies wholly within one standard deviation of the mean mixed layer depth [see Eq. (1) in Johnston and Rudnick (2009)]. This definition assumes that the transition layer thickness is related to the vertical displacements of the MLD and deeper isopycnals. These displacements have been found to be similar to the RMS displacement of a typical open ocean internal-wave spectrum (Johnston and Rudnick 2009; Munk 1981) which suggests that internal waves are heaving the MLD and deeper isopycnals into and out of contact with surface-intensified mixing, creating the transition layer defined here.

The first 5 days (days 261–266) of the period are characterized by relatively low winds, and an average buoyancy flux which is negative, that is, the ocean is gaining heat. During this early period the surface buoyancy flux shows a strong diurnal cycle. Although the mixed layer depth, obtained from the temperature profiles, is approximately constant with time, the depth of the OSBL implied by the thickness of the layer in which the dissipation rate is high, shows marked diurnal variations.

During the daytime, large values of TKE dissipation are generally limited to a layer near the surface, which is less than about 10 m deep. During the night, large values of TKE dissipation extend down to the stable layer at the base of the WML, and is consistent with the turbulence generated by the loss of buoyancy at the surface. This pattern of a shallow, stable OSBL during the day, followed by a deeper convective OSBL at night is repeated over several days, with the exception of day 263, when the depth of the OSBL remains elevated during the day coinciding with elevated winds and reduced buoyancy flux.

Between days 266 and 270, a significant storm passed through the area, with maximum winds speeds reaching ~20 m s^{−1} and significant wave heights of ~6 m. The beginning and end of the storm are marked by the dark gray vertical lines in Fig. 2. During the storm the sensible and latent heat fluxes at the surface increase, with an average buoyancy flux for the period of the storm which is negative, indicating cooling of the surface waters.

Figure 2c shows that the stratification at the base of the OSBL weakens during the storm. Temperature profiles obtained during the storm are shown in Fig. 3. During the storm the WML and the transition layer tend to cool, although there is significant variability between the profiles, which is probably associated with submesoscale variations, which are present in the area throughout the year (Thompson et al. 2016). Over the period the mixed layer depth increases from about 32 to 41 m while the depth of the base of the TL increases by about 4 m.

Just prior to the start of the storm (~day 266.5) there is a change in salinity of ~0.3 g kg^{−1} below the OSBL (not shown), that coincides with the temperature change at this depth (Fig. 2c), however during the storm the changes in salinity are generally small (~0.05 g kg^{−1}). At the start of the storm, Fig. 2c shows that there is a warming of the WML, which is accompanied by a reduction in the MLD.

The changes in temperature and salinity at the start of the storm may be due to advection, associated with horizontal changes in temperature, or changes in the position of glider relative to the horizontal gradients (Thompson et al. 2016). However, during the storm the data suggest that it is reasonable to consider that advective processes can be neglected, and that changes in the OSBL are primarily due to the surface forcing.

The track of the glider during the storm is shown in Fig. 1 and is consistent with the presence of inertial oscillations during the storm. The glider track shows clockwise rotations, which have a period of ~14.7 h, which is close to the local inertial period of 15.9 h at the latitude of the PAP site.

Figure 4 compares the velocity predicted by the bulk model (see appendix B) and observed currents from the ship’s ADCP. In the model it is assumed that the current below the OSBL is zero, and so there are no tidal or geostrophic currents. To isolate the wind-driven part of the current in the WML, Fig. 4 show the difference between the ADCP measured currents in the WML and at 49 m, the base of the transition layer.

The model predicts that inertial oscillations grow in amplitude from the start of the storm, and are superimposed on a mean wind-driven current. Toward the end of the storm the amplitude of the inertial oscillations is ~0.1 m s^{−1}.

The observed north–south component of the current (Fig. 4b) shows clear inertial oscillations from day 268, when the ship ADCP data are available. The amplitude of the observed oscillations are similar to those predicted by the model. The presence of oscillations is not as clear in the east–west component of the current. Because of the problems with the quality of the ADCP data during the storm, it is not clear whether the low amplitude of the oscillations in the east–west component of the current is real. With this caveat, the amplitude and timing of the inertial oscillations obtained from the model appear to be reasonable.

Figure 5 shows a time–depth cross section of the dissipation rate during the storm. The depth of the boundary between the high and low dissipation rate increases through the storm. However, in addition to the increase in the depth of the OSBL, there are also oscillations in the depth of the boundary superimposed on the overall increase.

The depth of the base of the OSBL, obtained from the bulk model, is shown by the dashed line in Fig. 5. The depth of the OSBL from the model increases steadily with time, due to entrainment. The magnitude of the deepening is consistent with the overall increase in the depth of the OSBL implied by the dissipation rate, although the model does not reproduce the higher-frequency variation in the depth of the WML shown by the dissipation rate. The high-frequency variations in the depth of the WML are probably not due to a stabilizing influence surface buoyancy flux during daytime. The effects of stabilizing surface fluxes on the depth of the OSBL are not included in the model, since in strong winds these effects should be small. In addition, the shoaling of the boundary layer due to a stabilizing surface buoyancy flux would not lead directly to changes in the depth of the stable WML determined from the temperature profiles, apparent in Figs. 2c and 2d.

Figure 6a shows the time series of the temperatures obtained from the glider, and the temperatures obtained from the Langmuir model, averaged over the depth of the OSBL. Between days 265, just before the start of the storm, and day 270, at the end of the storm, the upper ocean cools by about 1°C. However, the cooling rate is not constant with time, and during latter part of day 267, for example, the near-surface temperature actually increases.

From day 268, until the end of the storm, the change in the temperature of the OSBL obtained from the bulk model is similar to the observed change. However, before day 268 the model does not reproduce the variation in the temperature. For example, the rapid cooling observed at the start is not reproduced by the model, and the model does not reproduce the local minimum in the observed temperature on day 267. These variations in the observed temperature are consistent with the presence of submesoscale variations in temperature (Thompson et al. 2016; Whitt and Taylor 2017).

Figure 6b shows an expanded view of the period from day 268 to the end of the storm, when the winds are strongest. During this period the fluctuations in the observed temperatures about the general cooling trend are small, making this a good period to evaluate the heat budget of the OSBL. The observed cooling, estimated from a linear fit to the observed temperatures, is −0.45° ± 0.05°C (shown by the red line in Fig. 6b). The cooling predicted by the model is −0.38°C. The overall heat budget of the OSBL that is implied by the model is therefore consistent with the observed budget.

The cooling in the model is due to the surface heat flux and the entrainment flux. The entrainment flux obtained from the model is −250 W m^{−2}. This is in reasonable agreement with an entrainment flux of about −285 ± 65 W m^{−2}, estimated from the observed cooling, assuming the cooling is only due to the surface and entrainment fluxes. In the Langmuir model it is assumed that entrainment is due to a combination of convective and Langmuir turbulence. During this period, the entrainment flux in the model is primarily due to the Langmuir turbulence, which contributes about −200 W m^{−2} to the total entrainment. The entrainment flux due to convective turbulence is about −60 W m^{−2}.

The Langmuir model, described in appendix B, appears to be able to reproduce the evolution of important features of the OSBL, such as its thickness and cooling, during the storm. The Langmuir model assumes that entrainment is associated with Langmuir turbulence. Results obtained by assuming shear and convective turbulence is responsible for entrainment will be considered in section 5.

## 4. The transition layer

### a. The storm period

Figures 7a and 7b show the dissipation profiles from the glider for two periods, days 266.25–267.8 when the winds are increasing, and days 268.25–269.8 when the winds were strongest. The glider profiles have been scaled by $w*L3/hml$, where $w*L=\u2061(u*w2Us0)1/3$ is the velocity scale for Langmuir turbulence and *h*_{ml} is the depth of the WML (Grant and Belcher 2009). The structure for the OSBL from the dissipation profiles is consistent with the temperature profiles shown in Fig. 3, in that both sets of profiles show that the OSBL has two layers, the well-mixed layer and the stratified transition layer. The thickness of the transition layer appears to decrease between the two periods, due to the deepening of the WML, from 31.5 to 36 m, while the depth of the base of the transition layer remains at about 47 m.

A profile of the dissipation rate from one of the LES of Langmuir turbulence used in Grant and Belcher (2009) is also shown in Figs. 7a and 7b. In the WML the LES profile is in reasonable agreement with the observed profiles. The LES dissipation rate profile decreases more rapidly with depth in the transition layer than the observed dissipation, although this difference is less marked in the second period, due to the reduction in the thickness of the transition layer. The reduction in the thickness of the transition layer happens because, while the base of the well-mixed layer deepens, the depth of the base of the transition layer remains approximately constant. However, the thickness of the transition layer remains larger than the thickness of the transition layer in the LES profiles. Below the transition layer the observed dissipation rates are much larger than from the LES. This is because the dissipation rates below the transition layer in the real ocean are generally less than the noise level of the microstructure probe.

Skyllingstad et al. (2000) and Grant and Belcher (2011) have used LES to study the development of shear layers at the base of the WML. Skyllingstad et al. (2000) considered two cases. In the first, the surface stress rotated at the inertial period, and remained aligned with the current direction. Large shears developed across the TL, which increased in thickness. In the second case, the direction of the surface stress was kept constant, and the direction of the currents in the OSBL rotated relative to the surface stress. The shear and production of TKE at the base of the WML were much smaller than in the first case. The storm considered in this study corresponds to the second (nonresonant) case in Skyllingstad et al. (2000), since the winds associated with the present storm did not rotate at the inertial frequency.

To make a link between the inertial shear and the turbulence in the transition layer, a comparison is made with a simple parameterization of the maximum dissipation rate due to shear production at the base of the WML and the dissipation rates from the glider. This parameterization is tuned using the LES (to determine the coefficients *α* and *β*), which is described in appendix A, and is given by

where *D* is the maximum in the dissipation rate, ⟨**U**⟩_{ml} and ⟨**V**⟩_{ml} are the components of the current parallel and perpendicular to the direction of the surface stress, averaged over the depth of the WML, **U**_{ext} and **V**_{ext} are the current components below the OSBL, $u*w$ is the water-side surface friction velocity, and *h*_{bl} is the depth of the boundary layer. The surface friction velocity is defined by $u*w2=\u2212u\u2032w0\u2032\xaf$, where $u\u2032w0\u2032\xaf$ is the surface value of the momentum flux, *U*_{s0} is the surface Stokes drift, *δ* is the Stokes penetration depth, and *f* is the Coriolis parameter.

The first term in the brackets after the exponential function, on the right-hand side of Eq. (5) represents the dissipation due to shear production by the current shear in the direction of the surface stress, and is the same as the parameterization given by Grant and Belcher (2011) for the resonant conditions. The second term represents the shear production that is associated with the current shear perpendicular to the surface stress. Although $v\u2032w\u2032\xaf$ at the surface is zero, it increases rapidly with depth, and reaches maximum value just below the surface. The gradient of $v\u2032w\u2032\xaf$ above the maximum is assumed to be balanced by the Stokes–Coriolis force (Polton et al. 2005). The nonzero value of $v\u2032w\u2032\xaf$ within the OSBL leads to the production of TKE from the lateral shear at the base of the WML. Each of the components of the shear production must be greater than zero, that is, the shear terms in the TKE budget do not act as a sinks for the TKE.

Figures 8a–d show a comparison of time series of observed and predicted dissipation rate at different depths relative to the base of the WML. Figure 8a shows the dissipation rate 5 m above the base of the WML, that is, within the lower part of the WML. The dissipation rate at this depth gradually increases with time, as the surface wind increases. The curve in Fig. 8a shows the dissipation rate calculated from,

where the first term on the right-hand side is the dissipation rate due to Langmuir turbulence at the base of the WML. Grant and Belcher (2009) show for Langmuir turbulence the dissipation rate around the base of the WML is $~0.1\u2009w*L3/hml$, and the lower value for the coefficient in Eq. (6) reflects the rapid variations of the dissipation rate with depth around the base of the WML. The second term is the dissipation rate due to the convective turbulence, which is assumed to be constant with depth (Lombardo and Gregg 1989). There may be periodic variations in the dissipation rate at this level, but if they are present they are not clear.

The magnitude of the dissipation rate in the TL decreases with increasing depth (see Fig. 7). In addition, in the layer 5–10 m below the base of the WML, the magnitude of the dissipation rate shows clear oscillations in time, the amplitude of the oscillations decreasing with depth. The black curves in Figs. 8b–d show the dissipation rates obtained from Eq. (5), assuming that the thickness of transition layer is 10 m. The decrease in the amplitude of the oscillation in the dissipation rate with depth below the WML, is consistent with the amplitude of the oscillations going to zero around the base of the TL. Since Eq. (5) gives the dissipation rate at the base of the WML, to capture the decrease in the dissipation rate within the TL, the values obtained from Eq. (5) have been multiplied by 1.0, 0.5, and 0.25 in Figs. 8b–d. The magnitude of these factors is reasonable given that the dissipation rate tend to zero at the base of the TL, but having to use them means that a precise value of the coefficient in Eq. (5) cannot be assessed using the data. A more complete parameterization of the dissipation rate due to shear production in the TL would specify the depth dependence of the dissipation rate in the TL.

From the start of the storm, until the middle of day 268, the predicted dissipation rates from Eq. (5) are in good agreement with the observations. At 5 m below the base of the WML, the observations show some high dissipation rates, which may reflect the presence of the relative large dissipation rates just above this depth. After day 268.05, the dissipation rates calculated from Eq. (5) are about a factor of 2 larger than the observed dissipation rates, as shown by the dotted curves, which show results from Eq. (5), multiplied by the factors given above and divided by 2. The dotted curves match the observations during this period.

The bulk model does not give any information on the structure of the transition layer, and in estimating the dissipation rates from the parameterization the thickness of the transition layer has been assumed to be a constant 10 m. However, Figs. 7a and 7b show that the thickness of the transition layer decreases through the storm, and this may explain why, when matched to the dissipation rate before day 268.5, the parameterization overestimates the dissipation rate after day 268.5. A more sophisticated model and parameterization would be needed to predict the evolution of the structure of the TL.

Grant and Belcher (2011) found that the magnitude of the buoyancy flux at the base of the WML, due to the shear production of turbulence, was ~33% of the dissipation rate obtained from Eq. (5). Using this with Eq. (5) suggests that the heat flux at the base of the WML due to shear production between days 268.25 and 269.8 is ~−30 W m^{−2}. This is about 12% of the entrainment flux attributed to convective and Langmuir turbulence in the Langmuir model. This implies that for this storm the effects of shear turbulence on the evolution of the OSBL were small.

This analysis shows (i) that the evolution of the dissipation rate within the TL follows a different behavior to the evolution of the dissipation rate within the ML, (ii) that the oscillations in the dissipation rate are coherent through the TL, and (iii) that the period of the oscillations is close to the local inertial period, about 15.9 h.

### b. Post-storm shear spikes

Burchard and Rippeth (2009) observed that enhanced turbulence was linked to alignment of the shear across the base of the OSBL with the direction of the surface wind. They called these periods of enhanced shear, shear spikes. The observations used by Burchard and Rippeth (2009) were obtained on the continental shelf, with a water depth of ~100 m, where the currents were affected by tides. Brannigan et al. (2013) showed that shear spikes were present in the open ocean, although were unable to show that they were associated with periods of enhanced turbulent mixing. The models of Burchard and Rippeth (2009) and Brannigan et al. (2013) predict the occurrence of periods of enhanced shear, and so increased likelihood of shear instability. However, they do not predict the dissipation rate associated with the shear spikes.

Rumyantseva et al. (2015) observed similar shear spikes, with enhanced turbulence, at the PAP site on day 271, after the storm and shortly after the microstructure glider had been recovered. The surface winds were much lighter on day 271 than during the storm, but the currents continued to show the large amplitude inertial oscillations generated by the storm. Rumyantseva et al. (2015) measured the dissipation rate using an MSS90 microstructure profiler, and found that dissipation rates increased by over an order of magnitude during the shear spikes, reaching a magnitude of ~10^{−7} W kg^{−1}.

The model simulation was continued to cover the period of the observations of Rumyantseva et al. (2015) (not shown). The dissipation rates estimated from Eq. (5) for day 271 were much smaller than the observed dissipation rates during the shear spikes. The temperature–depth cross section given in Rumyantseva et al. (2015) suggests that the thickness of the pycnocline is smaller than it was during the storm, and that the temperature change across the pycnocline is smaller. The appearance of spikes in the dissipation rate during this period may be due to these changes, although they do not explain why Eq. (5) does not apply during this period. This point will be considered further in section 5.

## 5. Discussion

We have presented measurements of the dissipation rate in the northeastern Atlantic which were obtained using a microstructure glider. During the period when the glider was deployed a storm passed over the area, and the data from the microstructure glider showed that there were oscillations in the dissipation rate in the transition layer at the base of the OSBL. The period of the oscillations was close to the local inertial period. This study has compared the behavior of the dissipation rate in the transition layer with a simple parameterization of the dissipation associated with the production of TKE due to inertial oscillations in the current shear. The parameterization was derived using results obtained from LES.

Reliable measurements of the current shear were not available, and to determine the current shear needed by the parameterization, a simple bulk model (described in appendix B) was used (Niiler and Kraus 1977). This model included a parameterization of entrainment, which for the results shown in the previous sections assumed that turbulence in the OSBL was due to a combination of convective and Langmuir turbulence. In addition to predicting the currents during the storm, the model also predicted the evolution of the temperature and the depth of the base of the OSBL, which could be compared with the observations.

The change in the depth of the OSBL base, and the cooling of the OSBL obtained from the model were in reasonable agreement with the observed changes (Figs. 5 and 6). The agreement suggests that the assumption that Langmuir turbulence was present in the OSBL, and the parameterization of entrainment by Langmuir turbulence obtained from LES are reasonable.

However, variations in the near-surface temperature, which are probably due to the presence of submesoscale variability (Thompson et al. 2016), make it difficult to conclude from this comparison that the Langmuir turbulence must be present. In particular, the observed cooling of the OSBL may have been affected by advection, and it is possible that the observations are consistent with other assumptions about entrainment. It is useful to compare the results from the Langmuir model with results obtained from a model in which entrainment is assumed to be due to conventional shear driven turbulence.

Over the last 1.5 days of the storm the cooling from the Langmuir model is −0.38°C, which is similar to the observed cooling of −0.45° ± 0.05°C. The model suggests that entrainment is a significant term in the heat budget of the OSBL, and the agreement between the modeled and observed cooling depends on the parameterization of entrainment.

The parameterization for entrainment can be changed for one in which entrainment is assumed to be driven by convective and conventional shear turbulence (see appendix B). Using this model, the cooling obtained over the same period is −0.31°C. The reduction in the cooling is due to the parameterized entrainment flux being smaller than that in the Langmuir model. The entrainment flux due to shear turbulence, obtained from the shear model, is ~−116 W m^{−2} compared to ~−250 W m^{−2} from the Langmuir model. The cooling from the shear model is just about consistent with the observations, particularly if the observed cooling is influenced by more than just the surface and entrainment fluxes. In addition to the uncertainties associated with the observations, the constants used in the parameterization are derived from LES, and may differ from values that might be obtained from direct observations. From this comparison it is not possible to conclude that the Langmuir model is better than the shear model, although the Langmuir model gives reasonable results.

The change in the thickness of the OSBL during the storm obtained from the models also depends on the parameterization of entrainment. For the Langmuir model the change in the thickness of the OSBL over the storm is about 8.8 m, while for the shear model the change is 5.4 m. Figure 5 shows that there is significant variability in the depth of the OSBL base, in addition to the overall increase in the thickness. This variability may be associated with the submesoscale variability in the area of the observations. The presence of this variability in the thickness of the OSBL in the observations means it is not possible to conclude that the Langmuir model is better than the shear model, although again the results from the Langmuir model are reasonable.

What might be needed to come to a more definite conclusion? The differences between the two models used here increase with time, and in a storm of longer duration it is possible that the differences in the cooling and the change in the depth of the OSBL might become large enough for a more definite conclusion to be reached. With glider technology it should be possible to obtain more data during storms to help confirm the general presence of Langmuir turbulence in such situations, and the usefulness of LES in developing parameterizations, through studies such as this.

The agreement between the nondimensional dissipation rates within the OSBL with the results from LES of Langmuir turbulence, given by Grant and Belcher (2009), cannot be used as support for the presence of Langmuir turbulence during the storm. The reason is that, when scaled with $u*w3/hml$, Grant and Belcher (2009) showed that when the Langmuir number is ≈ 0.3, the nondimensional profiles agree with profiles of LES of conventional shear turbulence. Sutherland et al. (2014) have presented observations which suggest that the dependence of the nondimensional dissipation rate on the Langmuir number is consistent with that found by Grant and Belcher (2009), but the variation in Langmuir number in the present data is not sufficient to confirm this.

The parameterization of the dissipation rate given in Eq. (5) was obtained by assuming the base of the TL corresponds to the base of the OSBL, and that the profiles of $u\u2032w\u2032\xaf$ and $v\u2032w\u2032\xaf$ go to zero at the base of the OSBL/TL. The dissipation rate in the TL varies with time because of the interaction between the boundary layer stresses and the time varying shear across the TL due to the inertial oscillations of the currents within the OSBL. In the LES, which were used to develop the parameterization of the dissipation rate [Eq. (5)], the bulk Richardson number for the TL was between 0.2 and 0.4 (Grant and Belcher 2011), and the thickness of the shear layer increased with time. However, during the storm, the bulk Richardson number of the TL was ~4 when the turbulence dissipation was a maximum, and rather than increasing with time, the thickness of the transition layer decreased with time (see Fig. 7). Despite these differences in the stability of the TL, the comparison shown in Fig. 8 suggests that the Eq. (5) is a reasonable parameterization, although it is not obvious why this should be.

A possible reason why Eq. (5) works during the storm is that the TL is not an isolated shear layer, but is connected to the well-mixed region of the OSBL through the transport of turbulent kinetic energy into the upper part of the TL. Studies of entrainment in the sheared convective atmospheric boundary layer show that shear production of turbulence in the inversion can occur for gradient Richardson number significantly above 1/3 (Haghshenas and Mellado 2019), similar to the situation during the storm. Haghshenas and Mellado (2019) found that as the shear increases, the gradient Richardson number tends to a value of about 1/3, similar to the values in the LES (Grant and Belcher 2011). While the present situation is not directly comparable to entrainment in the atmospheric boundary layer, the interaction between the well-mixed layer and the transition layer through the transport of TKE may help explain why the results obtained with Eq. (5) are reasonable. Further studies are needed to improve our understanding of processes in the transition layer.

In this study the parameterization for the dissipation in the TL has only been used diagnostically, and the effects of the shear generated turbulence were not included in the bulk model. This is reasonable since the shear production of turbulent kinetic energy during this storm was relatively low, but in a more complete model, in which the effects of the shear production of TKE are included, a representation of the evolution of the transition layer may be possible. In particular, the evolution of the depth of the transition layer, which was simply specified in the present study, would need to be modeled.

The microstructure glider was recovered after the end of the storm, but further measurements of the dissipation were obtained after the storm using a profiler deployed from the *Discovery* (Rumyantseva et al. 2015), finding evidence of shear spikes (Burchard and Rippeth 2009; Lenn et al. 2011; Brannigan et al. 2013; Lincoln et al. 2016). Shear spikes are associated with the generation of turbulence within the transition layer, as a result of surface wind and current shear alignment which produces enhanced shear. The dissipation rate during the shear spikes was ≈10^{−7} W kg^{−1}, which is similar to the peak dissipation rates due to shear production that were observed in the TL during the storm. However, during the post-storm period the surface winds were much lighter than winds during the storm, and the dissipation rates implied by the present model are negligible, due to the dependence of the dissipation rate on $fhbl/u*w$.

The rate of work by the surface stress acting on the inertial currents in the OSBL is $u*w2U/hbl$, where *U* is the current parallel to the surface stress, averaged over the depth of the OSBL, (Grant and Belcher 2011). This can be thought of as the divergence of a flux of mean kinetic energy, where the surface flux is $u*w2U$. For the post-storm period, currents from the model give $u*w2U/hbl\u224810\u22127$ W kg^{−1}. The coincidence in the magnitude of the dissipation rate during the shear spikes and the rate at which work is being done by the surface stress acting on the inertial currents suggests that the turbulence associated with shear spikes arises directly from the breakdown of the shear at the base of the OSBL. Since the turbulence is assumed to occur because of the work done on the mean flow by the surface stress, there is no contribution to the production of TKE from the component of the current that is perpendicular to the surface stress. This is in contrast to the parameterization given in Eq. (5) where the production of TKE by the lateral shear is assumed to occur, as the steady-state momentum balance of the OSBL implies $v\u2032w\u2032\xaf$ is not zero below the surface. The changes in the properties of the pycnocline that occurred after the storm, which should have reduced the bulk Richardson number, would make the generation of turbulence from the simple breakdown of the shear possible.

The data presented in this study has been analyzed using a one-dimensional framework, and any effects due to submesoscale processes have been neglected. Whitt and Taylor (2017) have recently presented results from a large-domain LES (horizontal domain 1.9 km × 1.9 km) of the storm in this study. By coincidence their domain size was comparable to the diameter of the circular path taken by the glider during the storm (see Fig. 1a).

This simulation shows submesoscale features, with scales of order 1 km, develop during the storm, and help to maintain stable stratification within the mixed layer. The variability in the temperatures measured by the glider, shown in Fig. 6, may be due to this submesoscale variability. Whitt and Taylor (2017) also found that turbulence levels showed significant horizontal variability that was related to the submesoscale variability in the stratification. Given these results, it is reasonable to ask if the one-dimensional approach used here is valid.

The results of this study suggest, that despite the presence of submesoscale variability, the one-dimensional assumption is reasonable. Dissipation rates within the bulk of the OSBL are consistent with scaling results from LES. The bulk model produced reasonable estimates for the cooling of the OSBL, and the increase in the thickness of the OSBL. While the presence of submesoscale variability may have made it difficult to decide which form of turbulence was present in the WML (Langmuir or shear), the evolution of the OSBL could be described reasonably well using a one-dimensional framework.

## 6. Summary

This paper has used observations from a microstructure glider, together with a simple bulk model of the OSBL and dissipation in the transition layer, to help understand the turbulence and evolution of the OSBL during a storm. The key results from this study are as follows:

The OSBL has a two layer structure, a well-mixed layer separated from the deeper water by a stratified transition layer. The flow in the transition layer is turbulent, with the dissipation rate showing periodic variations, with a period close to the local inertial period (~ 15.9 h).

A parameterization of the dissipation rate due to shear turbulence was developed using results from LES. The dissipation rates obtained from the parameterization were in reasonable agreement with the dissipation rates obtained by the microstructure glider in the TL. The Richardson number for the TL was about 4, suggesting that it was too stable for shear turbulence to develop. It is possible that the transport of TKE from the WML into the TL plays a role.

The evolution of the thickness of the OSBL and its heat budget were obtained using a simple bulk model in which entrainment was assumed to be due to Langmuir turbulence. The parameterization of entrainment due to Langmuir turbulence was obtained from LES. The model results suggest that cooling of the OSBL due to entrainment was significant. Although the observations are consistent with this model, and the presence of Langmuir turbulence, the duration of the storm was not long enough conclude that this model was better than one in which entrainment is assumed to be due to conventional shear turbulence.

The parameterization of the dissipation rate developed in this study did not predict the occurrence of dissipation due to shear spikes after the storm. It is not clear why, but it was noted that the observed dissipation rate during the shear spikes was comparable to the rate of work done by the surface stress on the mean currents in the OSBL, and that changes to the pycnocline probably made the breakdown of the shear likely. However, further work is needed to understand the generation of turbulence associated with shear spikes.

The parameterizations for entrainment in the bulk model were obtained from LES. Since the constants in the bulk model were obtained from the LES it was not necessary to tune the model to other observations, and the results from the study provide an example of the usefulness of LES in developing parameterizations.

## Acknowledgments

We are grateful to the engineers and scientists at NMFSS, captain and crew of the RRS *Discovery*, and numerous scientists and technicians that helped during deployment and recovery of the glider and wave buoy and aboard the RSS cruise D381b. All data are archived at the British Oceanographic Data Centre. This research was funded by grants from the Natural Environmental Research Council OSMOSIS (NE/1020083/1) and FASTNEt (NE/1030224/1).

### APPENDIX A

#### Parameterization of Shear Production at the Base of the OSBL

In this appendix a parameterization of the shear production of TKE at the base of the OSBL is developed. The parameterization is based on the work of Grant and Belcher (2011), which is extended to include the effects of inertial oscillations in the OSBL. Grant and Belcher (2011) only considered the case where the direction of the surface stress and the currents in the OSBL were aligned and constant in time (which was approximated by setting the Coriolis parameter to zero). They developed the following parameterization for the generation of turbulence kinetic energy (TKE) by the shear,

where *S*_{U} is the shear production, **U** is the component of the current in the direction of the surface stress, ∂**U**/∂*z*|_{ml} is the shear at the base of the well-mixed layer, $u\u2032w\u2032\xafml$ is the turbulent momentum flux at the base of the well-mixed layer, $u*w$ is the surface friction velocity of water, ⟨**U**⟩_{ml} is the average velocity in the well-mixed layer, *h*_{bl} is the depth of the boundary layer, and **U**_{ext} is the current velocity at the base of the boundary layer. The surface friction velocity is defined as $u*w2=\u2212u\u2032w0\u2032\xaf$, where $u\u2032w0\u2032\xaf$ is the surface value of the momentum flux.

When the currents and surface stress rotate, the total shear production is

where *S* is the total shear production, **V** is the component of the current perpendicular to the surface stress, $v\u2032wml\u2032\xaf$ is the lateral component of the momentum flux at the base of the well-mixed layer (and $u\u2032wml\u2032\xaf$ is the component in the direction of the surface momentum flux), and ∂**V**/∂*z*|_{ml} is the shear in the lateral component of the current, at the base of the well-mixed layer.

Following Grant and Belcher (2011), $u\u2032wml\u2032\xaf$ is parameterized as $\u2212u*w2\u2061(1\u2212hml/hbl)$, where *h*_{ml} is the depth of the mixed layer (which for the LES results is defined as the level of the minimum buoyancy flux), and ∂**U**/∂*z*|_{ml} ∝ (⟨**U**⟩_{ml} − **U**_{ext})/(*h*_{bl} − *h*_{ml}), so that the first term of Eq. (A2) is given by Eq. (A1).

At the surface, $v\u2032w0\u2032\xaf=0$, but within the WML $v\u2032w\u2032\xaf$ is positive (in the Northern Hemisphere), with a maximum in the lower part of the OSBL (Zikanov et al. 2003). The maximum value of $v\u2032w\u2032\xaf$ can be estimated by considering the balance between the Coriolis force on the near-surface drift and stress gradient. The magnitude of the near-surface drift is about 3% of the wind speed (Wu 1983), and in the open ocean is mainly due to the Stokes drift that is associated with the surface waves (Wu 1983). Taking the depth of the maximum in $v\u2032w\u2032\xaf$ to be ~*δ*, where *δ* is the Stokes penetration depth, the balance between the stress gradient and the Stokes–Coriolis force gives,

where *f* is the Coriolis parameter and *u*_{s0} is the surface Stokes drift.

Using Eq. (A3) the shear production associated with the velocity component perpendicular to the surface stress can be parameterized as

where ⟨**V**⟩_{ml} is the average of **V** over the well-mixed layer, **V**_{ext} is **V** at the base of the OSBL, Δ*h* is the thickness of the pycnocline, and *β* is a coefficient.

The total shear production is the sum of *S*_{u} from Eq. (A1) and *S*_{υ} from Eq. (A4) and is given by

The results from the LES are consistent with the shear production associated with each of the velocity components having to be greater than zero, which is represented in Eq. (A5) by the MAX functions.

Figure A1a shows the time series of maximum shear production from the LES, compared to the parameterized shear production determined from Eq. (A5). The time in Fig. A1 has been normalized by the inertial period, *T*_{I} = 2*π*/*f*. The LES used to obtain the estimates of the shear production was the same as the simulations described in Grant and Belcher (2011), but with the Coriolis parameter set to 0.25 × 10^{−4} s^{−1}. The averages over the well-mixed layer, the thicknesses of the well-mixed layer and the shear layer that are needed to calculate the shear production from Eq. (A5) were also determined from the LES for this comparison.

For the first 0.4*T*_{I} the shear production from the LES is large and approximately constant. From 0.4*T*_{I} to 0.7*T*_{I} the shear production decreases, becoming constant after 0.7*T*_{I}. Before 0.7*T*_{I}, reasonable agreement between the shear production obtained from Eq. (A5), and the shear production from the LES, is obtained with *α* = 0.2 and *β* = 1.5.

Before 0.7*T*_{I} the time variation on the shear production is associated with the rotation of the inertial current with respect to the surface stress. After 0.7*T*_{I} the shear production associated with the inertial shear is zero, and the shear production is due to current shear in the well mixed layer that is associated with Langmuir turbulence (Grant and Belcher 2009).

The value of *α* = 0.2 in the parameterization of the shear production is smaller than *α* = 0.4 obtained by Grant and Belcher (2011), and suggests that *α* is a function of a nondimensional parameter. The most obvious candidate for this nondimensional parameter is $fhbl/u*w$, the ratio of the depth of the OSBL to the Ekman depth $u*w/f$.

Figure A1b shows the values of *α* obtained by Grant and Belcher (2011), the present value and the value obtained from a third LES with *f* = 0.5 × 10^{−4} s^{−1}, as a function of $fhbl/u*w$. The parameter *α* decreases as $fhbl/u*w$ increases. The curve in Fig. 6 shows $\alpha =0.4\u2009exp\u2061(\u22124.5\u2009fhbl/u*w)$, and is used as an approximate parameterization for *α* in the main text.

The observations used in the present study are of the dissipation rate *ε* rather than shear production. Grant and Belcher (2011) found that the dissipation was about 75% of the shear production. Assuming that this holds when the Coriolis parameter is not zero, the parameterization of the dissipation rate is

### APPENDIX B

#### The Bulk Model

In this appendix we describe the simple bulk model of the OSBL that is used in the body of the paper to understand the evolution of the OSBL shown in the observations. The inertial shear needed to estimate the dissipation is obtained using a bulk model of the OSBL. Following Niiler and Kraus (1977), this model predicts the time evolution of the buoyancy and components of the current averaged over the OSBL. From the observations, the dissipation rate at the base of the OSBL is small, and so the model does not include the effects of shear turbulence on the evolution of the OSBL. The equations for the currents and buoyancy are

where $u\u2032w0\u2032\xaf$, $v\u2032w0\u2032\xaf$, and $w\u2032b0\u2032\xaf$ are the surface momentum and buoyancy fluxes and $w\u2032bent\u2032\xaf$ is the buoyancy flux due to entrainment. In Eqs. (B1) and (B2), the components of the current and turbulent fluxes are relative to a fixed, geographic frame.

Following Grant and Belcher (2009), the entrainment buoyancy flux is parameterized as

where the first term on the right-hand side is the entrainment flux associated with the forcing by the surface buoyancy flux and the second term represents entrainment due to Langmuir turbulence. In Eq. (B4) it has been assumed that the total entrainment due to the combination of convective and Langmuir turbulence is just the sum of the individual entrainment fluxes.

From the results in Grant and Belcher (2009), the entrainment flux due to shear turbulence can be parameterized as

The equation for the depth of the base of the OSBL is

where Δ*B* = ⟨*B*⟩_{ml} − *B*_{ext}.

Equation (B6) is describes the evolution of the depth of an entraining boundary layer in mixed-layer models (Niiler and Kraus 1977). The effects of shear turbulence are expected to be small for this storm, and so the effects of current shear on the depth of the OSBL have not been included in Eqs. (B4) and (B6). Given the strong winds during the storm, the shortwave irradiance was not large enough to lead to the formation of a shallow, stable boundary layer, and the model only considers the evolution of the depth of the OSBL due to entrainment.

## REFERENCES

*Dynamics of the Upper Ocean*. Cambridge University Press, 336 pp.

_{2}in temperate seasonally stratified shelf seas

## Footnotes

^{f}

These authors contributed equally to this work.

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