Transition to turbulence in wind-drift layers

The evolution of initially-surface-concentrated dye was observed with a Laser-Induced Fluorescence (LIF) system. Images were acquired with a CCD camera (Jai TM4200CL, 2048 × 2048 pixels) equipped with an 85 mm Canon EF lens focused at the air-water interface. Illumination was provided by a thin 3 mm thick laser light sheet generated by a pulsed dual-head Nd-Yag laser (New Wave Research, 120 mJ/pulse, 3–5 ns pulse duration). The laser light illuminated a thin layer of fluorescent dye carefully applied to the water surface prior to each experiment. Observations were conducted with the light sheet in both along-wind and transverse directions. The LIF camera collected images at a 7.2Hz frame rate and with a field of view of 11.6 x 11.6 cm in the along-wind configuration, and 13.9 x 13.9 cm in the transverse direction. An edge detection algorithm based on local variations of image intensity gradients, computed by kernel convolution was used to identify the location of the surface in the LIF images (Buckley & Veron, 2017).

In addition to LIF, we employed Thermal Marking Velocimetry (TMV), as developed by Veron & Melville (2001) and Veron et al. (2008), to measure the surface velocity by tracking laser-generated Lagrangian heat markers in thermal imagery of water surfaces. In the present experiment, infrared images of the surface were captured by a 14-bit, $640\times 512$ quantum well infrared photodetector (QWIP — 8.0-9.2$\mu$m) FLIR SC6000 infrared camera operated at a 43.2 Hz frame rate, with an integration time of 10 ms, and a stated rms noise level below 35 mK. After image correction to account for the slightly off-vertical viewing angle of the imager, the resulting image sizes were $24.6\times 24.6$ cm.

The infrared imager is sensitive enough to detect minute, turbulent temperature variations in the surface thermal skin layer (Jessup et al., 1997; Zappa et al., 2001; Veron & Melville, 2001; Sutherland & Melville, 2013). It thus easily detects active weakly-heated markers, generated by a 60 W air-cooled $\mathrm{CO_{2}}$ laser (Synrad Firestar T60) equipped with an industrial marking head (Synrad FH Index) and two servo-controlled scanning mirrors programmed to lay down a pattern of 16 spots with 0.8 cm diameter and at a frequency of 1.8 Hz.

The spatially averaged surface velocity was estimated by tracking the geometric centroid of these Lagrangian heat markers for approximately one second. Both Gaussian interpolation (which has sub-pixel resolution due to the Gaussian pattern of the laser beam) and a standard cross-correlation technique yielded similar estimates for the surface velocity.

Figure 1 summarizes the experimental results. Figure 1a shows a time series of estimated wave steepness, and figure 1b shows the measured average surface velocity using TMV. The thick gray line in figure 1 plots $At$, where $A=1\,\mathrm{cm\,s^{-1}}$ showing that the surface current increases linearly in time. (The time axis for laboratory measurements is adjusted to meet this line, which constitutes a definition of “$t=0$”.) Following Veron & Melville (2001), figures 1a–b divide the development of the waves and currents into four stages:

1. 1.

   Viscous acceleration, $t=0$–$16$ s. In the first stage, viscous stress between the accelerating wind and water accelerates a shallow, laminar, viscous wind-drift layer.

2. 2.

   Wave-catalyzed “Langmuir” shear instability, $t=16$–$18$ s. At $t\approx 16$ s, detectable capillary ripples appear. A wave-catalyzed shear instability — which obey the same dynamics as “Langmuir circulation”, which often refers to much larger scale motions in the ocean surface boundary layer Craik & Leibovich (1976) — immediately starts to develop and grow in the wind-drift layer.

3. 3.

   Self-sharpening, $t=18$–$20$ s. When the shear instability reaches finite amplitude, nonlinear amplification due to perturbation self-advection sharpens the instability features into narrow jets and downwelling plumes.

4. 4.

   Langmuir turbulence, $t>20$ s. The self-sharpened circulations develop significant three-dimensional characteristics and transition to fully developed Langmuir turbulence.

As the wind starts to accelerate, viscous stress across the air-water interface drives a laminar wind-drift current. The thick gray line in figure 1 indicates that the average surface velocity nearly obeys

where $A\approx 1\,\mathrm{cm\,s^{-2}}$. Veron & Melville (2001) point out that the viscous stress consistent with linear surface current acceleration is

where $\boldsymbol{\tau}$ is the downwards kinematic stress across the air-water interface, $\boldsymbol{\hat{x}}$ is the along-wind direction ($\boldsymbol{\hat{y}}$ and $\boldsymbol{\hat{z}}$ are the cross-wind and vertical directions), and $\alpha\approx 1.2\times 10^{-5}\,\mathrm{m^{2}\,s^{-5/2}}$ produces $A=\alpha\sqrt{4/\pi\nu}\approx 1\,\mathrm{cm\,s^{-1}}$ given the kinematic viscosity of water, $\nu=1.05\times 10^{-6}\,\mathrm{m^{2}\,s^{-1}}$. The laminar viscous solution beneath (2) is (Veron & Melville, 2001)

Viscous acceleration continues until gravity-capillary ripples appear at the air-water interface when $t\approx\tilde{t}\equiv 16$ seconds and thus $U_{0}(t)\approx\tilde{U}_{0}\equiv 16\,\mathrm{cm\,s^{-1}}$.

As soon as ripples appear on the water surface, a second, slower, non-propagating instability begins to grow. Remarkably, this second instability is catalyzed by and therefore requires the presence of capillary ripples: for example Veron & Melville (2001) show that instability and turbulence are suppressed if ripple generation is suppressed by layering a surfactant on the water surface.

Interestingly, ripple inception is not explained by linear theory (Young & Wolfe, 2014). We neverthless assume that the ensuing dynamics are described by the weakly-nonlinear and wave-averaged “Craik-Leibovich” momentum equation. This ansatz is supported by the data in figure 1a, which plots an estimate of the wave steepness $\epsilon\equiv ak$, where $a$ is the surface wave amplitude (readily measured in the lab) and $k\approx 2\pi/3\,\mathrm{cm^{-1}}$–$2\pi/5\,\mathrm{cm^{-1}}$ is the gravity-capillary wavenumber for ripples with wavelengths between 3 and 5 cm (slightly larger than the minimum gravity-capillary wavelength 1.7 cm). Apparently, by the time ripples reach detectable amplitudes they have small slopes with $\epsilon<0.1$.

The wave-averaged Craik-Leibovich equation (Craik & Leibovich, 1976) formulated in terms of the Lagrangian-mean momentum $\boldsymbol{u}^{\mathrm{L}}$ of the wind-drift layer is

where $\boldsymbol{u}^{\mathrm{S}}$ is the Stokes drift of the field of capillary ripples and $p^{\mathrm{E}}$ is the Eulerian-mean pressure. The asymptotic derivation of the CL equation (4) requires $\epsilon\ll 1$. We require $\boldsymbol{u}^{\mathrm{L}}$ to be divergence-free (Vanneste & Young, 2022),

The Stokes drift associated with monochromatic capillary ripples propagating in the along-wind direction $\boldsymbol{\hat{x}}$ is

is the phase speed of gravity-capillary waves in deep water with gravitational acceleration $g=9.81\,\mathrm{m\,s^{-2}}$ and surface tension $\gamma=7.2\times 10^{-5}\,\mathrm{m^{3}\,s^{-2}}$. In all cases considered here, the Stokes drift (6) is minuscule compared to the mean current $\boldsymbol{u}^{\mathrm{L}}\sim U$.

We invoke the power method (Constantinou, 2015) to investigate the stability of the wind-drift layer immediately following ripple generation. In our implementation of the power method we extract the fastest growing mode perturbing the steady shear flow $\tilde{U}(z)$, such that

where $\tilde{U}(z)\equiv U(z,\tilde{t})$ represents the wind-drift profile “frozen” at $\tilde{t}=16$ s and $\boldsymbol{u}=(u,v,w)$ is the perturbation velocity. Inserting (7) into (4)–(5), introducing a streamfunction $\psi$ with the convention $(v,w)=(-\psi_{z},\psi_{y})$, and neglecting terms that depend only on the mean flow $U$ or $u^{S}$ yields the two-dimensional system

where $\triangle\psi=w_{y}-v_{z}$ is the $x$-component of the perturbation vorticity, $\mathrm{J}(a,b)=a_{y}b_{z}-a_{z}b_{y}$ is the Jacobian operator, and $\Omega=\tilde{U}_{z}-u^{\mathrm{S}}_{z}$ is the Eulerian-mean shear — or, as we prefer, the mean, total cross-wind vorticity $\bnabla\times\left(\tilde{U}\boldsymbol{\hat{x}}-\boldsymbol{u}^{\mathrm{S}}\right)=\Omega\boldsymbol{\hat{y}}$. The power method isolates the fastest growing linear modes of (8)–(9) by iterating over short integrations of the fully nonlinear CL equations (4)–(5) from $t=\tilde{t}$ to $t=\tilde{t}+\Delta t$ and repeatedly downscaling $\boldsymbol{u}$ via (7) to ensure each integration has essentially linear dynamics.

Figure 2 shows the structure of the most unstable mode for $\epsilon=0.1$. The mode is necessarily sinusoidal in the horizontal because (8)-(9) is horizontally-uniform and has alternating upwelling and downwelling regions, with a much larger streamwise component $u^{\prime}$ than cross-stream and vertical components $v^{\prime}$ and $w^{\prime}$. Figure 3 shows the results of a parameter sweep from $\epsilon=0.04$ to $\epsilon=0.3$, illustrating how increasing $\epsilon$ leads to faster growth rates and decreasing wavelengths of the fastest growing mode. We recover Veron & Melville (2001)’s experimental result that shear instability does not occur without surface waves.

Note that the kinetic energy source for growing perturbations is the mean shear $\tilde{U}(z)$ and there is no energy exchange between perturbations and the surface wave field within the context of the wave-averaged equations (4)–(5). To see this, consider that (4)–(5) conserves total kinetic energy $\int\tfrac{1}{2}|\boldsymbol{u}^{\mathrm{L}}|^{2}\,\mathrm{d}V$ when $\nu=0$ and $\partial_{t}\boldsymbol{u}^{\mathrm{S}}=0$ (viscous stress is negligible). This is why we characterize the shear instability as “wave-catalyzed”: while the presence of waves is necessary for instability, and while the instability growth rate is strongly affected by wave amplitude, the kinetic energy of the growing perturbation is derived solely from the mean shear.

When the wave-catalyzed shear instability reaches finite amplitude, it begins to self-sharpen, producing narrow along-wind jets and downwelling plumes. The sharpening — but still two-dimensional — plumes then transport a measureable amount of mean momentum downwards before becoming unstable to three-dimensional perturbations and thereby transitioning to fully-developed turbulence. This nonlinear sharpening and depletion of the average near-surface momentum occurs between 17.5–20 seconds, as evidenced by the red shaded region in figure 1b.

To simulate the second wave-catalyzed instability through finite amplitude and toward transition to turbulence, we propose a simplified model based on the wave-averaged equations (4)–(5) with two main components representing (i) the capillary ripples and (ii) the initial condition at $\tilde{t}=16$ seconds. We model the evolving capillary ripples as steady, monochromatic surface waves with wavenumber $k=2\pi/3\,\mathrm{cm^{-1}}$. We model the condition of the wave tank at $t=16$ seconds with

where $\mathbb{U}^{\prime}$ is an initial noise amplitude and $\boldsymbol{\Xi}$ is a vector whose components are normally-distributed random numbers with zero mean and unit variance.

Through experimentation, we find that the instability and transition to turbulence are only weakly sensitive to the wavenumber $k$. The tuning parameters of our model are therefore (i) the amplitude of the initial perturbation $\mathbb{U}^{\prime}$, and (ii) the wave steepness $\epsilon$. We discuss possible interpretations for $\mathbb{U}^{\prime}$ in section 4.

We also simulate the evolution of dye concentration $d$ via

with molecular diffusivity $\kappa=10^{-7}\,\mathrm{m^{2}\,s^{-1}}$, the smallest we can reasonably afford computationally. (The correspondance between $\theta$ and rhodamine is imperfect because the molecular diffusivity of rhodamine is $\kappa=10^{-9}\,\mathrm{m^{2}\,s^{-1}}$.) We initialize $\theta$ with a $\delta$-function at the surface. We integrate (4)–(5) and (11) given (6) using Oceananigans (see Ramadhan et al. 2020 and Wagner et al. 2021) with a second-order staggered finite volume method on a single A100 GPU in a $10\times 10\times 5$ cm domain in $(x,y,z)$ with 0.13 mm regular spacing in $x,y$ and variable spacing in $z$ with $\min(\Delta z)\approx 0.26\,\mathrm{mm}$, corresponding to $768\times 768\times 512$ finite volume cells.

The results are visualized in figure 4. Figure 4a compares the average surface velocity diagnosed from the simulation with the laboratory measurements presented in figure 1b. Figure 4b plots the maximum absolute vertical velocity, $\max|w^{L}|$. Figures 4c–l compare the simulated dye concentration on $(y,z)$ slices with LIF measurements of rhodamine from the laboratory experiment, showing the particularly good qualitative agreement between simulated and measured dye features. Visualizations are shown at $t=18.1$, 19.3, 20.0, 20.7, and 21.7 seconds.At $t=18.1$ s (figures 4c and d), the sharpened plumes have only just started advect appreciable amounts of dye. At $t=19.3$ seconds (figures 4 e and f), the plumes are beginning to roll up into two-dimensional mushroom-like structures. (Note the small, unexplained discrepancy between simulated and measured average surface velocities around $18.1<t<19.3$.) At $t=20.0$ seconds and $t=20.7$ (figures 4g, h, i, and j) three-dimensionalization and transition to turbulence are underway. At $t=21.7$ both the simulated and measured dye concentrations appear to be mixed by three-dimensional turbulence.

Figure 5 illustrates the sensitivity of the wave-averaged model to amplitude of the specified surface ripples and to the amplitude of the random initial perturbation. Figures 5a and b plot the surface-averaged along-wind velocity $u$ and maximum vertical velocity $w$ for two wave amplitudes $\epsilon=0.1$ and $\epsilon=0.11$, and for two initial perturbaiton ampltidues $\mathbb{U}^{\prime}=5\,\mathrm{cm\,s^{-1}}$ and $\mathbb{U}^{\prime}=10\,\mathrm{cm\,s^{-1}}$. The dependence of the maximum vertical velocity is the most evocative: doubling the initial perturbation shortens the self-sharpening phase (in which the maximum vertical velocity in figure 5b flattens before increasing sharply during the transition to turbulence) by a factor of five. Of the three cases plotted in figure 5, only $\epsilon=0.11$ and $\mathbb{U}^{\prime}=5\,\mathrm{cm\,s^{-1}}$ yield the satisfying agreement depicted in figure 4.

Following three-dimensionalization, momentum and dye are rapidly mixed to depth. Figure 6 visualizes (a) the $x$-momentum and (b) dye concentration at $t=23.4$ seconds, showing how the flow is organized into narrow along-wind streaks and broader downwelling regions — classic characteristics of Langmuir turbulence (Sullivan & McWilliams, 2010).

In figure 7, we compare the rate which dye is mixed to depth in the simulations versus measured by LIF in the lab. Figure 7a shows the simulated horizontally-averaged tracer concentration in the depth-time $(z,t)$ plane, while figure 7b shows a corresponding laboratory measurements extracted from LIF measurements in the $(y,z)$-plane. In figures 7a and b, a light blue line shows the height $z_{99}(t)$ defined as the level above which 99% of the simulated tracer concentration resides,

Using $z_{99}(t)$ to compare the tracer mixing rates exhibited in figures 7a and b, we conclude that the simulations provided a qualitatively accurate prediction of dye mixing rates. If anything, the simulation overpredicts the dye mixing rate — but the data probably does not warrant more than broad qualitative conclusions.