Several features of Langmuir turbulence remain unquantified despite its potentially large impacts on ocean surface mixing. For example, its vertical velocity variance, expected to be proportional to based on numerical simulations, was proportional to in recent field observations, where is the friction velocity and is surface Stokes velocity. To investigate unquantified features of Langmuir turbulence, we conducted a field experiment around a marine observation tower in a shallow sea off the southern coast of Japan in early winter when winds and waves (often swells) were often misaligned. Coherent structures similar to Langmuir cells were successfully identified in the horizontal and vertical structures of turbulent flows measured with upward- and horizontally looking acoustic Doppler current profilers (ADCPs). ADCPs and several anemometers attached at the tower showed that turbulent vertical velocity variance was large when the Langmuir number and Hoenikker number (; where B is surface buoyancy flux and H is the water depth) were both small and that the orientation of the cells was generally aligned in the direction of Lagrangian current shear. These results agree well with the previous numerical results. As in the previous observations, however, the vertical velocity variance appeared to be proportional to . In our experiment, this curious feature was explained by compensatory effects between waves and convection. Misaligned wind with waves also seems to characterize the observed Langmuir turbulence, though further quantitative analysis is required to confirm this result.
Turbulent mixing in the ocean surface boundary layer (OSBL), a key process for the transfer of momentum, heat, and gases between the ocean and the atmosphere, is driven by winds (shear-driven turbulence), destabilizing buoyancy flux (convective turbulence), and surface waves (wave breaking and Langmuir turbulence) (e.g., Thorpe 2007). Among these factors, Langmuir turbulence has garnered increased attention, as several features of Langmuir turbulence remain unquantified despite its potentially large impacts on ocean surface mixing, sea surface temperature, and consequently Earth’s climate. An example of this impact is found in the Southern Ocean where mixed layer depth simulated in ocean and coupled general circulation models is significantly improved by proper parameterization of the effects of Langmuir turbulence (e.g., Belcher et al. 2012; Li et al. 2016).
Typical features of Langmuir circulations (LCs), revealed by early field observations (e.g., Weller et al. 1985), were successfully simulated by the large-eddy simulations (LESs) with a vortex force (representing effects of given surface waves on wave-filtered flows; e.g., Skyllingstad and Denbo 1995; McWilliams et al. 1997). Since then, LESs have been used to investigate Langmuir turbulence (e.g., Noh et al. 2004; Grant and Belcher 2009). Turbulent kinetic energy (TKE) analysis of such simulated flows shows that when Langmuir turbulence is dominant, production of TKE due to surface waves (Stokes production) dominates over production due to winds (shear production) and convection (buoyancy production; Polton and Belcher 2007; Grant and Belcher 2009). The scales of vertically integrated shear, buoyancy, and Stokes production ( for shear production, for buoyancy production, and for Stokes production, where is the water-side friction velocity, is surface Stokes drift velocity, B is surface buoyancy flux, and H is the mixing layer depth) are used to derive nondimensional parameters representing the relative importance of the three production terms, and , where La is the (turbulent) Langmuir number (e.g., McWilliams et al. 1997) showing a ratio of shear production to Stokes production, and Ho is the number that characterizes relative contribution of buoyancy production to Stokes production (e.g., Belcher et al. 2012). Here we refer to this number as the Hoenikker number, though the e-folding depth of waves instead of water depth H is used in its original definition (e.g., Li and Garrett 1995). Langmuir turbulence is expected to dominate when La and Ho are both small. Belcher et al. (2012) estimated the frequency of Langmuir turbulence using climatological values of La and showed that Langmuir turbulence is often dominant over the global oceans. The turbulence velocity scale is also obtained from the scales of the productions defined above (because the magnitude of the production is roughly proportional to the velocity cubed): is the scale of shear-driven (wind driven) turbulence, is the scale of convective turbulence, and is the scale of Langmuir turbulence (e.g., Li and Garrett 1995; Grant and Belcher 2009).
In contrast to numerical simulations, field measurements of ocean surface turbulence are not many because quantitative measurements of the flow in the OSBL are quite difficult to obtain (e.g., D’Asaro 2001). LCs in the field are affected by more factors and hence are more complicated than those in numerical simulations where forcing and environmental conditions were usually simplified for better understanding of the essential dynamics of LCs. Consequently, field studies did not provide a consistent view of Langmuir turbulence in the field. For example, Smith (1998) reported that surface velocity variances were proportional to . In recent studies investigating vertical velocity variance (VVV) in the field (D’Asaro 2001; Tseng and D’Asaro 2004; Gargett and Grosch 2014), on the other hand, the observed VVV was found to be proportional to . The observed proportionality factor between VVV and in these observations was larger than that expected from shear-driven turbulence, indicating the influence of surface waves and/or buoyancy. Because of weak buoyancy forcing and poor correlations between VVV and the convective scale [; e.g., D’Asaro 2001], was considered as a possible explanation for the apparent proportionality between VVV and (e.g., D’Asaro 2014; Gargett and Grosch 2014). Note, however, that is equivalent to constant La, which is not always true in the oceans, as will be shown in the present study. On the other hand, Gerbi et al. (2009) made quantitative analysis of observed turbulence in the OSBL and found less significant impacts of Langmuir turbulence. These field studies are apparently curious in that they are inconsistent with the above numerical simulations that suggest frequent and global occurrence of Langmuir turbulence (e.g., Belcher et al. 2012) whose VVV is proportional to (e.g., Li and Garrett 1995; Grant and Belcher 2009).
Several studies have recently been made to fill these gaps between field measurements and numerical simulations. For example, Harcourt and D’Asaro (2008) performed LESs to investigate effects of wave parameters on Langmuir turbulence and found that the vertical decay scale of the Stokes drift and wave age affect VVV of the Langmuir turbulence. Kukulka et al. (2011) found that crosswind tidal currents modulate the Langmuir turbulence in a shallow sea. Van Roekel et al. (2012) performed LESs and found that wind misaligned with waves affects VVV of the Langmuir turbulence. Nevertheless, there still remains the gap between the observed VVV and the simulated VVV. Further investigations, especially field experiments, are still necessary to quantify the Langmuir turbulence and its possible roles in the ocean surface mixing.
To investigate unquantified features of Langmuir turbulence, we conducted a field experiment around a marine observation tower in a shallow sea off the southern coast of Japan in early winter when the dominant wind was often misaligned with waves and measured turbulent flows in the OSBL using upward- and horizontally looking ADCPs. Simultaneous measurements of surface waves and atmospheric fluxes (momentum and heat fluxes) at the tower allow us to quantify the response of ocean surface turbulence to winds, waves, and the buoyancy flux. Section 2 of this paper describes our field experiment and the data analysis. Section 3 describes the observed features of ocean surface turbulence. In this section, the dependence of VVV on La and Ho is investigated to fill the gap between observed and simulated VVVs of the Langmuir turbulence. The orientation of Langmuir cells is also identified and compared with previous studies. The TKE budget is discussed in section 4, and section 5 presents concluding remarks.
2. Field experiment and data analysis
a. Overview of field experiment
The field experiment was conducted in November and December in 2014 and 2015 at the mouth of the Tanabe Bay, Japan, facing the western North Pacific (Fig. 1). The observation tower of Shirahama Oceanographic Observatory (http://rcfcd.dpri.kyoto-u.ac.jp/frs/shirahama/En_index.html) stands at 33.71°N, 135.33°E, on top of a small underwater seamount rising from a relatively flat seabed of about 30-m depth. The top of the seamount is at about 10-m depth. Both wind sea and swell usually propagate from the west in this experimental site. The vertical decay scale of the velocities associated with the wind sea of a 5-s wave period (6.2 m) is smaller than the depth at the top of the seamount, while the vertical decay scale of the swell of a 10-s wave period (25 m) is similar to the water depth. Tidal range at this observation tower is ±1 m relative to the mean sea level (MSL) at the tower.
The tower is equipped with several meteorological sensors (http://rcfcd.dpri.kyoto-u.ac.jp/frs/shirahama/En_tower.html). A supersonic anemometer (Sonic Corporation, SAT-550) is attached at 23-m height above the MSL, recording wind velocities (three components) as well as air temperature at 10-Hz sampling interval. Downward shortwave radiation is measured at 15-m height above the MSL. A thermometer, hydrometer, and barometer are also attached at 12.5-m height to record hourly air temperature, humidity, and pressure, respectively. A microwave radar measures wave heights, while an infrared radiometer measures sea skin temperature. Below the sea surface, water temperature is measured at 5- and 10-m depth. During our experiment, water temperatures at these depth levels were nearly identical; the preexisting stratification was neutral.
In addition to the above routinely operating sensors, a total of three acoustic Doppler current profilers (ADCPs) were deployed during our field experiment. One was an upward-looking ADCP (RDI, Sentinel V; 1000 kHz) deployed at the sea floor of about 10-m depth near the tower (Fig. 1). It measured vertical profiles of velocities at 0.3-m depth intervals. The transducer of this upward-looking ADCP (UADCP) was set at about 1.5-m height above the sea floor, resulting in a measurement range from 8- to 2-m depth below the MSL. This UADCP has five transducers (one vertical and four slanted transducers) and operated for 20 min of every hour (burst mode) at a 3-Hz sampling interval from 1100 LT 11 November to 0820 LT 10 December 2014 and from 1100 LT 20 November to 0820 LT 25 December 2015. Hereafter, each burst was labeled by the elapsed time (hour) since the first burst in each year (1100 LT 11 November 2014 and 1100 LT 20 November 2015). During the burst, the data from each single ping were recorded.
The other two ADCPs are horizontally looking ADCPs (RDI, H-ADCP; 600 kHz) attached to the eastern side of a pier of the tower at 3-m depth below the MSL. These two horizontally looking ADCPs (HADCPs) have a total of six transducers and were oriented at right angles to cover ±70° azimuth range from the east; the beam direction of each transducer is 20°, 45°, 70°, 110°, 135°, 160° (clockwise from the north), respectively. The two HADCPs were synchronized to each other and pinged at a 0.75-Hz sampling interval. The HADCPs were designed to operate for the same 20 min as the UADCP, and the HADCPs began collecting data almost simultaneously with the UADCP in both 2014 and 2015. However, in 2014, HADCPs stopped pinging at 1520 LT 26 November. In 2015, HADCPs stopped pinging at 1602 LT 5 December, restarted pinging later that same day at 1824 LT, and stopped completely at 0444 LT 19 December. Consequently, in 2014, UADCP and HADCPs operated simultaneously from burst 0 (1100–1120 LT 11 November) to burst 360 (1500–1520 LT 26 November), but only UADCP operated after burst 361. In 2015, on the other hand, they operated simultaneously from burst 0 (1100–1120 LT 20 November) to burst 364 (1500–1520 LT 5 December), while they operated at different times of each hour from burst 367 [1800–1820 (UADCP) and 1824–1844 LT (HADCPs) 5 December] to burst 689 [0400–0420 (UADCP) and 0424–0444 LT (HADCPs) 19 December], and only UADCP operated after burst 690. Radial coverage of HADCP was ≤80 m with a 1-m range interval. The data from each single ping were recorded.
b. Meteorological fluxes
The vertical momentum flux in the atmospheric boundary layer [ (, ), where the overbar represents 20-min averages and the prime represents an anomaly from the average, = 1.2 kg m−3 is air density, and , and are the zonal, meridional, and vertical wind velocity, respectively] and the sensible heat flux (, where J kg−1 K−1 is specific heat capacity of air, and is air temperature) were estimated using the eddy-correlation method from the 10-Hz wind velocity and air temperature data. To reduce noise effects due to a malfunction of the logger system, a five-points running mean was applied to hourly time series of the vertical momentum flux and the sensible heat flux. Then, the heat transfer coefficient was estimated from the calculated sensible heat flux , averaged sea surface minus air temperatures , and averaged wind speed over the 20-min periods. When the temperature difference is smaller than 0.5 K, the estimates were very scattered and were determined to be less reliable. These unreliable values, in addition to unrealistically large () and small () values, were regarded as erroneous and were replaced by linearly interpolated values. These unreliable or unrealistic values accounted for 17% of the total. After smoothing, the mean and standard deviation of the estimated were 2.88 × 10−3 and 1.36 × 10−3, respectively. Assuming that latent and sensible heat have the same transfer coefficients , the latent heat flux was calculated using the bulk formula of Kara et al. (2000) with averaged humidity. The formulation described in Rosati and Miyakoda (1988) was used to estimate the longwave radiation flux. The cloud cover rates used in this formulation were estimated by dividing the hourly shortwave radiation (which is reduced in magnitude by clouds) by the largest shortwave radiation measured at the corresponding hour in November and December in 2014 and 2015 (which is assumed to represent shortwave radiation under clear sky). The net heat flux , the sum of the shortwave and longwave radiations and the sensible and latent heat fluxes, was then converted to a buoyancy flux , where α is the thermal expansion rate of seawater, and =3900 J kg−1 K−1 is the specific heat capacity of seawater.
c. Mean flow
In this study, “mean” refers to 20-min averages. Vertical profiles of horizontal mean flows were estimated using UADCP in a standard manner, using the four slanted radial velocities. [The fifth radial (vertical) velocity was not used.] The nominal error of the mean velocity was 2 × 10−3 m s−1. Note that in our study area, fish motion often contaminated the radial velocity estimation. We could not completely remove this contamination using other measured data such as echo intensity, because the echo intensity from the injected clouds of buoyant bubbles into the surface ocean was often as large as the fish-intensified echo intensity. Therefore, we checked the vertical profiles of UADCP’s echo intensity by eye and discarded the burst data if fish contamination was evident in the burst.
Horizontal distributions of horizontal mean flows at 3-m depth were also estimated using the radial velocities of the HADCPs. For this estimate, a set of two radial velocities, and where is the HADCP’s radial velocity, r is radial distance from the HADCP, and is azimuth angle of the radials, were used to estimate the horizontal velocity vector at [r, ] in the polar coordinate. A pair of and with were selected to estimate the horizontal velocity vector. Consequently, the vectors were obtained in a total of seven radial directions covering ±45° from the east. The nominal error in the horizontal velocity estimated from a single HADCP is 1.4 × 10−3 m s−1.
d. Turbulent velocity variances and dissipation rate of TKE
Turbulent flows as well as orbital motions of surface waves were dominant in the velocity field measured with the ADCPs. To extract the turbulent flow component, a low-pass filter was applied to the radial velocity time series. Based on the power spectral density of the velocity (not shown), we set the cutoff frequency to 1/30 s−1. Manufacturer-provided software (RDI’s ReadyV) reports that the nominal error in the horizontal velocity of the UADCP after averaging over 30 s was 1.3 × 10−2 m s−1. To reduce this relatively large nominal error due to short bin length (0.3 m), we performed a three-point along-beam running mean to the low-pass filtered radial velocities. The error in the radial velocity is expected to be about 3.9 × 10−3 m s−1. The variance method (Lohrmann et al. 1990; Lu and Lueck 1999) was then applied to the radial velocities of the UADCP to estimate the vertical profiles of the turbulent velocity variances (, , ) and momentum fluxes (). Measurement errors on the order of 1 × 10−2 m s−1 in the radial velocity likely contaminated the estimated turbulent quantities. In fact, application of this method to our data sometimes resulted in negative and/or . In the present study, the estimates with negative or were considered unsuccessful, and these data were discarded as erroneous. The UADCP was tilted from the vertical by ≤6°. This tilt introduced an additional error to the momentum fluxes (e.g., Lu and Lueck 1999). A detailed procedure of the tilt correction is described in the appendix where error in the estimated momentum flux due to the tilt of the UADCP is estimated to be less than 10%.
Vertical profiles of the TKE dissipation rate were estimated using the spectral method (McMillan et al. 2016),
where is the transverse Kolmogorov constant, is the horizontal velocity speed, is the power spectral density (PSD) of the fifth (vertical) radial velocity in the inertial subrange, and is the noise level in the PSD. Only the fifth radial velocity was used in estimating ε in this study because additional use of other radial velocities resulted in more frequent unsuccessful estimates due to fish contamination. Figure 2 shows examples of the PSDs. The inertial subrange where the PSD follows the Kolmogorov scaling (), a range of surface waves, and a range contaminated by noise are clearly identified.
The frequency ranges used to estimate in the inertial subrange and in the noise range were set as 5 × 10−3 s−1 5 × 10−2 s−1 (referred to as range) and 0.8 s−1 (referred to as range), respectively (where f is frequency). After the frequency ranges were converted to wavenumber ranges using mean horizontal velocity at each depth [e.g., ], we applied the regression analysis to the measured PSD in the range. The power of the least squares regression line was used as a measure of estimate quality; if the power was within ±50% of the Kolmogorov power (−5/3), then the data were considered reliable. Otherwise, the data quality was considered low and the data were discarded. The observed noise-subtracted PSD  was then fitted to the Kolmogorov scaling () over the range in the least squares sense, and this fitted noise-subtracted PSD (that is proportional to ) was inserted in Eq. (1) to estimate the dissipation rate .
e. Wave height, wave spectra, and Stokes drift velocity
Hourly directional spectra of surface waves were estimated from the observed velocities, water pressures, and sea surface heights measured with the UADCP using the manufacturer-provided software (RDI’s Waves Mon). Significant wave heights from the UADCP agreed well with significant wave heights measured with the microwave radar of the tower (correlation coefficient = 0.97 and root-mean-square difference = 0.10 m), indicating that the UADCP wave measurements were of good quality. Following Kenyon (1969), hourly profiles of Stokes drift velocity were estimated from the hourly directional wave spectra estimated with the UADCP as
where ω (rad s−1) is wave angular frequency, θ is wave direction, is wavenumber vector, , g = 9.8 m s−2 is the gravitational acceleration, and H = 10 m is water depth. Here and were set to 2π/30 rad s−1 and 2π/2.5 rad s−1, respectively. The higher-frequency waves than were not considered in estimating , assuming that such shallower waves have less impact on deep Langmuir turbulence captured by our ADCPs (e.g., Gargett and Grosch 2014).
a. Overview of winds, waves, buoyancy flux, and turbulence
Time series of hourly winds (speed and direction), waves (significant wave height and surface Stokes drift speed and direction) estimated from UADCP, and the net surface heat flux as well as water temperature are displayed in Fig. 3. In both 2014 and 2015, strong winds occurred at about 4–7-day intervals (Figs. 3a,g). Strong winds were mostly associated with the northerly or northwesterly monsoon, a typical weather pattern for this season in this area. The strong southerly wind observed on 11 December 2015 (around burst 500) was accompanied by the development of an atmospheric low pressure system that traveled eastward along the southern coast of Japan. These strong winds were followed by large significant wave heights (≤2.8 m) and Stokes drift velocities (≤0.91 m s−1; Figs. 3b,h). In our field experiment, westerly (eastward) waves and Stokes drift were dominant because the western side of the bay is open to the North Pacific where the waves developed. As will be shown later, these waves were often swell. Thus, winds and waves were often misaligned. The net heat flux (Figs. 3c,i) was almost always positive (sea surface cooling), although large diurnal variations made the net heat flux negative (sea surface heating) for a short period around noon. Because strong winds induce large latent heat fluxes, the surface cooling was partly correlated to the winds. This correlation was also evident in time series of La and Ho (Figs. 3d,j). The VVV () appeared generally large when La and Ho were both small (Figs. 3e,k). The relative contributions of winds, waves, and the surface buoyancy flux to VVV, however, need to be investigated, and this will be done in section 3c.
b. Typical events
In this subsection, five events are selected from the 2015 data to illustrate the typical features of the turbulent flows measured in our field experiment.
Figure 4 shows the directional wave spectrum, vertical profiles of the Eulerian horizontal velocity, and the Stokes drift velocity as well as vertical shear of these velocities (and the turbulent Reynolds stress described in the next section) and a horizontal distribution of horizontal velocities at 3-m depth measured in burst 319 (1800–1820 LT 3 December 2015; labeled L1 in Fig. 3). The burst-averaged wind was 13.7 m s−1 in speed and north-northwesterly. The significant wave height was 1.53 m, and the wave direction was predominantly westerly. The surface Stokes drift velocity was from the west-northwest and was 0.40 m s−1 in speed. The direction of the Stokes drift corresponded to that of the largest wave component. The net surface heat flux was large and positive (561 W m−2; sea surface cooling). The Eulerian current at 3-m depth was from the north-northeast and was 0.16 m s−1 in speed (measured by the UADCP). This Eulerian current seems to be a wind-induced coastal current because similar currents were observed for similar northwesterly winds as shown later (Figs. 6, 8). The current was vertically uniform despite strong winds, except very near the surface where the shear was downwind. This indicates that the turbulence was nonlocal, that is, wave driven and/or buoyancy driven (not wind driven).
Figure 4 also shows the time–depth variations of the vertical velocity (fifth radial velocity measured with the UADCP) and the time–horizontal range variations of the radial velocities (measured with the HADCP) during burst 319. These velocities were low-pass filtered (turbulent component). The corresponding (but unfiltered) echo intensities are also shown in Fig. 4. Vertical velocities sometimes exceeded 0.1 m s−1. Strong downdrafts were accompanied by downward intrusions of the high echo intensity, probably associated with clouds of bubbles generated by surface wave breaking. Note that coherent structures are evident in the horizontal velocity and echo intensity measured with the HADCP. (Incoherent isolated large velocities and echo intensities apparent in 10–60-m range were due to fish motion.) This suggests that the turbulence was organized into rolls like Langmuir cells.
To quantify the horizontal orientation of the rolls, the observed coherent structure was fitted with a monochromatic one-dimensional sinusoidal waveform advected by the observed background mean horizontal flow. This waveform was selected to estimate a single orientation of the rolls in each burst in the simplest way. Given the wavelength L and the direction of the wavenumber vector , where θ is a direction of the roll axis, the structure (horizontal pattern) S advected by the horizontal mean current will be detected in the radial velocity of the i’s beam of the HADCP as
where r is distance from the HADCP, t is time, is the direction of the i’s radial beam, and δ is phase. Here we considered that the θ and L that provide the largest correlation coefficient between the filtered echo intensity and for correspond to the orientation and width of the observed roll, respectively. An example is shown in Fig. 5, where the filtered echo intensity measured with the HADCP and the corresponding best-fitted monochromatic wave structure are displayed. In this case, = 95° and m provides the largest correlation coefficient (0.15). Direction of beam 6 (3) corresponded roughly to the direction perpendicular (parallel) to the roll axis that was advected by the background Eulerian current in the direction of beam 6. Thus, more (less) evident coherent structure was found along beam 6 (3). The apparent mismatch between and for beam 3 and 4 was also due to irregular distribution of the observed structure and inhomogeneous horizontal velocity field. Note also that the wavelength estimation was not stable; the hourly changes of the estimated length were unreasonably large probably because of less uniform width of the rolls as shown in Fig. 5. For this reason, only the orientation of the rolls θ is considered in the following analysis.
Figure 4 also shows the estimated orientation of the observed roll. The orientation fell between the wind and the wave directions. Previous studies (e.g., Gnanadesikan and Weller 1995; Polonichko 1997; Van Roekel et al. 2012) suggested that rolls of LCs are aligned with the shear vector of Lagrangian (Eulerian plus Stokes) current, whose orientation θ was calculated by Van Roekel et al. (2012) as
where (, ) and (, ) are a difference of the wind-induced logarithmic boundary layer flow and Stokes drift between 0.12- and 8-m depths, respectively. The overall agreement between the observed and predicted orientations demonstrates that the observed rolls were primarily aligned with the Lagrangian current shear as were Langmuir cells. These observed features indicate that the turbulence in burst 319 was associated with LCs (Langmuir turbulence).
Figure 6 shows another event measured in burst 329 (0400–0420 LT 4 December 2015; labeled L2 in Fig. 3), which occurred 10 h after the event described above. At this time, the northwesterly wind reached 19.4 m s−1 in speed, and the surface Stokes drift velocity was further enhanced (0.57 m s−1) and farther eastward. The Eulerian current was still from the north-northeast but was slightly intensified (0.19 m s−1). The net surface heat flux had also intensified (778 W m−2). VVV had intensified, although the coherent structure was quite similar to that in burst 319. The orientation of the observed coherent structure was again similar overall to the orientation predicted from winds and Stokes drift.
Figure 7 shows another example of turbulence observed during 0200–0220 LT (0224–0244 LT for the HADCP) 11 December 2015 (burst 495), labeled L3 in Fig. 3. Though the UADCP and the HADCP were not operated simultaneously at this time, the small time lag allows us to neglect this effect. On this day, an atmospheric low pressure system traveling along the southern coast of Japan induced intense winds with large temporal variability (Fig. 3). The wind at this time was south-southeasterly, while the surface waves and the mean Eulerian current were from the southwest. The wind speed was high (16.7 m s−1) and the surface Stokes drift velocity was large (0.30 m s−1). The Eulerian mean current speed was 0.11 m s−1. The negative surface heat flux (surface heating; −323 W m−2) excludes the possibility of vertical convection at this time. Nevertheless, the Eulerian velocity was well mixed in the vertical. Turbulent velocities and echo intensity showed coherent structures (rolls) of turbulent flows advected northeastward by the mean Eulerian current. The orientation of the observed rolls was in good agreement with the predicted orientation of Langmuir cells.
Figure 8 shows the flows measured from 2200 to 2220 LT 26 November 2015 (burst 155; labeled S in Fig. 3). At this time, wind was 4.87 m s−1 in speed and north-northwesterly, while the Stokes drift was 0.24 m s−1 and westerly. The surface Eulerian current was 0.13 m s−1 in speed and was from the northeast. The net heat flux was 558 W m−2. The observed features were mostly similar to those in burst 319, except that the Eulerian velocity had large vertical shear throughout the water column. The orientation of the shear velocity vector was down–Eulerian flow rather than downwind, except very near the surface. These results suggest that the turbulence measured by the UADCP was shear-driven turbulence, probably originating from bottom stress. A slight difference in the Eulerian current direction between burst 155 (from the east-northeast at greater depths along which water depth changes rapidly) and 319 and 329 (from the north-northeast along which water depth changes smoothly) seems responsible for this large difference in the bottom-stress contribution to the observed turbulence. In fact, large vertical shear in the Eulerian velocity was often found when the Eulerian currents were from the east-northeast (not shown). Further quantitative investigations, with more specific measurements of local bottom topography and near-bottom flows as well as realistic numerical simulations, are necessary for more comprehensive understanding of the bottom-stress contribution and will be our future study. Note, however, that coherent structures of turbulent flows were also found in the HADCP-measured radial velocities and echo intensities at 3-m depth away from the tower where the water depth is large (≥20 m), and hence the effect of bottom stress is expected to be small. Because the observed orientation of the roll was again similar to the predicted orientation, we considered that turbulence similar to that in bursts 319, 329, and 495 occurred away from the tower in this burst (155). Thus, the shear-driven turbulence caused by bottom stress is considered to have dominated locally at the UADCP deployment location. Note that the bottom-stress-driven turbulence was only occasionally dominant (shown later in Fig. 14).
The last event shown in Fig. 9 was obtained from 1800 to 1820 LT 30 November 2015 (burst 247; labeled C in Fig. 3). At this time, wind speed (4.87 m s−1), waves (significant wave height was 0.166 m and the surface Stokes velocity was 0.015 m s−1), and the Eulerian current (0.066 m s−1 in speed) were all small, while the net heat flux was large and positive (865 W m−2). Turbulence was weak, but was captured by the UADCP (Fig. 2b). VVV was 1.2 × 10−4 m2 s−2 = 0.37W*, close to the scale of convective turbulence (0.3). The weak turbulence is thus considered buoyancy driven (i.e., convection).
c. Dependence of turbulent intensity on La and Ho
The previous subsection shows that the observed turbulence was driven by waves, winds, and surface buoyancy flux (and occasionally by the bottom stress). To quantify the relative contributions of these surface external forcings to the observed turbulence, that is, the dependence of the observed turbulence on these external forcings, vertical averages of VVV () and ε observed in 2014 and 2015 are plotted in La–Ho space (Fig. 10). Here VVV and ε in each burst were first averaged over the observed depths and then gridded within a box in La–Ho space with dimensions of . (Results described below are less sensitive to the box dimension.) Larger (smaller) La represents a larger (smaller) contribution of winds relative to waves, while larger (smaller) Ho means a larger (smaller) contribution of sea surface cooling relative to waves. Note that VVVs under sea surface cooling (positive buoyancy flux) were analyzed in this section because the number of VVV data under sea surface heating was too small as a result of short periods of sea surface heating. Note also that the bottom-stress-driven turbulence (as observed in burst 155) was excluded from this analysis by discarding the burst data in which the vertically averaged vertical shear of the Eulerian mean velocity below 2.3 m exceeded 0.015 s−1. (This threshold value was determined by visual inspection.)
Figure 10 shows two important features revealed in our field experiment. First, large VVV and ε are found at smaller La and Ho, when waves are expected to be dominant in driving turbulence. This is consistent with the fact that coherent structures similar to Langmuir cells were often observed in our field experiment. This result strongly suggests that the strong turbulence observed in our field experiment was Langmuir turbulence. Second, the distribution of the observed VVV in La–Ho space is along the positive diagonal; large (small) La is generally accompanied by large (small) Ho (see also Figs. 3d,j). This is because strong (weak) winds are frequently accompanied by strong (weak) sea surface cooling.
Previous observations (e.g., D’Asaro 2001; Tseng and D’Asaro 2004; Gargett and Grosch 2014) found that was proportional to , the scale of the shear-driven (wind driven) turbulence. Proportionality between the friction velocity and the Stokes drift velocity , that is, constant , was considered a possible explanation for this proportionality. This proportionality between and , as well as and , was also found in our field experiment (Figs. 11a,b), although scatter was larger than in the previous observations. Note, however, that La was not constant in our data. In fact, the proportionality factor between and (85) is much larger than that expected from the equilibrium wind sea (11–15) calculated from the Pierson–Moskowitz spectrum (e.g., Gargett and Grosch 2014). Larger values of than that expected from the equilibrium wind sea indicates that swells were dominant in our field. Correlations between and (0.13; 95% confidence interval being 0.03–0.23) or (0.31; 95% confidence interval being 0.21–0.40) were found to be small in our data (e.g., Fig. 11c), as D’Asaro (2001) reported.
To compare the observed dependence of VVV and ε on La and Ho with those simulated in the previous numerical studies (e.g., Li et al. 2005; Grant and Belcher 2009), they were normalized by and , respectively (Figs. 10c,d). The normalized VVV was larger for smaller La and larger Ho; it is never constant if plotted on La–Ho space. More noteworthy is that the observed dependence of the normalized VVV on La and Ho was in qualitative and quantitative agreement with the simulated dependence (e.g., Li et al. 2005). The isolines of the normalized VVVs in La–Ho space run in the positive diagonal direction over the observed ranges of La and Ho, and magnitudes of the normalized VVVs (0.95–16) were similar to those values (0.5–20) in the simulation of Li et al. (2005).
The reason for the apparent proportionality between and (Fig. 11a) in our data is thus explained primarily by the fact that La is generally large when Ho is large; La and Ho are distributed along the positive diagonal in La–Ho space, while VVV/ does not change significantly along the diagonal. The diagonal distribution is due to the high correlation between surface cooling (primarily latent heat flux) and wind speed. As La increases (the wave effect reduces), Ho increases (the convective effect increases). In other words, the wave effect is compensated for by the convective effect. Consequently, the normalized VVV does not appear to change significantly with La if the Ho dependence is not considered. This compensation effect may partly explain the apparent proportionality between VVV and found by D’Asaro (2001) and Tseng and D’Asaro (2004), who conducted field experiments in early winter (when correlation between La and Ho was expected) but did not focus on the compensation effect.
To see the dependence of VVVs on La more quantitatively, the normalized VVVs were plotted as a function of La (Fig. 12). It was found that when Ho < 0.1, that is, VVVs were less affected by the surface buoyancy flux, the power of the dependence on La was −2.0 (95% confidence interval was from −2.3 to −1.6), slightly larger in magnitude than the power (−4/3) obtained in the LES of Grant and Belcher (2009). The power of the observed dependence of the normalized VVVs at La < 0.3 was further large (−2.4). Therefore, other factors than considered in Grant and Belcher (2009), where a monochromatic wave aligned with wind was considered, influenced VVVs observed in our field experiment. Among several factors that can cause this slight difference, misalignment of wind and wave directions seemed one possible cause. Van Roekel et al. (2012) suggested that the Langmuir turbulence under the misaligned wind and wave conditions is adapted to that under the aligned condition if projected Langmuir number is used instead of La, where , , and are the direction of Lagrangian current, Stokes drift, and wind, respectively. The power of the observed dependence of the normalized VVVs on Lap at Ho < 0.1 was −1.2 (95% confidence interval was from −1.5 to −0.9; Fig. 12b), indicating some misalignment effects on the observed VVVs. However the power of the dependence on Lap was still large (−2.1) at Lap < 0.3. Thus, other factors likely affected the observed VVVs. Harcourt and D’Asaro (2008) and Kukulka and Harcourt (2017) showed that vertical decay scale D of the Stokes drift changes VVVs of the Langmuir turbulence if where H is the mixing layer depth in their simulations. We estimated D as the e-folding depth scale of the Stokes depth profile by the least squares fit from the surface to 8-m depth and calculated with H being the water depth (≤30 m) and found that was mostly larger than 0.1 (not shown). Thus, the vertical depth scale of the Stokes drift had negligibly small effects on the observed VVVs. On the other hand, Harcourt and D’Asaro (2008) also showed that VVVs of the Langmuir turbulence increase with wave age . The wave age, estimated as and averaged in La–Ho space (Fig. 10e), was larger where the normalized VVVs were larger. The range of the observed wave age (0.7–1.5) combined with the results of Harcourt and D’Asaro (2008) indicates that VVVs may have been magnified by a factor of 2–3 because larger wave age frequently occurred at smaller La and larger Ho.
d. Orientation of Langmuir cells
Figure 3 shows the time series of the observed orientation and predicted orientation of Langmuir cells, as well as wind and Stokes drift directions, obtained in the burst with La < 0.3 (i.e., Langmuir turbulence was expected to be dominant). In this figure, orientation was measured clockwise from the north and was defined from 45° to 225°. (No distinction was made between θ and θ + 180°.) Figure 13 shows the frequency distribution of and as well as wind and Stokes directions. Overall, was found between wind and wave directions. The largest frequency of and almost coincide with each other, indicating that orientation of Langmuir cells is overall along the orientation of the Lagrangian current shear. A wider spread of than was also found. This spread may be related partly to temporal change in winds and waves and partly to uncertainty of the present direction finding method using Eq. (3).
4. Discussion: TKE budget
To further investigate the source of the observed turbulence, the TKE budget was analyzed based on the TKE tendency equation:
where is the horizontal Eulerian velocity vector, is the Stokes drift velocity vector, b is buoyancy, and p is pressure. The first term of the right-hand side represents the production rate of TKE by sheared Eulerian mean flow (shear production), the second term is the production rate by wave–current interaction (Stokes production), the third term is the conversion rate from potential energy into TKE (buoyancy production), the fourth term is the vertical transport of TKE, and the last term is the TKE dissipation rate. In the present study, , , , and were burst-averaged quantity at each depth that were measured/estimated with the UADCP. The estimated Reynolds stress values are shown inFigs. 4 and 6–9. See the appendix for details of the estimation procedure. Thus, we evaluated the shear and Stokes productions and the dissipation of TKE in each burst at each depth. Vertical profiles of buoyancy production were not measured, but the vertically uniform water temperature allows us to assume , where B is the (measured) surface buoyancy flux and H is the water depth. The TKE transport term is another unmeasured quantity. This term represents downward transport of high TKE generated by waves near the surface (e.g., Grant and Belcher 2009) and is expected to be negative near the surface and positive at greater depths.
Note that the present TKE budget analysis is more qualitative than quantitative, because the Reynolds stress estimate has 10% error ( appendix). We believe, however, that this analysis still reveals important features of the observed Langmuir turbulence, as described below.
Figure 14 shows the vertical profiles of shear production, Stokes production, buoyancy production, and dissipation during our field experiment in 2014 and 2015. In most cases, shear production and Stokes production dominated over buoyancy production. Dissipation was large when the productions were large. To illustrate relative contribution of each production terms to the observed Langmuir turbulence in more detail, the production rates and the dissipation rate in each burst were first normalized by , the scale for the Langmuir turbulence (e.g., Grant and Belcher 2009), at the corresponding time and then averaged over the bursts in which Ho was less than 0.1 and the bottom-stress effects were small, that is, when pure Langmuir turbulence was expected (Fig. 15). A larger (smaller) total (shear plus Stokes plus buoyancy) production term than the dissipation term at shallower (greater) depths was qualitatively consistent with the unestimated transport term that should be negative (positive) at shallower (greater) depths. Relative contribution of the vertically integrated Stokes, shear, and buoyancy production rates was 59%, 36%, and 5%, respectively. The largest contribution by the Stokes production was consistent with the Langmuir turbulence. Magnitudes of the normalized production and dissipation terms were not largely different from those in LESs (e.g., Polton and Belcher 2007; Grant and Belcher 2009). Note, however, that the estimated shear production term was positive and nonnegligible over the measurement depth range. This feature was different from most of the previous LESs in which shear production was negligibly small. In our field experiment, a larger shear production rate than Stokes production rate was often associated with the Reynolds stress directed to a more upwind direction than up-Stokes direction as found in bursts 329 and 495 (Figs. 6, 7). This does not seem curious in view of the LES results of Van Roekel et al. (2012), who investigated Langmuir turbulence under wind misaligned with waves and found that Reynolds stress is in up-Stokes direction in the growing stage of the LCs but is in upwind direction in its later mature stage, while structure of the LCs is left unchanged. Thus, the present field experiment suggests nonnegligible contribution of shear production due to LCs under wind misaligned with waves to the ocean surface turbulence.
5. Concluding remarks
In the present study, turbulence in the ocean surface boundary layer was measured, and its response to external forcings (winds, waves, and buoyancy flux) was investigated through a field experiment conducted around the marine observation tower off the southern coast of Japan. Turbulent velocities were measured with upward-looking and horizontally looking ADCPs, while atmospheric fluxes were measured at the observation tower. Surface waves were measured with the upward-looking ADCP. Observed coherent turbulence structures (roll structure) and larger vertical velocity variances (VVV) at small La and Ho suggest that the observed turbulence was Langmuir turbulence. In contrast to previous observations (D’Asaro 2001; Tseng and D’Asaro 2004; Gargett and Grosch 2014), VVV was found to depend on both La and Ho, as reported in the previous numerical simulations (e.g., Li and Garrett 1995; Harcourt and D’Asaro 2008; Grant and Belcher 2009). We also found as in the previous observations, if dependence on Ho was not taken into account. In the present experiment, this feature is explained by or compensatory effects between surface waves and convection; as La increases (the wave effect reduces), Ho increases (the convective effect increases). Because turbulent heat flux is often dominant in surface heat fluxes, can be expected in other regions. Thus, we suggest that , in addition to , can contribute to apparent .
The TKE budget analysis suggests that shear production plays a nonnegligible role in maintaining the high TKE of the Langmuir turbulence. Misalignment of winds and waves seems to be a key factor in this nonnegligible contribution. However, the quality of the measurements was not high enough to evaluate these processes in further detail. The easiest way to increase the measurement quality is to increase the ping rates of the ADCP up to, say, 16 Hz, which was not possible at the time of our field experiment but is now possible for the recent ADCP with a fifth beam. In future studies, we will obtain higher-quality measurements that will contribute to a better understanding of sea surface mixing not only in coastal regions but also over the ocean, where the misalignment is expected to be nonnegligible (Li et al. 2016). For comprehensive understanding of the present results, effects of bottom topography, such as the bottom-stress-induced turbulence and the modulation of swell by the seamount and its impact on LCs, need to be investigated. These will be done in our future study.
We express our sincere thanks to two anonymous reviewers who gave us instructive comments to improve our manuscript. We also express our hearty thanks to Mr. Yoshikazu Yamamoto for his support in our field experiment. This work was supported by JSPS KAKENHI Grant numbers JP25287123 and JP15H05824. The data used in figures of this manuscript are available from authors upon request or online (https://fsv.iimc.kyoto-u.ac.jp/public/ZkYoAAMcdknAS3ABOIhh0RXgLPwyAgrRX06gCRGxpY0F).
UADCP Tilt Correction
Denoting () as velocity components in the UADCP’s Cartesian coordinate [positive () is in the direction from transducer 2 to 1 (from 4 to 3) and positive is from the tail to the head of the UADCP], the turbulent quantities are calculated from the UADCP’s radial velocities () as
where = 25° is the slanted beam angle. Since the UADCP was tilted from the vertical by ≤6° during our field experiment, this tilt needs to be corrected. With () being velocity components of the corrected ADCP coordinate (say, is vertical velocity), the corrected turbulent quantities are written as (e.g., Lu and Lueck 1999)
where and are the pitch and roll angles, respectively. In the above, was not measured and needs to be estimated.
In this study, we indirectly estimated using the radial velocities of the UADCP and the HADCP. First, the radial velocities of the HADCP (where r is distance from the HADCP and represents the direction of the ith beam) were used to estimate at 3-m depth. Using the relation
where is measured with respect to the direction of the upward-looking ADCP, the following expressions are obtained:
where spatial homogeneity was assumed in , , and . The above expressions were used to estimate as well as and from the HADCP’s radial velocities . Here values at 5 35 m were used. Validity of the spatial homogeneity assumption in , , and was expected to be relatively high in this short range. Data with large fish-motion-induced errors were excluded by visual inspection.
Figure A1a shows scatterplots between (red) and (blue) measured with HADCP and UADCP at 3-m depth. Though the scatter is large and the slope of the regression line is larger than unity, the overall correspondence between the magnitudes of the two ADCP’s velocity variances were found, indicating fairly good estimation of and , and probably .
Next, two neighboring radial velocities of the UADCP [e.g., and where r is distance from UADCP] were used to estimate the covariance between diagonal velocities (e.g., u and υ) at points separated by (where = 25° is the slanted beam angle of UADCP) as
where represents the correlation between u and υ measured at points separated by . Scatterplots between and and between and show that magnitudes of and were slightly smaller than magnitudes of and (Fig. A1b). However, the order of their magnitudes remains unchanged. Scatterplots between measured at 3-m depth by UADCP and measured by HADCP also show overall similarity in magnitude. These results indicate that and that and (Fig. A1b). As a result of , estimation errors in and were thus found to be less than 10%.
Finally, velocity components in the Earth coordinate system (where u, υ, and w are eastward, northward, and upward velocity) were calculated using the UADCP heading.
Current affiliation: Hydro Technology Institute Co., Ltd., Osaka, Japan.