2. Defining the surface layer in FST

The defining property of the surface layer is the heterogeneous mixture of two fluids, here air and water. At any instant, this mixture can be described by a binary colour function field:

(2.1) \begin{equation} c(\boldsymbol{x},t) \equiv \begin{cases} 0 & \text{if}\quad \boldsymbol{x}\in \text{air} \\ 1 & \text{if}\quad \boldsymbol{x}\in \text{water}. \end{cases} \end{equation}

Our interest is the average behaviour of the surface layer, and we consider the simplest case where FST is horizontally and temporally homogeneous. We use an average along the homogeneous dimensions,

(2.2) \begin{equation} {\overline {\boldsymbol{\cdot }}} \equiv \frac {\iiint \boldsymbol{\cdot }\,\mathrm{d} x \,\mathrm{d} y\, \mathrm{d} t}{\iiint \mathrm{d} x\, \mathrm{d} y \,\mathrm{d} t}, \end{equation}

to define the intermittency factor

(2.3) \begin{equation} \gamma (z) \equiv {\overline {c(\boldsymbol{x},t)}}, \end{equation}

which here is only a function of depth $z$ . This intermittency factor gives the average portion of volume occupied by water, with $\gamma =1$ corresponding to all water and $\gamma =0$ to all air (Brocchini & Peregrine Reference Brocchini and Peregrine2001b ). The opposite, $1-\gamma$ , is the void fraction.

Because the free surface is highly distorted, we first need a definition of the mean free surface. While alternative definitions may be more appropriate for wave-driven flows (Brocchini & Peregrine Reference Brocchini and Peregrine1996, Reference Brocchini and Peregrine2001b ), for FST we find a simple definition is sufficient. We define the mean free-surface height $\bar {\eta }$ as the height where air and water are evenly mixed:

(2.4) \begin{equation} \gamma (z=\bar {\eta }) = 0.5. \end{equation}

We note that, because FST entrains air bubbles deep beneath the free surface, $\bar {\eta }$ will be higher than the location of the quiescent free surface. In this way, $\bar {\eta }$ is a measure of bubble entrainment, which is the focus of Gaylo & Yue (Reference Gaylo and Yue2025). Our interest here is the turbulence near the free surface rather than bubbles deep beneath, so we focus on the depth relative to $\bar {\eta }$ , $z-\bar {\eta }$ .

Because the defining characteristic of the surface layer is mixing of air and water, it is common to choose $\gamma _{\textit{max}}$ and $\gamma _{\textit{min}}$ , and define the surface layer as the depths $z$ where $\gamma _{\textit{max}}\gt \gamma (z)\gt \gamma _{\textit{min}}$ (Brocchini & Peregrine Reference Brocchini and Peregrine2001b ). Hendrickson & Yue (Reference Hendrickson and Yue2019) express this in terms of density $\rho _{\textit{max}}\gt {\overline {\rho }}(z)\gt \rho _{\textit{min}}$ , which is equivalent. The choice of these limits is somewhat arbitrary. To justify a choice, Brocchini & Peregrine (Reference Brocchini and Peregrine2001b ) start by considering a single-valued free surface where instantaneous free-surface height $\eta (x,y,t)$ follows a Gaussian distribution with mean $\bar {\eta }$ and variance $\sigma ^2$ . It can be shown that this gives an intermittency factor

(2.5) \begin{equation} \gamma _{\textit{Gaussian}}(z) = \frac {1}{2} \left [1 - {\textrm {erf}}\left (\frac {z-\bar {\eta }}{\sqrt {2} \sigma }\right )\right ]\!. \end{equation}

Using $z-\bar {\eta }\in [-3 \sigma , 3 \sigma ]$ to define the surface layer gives $\gamma _{\textit{min}}\approx 0.001$ and $\gamma _{\textit{max}}\approx 0.999$ . In § 4 we show that FST produces a non-negligible volume of bubbles deep beneath the free surface, causing a non-Gaussian distribution, where $\gamma$ approaches $1$ very slowly. Due to these bubbles, we find $\gamma _{\textit{max}}\approx 0.999$ overestimates the depth of the surface layer. More generally, a challenge with using $\gamma _{\textit{max}}$ and $\gamma _{\textit{min}}$ is that the resulting definition of the surface layer depends on the tail behaviour of $\gamma$ , which may be only indirectly linked to the behaviour close to the mean free surface.

As a more robust definition of the surface layer, we use the behaviour of $\gamma$ at the mean free surface rather than the tail behaviour far from the mean free surface. To do this, we define a free-surface thickness in terms of the vertical derivative of $\gamma$ evaluated at the mean free surface $\bar {\eta }$ :

(2.6) \begin{equation} \delta _s \equiv C_{\mathrm{d}\gamma } \left (\left .\frac {\mathrm{d} \gamma }{\mathrm{d} z}\right |_{z=\bar {\eta }}\right )^{-1} \!. \end{equation}

Asserting that the surface layer is centred about $\bar {\eta }$ , depths $z-\bar {\eta }\in [- \delta _s/2, \delta _s/2]$ define the surface layer. For characterising depth relative to the surface layer, we define a non-dimensionalised depth,

(2.7) \begin{equation} z^* \equiv \frac {z-\bar {\eta }}{\delta _s}, \end{equation}

in which case the surface layer is non-dimensional depths $z^*\in [-0.5, 0.5]$ .

What remains in (2.6) is to determine the constant $C_{\mathrm{d}\gamma }$ such that $\delta _s$ does a good job in characterising the surface layer. To start, we choose the constant $C_{\mathrm{d}\gamma } = 6/\sqrt {2\pi }$ so that in the case where $\gamma =\gamma _{\textit{Gaussian}}$ , $\delta _s=6 \sigma$ and our definition of the surface layer is equivalent to $z-\bar {\eta }\in [-3 \sigma , 3 \sigma ]$ from Brocchini & Peregrine (Reference Brocchini and Peregrine2001b ). While initially somewhat arbitrary, it turns out that $C_{\mathrm{d}\gamma } = 6/\sqrt {2\pi }$ is an appropriate choice. Using DNS of FST at a wide range of turbulent Froude numbers ${\textit {Fr}}_T^2$ , we show in §§ 4 and 5 that the resulting $\delta _s$ indeed characterises the surface layer. Specifically, scaling by $z^*$ and turbulence properties measured at the bottom of the surface layer, $z^*=-0.5$ , collapses measurements in the surface layer across a wide range of ${\textit {Fr}}_T^2$ ; and $z^*=-0.5$ is appropriate to separate the surface layer where turbulence is directly affected by the free surface from the region beneath where turbulence is only indirectly affected. These two insights suggest the possibility of distinct turbulence closure modelling within the surface layer, and in § 6 we show how $\delta _s$ can be predicted a priori.

3. Direct numerical simulation of statistically steady FST

To study FST, we pick a flow which isolates FST from other free-surface mechanisms (e.g. waves) and is homogeneous in the horizontal and temporal dimensions. This flow is a forced free-surface turbulence, where isotropic turbulence is continuously forced deep beneath the surface to obtain statistically steady turbulence at the surface. Simulations of this flow were first performed by Guo & Shen (Reference Guo and Shen2009) for ${\textit {Fr}}_T^2\ll 1$ , and have recently been extended to moderate and large ${\textit {Fr}}_T^2$ (Calado & Balaras Reference Calado and Balaras2025; Gaylo & Yue Reference Gaylo and Yue2025). In this work, we analyse the simulations from Gaylo & Yue (Reference Gaylo and Yue2025) and perform additional simulations to obtain even larger ${\textit {Fr}}_T^2$ (see table 1). For completeness, we summarise the simulations from Gaylo & Yue (Reference Gaylo and Yue2025) and note the small change needed for the new, larger- ${\textit {Fr}}_T^2$ simulations.

Table 1.

For simulations from Gaylo & Yue (Reference Gaylo and Yue2025) (GY25), the labels used in that paper are provided. All values are non-dimensionalised by the characteristic turbulent scales in the forcing region. Turbulence properties $u_{\textit{rms}}$ and $\varepsilon$ are measured at the bottom of the surface layer ( $z^* = -0.5$ ). The characteristic length scale $\ell = u_{\textit{rms}}^3/\varepsilon$ is used to calculate the near-surface turbulent Froude number ${\textit {Fr}}^2_T =(u_{\textit{rms}}^2/\ell ) / {g}$ and the turbulent Reynolds number $\textit{Re}_T = (u_{\textit{rms}} \ell ) / \nu _w$ .

We perform DNS by solving the three-dimensional incompressible Navier–Stokes equation for two immiscible phases:

(3.1a) \begin{equation} \frac {\partial \rho \boldsymbol{u}}{\partial t} + \boldsymbol{u} {\boldsymbol{\cdot }} {\boldsymbol{\nabla }} \left (\rho \boldsymbol{u}\right ) = - {{\boldsymbol{\nabla }} p} + {{\boldsymbol{\nabla }}{\boldsymbol{\cdot }} \boldsymbol{\tau }} - \rho {g} \hat {\boldsymbol{e}}_z + \rho \boldsymbol{f}, \end{equation}

(3.1b) \begin{equation} {{\boldsymbol{\nabla }}{\boldsymbol{\cdot }} \boldsymbol{u}} = 0, \end{equation}

(3.1c) \begin{equation} \frac {\partial c}{\partial t} + \boldsymbol{u} {\boldsymbol{\cdot }} {\boldsymbol{\nabla }c} = 0 , \end{equation}

where $\boldsymbol{\tau } = \mu ({{\boldsymbol{\nabla }} \boldsymbol{u}}+{{\boldsymbol{\nabla }} \boldsymbol{u}}^{\mathrm{T}} )$ is the viscous stress tensor. Numerically, we solve these equations on a uniform Cartesian grid using a second-order finite-volume scheme (Campbell Reference Campbell2014; Yu et al. Reference Yu, Hendrickson, Campbell and Yue2019). In (3.1a ) all terms but pressure are treated explicitly, and the pressure projection method is used to solve for pressure and satisfy (3.1b ). For (3.1c ), we use the conservative volume-of-fluid method, which maintains a sharp interface between air and water and ensures volume conservation of both to machine precision (Weymouth & Yue Reference Weymouth and Yue2010). Density is given by

(3.2) \begin{equation} \rho \equiv \rho _a c + (1-c) \rho _w, \end{equation}

and similar for viscosity $\mu$ . We use realistic air–water values of $\rho _{a}/\rho _{w}=0.00123$ and $\mu _{a}/\mu _{w}=0.0159$ . To focus on the scaling with ${\textit {Fr}}_T^2$ rather than Weber number $\textit {We}$ , we neglect surface tension, in effect setting ${\textit {We}}=\infty$ . While this means that the behaviour of small, sub-Hinze-scale bubbles is not resolved (Hinze Reference Hinze1955; Ruth et al. Reference Ruth, Aiyer, Rivière, Perrard and Deike2022), studies of this flow including surface tension have confirmed that its effect on turbulence and large bubbles is small for large Froude and Weber numbers (Calado & Balaras Reference Calado and Balaras2025; Gaylo & Yue Reference Gaylo and Yue2025).

To obtain statistically steady turbulence, energy is added deep beneath the free surface through linear forcing (Lundgren Reference Lundgren2003; Rosales & Meneveau 2005), where $\boldsymbol{f} = A \boldsymbol{u}' F_b(z_c)$ and

(3.3) \begin{equation} F_b(z_c) = \begin{cases} 0 & z_c \gt l_b + l_d \\ \dfrac{1}{2} \left (1 - \cos {\left [\dfrac{\pi }{l_d}(z_c - l_b - l_d)\right ]}\right ) & z_c\le l_b + l_d \\ 1 & z_c \le l_b \\ \end{cases} \end{equation}

keeps the forcing primarily in a forcing region of height $2 l_b = 3 \pi$ centred around $z_c$ . Using the method described by Rosales & Meneveau (2005), the value of $A$ is chosen dynamically to give characteristic turbulent scales of unity ( $\varepsilon$ and $u_{\textit{rms}}=1$ ) within this forcing region. There are damping regions $l_d=\pi /2$ on either side of the forcing region, then a free length $l_f=\pi /2$ before the bottom of the domain or the quiescent free surface at $z=0$ . Gaylo & Yue (Reference Gaylo and Yue2025) use an air gap of height $\pi$ above the quiescent free surface. For the larger- ${\textit {Fr}}_T^2$ simulations performed in this work ( ${g}^{-1}\gt 2$ ), the air gap was increased to $1.75 \pi$ to avoid the top boundary affecting the turbulence. The horizontal domain length is $2 \pi$ . Gaylo & Yue (Reference Gaylo and Yue2025) show that a $256^2\times 768$ grid resolves the Kolmogorov microscale ( $\eta _T / \Delta \approx 1.8$ near the surface), and our increased air gap makes the grid $256^2\times 864$ . After the simulation has reached a statistically steady state, 400 evenly spaced samples over 128 eddy turnover times $u_{\textit{rms}}^2/\varepsilon$ (based on turbulence in the forcing region, $\varepsilon$ and $u_{\textit{rms}}=1$ ) are used for the temporal averaging in (2.2).

Our interest is the turbulence levels near the surface, so it is useful to define non-dimensional numbers based on the near-surface turbulence. At $z^*=-0.5$ , we measure the turbulence kinetic energy

(3.4) \begin{equation} {\overline {\rho }} \tilde {k} \equiv {\overline {\frac {1}{2} \rho \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}}}, \end{equation}

which gives a characteristic velocity $u_{\textit{rms}} = \sqrt {2\tilde {k}/3}$ , and the dissipation rate

(3.5) \begin{equation} {\overline {\rho }} \varepsilon \equiv {\overline {\tau _{\textit{ij}}\, \partial u_i/\partial x_{\!j}}}. \end{equation}

These measurements are shown in table 1, and are used to define near-surface non-dimensional numbers. Most relevant here is turbulent Froude number ${\textit {Fr}}^2_T =\varepsilon /u_{\textit{rms}} {g}$ . We note that using the turbulence levels at $z^*=-0.5$ gives only slightly different non-dimensional numbers from those using a near-surface water-phase average (Gaylo & Yue Reference Gaylo and Yue2025, (3.4)). In § 5 we show that the turbulence levels at $z^*=-0.5$ best characterise the turbulence within the surface layer.

4. Intermittency in the surface layer

From the DNS results, we start by analysing the intermittency factor $\gamma (z)$ . Figure 1 shows $\gamma (z)$ for the range of ${\textit {Fr}}_T^2$ , before and after scaling relative depth $z-\bar {\eta }$ by $\delta _s$ . Given that $\delta _s$ is defined by the derivative of $\gamma$ at $z=\bar {\eta }$ , the distributions necessarily collapse at the mean free surface. However, we also see that scaling by $\delta _s$ also does a good job in collapsing the distribution of $\gamma (z)$ throughout the surface layer. For $z^*\lt -0.5$ we do see some small ${\textit {Fr}}_T^2$ -dependent differences, and in § 4.1 we confirm that these are almost entirely the result of ${\textit {Fr}}_T^2$ -dependent bubble entrainment rates. In § 4.2 we investigate the shape of $\gamma (z)$ and show that it is closer to a logistic distribution than to the Gaussian distribution discussed by Brocchini & Peregrine (Reference Brocchini and Peregrine2001b ).

Figure 1.

Intermittency factor $\gamma$ as a function of (a) height relative to the mean free surface and (b) height scaled by surface-layer thickness. Solid curves indicate the mean free surface $\bar {\eta }$ and in (b) dashed lines show the extent of the surface layer.

4.1. Separating out the effects of bubbles and droplets

The intermittency factor $\gamma$ includes the effects of three different types of features: the free surface, bubbles and droplets. To separate these, we start by applying a connected component-labelling algorithm (Hendrickson, Weymouth & Yue Reference Hendrickson, Weymouth and Yue2020) to identify contiguous regions of air that are disconnected from the bulk region of air above the free surface (bubbles) and contiguous regions of water that are disconnected from the bulk region of water below the free surface (droplets). From this post-processing, we define two new binary colour functions:

(4.1a) \begin{equation} c_{\!D}(\boldsymbol{x},t) \equiv \begin{cases} 1 & \text{if } \boldsymbol{x}\in \text{droplet of water} \\ 0 & \text{otherwise}; \end{cases} \end{equation}

(4.1b) \begin{equation} c_{\!B}(\boldsymbol{x},t) \equiv \begin{cases} -1 & \text{if } \boldsymbol{x}\in \text{bubble of air} \\ 0 & \text{otherwise}. \end{cases} \end{equation}

Subtracting these from the original colour function, we obtain a third binary colour function:

(4.1c) \begin{equation} c_0 \equiv c - (c_{\!D} + c_{\!B}), \end{equation}

which describes the mixture of air and water if all bubbles were replaced with water and all droplets replaced with air. This $c_0$ isolates the direct effects of the free surface. Figure 2 shows a representative example of splitting the original colour function $c$ into these three colour functions.

Figure 2.

Two-dimensional slice over $z^{*}\in [-1,+1]$ from the ${\textit {Fr}}_T^2=0.13$ simulation (GY25, a.iv), showing: (a) the original colour function $c$ ; (b) the free-surface colour function $c_0$ ; (c) the droplet colour function $c_{\!D}$ ; and (d) the bubble colour function $1+c_{\!B}$ . Colour scale from $0$ (white) to $1$ (blue).

Figure 3.

The three subcomponents of intermittency factor $\gamma$ across a range of ${\textit {Fr}}_T^2$ (see figure 1 a for colour legend). Here $\gamma _0$ describes the free surface, $\gamma _{\!B}$ bubbles and $\gamma _{\!D}$ droplets. Note the difference in magnitude of the horizontal axes.

Qualitatively, we note some key differences between figure 2 from DNS (see also the rendering from Gaylo & Yue (Reference Gaylo and Yue2025, figure 3)) and the notional sketch of the surface layer from Brocchini & Peregrine (Reference Brocchini and Peregrine2001b , figure 1). Even for ${\textit {Fr}}_T^2$ as high as $0.26$ , DNS shows that an intact free surface is still easily identifiable. This observation is consistent with experimental work on open-channel flow. While the surface layer appears to the human eye like a uniform mixture of air and water, high-speed photography shows the free surface, while highly distorted, is still mostly intact (Killen Reference Killen1968; Wilhelms & Gulliver Reference Wilhelms and Gulliver2005).

To quantify the air–water mixture in the surface layer, we apply averaging like in (2.3) to each of the three new colour functions to obtain $\gamma _0$ , $\gamma _{\!B}$ and $\gamma _{\!D}$ , shown in figure 3. These split the intermittency factor into the contributions from each of the three phenomena:

(4.2) \begin{equation} \gamma =\gamma _0+\gamma _{\!B}+\gamma _{\!D}. \end{equation}

By rough order of magnitude, we see $\gamma _0$ is 100 times larger than $\gamma _{\!B}$ , which is 100 times larger than $\gamma _{\!D}$ . Consistent with the qualitative observations, we see that the effects of the free surface ( $\gamma _0$ ) dominate the surface layer. For droplets, we see a peak in $\gamma _{\!D}$ within the upper surface layer ( $z^*\approx 0.2$ ) and then a rapid decay with increasing height. The behaviour for bubbles is very different: $-\gamma _{\!B}$ increases with increasing depth throughout the surface layer, and then only very slowly decays with further increasing depth. Because $\rho _a/\rho _w \ll 1$ , the mass of a bubble is negligible compared with the water surrounding it, meaning that bubble moment is dominated by turbulence in the water rather than the bubble’s momentum. For these ${\textit {Fr}}_T^2$ where turbulence overcomes the bubble’s buoyancy, this means turbulence in the water advects the bubble deep beneath the surface. As expected, $-\gamma _{\!B}$ increases with ${\textit {Fr}}_T^2$ , as bubble entrainment scales like ${\textit {Fr}}_T^6$ (Gaylo & Yue Reference Gaylo and Yue2025).

Focusing on the free-surface effects, we see in figure 3 that the distribution of $\gamma _0$ collapses very well when scaling by $z^*$ . Because $\gamma _{\!D}$ is so small, droplet effects are negligible. The small- ${\textit {Fr}}_T^2$ dependence of $\gamma$ for $z^*\lt 0$ in figure 1(b) is confirmed to be due to the effects of bubbles. For significantly larger bubble void fractions (i.e. larger ${\textit {Fr}}_T^2$ ) one may need to calculate $\delta _s$ using $\gamma _0$ rather than $\gamma$ in (2.6); however, for these simulations the difference is negligible.

4.2. Describing the distribution of intermittency

We now consider the distribution of the intermittency factor. To highlight the tail behaviour of the distribution, figure 4 shows the derivative of intermittency with $z^*$ . Focusing first on $z^*\gt 0$ , we see tail behaviour like $\mathrm{d} \gamma /\mathrm{d} z\propto - \exp [\lvert z\rvert ]$ . This is opposed to $\mathrm{d} \gamma /\mathrm{d} z\propto - \exp [z^2]$ predicted by (2.5). Instead of Gaussian, this tail behaviour is characteristic of a Logistic distribution:

(4.3) \begin{equation} \gamma _{\textit{Logistic}}(z) = \frac {1}{2} \left [1 - \tanh \left (\frac {z-\bar {\eta }}{s}\right )\right ]\!, \end{equation}

where $s$ is a scale parameter proportional to the standard deviation $\sigma$ . The behaviour of $\gamma _0$ (figure 4 b) is the same as $\gamma$ (figure 4 a) for $z^*\gt 0$ , demonstrating that this logistic behaviour is the result of free-surface phenomenon rather than droplets. We can compare this result with experiments by Ruth & Coletti (Reference Ruth and Coletti2024), who studied a similar horizontally and temporally steady FST forced from below, but at ${\textit {Fr}}_T^2\ll 1$ where the free surface is not broken. For such a single-valued surface, $-\mathrm{d} \gamma /\mathrm{d} z$ is equivalent to the probability distribution function of free-surface deformations. Their results also appear to follow $\propto \exp [\lvert z\rvert ]$ (Ruth & Coletti Reference Ruth and Coletti2024, figure 5a). A physical interpretation of the logistic distribution of $\gamma$ comes from Waclawczyk (Reference Waclawczyk2021), who noted analogies between the effects of macroscopic turbulence and mesoscopic thermal fluctuations on intermittency at the respective scales.

Figure 4.

Distribution of (a) total intermittency and (b) intermittency excluding droplets and bubbles across a range of ${\textit {Fr}}_T^2$ (see figure 1 a for colour legend). Distributions are compared with $\gamma _{\textit{Gaussian}}$ (small-dashed curves) from (2.5) and $\gamma _{\textit{Logistic}}$ (black solid curves) from (4.3). Vertical dashed lines show the extent of the surface layer.

For $z^*\lt 0$ , the distribution of intermittency is more complicated. Within the lower surface layer ( $z^*\in [-0.5, 0]$ ) it roughly follows a Gaussian distribution; however, beneath the surface layer ( $z^*\lt -0.5$ ) there is a clear departure and the tail is much larger than predicted by (2.5). While removing the effect of bubbles with $\gamma _0$ reduces the tail, there is still a slow decay compared with (2.5). The logistic distribution does a better job in modelling the tail behaviour, but still does not fully capture it. Physically, the complexity of the tail behaviour can be explained by examining $c_0$ in figure 2(a). Even after removing bubbles, there are still pockets of air which, in the two-dimensional slice of the domain, appear separate from the bulk air above the free surface. By moving the location of the two-dimensional slice, we confirm that these are long filaments of air that reach deep beneath the free surface while still being connected to the bulk air above. Long filaments are also evident in the three-dimensional rendering from Gaylo & Yue (Reference Gaylo and Yue2025, figure 3). The filaments are presumably the result of bubbles that recently reconnected to the surface (degassing) or are soon to disconnect from the surface to form a bubble (entrainment), making it rather unclear if filaments should be treated like free-surface phenomena (as done here) or bubble phenomena.

5. Turbulence in the surface layer

We now examine the turbulence in and around the surface layer. First, § 5.1 examines the isotropy, confirming that, as noted before (Yu et al. Reference Yu, Hendrickson, Campbell and Yue2019; Gaylo & Yue Reference Gaylo and Yue2025), for sufficient ${\textit {Fr}}_T^2$ near-surface turbulence is nearly isotropic. In § 5.2 we examine the terms in the averaged vertical momentum ( $\overline {\rho w}$ ) conservation equation and in § 5.3 the terms in the averaged kinetic energy ( ${\overline {\rho }}\tilde {k}$ ) conservation equation. Throughout this analysis we show that relevant properties collapse across a wide range of ${\textit {Fr}}_T^2$ when properly scaled against $z^*$ , and that, partially for strong FST ( ${\textit {Fr}}_T^2\gt 0.1$ ), the free surface is not felt by the turbulence outside of the surface layer, with $z^*=-0.5$ being the transition point. Later in § 6 we discuss how these insights could be applied to RANS turbulence modelling.

5.1. Isotropy

For weak FST ( ${\textit {Fr}}_T^2\ll 1$ ), the restoring force of gravity prevents large deformations of the free surface, suppressing vertical fluctuations. This blockage effect creates highly anisotropic turbulence (Shen et al. Reference Shen, Zhang, Yue and Triantafyllou1999; Guo & Shen Reference Guo and Shen2010; Ruth & Coletti Reference Ruth and Coletti2024). However, as ${\textit {Fr}}_T^2$ increases the strength of the turbulence increases relative to gravity. Yu et al. (Reference Yu, Hendrickson, Campbell and Yue2019) show that for sufficiently large ${\textit {Fr}}_T^2$ the blockage effectively disappears and the near-surface turbulence is nearly isotropic. Calado & Balaras (Reference Calado and Balaras2025) show that for moderate ${\textit {Fr}}_T^2$ vertical fluctuations are partially suppressed. Recently, Gaylo & Yue (Reference Gaylo and Yue2025) proposed ${\textit {Fr}}_T^2=0.1$ as the approximate transition between moderate FST where vertical fluctuations are partially suppressed and strong FST where the turbulence is nearly isotropic.

Figure 5.

(a,b) Isotropy metric $J$ and (c) relative contributions of vertical fluctuations $\phi$ across a range of ${\textit {Fr}}_T^2$ (see figure 1 a for colour legend).

Figure 6.

Minimum value of the isotropy metric $J$ in the lower surface layer ( $z^*\in [-0.5,0]$ ) as a function of ${\textit {Fr}}_T^2$ , with the approximate transition from moderate to strong FST indicated by the vertical dashed line.

To quantify isotropy, we consider the isotropy metric

(5.1) \begin{equation} J \equiv 1 - 9 \left (\frac {1}{2} b_{\textit{ij}} b_{\textit{ij}} - b_{\textit{ij}} b_{\textit{jk}} b_{\textit{ki}} \right )\!, \end{equation}

where $b_{\textit{ij}}\equiv {\overline {u_i u_{\!j}}}/{\overline {u_k u_k}} - \delta _{\textit{ij}}/3$ is the anisotropy tensor. Here $J=1$ corresponds to perfectly isotropic turbulence and $J=0$ corresponds to either one- or two-component turbulence (Pope Reference Pope2000). Figures 5(a) and 5(b) show $J$ as a function of depth, unscaled and scaled by $\delta _s$ , respectively. Scaling by $\delta _s$ does a good job in collapsing the behaviour of $J$ with depth, and across all ${\textit {Fr}}_T^2$ that there is a local minimum in isotropy in the lower surface layer, generally at $z^*\approx -0.25$ . Figure 6 shows the values of $J$ at these local minimums. As first suggested by Gaylo & Yue (Reference Gaylo and Yue2025), figure 6 confirms a clear transition from moderate FST to strong FST at ${\textit {Fr}}_T^2\approx 0.1$ . Returning to figure 5(b), we highlight that for strong FST there is not a significant decrease in isotropy until near the bottom of the surface layer ( $z^*\lt -0.5$ ). This is the first indication that, for strong FST, turbulence beneath the surface layer does not ‘feel’ the presence of the free surface.

The cause of anisotropy in FST is the effect of gravity on vertical fluctuations, so, to understand $J$ , it is useful to consider the relative magnitude of vertical fluctuations,

(5.2) \begin{equation} \phi \equiv 3{\overline {ww}}/ {\overline {u_k u_k}}, \end{equation}

shown in figure 5(c). For moderate FST, we see $\phi \lt 1$ below and throughout the surface layer, indicating that vertical fluctuations are suppressed. Focusing on strong FST, we see that in the lower surface layer ( $z^*\in [-0.5,0]$ ) vertical fluctuations are suppressed, but in the upper surface layer ( $z^*\in [0,0.5]$ ) vertical fluctuations are in fact amplified ( $\phi \gt 1$ ). This helps explain the shape of $J$ in strong FST. By symmetry, our FST flow statistics will satisfy ${\overline {uu}}={\overline {vv}}$ and ${\overline {uv}}=0$ , from which one can show that

(5.3) \begin{equation} J = \frac {1}{4}\, \phi \, (3-\phi )^2. \end{equation}

We note that, due to finite sample sizes, measurements will satisfy (5.3) approximately, rather than exactly. Figure 5(c) shows that the local minimum in $J$ at $z^*\approx -0.25$ is because this is where vertical fluctuations are most suppressed. In the upper surface layer, figure 5(c) shows that amplification of vertical fluctuations is strongest at $z^*\approx 0.3$ . While less pronounced than the local minimum at $z^*\approx -0.25$ , we do see for the larger ${\textit {Fr}}_T^2$ in figure 5(b) evidence of a local minimum in $J$ at $z^*\approx 0.3$ . For $\phi$ to go from ${\lt } 1$ in the lower surface layer to ${\gt } 1$ in the upper surface layer, it must pass through $\phi =1$ around the mean free surface. This explains why near-perfect isotropy ( $J\approx 1$ ) is observed at the mean free surface ( $z^*=0$ ) for strong FST in figure 5(b).

5.2. Momentum

We now examine the vertical momentum equation. We define density fluctuations $\rho '\equiv \rho - {\overline {\rho }}$ , in which case the vertical component of (3.1a ) is

(5.4) \begin{equation} \frac {\partial \rho w}{\partial t} = - {\boldsymbol{\nabla }}{\boldsymbol{\cdot }} \left (\boldsymbol{u} \rho w\right ) - \partial p_d/\partial z + \partial \tau _{i3}/\partial x_i + \rho ' {g}, \end{equation}

where dynamic pressure $p_d$ is the total pressure minus a reference pressure $p_0$ and hydrostatic pressure:

(5.5) \begin{equation} p_d \equiv p - p_0 - {g} \int ^\infty _z {\overline {\rho }} \mathrm{d} z'. \end{equation}

Averaging (5.4) and splitting Reynolds stress into two terms ${\overline {\rho ww}} = {\overline {\rho }} \,{\overline {ww}} + {\overline {\rho ' ww}}$ , we are left with the derivative of four momentum flux terms:

(5.6) \begin{equation} 0 = - \frac {\partial }{\partial z} \left ({\overline {\rho }}\,{\overline {ww}} + {\overline {\rho ' ww}} + {\overline {p_d}} - {\overline {\tau _{33}}} \right )\!. \end{equation}

The appropriate choice of $p_0$ makes the sum of momentum flux terms zero.

Figure 7.

Momentum flux across a range of ${\textit {Fr}}_T^2$ (see figure 1 a for colour legend), normalised by mean density $\overline {\rho }$ measured for each $z^*$ and characteristic velocity $u_{\textit{rms}}$ measured at $z^*=-0.5$ (see table 1). For the viscous term, we further multiply by $\textit{Re}_T$ to bring it to a similar scale.

In figure 7 we plot (per unit mass) the mean-density component of Reynolds stress $\overline {ww}$ , the fluctuating-density component of Reynolds stress ${\overline {\rho 'ww}}/{\overline {\rho }}$ and the viscous stress $-{\overline {\tau _{33}}}/{\overline {\rho }}$ . Pressure stress ${\overline {p_d}}/{\overline {\rho }}$ cancels the sum of these three stresses. We see that plotting against $z^*$ and normalising the stresses by $u_{\textit{rms}}^2$ measured at the bottom of the surface layer ( $z^*=-0.5$ ) does a reasonable job in collapsing the results, and makes clear the distinct behaviour between strong and moderate FST. Beneath the surface layer ( $z^*\lt -0.5$ ), there is no contribution from the fluctuating-density component of Reynolds stress or viscous stress. Within the surface layer, viscous stress is of order $\textit{Re}_T^{-1}$ , meaning it is negligible in this DNS and certainty negligible at larger $\textit{Re}_T$ typical of real-world FST. For strong FST, we see that $\overline {ww}$ slightly decreases in the lower surface layer as vertical fluctuations are suppressed, then slightly increases in the upper surface layer as vertical fluctuations are amplified, consistent with figure 5(c).

We now focus on the fluctuating-density component of Reynolds stress $\overline {\rho 'ww}$ . Within the lower surface layer this term becomes slightly negative, then increases through the upper surface layer. As with intermittency in § 4.1, we can split this term into separate contributions from bubbles, droplets and the free surface. Recalling (3.2), we can write density fluctuations as $\rho ' = \Delta \rho (c - \gamma )$ , where $\Delta \rho =\rho _w-\rho _a$ . Thus, the density fluctuations caused by the free surface can be obtained by $\rho '_0 \equiv \Delta \rho (c_0 - \gamma _0)$ , and similar for $\rho '_{\!B}$ describing fluctuations caused by bubbles and $\rho '_{\!D}$ describing fluctuations caused by droplets. This allows us to decompose $\overline {\rho 'ww}$ into these three terms:

(5.7) \begin{equation} {\overline {\rho ' ww}} = {\overline {\rho _0' ww}} + {\overline {\rho _{\!B}' ww}} + {\overline {\rho _{\!D}' ww}}, \end{equation}

shown in figure 8. We see that, by at least two orders of magnitude, the free-surface effect ( $\overline {\rho _0'ww}$ ) dominates bubble or droplet effects. As one might expect, droplets only have a (small) effect above the mean free surface ( $z^*=0$ ). Bubbles, on the other hand, never have a significant effect through $\overline {\rho _{\!B}'ww}$ , even for larger ${\textit {Fr}}_T^2$ where the void fraction is larger. To explain these features, we note that $-{\overline {\rho '_{\!B} ww}}$ can be interpreted as the correlation between the presence of a bubble ( $\rho '_{\!B}\lt 0$ ) and large vertical velocities ( $ww\gt {\overline {ww}}$ ). Zero correlation is what one would obtain if bubbles were treated as passive particles, showing that, for the ${\textit {Fr}}_T^2$ considered here, bubbles are primarily advected by turbulence, with buoyancy having a smaller contribution.

Figure 8.

Components of $\overline {\rho ' ww}$ from (5.7), normalised by mean density ${\overline {\rho }}(z^*)$ and characteristic velocity $u_{\textit{rms}}$ measured at $z^*=-0.5$ (see table 1) across a range of ${\textit {Fr}}_T^2$ (see figure 1 a for colour legend).

5.3. Turbulence kinetic energy

We now examine the turbulence kinetic energy equation. Taking the inner product of (3.1a ) and $\boldsymbol{u}$ :

(5.8) \begin{equation} \frac {\partial }{\partial t}\left (\frac {1}{2} \rho \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}\right ) = - {{ \boldsymbol{\nabla }}{ \boldsymbol{\cdot }} \left [\boldsymbol{u} \left (\frac {1}{2} \rho \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}\right ) + \boldsymbol{u} p - \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\tau }\right ]} - \boldsymbol{\tau } \boldsymbol{:} {{ \boldsymbol{\nabla }} \boldsymbol{u}} - {g} \rho w. \end{equation}

To obtain the governing equation for ${\overline {\rho }} \tilde {k}$ , we take the average. Noting that the averaged conservation of mass equation gives ${\overline {\rho w}}=0$ , the gravity term cancels and we have

(5.9) \begin{equation} 0=- \frac {\partial }{\partial z} \left (\frac {1}{2}{\overline {\rho w u_i u_i}} + {\overline {w p}} - {\overline {u_i \tau _{i3}}}\right ) - {\overline {\rho }} \varepsilon. \end{equation}

In § 5.3.1 we compare the evolution of $\tilde {k}$ and $\varepsilon$ between the simulations with the free surface and single-phase simulation to show, especially for strong FST, that free-surface effects are restricted to within the surface layer. In § 5.3.2 we examine $\tilde {k}$ and $\varepsilon$ within the surface layer. In § 5.3.3 we focus on the energy flux terms in (5.9). Predicting the flux of energy into the surface layer will form an important part of the turbulence model discussed in § 6.

5.3.1. Extent of the effect of the free surface on turbulence

Recall that our simulation set-up injects turbulence kinetic energy $\tilde {k}$ in a region deep beneath the free surface, which then diffuses towards the surface, while some energy is dissipated by $\varepsilon$ . Because of this set-up, turbulence levels vary with depth, even in the absence of free-surface effects (Guo & Shen Reference Guo and Shen2009). To elucidate the effect of the free surface on $\tilde {k}$ and $\varepsilon$ , we first perform a new DNS without any free surface. This simulation has the same set-up as the air–water simulations described in § 3, except the entire domain is filled with water. We find that placing the upper boundary at $z=\pi$ (leading to the same domain size as that of the air–water simulations with an air gap of $\pi$ ) is sufficient such that the upper boundary has negligible effect on the turbulence at $z=0$ (the location of the quiescent free surface in the other simulations).

Figure 9.

(a,b) Turbulence kinetic energy $\tilde {k}$ and (c,d) dissipation rate $\varepsilon$ across a range of ${\textit {Fr}}_T^2$ (see figure 1 a for colour legend), compared with a simulation without a free surface (solid black line).

Figure 9 compares $\tilde {k}$ and $\varepsilon$ from the simulations with a free surface and the new simulation without a free surface. Turbulence values are normalised by the value at $z=-\pi /2$ , the depth above which there is no forcing and the flow evolves naturally. Focusing first on figures 9(a) and 9(c), in the absence of the free surface we see the expected exponential decay of both $\tilde {k}$ and $\varepsilon$ with decreasing depth (Guo & Shen Reference Guo and Shen2009). This is also true for simulations with the free surface at deeper depths; however, at shallower depths both $\tilde {k}$ and $\varepsilon$ depart and are much larger than in the absence of the free surface.

To quantify the amplification of $\tilde {k}$ and $\varepsilon$ due to the free surface, in figures 9(b) and 9(d) we compare the value of $\tilde {k}$ and $\varepsilon$ in simulations with the free surface to the value at the same $z$ in the simulation without the free surface. We then normalise and plot in terms of $z^*$ . While there are some small variations beneath, we see that, especially for strong FST, significant amplification of $\tilde {k}$ and $\varepsilon$ only happens within and above the surface layer ( $z^*\gt -0.5$ ). These results demonstrate that the free surface does not directly affect the turbulence beneath the surface layer, and that $z^*=-0.5$ is an appropriate definition of the lower limit of the surface layer. This would suggest that outside of the surface layer standard turbulence closure models are applicable; however, within the surface layer a model must account for the free-surface effects, namely an increase of $\tilde {k}$ by as much as a factor of 5 and $\varepsilon$ by as much as a factor of 10 relative to flow without a free surface, as seen in figures 9(b) and 9(d), respectively.

5.3.2. Turbulence within the surface layer

We now consider the turbulence levels within the surface layer, shown in figure 10. For the turbulence kinetic energy per unit mass $\tilde {k}$ , we see that there is relatively little variation within the surface layer, especially for the lower surface layer ( $z^*\in [-0.5,0]$ ); throughout the surface layer, $\tilde {k}$ does not depart significantly from the value at the bottom of the surface layer ( $z^*=-0.5$ ). Furthermore, due to $\rho _a\ll \rho _w$ the small variation in the upper surface layer ( $z^*\in [0, 0.5]$ ) is negligible when one considers the turbulence kinetic energy per unit volume, ${\overline {\rho }} \tilde {k}$ .

Figure 10.

Turbulence kinetic energy $\tilde {k}$ and dissipation rate $\varepsilon$ , per unit volume and per unit mass across a range of ${\textit {Fr}}_T^2$ (see figure 1 a for colour legend).

For the turbulent dissipation rate per unit mass $\varepsilon$ , the variations in the surface layer are slightly larger than for $\tilde {k}$ , but still $\varepsilon$ measured at the bottom of the surface layer does a good job in characterising dissipation rates within the surface layer, especially the lower surface layer ( $z^*\in [-0.5,0]$ ). Because $\rho _a\ll \rho _w$ the majority of the dissipation of energy happens in the lower surface layer. For turbulent dissipation rate per unit volume, ${\overline {\rho }} \varepsilon$ , we see a strong collapse for ${\textit {Fr}}_T^2\gt 0.2$ , and note that the transition towards this behaviour begins around ${\textit {Fr}}_T^2=0.1$ .

5.3.3. Energy flux into the surface layer

We now investigate the three energy flux terms in (5.9), shown in figure 11. Similar to the momentum flux terms, we plot against $z^*$ and normalise the stresses by $u_{\textit{rms}}^3$ measured at the bottom of the surface layer ( $z^*=-0.5$ ). First, we note that, as expected, the viscous diffusion term $\overline {u_i \tau _{i3}}$ is of order $\textit{Re}_T^{-1}$ , meaning it is negligible in this DNS and certainly negligible at larger $\textit{Re}_T$ typical of real-world FST. For the remaining terms, we see turbulent diffusion $\overline {\rho w u_i u_i}$ moves energy upward into the surface layer, and pressure diffusion $\overline {w p}$ partially counteracts this, moving some energy out of the surface layer. Of note, despite scaling by $z^*$ and $u_{\textit{rms}}$ , we still see ${\textit {Fr}}_T^2$ dependence in these energy flux terms.

Figure 11.

Turbulence kinetic energy flux terms across a range of ${\textit {Fr}}_T^2$ (see figure 1 a for colour legend), normalised by water density $\rho _w$ and characteristic velocity $u_{\textit{rms}}$ measured at $z^*=-0.5$ (see table 1). For the viscous term, we further multiply by $\textit{Re}_T$ to bring it to a similar scale.

While we have shown that the turbulence beneath the surface layer does not directly feel the free surface, the increased energy dissipation within the surface layer is still important to the overall energy budget of the flow. Thus, the energy flux into the surface layer (i.e. the net work done on the surface layer by the turbulence beneath) is an important quantity for any closure model. This net work is given by the sum of turbulent, pressure and viscous diffusion measured at the bottom of the surface layer:

(5.10) \begin{equation} W = \left [\frac {1}{2}{\overline {\rho w u_i u_i}} + {\overline {w p}} - {\overline {u_i \tau _{i3}}}\right ]_{z^* = -0.5}. \end{equation}

Because the viscous term is negligible in this DNS and certainly negligible at larger $\textit{Re}_T$ typical of real-world FST, we consider

(5.11) \begin{equation} W \approx \left [\frac {1}{2}{\overline {\rho w u_i u_i}} + {\overline {w p}}\right ]_{z^* = -0.5}. \end{equation}

In figure 12 we show the net work non-dimensionalised by $\rho _w$ and $u_{\textit{rms}}$ measured at $z^*=-0.5$ . We start by considering a linear scaling with ${\textit {Fr}}_T^2$ . Linear regression gives

(5.12) \begin{equation} W / ( \rho _w u_{\textit{rms}}^3) = C_{W} \, {\textit {Fr}}^2_T, \end{equation}

where the 95 % confidence interval is $C_W = 4.6\pm 0.8$ . To avoid overfitting we stick to this linear model, but note that the relatively low $R^2=0.641$ and some indication of a split at ${\textit {Fr}}_T^2=0.1$ suggest some improvement could be possible with additional simulations. Dimensionalising, we see how the work done on the surface layer depends on gravity and the turbulence ( $\varepsilon$ and $u_{\textit{rms}}$ ) measured at the bottom of the surface layer:

(5.13) \begin{equation} W /\rho _w = C_{W} \, u_{\textit{rms}}^2 \varepsilon {g}^{-1}. \end{equation}

While this described the ${\textit {Fr}}_T^2\in [0.03,0.3]$ we consider, we note that (5.13) is not expected to apply to weak FST ( ${\textit {Fr}}_T^2\ll 1$ ), where near-surface turbulence is strongly affected by viscosity (Shen et al. Reference Shen, Zhang, Yue and Triantafyllou1999).

Figure 12.

Work done on the surface layer by turbulence beneath (5.11) as a function of ${\textit {Fr}}_T^2$ . The solid line shows the linear fit to (5.12) with $C_W=4.6$ ( $R^2=0.641$ ).

6. Surface-layer boundary conditions for RANS

We have shown that, particularly for strong FST, there is a clear difference between the turbulence beneath the surface layer and the turbulence within the surface layer, with $z^*=-0.5$ doing a good job in separating the two. Because the turbulence beneath the surface layer is not directly affected by the free surface, standard RANS closure models will likely still be appropriate there. One approach then is to model the (indirect) effects of the surface layer through a boundary condition applied at $z^*=-0.5$ (dimensionalised, $z=\bar {\eta }- \delta _s/2$ ), analogous to wall models for boundary layers.

We first need $\delta _s$ to know where to apply this boundary condition. So far, we have obtained the surface-layer thickness $\delta _s$ a posteriori using measurements of $\gamma (z)$ , as described by (2.6). For a predictive model appropriate to RANS, we need to be able to predict $\delta _s$ based on the turbulence levels $\varepsilon$ and $u_{\textit{rms}}^2=2 \tilde {k}/3$ . Figure 13 shows the scaling of $\delta _s$ , normalised by the characteristic turbulence length scale $\ell = u_{\textit{rms}}^3/\varepsilon$ , as a function of ${\textit {Fr}}_T^2$ . We see a strong linear relationship:

(6.1) \begin{equation} \delta _s/\ell = C_\delta \, {\textit {Fr}}^2_T, \end{equation}

where the 95 % confidence interval is $C_\delta = 11.1\pm 1.1$ . Dimensionalising:

(6.2) \begin{equation} \delta _s = C_\delta \, u_{\textit{rms}}^2 {g}^{-1}. \end{equation}

This is the scaling one would expect from a simple wave argument: that energy is roughly evenly distributed between gravitational potential energy ( $\propto \rho _w {g} \delta _s$ ) and kinetic energy ( $\propto \rho _w u_{\textit{rms}}^2$ ). We note that $u_{\textit{rms}}$ is measured at $z^*=-0.5$ , which is a function of $\delta _s$ . This makes (6.2) an implicit equation, but because $u_{\textit{rms}}$ changes smoothly with depth (see $\tilde {k}$ in figure 10), solving (6.2) should not be a challenge.

Figure 13.

Surface-layer thickness $\delta _s$ normalised by the characteristic length scale $\ell = u_{\textit{rms}}^3/\varepsilon$ measured at $z^*=-0.5$ , as a function of ${\textit {Fr}}_T^2$ . The solid line shows the linear fit to (6.1) with $C_\delta =11.1$ ( $R^2=0.909$ ).

For a $k$ – $\varepsilon$ RANS model, we need to apply a boundary condition to the turbulence kinetic energy transport equation. This boundary condition describes the flux of turbulence kinetic energy across the surface-layer boundary at $z^*=-0.5$ . Brocchini (Reference Brocchini2002) derived the form of this boundary condition for free-surface flows in general, which could include turbulent energy production by mean shear within the surface layer. In our horizontally and temporally homogeneous simulations, energy flux is only driven by turbulent dissipation within the surface layer. This flux of $\tilde {k}$ into the surface layer (per unit area) is given by (5.13). Using (6.2), we can rewrite (5.13) as

(6.3) \begin{equation} W / \rho _w = \left (C_W/C_\delta \right ) \varepsilon \, \delta _s. \end{equation}

While Brocchini (Reference Brocchini2002) considers $k={\overline {\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}}}/2$ rather than $\tilde {k} = {\overline {\rho \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}}}/2{\overline {\rho }}$ , beneath the surface layer $\gamma \approx 1$ and the difference should be negligible. Therefore, this $C_W/C_\delta = 0.41$ is equivalent to $C_{{diss}}$ discussed by Brocchini (Reference Brocchini2002, (34)).

Finally, while we calculate the energy flux across the surface-layer boundary directly, if we assume that dissipation in air is negligible (vacuum assumption), it could also be determined based on the integral of dissipation (Brocchini Reference Brocchini2002; Iravani et al. Reference Iravani, Badiei and Brocchini2020):

(6.4) \begin{equation} W / \rho _w \approx \delta _s \int _{-0.5}^{0.5} \varepsilon \gamma\, \mathrm{d} z^*. \end{equation}

Recalling figure 10, as a first approximation one could model the dissipation rate per unit mass in the surface layer as a constant equal to $\varepsilon$ measured at $z^*=-0.5$ . Using any distribution of $\gamma$ where $\gamma -0.5$ is symmetric about $\bar {\eta }$ (this includes both the Gaussian (2.5) and logistic (4.3) distributions), we obtain $W/\rho _w \approx 0.5 \varepsilon \, \delta _s$ , or $C_W/C_\delta \approx 0.5$ . For comparison, we also consider dissipation decreasing linearly to zero at the top of the surface layer, as assumed by Iravani et al. (Reference Iravani, Badiei and Brocchini2020). In that case, (6.4) gives $C_W/C_\delta \approx 0.361$ if one uses (2.5) for $\gamma$ and $C_W/C_\delta \approx 0.329$ if one uses (4.3) for $\gamma$ . Thus, our directly measured $C_W/C_\delta = 0.41 \approx 0.5$ is consistent with the simplified description of dissipation being constant in the surface layer, although the difference if one assumed linearly decreasing dissipation is not too large.