Langmuir Turbulence Controls on Observed Diurnal Warm Layer Depths

The turbulent ocean surface boundary layer (OSBL) shoals during daytime solar surface heating, developing a diurnal warm layer (DWL). The DWL significantly influences OSBL dynamics by trapping momentum and heat in a shallow near‐surface layer. Therefore, DWL depth is critical for understanding OSBL transport and ocean‐atmosphere coupling. A great challenge for determining DWL depth is considering wave‐driven Langmuir turbulence (LT), which increases vertical transport. This study investigates observations with moderate wind speeds (4–7 m/s at 10 m height) and swell waves for which breaking wave effects are less pronounced. By employing turbulence‐resolving large eddy simulation experiments that cover observed wind, wave, and heating conditions based on the wave‐averaged Craik‐Lebovich equation, we develop a DWL depth scaling unifying previous approaches. This scaling closely agrees with observed DWL depths from a year‐long mooring deployment in the subtropical North Atlantic, demonstrating the critical role of LT in determining DWL depth and OSBL dynamics.

• Based on observations and large eddy simulation (LES), we examine diurnal warm layer (DWL) depths for moderate wind conditions with swell • Based on LES results, we develop a unified scaling for the minimum DWL depth considering shear and wave-driven Langmuir turbulence (LT) • Observed DWL depths are consistent with LES results that include LT indicating the importance of wave-driven turbulence on DWL dynamics 2 of 11 the OSBL turbulence (Webster et al., 1996;Weller & Anderson, 1996). By using novel instruments, recent studies provide detailed observations of turbulence characteristics during the DWL developments (Bogdanoff, 2017;Hughes et al., 2020aHughes et al., , 2021Moulin et al., 2018;Sutherland et al., 2016). Decaying near-surface turbulence dissipation rate ϵ is observed 1-2 hr after heating onset, indicating weakening turbulence due to OSBL stratification (Moulin et al., 2018). Because the trapped momentum in the DWL greatly enhances shear instabilities, the decayed ϵ is often followed by a rapid growth during the later heating period (Moulin et al., 2018), even exceeding the initial ϵ close to the ocean surface (Bogdanoff, 2017;Sutherland et al., 2016). Measured turbulent heat fluxes indicate that turbulent mixing transports surface heat downwards (Hughes et al., 2020b), a process that depends significantly on wind speed (Hughes et al., 2021). Despite those recent insightful observations revealing DWL dynamics, wave effects are rarely measured and usually not included in previous observational analyses.
Wind-driven surface gravity waves influence OSBL mixing through wave breaking and LT (D'Asaro, 2014). Wave breaking injects turbulent kinetic energy (TKE) into the ocean surface layer, which quickly decays within one significant wave height (Terray et al., 1996). LT is caused by the interaction between surface gravity waves and wind-driven shear through the Craik-Leibovich instability mechanism (Craik & Leibovich, 1976), considered as a critical component for OSBL mixing (Belcher et al., 2012;Thorpe, 2004). The wave-induced Stokes drift tilts the Eulerian current vorticity into the wave propagation direction, forming coherent and alternating LT cells (Craik & Leibovich, 1976). This LT cell-structure leads to strong vertical mixing, deepening the OSBL and homogenizing upper ocean properties (Kukulka et al., 2009;McWilliams et al., 1997;Thorpe, 2004). A turbulence-resolving large eddy simulation (LES) method based on the filtered Craik-Leibovich equation is commonly used to understand LT dynamics (McWilliams et al., 1997;Skyllingstad & Denbo, 1995), which has successfully captured many important LT features observed in earlier studies such as enhanced vertical velocities (e.g., Gargett et al., 2004;Weller et al., 1985) and organized surface convergent regions (e.g., Farmer & Li, 1995;Plueddemann et al., 1996;Smith, 1992). Relatively recent LES studies comprehensively diagnose LT dynamics under heating conditions and indicate a weakly stratified DWL in the presence of LT (Min & Noh, 2004;Noh et al., 2009;Pearson et al., 2015). Kukulka et al. (2009) illustrated the important role of LT in stratified OSBL deepening by comparing LES results and observations. A LES study driven by observed wind and wave forcing reveals that swell waves may drive strong LT to inhibit the DWL formation, violating the traditional Monin-Obukhov similarity theory (Kukulka et al., 2013). In contrast, strong heating also prevents the formation of LT (Min & Noh, 2004). Pearson et al. (2015) scaled DWL depths with LT based on a suite of LES results, suggesting the crucial effects of LT on DWL dynamics. Nevertheless, a systematic comparison of LES and observed DWL characteristics has not been conducted, so observational evidence of LT effects on DWL dynamics is still required.
This study compares high-resolution LES results of DWL depth to year-long observations in moderate wind conditions with swell waves, revealing significant effects of LT on DWL deepening and demonstrating the dependence of DWL depth on LT. Our results provide rare observational evidence of LT effects under heating conditions, which is important for understanding OSBL mixing and improving ocean-atmosphere models.

Mooring Observation
This study utilizes mooring data from the first Salinity Processes in the Upper Ocean Regional Study (SPURS-I) field campaign (Farrar et al., 2015;Lindstrom et al., 2015) to understand the evolution of the DWL. The mooring was deployed from September 2012 to September 2013 in the subtropical North Atlantic (around 24.5°N, 38.0°W), which includes the measurements of (a) meteorological data every minute, (b) frequency-directional wave height spectra with maximum frequency of F max = 0.485 Hz every hour, and (c) temperature and salinity profiles every 5 min. The vertical resolutions of the ocean sensors range from 1 to 5 m in the upper 50 m with finest resolution near the ocean surface. Using the TOGA COARE 3.0b algorithm (Fairall et al., 2003), the bulk air-sea fluxes are calculated based on the measurements. The Stokes drift vector for LT diagnosis is derived by integrating radially over the observed wave height spectra and subsequently projecting into the horizontal coordinates (Kenyon, 1969). The contribution from unresolved high frequency (>F max ) waves is included by adding an analytical tail described in Zheng et al. (2021).
To focus on observed surface heating effects, a subset of SPURS-I observations is selected (170 diurnal cycles), which satisfies the following three criteria: (a) changes in vertically integrated heat content are approximately 3 of 11 due to surface heat fluxes, so that we focus on events with small lateral heat advection, (b) winds are relatively constant and net surface heating occurs for at least 6 hr, and (c) the observed DWL depths must exceed 4.4 m so that there are at least four density profile measurements within the DWL. For the selected subsamples, solar radiation accounts typically for 80% of the changes in OSBL heat content, and the standard deviation of wind speed is about 0.5 m/s, so that wind speeds are relatively constant. The selected subsamples are dominated by moderate winds (U 10 = 4 -7 m/s) and moderate LT activities (La t = 0.4 -0.6, Figure S1 in Supporting Information S1).
Here, = √ * ∕ (0) is the turbulent Langmuir number, which is inversely related to LT intensity (McWilliams et al., 1997;M. Li et al., 2005), where u * is the frictional velocity and u s (0) is the magnitude of surface Stokes drift. La t = 0.3 usually represents strong LT and La t > 0.7 indicates weak LT or shear-driven turbulence (ST). For the selected samples, wind and surface Stokes drift directions are roughly aligned with an average direction difference of 30°. The mean observed wave length of wind waves is about 40 m, and these wind waves primarily contribute to the Stokes drift driving LT. The surface heat flux usually exceeds 500 W/m 2 around local noon and sometimes is up to 800 W/m 2 , able to form obvious DWLs.
This study focuses on wave-driven LT because of the relatively moderate wind conditions with substantial swell waves. In our study, the wave age (C p /U 10 = 1.1 -5.5, where C p is the peak wave phase speed based on the observed peak wave frequency) substantially exceeds those of fully developed equilibrium wind seas, indicating relatively weak breaking wave effects, unlike in previous LES studies (Noh et al., 2009). The depth of breaking wave effects and associated enhanced TKE levels may scale with the significant wave height determined from the part of the wave spectrum that is actively driven by winds (Melville, 1996;Terray et al., 1996). An estimate for such a significant wave height can be obtained for fully developed wind equilibrium seas (e.g., Donelan et al., 1985;Pierson & Moskowitz, 1964) and ranges in our study from 0.1 to 2.1 m with an average of 0.7 m, which is much smaller than the observed DWL depths. Therefore, as a first step, we consider only LES results for shear and wave-driven LT for the comparison with observations.

LES Experiments
A turbulence-resolving LES model based on the Craik-Leibovich equation is employed to simulate the DWL responses (McWilliams et al., 1997). The LES method decomposes the variables (e.g., current velocity, temperature, and pressure) into the resolved components and the subgrid scale (Moeng, 1984;Sullivan et al., 1996). Several LES experiments are designed with constant and time-varying surface heat fluxes Q = I 0 + Q 0 , where I 0 and Q 0 represent the penetrative and non-penetrative heat fluxes, respectively.
To investigate the influences from wind, LT, heating, and Coriolis force, the constant heating experiments are driven by different combinations of surface heat flux (I 0 = 10-1000 W/m −2 and Q 0 = 0), friction velocity (u * = 0.0061, 0.0087 m/s), Coriolis parameter (f = 0.5 × 10 −4 , 1 × 10 −4 , 1.4 × 10 −4 1/s), and surface Stokes drift (u s (0) = 0 -0.096 m/s). The corresponding Monin-Obukhov length = 3 * ∕ 0 and turbulent Langmuir number La t range from L = 0.415 -147 m and La t = 0.3 to ∞, respectively. Note that most experiments are driven by constant wind forcing (u * = 0.0061 m/s) so that La t change is mainly controlled by wave states. The surface buoyancy flux is given by B 0 = αgI 0 /ρ 0 c p , where α is the thermal expansion coefficient of seawater; c p is the heat capacity; ρ 0 = 1024 kg/m 3 is the typical seawater density; and g = 9.81 m/s 2 is the gravitational acceleration. Each simulation first spins up for 12 hr without heating to develop turbulence, and then the heat flux is imposed for another 14 hr resulting in quasi-stationary turbulent fields.
For the diurnal heating conditions, two LES simulations with (La t = 0.3) and without LT are driven by the same friction velocity (u * = 0.0061 m/s) and diurnal heat flux cycle. During the diurnal heat flux cycle, the surface heat flux gradually increases to the heating peak of Q = 500 W/m 2 from the initial cooling of Q = −200 W/m 2 within 8 hr, followed by a decrease in surface heating to the initial cooling after another 8 hr (Figure 2a). Note that Q = Q 0 when Q < 0 and Q = I 0 when Q > 0. All transient simulations first spin up for 12 hr with Q = −200 W/m 2 before the diurnal cycle is applied.
Following Kukulka et al. (2013), the penetrative effect of solar short-wave radiation is described by, where z is the vertical Cartesian coordinate positive upward and z = 0 is at the ocean surface; I(z) is the depth-dependent net heat flux; and the constant coefficients are R = 0.62, ζ 1 = 1.5 m, and ζ 2 = 20 m (Paulson & 4 of 11 Simpson, 1977). For cooling conditions (Q < 0), the buoyancy flux is directly applied at the ocean surface. For all LES experiments, the initial temperature is set constant within an initial mixed layer depth H i = 50 m, below which the ocean has a constant stratification with a temperature gradient of 0.04°C/m. The Stokes drift of this study is described by monochromatic waves with a constant wave length λ = 60 m but different wave amplitudes (McWilliams et al., 1997).
To resolve small-scale turbulence due to stratification (Wang, 2001), the LES experiments deploy a high-resolution domain, with the along-wind, crosswind, and vertical directions having, respectively, domain sizes of [150,150,90] m and [500, 500, 720] evenly-spaced grid points. Diagnostic analysis indicates that the selected resolution resolves over 85% TKE and is smaller than the Ozimidov scale so that turbulence is well resolved. Furthermore, additional sensitivity tests indicate that a larger domain or a higher resolution has minor effects on the presented LES results.

Approaches to DWL Depth Estimation
Following previous LES work (e.g., Zilitinkevich et al., 2007), we define the OSBL depth H b as the depth with active turbulent transport so that H b is where the magnitude of the Reynolds stress decreases to 5% of its surface value. During maximum heating, H b is close to the depth of maximum temperature gradient. However, it is challenging to measure Reynolds stress in the ocean. In order to compare LES results with observations, this study also explores alternative estimates of DWL depths that can be readily obtained from density profile observations.
Traditionally, mixed layer depth H m is defined by the depth at which the density increases by Δρ from the surface value and is often applied to observations (de Boyer Montégut et al., 2004). Because heat flux substantially varies over a diurnal cycle, this study uses a relatively small density threshold of Δρ = 0.01 kg/m 3 to enhance the sensitivity of DWL detection. Note that the change of density during diurnal heating is dominated by temperature changes and that the salinity changes hardly affect the DWL evolution. However, our results indicate that H m does generally not well represent H b because it is sensitive to the shape of density profiles (discussed in Section 3.1). Hence, we use an alternative much more robust DWL depth estimate based on the potential energy anomaly ϕ (Reichl et al., 2022), denoted as H ϕ . By definition, ϕ is where ρ m is the average density between the surface and depth H ϕ . In this study, we use ϕ = 1.5 J/m 2 to be consistent with our LES H b as discussed in the following section.

Scaling of DWL Depths
To focus on the effects of waves and heating on DWL depths, this study investigates non-dimensional DWL depths normalized by Ekman layer depth scale H e = u * /f. Our LES results (Figure 1a) reveal that the normalized turbulent boundary layer depth H b /H e decreases with normalized surface buoyancy flux expressed by ̂0 = 0∕ ( 2 * ) = ∕ because the DWL is shallower under strong heating and weak winds. This decrease of DWL depth with increasing buoyancy flux is consistent with previous studies (e.g., M. Li et al., 2005;Pearson et al., 2015;Price et al., 1986). The enhanced vertical mixing caused by LT yields larger H b /H e for fixed ̂0 . For sufficiently large ̂0 , H b approaches the ST result possibly because strong heating leads to a very shallow DWL and suppresses LT, consistent with Min and Noh (2004).
We extend previous work (e.g., Noh & Choi, 2018;Pearson et al., 2015;Zilitinkevich & Baklanov, 2002;Zilitinkevich et al., 2007) to develop a scaling of H b for the OSBL forced by constant heat fluxes, wind stresses, and Stokes drifts (detailed in Appendix A). Our scaling is unified for ST and LT, in the form of For F LT = 1 (i.e., no wave forcing), the scaling is consistent with the scaling for the stratified atmospheric boundary layer (Zilitinkevich et al., 2007), which is also validated for the OSBL (Yoshikawa, 2015). Figure 1a reveals the quantitative agreements between scaled H b /H e from Equation 3 and LES results with a coefficient of 5 of 11 determination r 2 = 0.98, indicating the success of our scaling in capturing the effects of wind, wave, heating, and Coriolis force on DWL depths. Our scaling also captures the weaker LT effects with increased ̂0 and only displays slight mismatches to LES results for relatively large ̂0 . Note that Equation 3 only applies for the fully turbulent DWL (and this excludes very weak winds) and captures shear and wave-driven LT. Both LES results and scaling indicate that moderate LT (La t = 0.4-0.52), close to SPURS-I conditions, still leads to a substantial LT influence on DWL depth (Figure 1a).
These LES results can only be compared to measured DWL depths, H m and H ϕ , based on density observations, since H b requires direct estimates of subsurface turbulent stress, which remains an observational challenge. The LES comparison between different DWL depth estimates reveals that H ϕ , the mixed layer depth estimated from the potential energy anomaly (Equation 3), much more closely approximates the turbulent boundary layer depths H b than more traditional density difference mixed layer depths H m (Figures 1b and 1c). H m usually overestimates H b for small ̂0 and underestimates H b for large ̂0 (Figure 1b), because H m is very sensitive to the shape of the temperature profiles (Figure 1d). The temperature profiles are also significantly influenced by LT (Figure 1d), leading to the large bias of H m for representing H b . This bias suggests that the H m threshold criterion for estimating H b is not universally applicable. Diurnal heating significantly changes the stratification and the density profile shape of the OSBL within a few hours, making H m a poor estimate for H b .
By contrast, H ϕ is an estimate based on the vertical integration of the density profiles so that it is less sensitive to the details of the profile shape. Such a robust estimate of the OSBL depth is critical for OSBL depths estimated from coarser density profile observations. Moreover, H ϕ connects the estimate of DWL depths to the physical process of OSBL mixing because it is determined by the potential energy distance from its well-mixed condition (Reichl et al., 2022). As expected, H ϕ /H e (r 2 = 0.96, Figure 1c) is much closer to the scaled H b /H e than H m /H e (r 2 = 0.81, Figure 1b). This comparison illustrates the applicability of H ϕ to estimate observed DWL depths, which is critical for a direct comparison between LES results and observations.

Application of Scaled H b to Diurnal Heating
As our LES scaling assumes constant forcing, we first assess the applicability of scaled turbulent boundary layer depth H b for transient diurnal heating before comparing with observations. Based on LES results driven by 7 of 11 transient diurnal heating, mismatches between the stationary scaling and transient LES H b appear at the beginning and the end of the heating cycle when heating is relatively weak, because the scaling is applied instantly to the transient forcing and does not consider that mean and turbulent fields take time to adjust to the forcing. During the beginning of the heating cycle, H b first sustains its depth by the substantial turbulence below scaled H b and then rapidly shoals when turbulence sufficiently decays. Conversely, during the end of the heating cycle, surface forcing forms strong stratification, decelerating the DWL deepening rate and leading to shallower H b than the scaling would suggest. In contrast, around the heating peak, the relatively strong and constant heat flux results in a roughly stationary DWL, which is thus well captured by scaled H b for cases with and without LT (Figure 2). Meantime, the match near the heating peak also implies that the scaled H b serves as an effective minimum of the DWL depths during a diurnal cycle.
Consistent with the earlier discussion, potential energy anomaly-determined H ϕ is close to H b during most periods of the diurnal cycle for both ST and LT (Figure 2). These results suggest that H ϕ is an accurate estimate for H b and can be directly compared to scaled H b during peak heating. Meanwhile, the more traditional density difference-determined H m can substantially differ from H b during a diurnal cycle because the shape of the density profile evolves over time.

Observational Evidence of LT Effects
Observed temperature profiles are frequently characterized by the development of DWLs during daytime heating (Figures 3a and 3b). For each diurnal cycle, the observed H ϕ evolves with substantial shoaling only 2-4 hr after the onset of heating, similar to LES results, indicating that LES successfully reproduces the observed DWL evolution (Figures 3a and 3b). As expected from LES, scaled H b is also closest to observed H ϕ around the heating peak. LT strength is expected to be moderate in Figures 3a and 3b with observed La t = 0.4-0.55. Observed H ϕ is more consistent with the scaling incorporating LT around the heating peak and, at times, is substantially deeper than the scaling without LT.
As introduced in Section 2.1, we choose a subset of SPURS-I observations with small lateral heat advection and relatively stationary wind for a direct comparison with LES scaling. Illustrated by Figures 3c and 3d, the comparison of scaled H b and observed H ϕ for the SPURS-I subsamples further demonstrates that LT yields a deeper DWL. The DWL depths are compared within a specific time window, 2:00 UTC to 5:00 UTC, for each subsample to reduce transient effects. The time window is estimated from the transient LES results, roughly 1 hr before and 2 hr after the heating peak. As shown in Figures 3c and 3d, most observed H ϕ falls in the range of 4-10 m, indicating a shallow DWL and strong stratification. The scaled H b based on LES results with LT agrees significantly better with observations, consistent with the expectation for observations in moderate LT conditions. For H b > 10 m, the DWL depths are relatively scattered but still show a consistent trend between bin-average and regression fit. The regression fit indicates that DWL depths are underestimated by about 20% without LT, implying the important role of LT on the DWL evolution. Note that we did not tune any parameters to observations in this study but that the scaled H b is derived from LES results.
It can be challenging to disentangle shear, buoyancy, and wave forcing, since they are not always independent from one another. Developing waves often result in a tight relationship with winds (see Figure S1 in Supporting Information S1). To better differentiate between the implicit correlation in forcing conditions, we delimit a narrow range of wind speeds (4.5 < U 10 < 5.5 m/s, see box in Figure S1 in Supporting Information S1) and heat fluxes (600 < I 0 < 750 W/m 2 ), corresponding to 110 <̂0 < 170 , from the observation used in Figures 3c and 3d.
In this way, we expect the variations in non-dimensional DWL depth, H ϕ /H e , to be mainly due to the background wave conditions because wind and buoyancy forcing are chosen to be nearly constant. As shown in Figure 4, larger wave forcing (La t less than ∼0.6) results in 20%-40% deeper non-dimensional DWL depths, H ϕ /H e , under nearly constant wind and heating. Because ̂0 and La t are not significantly correlated, we attribute the variations in DWL depths to variable LT, rather than wind stress or heating conditions (Figure 4). The magnitude of H ϕ /H e and the dependence of H ϕ /H e on La t are also remarkably close to LES scaling. Overall, the comparison between LES results and observations demonstrates that DWL is substantially influenced by LT.

Conclusion
In  9 of 11 evidence of LT effects on the DWL, providing further evidence for the necessity of incorporating LT into upper ocean mixing process and ocean-atmosphere prediction models.

Appendix A: Stationary Scaling of H b Under Heating
Following the studies of atmosphere boundary layer depth estimate (Zilitinkevich & Baklanov, 2002;Zilitinkevich et al., 2007), the depth of the boundary layer is estimated through a weighted power mean of H b(N) for neutral conditions and H b(S) for the stratified conditions, where n = 2 is empirical factor. For shear-driven boundary layer, Zilitinkevich et al. (2007) suggested H b(N) ∼ H e and ( ) ∼ 2 * ∕( 0) 1∕2 ∼ ( ) 1∕2 . To take account LT effects, we modify u * with enhancement factor F LT . By comparing and extending previous studies (McWilliams & Sullivan, 2000;Noh & Choi, 2018;Pearson et al., 2015;Reichl et al., 2016), we find an empirical form of = ( 1 + −2 ) scales well with our LES results, where C L is empirical factor. The scaled H b is then expressed as, where C N = 0.93, C S = 0.825, and C L = 0.025, which are determined through a nonlinear regression based on the LES H b . In addition to LT, this scaled DWL depth also includes effects due to ST, Coriolis force, and surface heating. Note that our scaling does not explicitly include the effects of Stokes drift depth decay scale which should be included in a more complete parameterization (Harcourt & D'Asaro, 2008;Kukulka & Harcourt, 2017;Q. Li & Fox-Kemper, 2017).

Data Availability Statement
The SPURS-I mooring data can be accessed through the Woods Hole Oceanographic Institution's Upper Ocean Processes Group website (https://uop.whoi.edu/projects/SPURS/spurs1data.html). The numerical results used in this study are simulated by the National Center for Atmospheric Research Large Eddy Simulation model initially established by Moeng (1984) and operated with the numeric scheme of Sullivan et al. (1996).