3.1. Flow problem

The governing equations for the flow are the incompressible Navier–Stokes equations and the continuity equation:

(3.1) \begin{align} \frac {\partial u_i}{\partial t} + \frac {\partial (u_i u_j)}{\partial x_j} & = - \frac {\partial p}{\partial x_i} + \frac {1}{\textit {Re}} \frac {\partial ^2 u_i}{\partial x_j \partial x_j} + f_i \, , \end{align}

(3.2) \begin{align} \frac {\partial u_i}{\partial x_i} & = 0 \, , \end{align}

where $u_i$ ( $i = 1, 2, 3)$ are the instantaneous velocity components, $p$ is the dynamic pressure, i.e. the pressure deviation from the hydrostatic pressure, normalised by $\rho U^2$ , where $\rho$ is the density, and $U$ is the characteristic velocity scale. The term $f_i$ on the right-hand side is a depth-dependent forcing term that generates forced isotropic turbulence (defined below). The Reynolds number is $Re = U L / \nu$ , with characteristic length $L$ and kinematic viscosity $\nu$ . In our dimensionless simulation set-up, we follow the convention of numerical simulations of isotropic turbulence and set the horizontal periodic domain dimensionless length to $2\pi$ (so that the smallest dimensionless wavenumber has unit value). Fixing the dimensionless domain size implicitly sets the reference length $L$ . As the large-scale motions of forced isotropic turbulence are proportional to the domain size, $L$ may be considered to be tied to the length scale of the large-scale motions in turbulence, specifically the integral length scale, which is further discussed in § 3.2. The characteristic scale $U$ is associated with the forcing $f_i$ , where we fix the forcing magnitude in the dimensionless system, which is equivalent to setting a dimensionless turbulence time scale (details below).

The flow variables are computed in a three-dimensional domain that is periodic in the horizontal directions ( $x_1, x_2$ , i.e. $x, y$ ) and restricted above by a deformable free surface at $x_3 = z = -\eta$ , where $\eta = \eta (x,y,t)$ is the surface elevation and positive $z$ -direction is taken to be from the surface and down into the flow (with average depth $\bar {z} = 0$ at the surface). At the bottom boundary, a free-slip boundary condition is used to eliminate the viscous effects of the bottom boundary.

The flow is subjected to the kinematic and dynamic boundary conditions at the free surface. The kinematic boundary condition states that the upper boundary is material, so that particles on the surface remain on the surface from a Lagrangian point of view, that is, we can have no penetration of the free surface. This constraint results in the following evolution equation for $\eta$ :

(3.3) \begin{equation} \frac {\partial \eta }{\partial t} = -w - u\frac {\partial \eta }{\partial x} - v \frac {\partial \eta }{\partial y} , \quad \text {at } z=-\eta (x,y). \end{equation}

The dynamic boundary condition states the stress balance at the free surface:

(3.4) \begin{align} \boldsymbol {n}\boldsymbol {\cdot }\mathsf {\tau }\boldsymbol {\cdot }\boldsymbol {n}^{\textrm{T}} &= {We}^{-1}\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {n}, \end{align}

(3.5) \begin{align} \boldsymbol {t}_x\boldsymbol {\cdot }\mathsf {\tau }\boldsymbol {\cdot }\boldsymbol {n}^{\textrm{T}} &= 0, \end{align}

(3.6) \begin{align} \boldsymbol {t}_y\boldsymbol {\cdot }\mathsf {\tau }\boldsymbol {\cdot }\boldsymbol {n}^{\textrm{T}} &= 0, \end{align}

where $\mathsf {\tau }=-(p+z/Fr^2)\mathsf {I} + {\textit {Re}}^{-1}(\boldsymbol {\nabla } \boldsymbol {u}+\boldsymbol {\nabla } \boldsymbol {u}^{\textrm{T}})$ is the stress tensor, where $z/Fr^2$ is the dimensionless hydrostatic pressure, $Fr=U/\sqrt {gL}$ is the Froude number with $g$ being the gravitational acceleration, and $\mathsf {I}$ is the identity matrix; $\boldsymbol {n}$ , $\boldsymbol {t}_x$ and $\boldsymbol {t}_y$ denote the unit vectors that are normal and tangential to the water surface, where $\boldsymbol {n}$ points from air to water; ${We}^{-1}\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {n}$ is the pressure difference caused by surface tension $\sigma$ and the Weber number is defined as $We=\rho U^2 L/\sigma$ . Here, there are no normal or tangential stresses imposed by the air above the water surface. The water-side stress is balanced with the surface tension-induced stress.

Figure 2.

Illustration of the computational set-up for the isotropic turbulence interacting with a deformable free surface, including details on the structure of the free region. Regions and surface deformations not to scale.

The simulation set-up mimics the experiments of free-surface homogeneous turbulence, where turbulence, generated by an oscillating grid far beneath the water surface, diffuses towards the surface through self-interaction of turbulence vortices and viscous diffusion (e.g. Thompson & Turner Reference Thompson and Turner1975; Hopfinger & Toly Reference Hopfinger and Toly1976; Brumley & Jirka Reference Brumley and Jirka1987; McKenna & McGillis Reference McKenna and McGillis2004; Herlina & Jirka Reference Herlina and Jirka2008; Chiapponi et al. Reference Chiapponi, Longo and Tonelli2012; Lacassagne et al. Reference Lacassagne, El-Hajem, Morge, Simoens and Champagne2017), and more recent experiments, where the turbulence is generated by jets with zero net flow (e.g. Variano & Cowen Reference Variano and Cowen2013; Jamin et al. Reference Jamin, Berhanu and Falcon2024). In our simulations, turbulence is generated by a random linear forcing method (Lundgren Reference Lundgren2003; Rosales & Meneveau Reference Rosales and Meneveau2005), which energises isotropic turbulence through the body force term $f_i = a_0 F(z) u_i$ in the momentum equation (3.1). To ensure that the imposed forces do not interfere with the free surface, $f_i$ decays to zero near the surface, controlled by the depth-dependent function $F(z)$ . This function is set to $F = 1$ in the ‘centre region’ of the simulation domain, away from the surface ( $\bar {z}\in [\pi , 4\pi ]$ ). Towards the surface, $F(z)$ decreases smoothly over a ‘damping region’ ( $\bar {z}\in [0.5\pi , \pi ]$ ), and becomes zero in the ‘free region’ ( $\bar {z} \in [-\eta , 0.5\pi ]$ ); see figure 2 for domain illustration. The flow reaches equilibrium when the energy injected by $f_i$ balances dissipation. It can be shown that $a_0$ is inversely proportional to the ratio of turbulent kinetic energy to the dissipation rate (Rosales & Meneveau Reference Rosales and Meneveau2005), meaning that $a_0^{-1}$ represents a characteristic time scale of the forced turbulence. We set the dimensionless parameter $a_0=0.1$ , fixing this turbulence characteristic time scale. In this set-up, the centre region acts as a source of statistically steady isotropic turbulence and the intensity of isotropic turbulence decays over the damping and free region as the surface is approached. The isotropy of decaying turbulence is validated by Guo & Shen (Reference Guo and Shen2009). Additionally, to prevent interference from the bottom boundary, the simulation domain is vertically extended to $5\pi$ , i.e. the total domain size is $L_x \times L_y \times L_z= 2\pi \times 2\pi \times 5\pi$ , with $F(z)$ set symmetrically about $\bar {z}=2.5\pi$ .

The grid used in the simulations has $256\times 256\times 660$ points for the high-Reynolds-number cases and $128\times 128\times 348$ points for the low-Reynolds-number cases (case details in § 3.2). The grid is equidistant in $x$ - and $y$ -directions, and stretched in $z$ -direction to ensure a very fine grid near the free surface. The stretching enhances the resolution of the turbulence dynamic processes near the surface. A posteriori evaluation of the Kolmogorov scales indicates that the grid resolution meets the DNS criteria for resolving small-scale turbulent motions (Guo & Shen Reference Guo and Shen2010; Xuan & Shen Reference Xuan and Shen2022). The duration of the simulations are $160 T_\infty$ – $200 T_\infty$ (where $T_\infty$ is the eddy turnover time, defined below) after becoming statistically stationary. Sampling is done at intervals of $0.0127T_\infty$ – $0.0182T_\infty$ .

The governing equations are discretised on the computational grid using a pseudo-spectral method in the $x$ - and $y$ -directions, and a finite-difference method in the $z$ -direction. This hybrid discretisation leverages the strengths of each method, with the pseudo-spectral method providing high accuracy in horizontal planes and the finite-difference approach handling the grid clustering flexibly. To accurately model the interactions between near-surface turbulence and free-surface motions, the simulations are performed on a deforming curvilinear grid that fits the free surface. The free-surface elevation is updated by using a second-order Runge–Kutta method to integrate the kinematic boundary condition. Within each stage of the Runge–Kutta integration, the Navier–Stokes equations (3.1)–(3.2) are solved using a fractional-step method to maintain the incompressibility constraint. The pressure Poisson equation, required by the projection step in the fractional-step method, involves variable coefficients arising from the grid transformation and is solved by a fixed-point iteration method (Xuan & Shen Reference Xuan and Shen2019). Simulating with the surface-conforming curvilinear grid is computationally demanding, but ensures that the complex dynamics within the free-surface turbulence boundary layer are accurately captured. (For details and validations of the computational methods, boundary resolutions and algebraic transformation used in the grid transformation, see Xuan & Shen Reference Xuan and Shen2019).

3.2. Flow properties: length scales and non-dimensional groups

In this section, we review and define the pertinent non-dimensional groups and length scales we will later use in the analysis of turbulent flow beneath a free surface. Three dimensionless numbers characterise the free-surface flow: the Reynolds number, the Weber number and the Froude number. We have performed simulations for two Reynolds numbers, $Re = 2500$ and $Re = 1000$ . To take into consideration the effects of surface tension, $\sigma$ , we have simulated flows with three different Weber numbers, $We = 10, 20$ and $\infty$ . The Froude number is fixed at 0.1 for all simulations. As indicated by the free-surface boundary condition (3.4), the Weber number and Froude number influence the free-surface dynamics, which can in turn affect the underlying turbulence. The Froude number and Weber number measure the ability of gravity forces and surface tension to restore the integrity of the surface disturbed by turbulence, respectively. Since $We \gg 1$ and $Fr \ll 1$ , all simulations are performed for flows with gravity-dominated surface motion, where the gravity ensures the turbulence-induced deformation of the free surface remains small. The surface tension can influence small-scale surface features.

In the absence of mean flow, the turbulent Reynolds number is an appropriate measure. It is defined as $Re_\infty = {\tilde {u} l} / {\nu }$ , where $\tilde {u}$ is the ‘representative velocity’ (Tennekes & Lumley Reference Tennekes and Lumley1972, § 3.1), defined as $\tilde {u}^2 = {1}/{3} (\overline {u_iu_i})$ , and $l$ is the macroscale of the turbulence. Thus, $\tilde {u}^2$ equals two-thirds of the turbulent kinetic energy. The exact manner of averaging implied by the overbar is important because most statistical turbulence properties beneath a free surface vary with depth, a point we discuss below.

While the representative velocity is computed directly from flow variables, the macroscale is approximated (under the assumption of isotropic turbulence, verified in § 4.1) by $l \approx {\lambda _{\textrm{T}}} Re_\lambda / 15 = \tilde {u}^3/\epsilon$ , where $Re_\lambda = \tilde {u} \lambda _{\textrm{T}} / \nu$ is the Taylor–Reynolds number, and $\lambda _{\textrm{T}} = \tilde {u}\sqrt {15 \nu / \epsilon }$ is the Taylor microscale (Tennekes & Lumley Reference Tennekes and Lumley1972, p. 67). To compute the latter, the dissipation $\epsilon = 2 \nu \overline {s_{ij}s_{ji}}$ is computed from the stress components $s_{ij}=({\partial }_iu_j+{\partial }_ju_i)/2$ whence the Kolmogorov length scale, $L_{\textrm{K}} = (\nu ^3 / \epsilon )^{1/4}$ , can be calculated. A more widely used length scale than the macroscale is the integral scale of the turbulence, approximated by $L_\infty = l / 2$ (Tennekes & Lumley Reference Tennekes and Lumley1972, p. 273), whence we arrive at the form of turbulent Reynolds number most commonly found in the literature on free-surface flows: $Re_\infty = 2 \tilde {u}L_\infty / \nu$ . The integral length scale and an appropriate velocity scale yield the eddy turnover time, $T_\infty$ . Here, $T_\infty = L_\infty / \tilde {u}$ . Note that although the turbulent properties differ in the six simulations (see table 1), the integral length scale varies little. This is consistent with Rosales & Meneveau (Reference Rosales and Meneveau2005), who reported an average macroscale $l \approx 19\,\%$ of the domain width for isotropic turbulence generated by linear forcing.

Table 1.

Flow and turbulence properties. From left: case number, Reynolds number, Weber number, Froude number, turbulent Reynolds number, Taylor Reynolds number, integral length scale, Taylor length scale, Kolmogorov length scale, viscous layer thickness. All length scales are normalised by the characteristic length $L$ .

The presence of a free surface fundamentally changes the properties of the turbulence immediately beneath it. Nearest the surface, a thin viscous layer is formed, beneath which is found a thicker blockage layer (see e.g. Shen et al. Reference Shen, Zhang, Yue and Triantafyllou1999). The viscous layer has thickness ${L_\nu }\sim Re_\infty ^{-1/2}$ and occurs because the tangential viscous stress at the free surface must be zero. The blockage layer extends deeper, to a depth ${L_\beta }\sim \mathcal {O}(L_\infty )$ , and occurs because velocities must be tangential to the surface. Based on measurements by Brumley & Jirka (Reference Brumley and Jirka1987), the viscous layer thickness is conventionally taken to be $L_\nu = 2 \textrm{Re}_\infty ^{-1/2} L_\infty$ (see discussion by Calmet & Magnaudet Reference Calmet and Magnaudet2003). No single definition for $L_\beta$ has become standard, yet ${L_\beta }=L_\infty$ is often used for convenience (see e.g. Magnaudet Reference Magnaudet2003; Guo & Shen Reference Guo and Shen2010; Herlina & Wissink Reference Herlina and Wissink2014). The structure of the surface-influenced region is illustrated in figure 2. A more detailed discussion of the layered structure of free-surface flows may be found, e.g. from Shen et al. (Reference Shen, Zhang, Yue and Triantafyllou1999).

We note in passing that the literature is not consistent regarding the factor $2$ in the expression for $L_\nu$ . The reason seems to be that both the turbulent macroscale and the integral scale are in common use, which are only a factor $2$ apart and are occasionally conflated.

In the present study, non-dimensional groups and characteristic length scales which characterise the turbulence are calculated in the bulk. Recall that we consider only the free region in our analysis (see figure 2). Since there is no forcing in this region, the turbulence decays as it approaches the free surface, and hence, dissipation, root-mean-square velocity, etc. vary with depth (details in § 4). Thus, we have no single measure of bulk properties. To make a comparison between the different cases as clear as possible, all ‘bulk properties’ are measured at a single reference depth, where averages and root-mean-square values (implied by the overbars in the definitions for $\tilde {u}$ and $\epsilon$ ) are computed plane-wise, before being averaged in time to get better statistics. The selected reference depth is outside the blockage layer to ensure that the theory for isotropic turbulence is a reasonable approximation, yet close enough to the blockage-layer boundary to minimise turbulence decay between the reference plane and the near-surface region. In all cases, the reference depth is chosen as $\pi /3$ ( $ \approx 1.6L_\infty$ for case 1).

A summary of the most important flow properties of our simulated cases can be found in table 1.

4.1. Surface effects on the velocity and vorticity components

Figure 3.

Root-mean-square values of velocity (blue curves and bottom-of-panel abscissae) and vorticity (red curves, top-of-panel abscissae) for cases 1–6 in panels (a)–(f), respectively. Flow variables are plotted against averaged depth, $\bar {z}$ , normalised by the integral length scale ( $L_\infty$ ). The horizontal lines mark the viscous boundary layer limit (dashed) and the reference depth (dash-dotted).

Consider components of the velocity vector $\textbf{u}_{\textrm{rms}}$ and the vorticity field ${\boldsymbol \omega }_{\textrm{rms}}= (\nabla \times {\textbf{u}} )_{\textrm{rms}}$ as a function of depth, for cases 1–6 in figures 3(a)–3(f), respectively. The subscript ‘rms’ denotes root-mean-square taken plane-wise and averaged in time, i.e. $(\cdots )_{\textrm{rms}}= \overline {\langle (\cdots )^2\rangle ^{1/2}}$ , where $\left \langle \cdots \right \rangle =(L_xL_y)^{-1}\int (\cdots )\,\textrm{d} x\textrm{d} y$ is the horizontal-plane average (depending on $t$ and $\bar {z}$ ), and the overbar is average over time. Note that only a single horizontal component is included for each case in figure 3, as the flow is horizontally isotropic.

An average measure for ‘horizontal’ grid plane depth is used here and later, denoted $\bar {z}$ . In short, the zero-level of $\bar {z}$ is the surface elevation measurement, i.e. $z(x,y,\bar {z})|_{\bar {z}=0}=-\eta$ , and $\bar {z} \gt 0$ is the average distance from grid points $(x_i,y_j,z_k)$ on grid plane $k\gt 0$ to the surface points $(x_i,y_j,z_0)$ . In effect, $\bar {z}$ is used to approximate the vertically undulating grid as a vertically stretched Cartesian grid in the analysis. This convention enables straightforward comparison of flow variables at a grid plane on the undulating grid without interpolation. As the variations in surface elevation in a grid plane are small ( $\mathcal {O}(10^{-4} L_\infty )$ for all cases), this is an acceptable convention. (An alternative to $\bar {z}$ is to use a signed distance function for the depth, which necessitates interpolation for all measurement values, see Guo & Shen Reference Guo and Shen2010).

For all cases, we observe a decrease in $u_{z\textrm{,rms}}$ and an increase in $u_{x\textrm{,rms}}$ as the free surface is approached, which corresponds to a redistribution of kinetic energy from the vertical to the horizontal velocity components (Shen et al. Reference Shen, Zhang, Yue and Triantafyllou1999). In particular, the vertical component tends towards zero as $\bar {z}\rightarrow 0$ , as reflected in the kinematic boundary condition which restricts movement in the vertical direction for low-Froude-number flows with small surface deformation. This is known as the blockage effect (Hunt & Graham Reference Hunt and Graham1978; Perot & Moin Reference Perot and Moin1995; Shen et al. Reference Shen, Zhang, Yue and Triantafyllou1999), which forces the turbulence fluctuations to become increasingly two-dimensional as the surface is approached.

We recognise that the flow is roughly isotropic for $\bar {z} \gtrsim 1.5L_\infty$ , agreeing with the notion that the blockage layer thickness is of order $\mathcal {O}(L_\infty )$ . It also suggests that the selected depth for bulk measurements, $\pi /3$ ( $\approx 1.6 L_\infty$ for case 1) is appropriate, as it is outside the surface-influenced region.

The horizontal vorticity component $\omega _{x,\textrm{rms}}$ (red solid lines in figure 3) decreases rapidly through the viscous boundary layer as the free surface is approached. This development of the vorticity demonstrates the influence of the dynamic boundary condition, which implies that only the surface-normal vorticity component is non-zero at the free surface.

The viscous boundary layer thickness, $L_\nu$ defined in § 3.2, coincides well with the local maximum in horizontal vorticity for all cases, although the local vorticity maxima are less sharp for the lower Reynolds number cases. This supports the notion that the viscous sublayer thickness can be approximated by the depth of the shallowest maximum of the horizontal root-mean-square vorticity components (see, e.g. Calmet & Magnaudet Reference Magnaudet2003).

The results for all six cases are in line with the literature on turbulent free-surface flows, which ties the blockage layer to the kinematic boundary condition and the viscous layer to the dynamic boundary condition. We note that there is some discrepancy in the numerical values for the root-mean-square velocity components between the cases with different Weber numbers, particularly prominent when comparing the value of $\textit{u}_{\textrm{rms}}$ at $\bar {z} = 0$ for the cases with $Re = 1000, We = \infty , 20, 10$ (figure 3 d–f). Since the higher value is observed for the intermediate Weber number case, the effect can hardly be credited as a Weber number effect (even less so, as the opposite effect is observed for the cases with $Re = 2500$ ). The discrepancy is likely due to a few more large upwelling events in case 5 than in cases 4 and 6, which results in a slightly higher $\textit{u}_{\textrm{rms}}$ for case 5.

4.2. Distribution of vortex inclination angles through the surface-influenced layer

The free-surface motions in the flow examined in this study are small, due to the low Froude number. Consequently, the dynamic boundary condition of zero tangential shear stress corresponds to small horizontal vorticity, making $\omega _z$ the only significant vorticity component at the surface. For a more detailed picture of the vortices in the vicinity of the surface, we quantify the orientation of the vortices in the flow by computing the distribution of the inclination angle of the vorticity vector along horizontal planes. Shen et al. (Reference Shen, Zhang, Yue and Triantafyllou1999) reported vorticity orientation in the presence of a shear flow. The absence of a mean flow in our case allows us to isolate particular features related to the turbulent vortices in the vicinity of the free surface.

The inclination angle of the vorticity vector ${\boldsymbol \omega } = (\omega _x, \omega _y, \omega _z) \equiv ({\boldsymbol \omega }_\perp ,\omega _z)$ is measured as the angle of displacement from the vertical direction $z$ , computed by

(4.1) \begin{align} \theta = \arccos \left (\omega _z/\omega \right ), \end{align}

where $\omega = |{\boldsymbol \omega }|$ . We consider statistics of the inclination angle weighted by the normalised magnitude of the projection of $\boldsymbol \omega$ , similar to Moin & Kim (Reference Moin and Kim1985), i.e. we compute the statistics of $\kappa \theta$ , where the weighting factor is computed for each grid point by

(4.2) \begin{align} \kappa = \omega /\langle \omega \rangle . \end{align}

Figure 4.

Normalised histograms of weighted inclination angles, $\kappa \theta$ , at different depths $\bar {z}$ beneath the mean surface level for all $x,y,t$ , in multiples of the viscous layer thickness $L_\nu$ . The red line represents the distribution of weighted inclination angles for a completely isotropic state.

The distributions shown as histograms in figure 4 are results for case 1, computed for all time steps and all grid points in each grid plane at distances $\bar {z}$ beneath the mean surface level (similar figures for the other cases can be found in the supplementary material). Note that an isotropic state does not yield a uniform distribution of inclination angles $\theta$ , since the state space is larger for horizontal orientation than for vertical (e.g. $\theta =0$ and $\theta =\pi$ means $\omega _\perp =0$ , while $\theta =\pi /2$ means $\omega _z=0$ ). For reference, the distribution of inclination angles for a completely isotropic state is included in the bottom right panel (generated synthetically from $10^7$ vectors isotropically distributed in $\mathbb {R}^3$ ).

Nine values of $\bar {z}$ are chosen, starting at the surface and extending throughout the blockage layer to approximately one integral length scale into the flow (recall from table 1 that $L_\infty \approx 10.5{L_\nu }$ for case 1). Near the base of the blockage layer, the flow approaches isotropy, as reflected in the distribution $\theta$ in the bottom right frame. In all other frames, $\theta$  shows preferred distribution around $0$ , $\pi /2$ or $\pi$ , indicating that $\omega _z$ is either much smaller than or much larger than $\omega _\perp$ (for $\theta \approx \pi /2$ or $\theta \approx 0, \pi$ , respectively). Since the angle $\theta$ is the displacement from the vertical axis, a high probability around $\theta = \pi /2$ indicates a predominance of horizontal vorticity, whereas a high probability around $\theta =0, \pi$ indicates the dominance of vertical vorticity. We discern a gradual transition from vertical ( $\theta = 0, \pi$ ) to horizontal ( $\theta = \pi /2$ ) preferred orientation as one moves away from the free surface. The dynamic boundary condition forces all vorticity to be normal to the free surface at the surface, and hence the vertical orientation is strongly dominant there. Yet already in the middle of the viscous layer, the horizontal orientation is equally prominent, and it remains dominant in the upper part of the blockage layer before an essentially isotropic state is reached at $\bar {z}\approx 10{L_\nu }\approx L_\infty$ , roughly the deep boundary of the blockage layer. The results extend (and agree with) those of Shen et al. (Reference Shen, Zhang, Yue and Triantafyllou1999) wherein only the surface and bulk were considered. Much light is shed on the transition from vertical to horizontal preference when dimples and scars are next considered.

5.1. Identification of dimples and vortices

It was shown by Babiker et al. (Reference Babiker, Bjerkebæk, Xuan, Shen and Ellingsen2023) that a suitable set of SOC for dimples is that they are long-lived and approximately circular surface indentations, which complements the established understanding – a PCD in our terminology – that a true dimple is the imprint of a surface attached vortex situated where a vortex tube terminates at the free surface. With a simple computer vision method, the SOC was implemented, showing excellent accuracy and precision in identifying true dimples according to the PCD, the latter implemented using the $\lambda _2$ method for vortex identification. Because the overlap is so good, it makes no practical difference for our statistical analysis of sub-dimple flow whether dimples are identified from the surface (SOC) or interior (PCD), so since the connection is not in question here, we choose to use the latter, and henceforth use the term ‘dimple’ for what we in § 2 denoted ‘true dimple’. (This contrasts the later study of scars.)

The $\lambda _2$ method uses the three eigenvalues of the matrix $s_{ik}s_{kj} + w_{ik}w_{kj}$ , where $w_{ij} = ({\partial }_iu_j - {\partial }_ju_i) / 2$ is the vorticity tensor, $s_{ij}$ is the stress tensor (§ 3.2), and the eigenvalues $\lambda _i$ are ordered such that $\lambda _1\lt \lambda _2\lt \lambda _3$ . Fluid points which satisfy $\lambda _2\lt \lambda _{2,\textrm{th}}\leqslant 0$ are identified as belonging to a vortex core, where $\lambda _{2,\textrm{th}}$ is a threshold value which must be chosen carefully as it influences all our statistical analysis. After some consideration, we choose $\lambda _{2,\textrm{th}}=-6/\Omega _{\textrm{T}}$ , where $\Omega _{\textrm{T}} \equiv \tilde {u}/\lambda _{\textrm{T}}$ is a characteristic frequency. Detailed reasons are discussed in Appendix A, where we also compare the $\lambda _2$ method with an alternative vortex identification method and consider the effects of different $\lambda _2$ -thresholds on our main results. Although $\lambda _{2,\textrm{th}}$ does quantitatively affect the subsequent analysis, the main features of the statistics remain consistent regardless of $\lambda _{2,\textrm{th}}$ , allowing us to draw conclusions about the turbulent structures associated with dimples and vortices.

5.2. Detection and preliminary identification of scars

Formulating a surface-only criterion for identifying scars is not equally obvious. The concept of scars simply signifies visual patterns on the surface that, according to Brocchini & Peregrine (Reference Brocchini and Peregrine2001), ‘occur where flow on at least one side is downward causing a trough in the surface’. They described scars as cusp- or corner-like indentations, which appear as lines on the surface at the edge of upwelling boils and in regions adjacent to submerged horizontal vortices. Other studies have shown how pairs of counter-rotating vortices generate scars at the free surface: Sarpkaya & Suthon (Reference Sarpkaya and Suthon1991) demonstrated this experimentally, generating the counter-rotating vortex pair through a submerged slit (see also Dommermuth Reference Dommermuth1993; Sarpkaya Reference Sarpkaya1996). Ohring & Lugt (Reference Ohring and Lugt1991) explored the same phenomenon in simulations of laminar flow, focusing on how different Froude and Weber numbers affect the position of the scar relative to the submerged vortices.

In the present study, we observe submerged horizontal vortices at the edges of upwelling events, with a scar located directly above such a vortex (see figure 5). This is similar to the low-Froude-number cases shown by Ohring & Lugt (Reference Ohring and Lugt1991), where they saw a scar-like indentation directly above the vortex (as opposed to higher Froude numbers where the scars appear to the side of the vortex). Hence, we find it reasonable to hypothesise that scars are inextricably linked to the said vortex, so that scars, like dimples, are coupled to a specific sub-surface coherent turbulent structure which is large and long-lived. Such information can be used both to make inferences about sub-surface behaviour below the scars detected at the surface, and, reversely, to classify scars as a particular surface signature of a sub-surface vortex. To test the hypothesis, we commence by identifying scars from the surface elevation alone and consider conditional and unconditional statistics of the flow beneath.

Heuristically, we demand that the following are necessary conditions for a proper definition of a scar based on surface deformation only:

1. (i) scars are depressions in the surface;

2. (ii) scars adhere to the surface for some time, long compared with noise imprints on the surface;

3. (iii) scars are elongated structures, i.e. they have areas much smaller than the smallest circle which circumscribes them;

4. (iv) scars are advected approximately with their surrounding flow.

Point (iv) is included to distinguish scars from propagating or standing surface waves. It is not of practical interest for the data we consider here since we do not have waves present, but becomes relevant in a field setting, and is included for completeness. Surface features will not, in general, move with precisely the mean flow since vortices also propel themselves. In a spectral sense, imprints of turbulence on a mean flow will produce broadband signals approximately on the advection line of the flow, not on the wave dispersion curve (Luo et al. Reference Luo, Dolcetti, Stoesser and Tait2023; Bullee et al. Reference Bullee, Weichert, Nore, Li, Ellingsen and Hearst2024).

We identify scars at the surface as defined by the SOC proposed above by adapting the method of dimple detection using wavelets, presented by Babiker et al. (Reference Babiker, Bjerkebæk, Xuan, Shen and Ellingsen2023), to the purpose. The conditions of the SOC are met by: (i) using so-called ‘Mexican hat’ wavelets; (ii) filtering out structures with shorter lifetimes than five recorded frames in the simulations (approximately $0.06 T_\infty$ for case 1); (iii) filtering out structures whose eccentricity is lower than a threshold. As the detected structures change shape over time, we use running averages with a window size of five recorded frames to compute the eccentricity. Structures detected by the wavelets which have averaged eccentricity ${\lt } 0.85$ are discarded, and hence, we have reversed the condition that Babiker et al. (Reference Babiker, Bjerkebæk, Xuan, Shen and Ellingsen2023) used to detect the near-circular dimples. Note that to successfully separate scars from dimples and other structures by eccentricity, the eccentricity requirement must be applied to structures which cover at least three grid points. Structures of two grid points only have one possible configuration, which has eccentricity ${\gt } 0.85$ , and, hence, they will always be classified as scars (which is incorrect, as there are a large number of dimples with area of two pixels, as we will see in § 5.3). To avoid false positives, the eccentricity filtering therefore enforces scar area $A_S \geqslant 3$ , with the possibility of excluding some very small scars in the analysis. We consider this heuristically reasonable since ‘scars’ are conventionally understood to be larger than the smallest surface features. The condition of minimum area is not included in the SOC above however, as it is a practical necessity in our data analysis, not a general condition. A video illustrating the detection of dimples and scars may be found in the Supplementary Material.

Figure 5.

Snapshot of the surface and subsurface vortical structures, demonstrating the relation of vortical structures to dimples and scars at the surface: (a) surface elevation, together with (b) detected surface features, (c) sub-surface vortices and (d) both. Vortices are iso-surfaces of $\lambda _2$ in the sub-surface velocity field, coloured green and red for distinguishability.

Consider the snapshot of the surface elevation in a region where the method above identifies a scar, in figure 5(a). By visual inspection, a scar can be discerned as an elongated depression on the surface next to a near-circular dimple; these match the features detected with the method described above, as shown in figure 5(b). Applying the $\lambda _2$ criterion on the flow below the scar brings a horizontal vortex tube aligned with the scar to attention (figure 5 c). This is an example of the hypothesised ‘scar/submerged vortex’-coupling, a phenomenon we investigate statistically in § 5.4. The well-known surface-attached vortex terminating in a surface dimple is also shown. The mean orientation of these vortices is roughly parallel and orthogonal to the surface, respectively, a point we explore in detail in § 5.6.

5.3. Conditional vortex probabilities beneath dimples

We begin by considering surface-attached vortices, adopting a simple method: we numerically estimate the depth-dependent, conditional probability that a point ${{\boldsymbol r}}_{\bar {z}}=(x,y,\bar {z})$ lies inside a vortex given that the corresponding point ${\boldsymbol r}_\eta =(x,y,-\eta (x,y))$ directly above lies within a (true) dimple. We denote these conditional probabilities $f^{\textrm{V}}_{\eta }(\bar {z})=P[V({{\boldsymbol r}}_{\bar {z}})|V({\boldsymbol r}_\eta )]$ , with:

1. (i) $V({{\boldsymbol r}}_{\bar {z}})$ , the event of ${{\boldsymbol r}}_{\bar {z}}$ lying inside a vortex;

2. (ii) $V({\boldsymbol r}_\eta )$ , the event of ${\boldsymbol r}_\eta$ lying inside a dimple at the surface.

Consequently, $f^{\textrm{V}}(\bar {z})=P[V({{\boldsymbol r}}_{\bar {z}})]$ denotes the (unconditional) probability of being inside a vortex at depth $\bar {z}$ for any arbitrary horizontal position. Since we identify dimples from surface-attached vortices, $\lim _{\bar {z}\to 0}f^{\textrm{V}}_{\eta }(\bar {z})=1$ by construction (using dimples detected from the SOC gives a value slightly smaller than $1$ due to the occasional false positive identification, but otherwise changes none of the conclusions). The results for both conditional and unconditional probabilities are averaged over horizontal grid planes and over time. The depth measure $\bar {z}$ is used to account for the surface-adhering grid (see § 4.1). In computing $f^{\textrm{V}}_{\eta }(\bar {z})$ , each dimple is given the same weight regardless of its size and strength, i.e. only points directly beneath the centroid of the surface region wherein $\lambda _2\leqslant \lambda _{2,\textrm{th}}$ are considered (we come to dependency on dimple size shortly). The impact of the parameter $\lambda _{2,\textrm{th}}$ on $f^{\textrm{V}}_{\eta }(\bar {z})$ is discussed in Appendix A.2.

Figure 6 depicts the results for case 1, with the conditional vortex probability beneath dimples (solid line) plotted along with a fitted Gaussian function (defined as $(1-P_\infty )\exp [- (1/2) (\bar {z}/ L_\infty )^2/\sigma ^2]+P_\infty$ with tuneable constants $\sigma$ and $P_\infty$ ) for reference (dashed line), and unconditional probabilities (dash-dotted line). The Gaussian fit closely follows the conditional probability curve, $f^{\textrm{V}}_{\eta }(\bar {z})$ , particularly in the vicinity of the free surface. Moreover, $f^{\textrm{V}}_{\eta }(\bar {z})$ approaches the unconditional probability, $f^{\textrm{V}}(\bar {z})$ , in the outer part of the blockage layer ( $z \gtrsim L_\infty$ ) where the influence of the free surface becomes negligible. Note, however, that $f^{\textrm{V}}_{\eta }(\bar {z})$ appears to dip below $f^{\textrm{V}}(\bar {z})$ rather than directly overlap with it outside the blockage layer. We return to this when discussing the equivalent probabilities for scars below. We also note that $f^{\textrm{V}}_{\eta }(\bar {z})$ consistently follows a Gaussian fit, albeit with varying $\sigma$ and $P_\infty$ , for different $\lambda _{2,\textrm{th}}$ , as shown in figure 15(a) and discussed in Appendix A.2.

Figure 6.

Conditional probability $P[V({{\boldsymbol r}}_{\bar {z}})|V({\boldsymbol r}_\eta )]$ of being inside a vortex, given that the horizontal position lies directly beneath the centroid of the free-surface cross-section of a surface-attached vortex. Unconditional probability $P[V({{\boldsymbol r}}_{\bar {z}})]$ (black, dash-dotted line) and Gaussian fit (red, dashed line) for reference. The vertical dashed line marks the limit of the viscous boundary layer.

While the conditional probability $f^{\textrm{V}}_{\eta }(\eta )=1$ at the surface by construction, the unconditional probability $f^{\textrm{V}}(\bar {z})$ is nearly indistinguishable from zero at $\bar {z}\to 0$ in figure 6. The unconditional probability is effectively an area fraction of vortices in each grid plane. Hence, the low value of $f^{\textrm{V}}(\bar {z})$ at the surface reflects how only a small fraction of near-surface vortices attach to the surface while the majority do not penetrate the viscous layer.

The computations of conditional probabilities are repeated, considering the dependency on vortex surface area (or more precisely, dimple area, $A_D$ ) by binning surface cross-sections of vortex tubes identified by the $\lambda _2$ criterion. Results for the conditional probabilities along with binning properties are shown in figure 7.

Figure 7.

(a) Distribution of vortex count by dimple area ( $A_D$ ), the latter measured as the area of the cross-section of a surface-attached vortex at $\bar {z}=0$ and given in number of pixels and as scaled by Taylor microscale squared (a single pixel has an area of approximately $10^{-2} \lambda _T^2$ ). Blue bars represent the vortex count per area bin, black rectangles delimit the bins so that each one, except for the rightmost one, contains at least $5000$ dimple counts. (b) Curves for conditional probability $P[V({{\boldsymbol r}}_{\bar {z}})|V({\boldsymbol r}_\eta )]$ . Each curve represents one bin from panel (a) coloured by increased dimple area. The vertical dashed line marks the limit of the viscous boundary layer. The dash-dotted black line is the unconditional probability, $P[V({{\boldsymbol r}}_{\bar {z}})]$ . (c) Variance ( $\sigma ^2$ ) of approximated Gaussian curves of conditional probability data in panel (b) by the weighted average of bin area (blue x-markers) and linear fit using $l_1$ -norm (dashed red line). Horizontal error bars indicate the range of dimple area covered by bins that span multiple sizes.

In the computation of size distribution and binning, vortex sizes are computed and counted for each sampled time step, meaning that vortices that persist over multiple snapshots are counted once for each snapshot. Although no weighting is performed with respect to dimple size or depth, larger vortices tend to have longer lifetimes and thus produce more counts. Also, one and the same vortex can be counted into different bins, as its surface cross-sectional area changes over time. The distribution of vortex count by size (blue bars in figure 7 a) is roughly exponential for vortices smaller than approximately $80$ pixels and levels off for the infrequent larger vortices. This could be an indication that a critical size exists above which the trend changes, although the comparatively small number of large vortices makes such conclusions tentative.

To get statistically significant results when conditional probabilities are computed for each size span of vortices, we require each bin to contain at least $5000$ vortex counts, i.e. for each plotted curve of conditional probabilities in figure 7(b), the results are averaged over a minimum of 5000 instances (much more for the frequently occurring smaller vortices, where single dimple sizes are counted more than $10^4$ times, as observed for the dimples with $A_D \lesssim 10$ in figure 7 a). Where bins contain more than a single area, this is indicated with a black box in figure 7(a). This results in a total of $25$ bins and ensures that bin averages are reasonably well converged. The rightmost bin (spanning vortices with area 93– $178$ pixels) counts only 485 vortices and is therefore not included in figures 7(b) and 7(c).

Figure 7(b) depicts the vortex size-dependent conditional probabilities as a function of depth $\bar {z}$ . Considering single curves, we notice that all the probability curves are near-Gaussian, like the probabilities averaged over all vortex sizes in figure 6, while remarkably, the variance of the Gaussian fits, plotted in figure 7(c), increases approximately linearly with respect to the vortex size (except for the outlier for the bin that covers a broad span of the intermediate and large vortices). All vortices greater than approximately $6$ pixels extend well beyond the viscous layer (probability $\gtrsim 95\,\%$ ), and even among the single-pixel vortices, $66\,\%$ remain present directly beneath the (tiny) dimple in the limit of the viscous layer. Larger vortices penetrate much further, with the largest extending vertically into the lower part of the blockage layer. Nevertheless, all size-dependent conditional probability curves are below 0.42 when measured at depth $\bar {z} \geqslant 0.5 L_\infty$ , and fall off to the average volume fraction of vortices within the blockage layer. This result indicates that the vortical structures, whose imprints manifest as dimples at the surface, are either deflected by the surrounding turbulence, bent upwards towards the surface again, or weakened to such an extent that they are no longer counted as a vortex by the $\lambda _2$ criterion before they reach the limit of the surface-influenced region.

5.4. Conditional vortex probabilities beneath scars

We next compute probabilities as described in § 5.3, but now conditional of ${\boldsymbol r}_\eta$ lying inside an area of the surface identified as a scar with the SOC implementation in § 5.2. We calculate the conditional probability $f^{\textrm{S}}_{\eta }(\bar {z})=P[V({{\boldsymbol r}}_{\bar {z}})|S({\boldsymbol r}_\eta )]$ , where $S({\boldsymbol r}_\eta )$ is the event of ${\boldsymbol r}_\eta$ being inside a scar at the surface. Note that a slight modification is necessary to compute the conditional probabilities of the scars as compared with those of vortices: since scars are non-circular, using their centroids as measurement points may entail considering points lying outside the scars themselves for the computations of conditional probabilities; this happens for approximately $6\,\%$ of all scars for case 1 when centroids are used without modification. A more fitting measurement point for a scar is a combination of the centroid and its approximated centreline. The latter is computed by skeletonisation. This process reduces a two-dimensional pattern to a centreline of single-pixel thickness while preserving its topology (we use the function bwskel in MATLAB’s Image Processing Toolbox (release 2022b, Mathworks, USA) for this purpose. After skeletonisation, the measurement point is set to the point on the centreline that has the shortest Euclidean distance to the centroid. This final step serves to avoid overemphasising long, thin tails of scars, as may happen if the midpoint of the centreline is selected directly.

Figure 8 depicts the conditional vortex probabilities beneath scars for case 1. One immediately observes that the conditional probability for vortices beneath scars increases rapidly through the viscous layer, with maximum probability at the limit of this layer. In other words, virtually all structures that were identified as scars from the set of SOC have a vortex directly beneath them situated approximately at $z=L_\nu$ . Correspondingly, by considering the cumulative probability, i.e. the probability that there is at least one point within a vortex along the vertical line from the surface down to a level $\bar {z}$ , we see that it has reached approximately $94\,\%$ at $\bar {z}={L_\nu }$ . Starting from an arbitrary point, the same probability is only 3.3 %. On this evidence, we can safely conclude that a scar, as identified with the set of SOC in § 5.2, is a signature of a vortex situated around $\bar {z}\approx {L_\nu }$ . Note that this does not imply that we should expect $f^{\textrm{S}}_{\eta }(\bar {z})$ to equal 1 at $z = L_\nu$ , since our scar identification implementation includes a non-zero number of false positives and since the amplitude of $f^{\textrm{S}}_{\eta }(\bar {z})$ is affected by the choice of $\lambda _{2,\textrm{th}}$ . Yet for all tested $\lambda _{2,\textrm{th}}$ values, $f^{\textrm{S}}_{\eta }(\bar {z})$ reaches its peak at $z={L_\nu }$ (figure 15), indicating that the cores of scar-related vortices are indeed situated at the boundary of the viscous layer (further details in Appendix A.2).

Figure 8.

Conditional probability $P[V({{\boldsymbol r}}_{\bar {z}})|S({\boldsymbol r}_\eta )]$ of ${{\boldsymbol r}}_{\bar {z}}$ being inside a vortex, given that there is a scar on the surface directly above (blue, solid line) the corresponding cumulative probability (red, dashed line), and the unconditional probability $P[V({{\boldsymbol r}}_{\bar {z}})]$ (black, dash-dotted line), for reference. The vertical dashed line marks the limit of the viscous boundary layer.

Another observation from figure 8 is that unlike $f^{\textrm{V}}_{\eta }(\bar {z})$ , $f^{\textrm{S}}_{\eta }(\bar {z})$ does not approach $f^{\textrm{V}}(\bar {z})$ inside the depth range plotted here. This is not surprising given the intermittency of the turbulence and the distribution of scars: scars are numerous only during strong turbulent (upwelling) events (see supplementary material available at https://doi.org/10.1017/jfm.2025.072 for videos depicting the evolution of scars and dimples at the surface). Thus, the conditional statistics for scars are biased towards times with increased vortex concentration in the upper bulk. The amount of dimples present in the flow is also connected to upwelling events. Yet, unlike scars, which spatially and temporally appear in the immediate vicinity of the upwelling structures, dimples appear as upwelling boils, mature and die out (with an associated lowering of vortex concentration), and are subsequently spread out and advected away from the upwelling cores. In fact, a weak opposite effect on conditional probabilities for dimples can be seen in figure 6: $f^{\textrm{V}}_{\eta }(\bar {z})$ approaches $f^{\textrm{V}}(\bar {z})$ yet lies slightly below it in the bulk. These trends for $f^{\textrm{V}}_{\eta }(\bar {z})$ and $f^{\textrm{S}}_{\eta }(\bar {z})$ are even more prominent for the low-Reynolds-number cases (see § 5.5), where lower turbulence intensity overall increases the spatio-temporal dependencies of the instantaneous conditional probabilities.

Figure 9.

(a) Distribution of scar count by scar area ( $A_{\textrm{S}}$ ), where the latter is the surface area covered by a detected scar in number of pixels and as scaled by Taylor microscale. Blue bars represent the scar count for each pixel size. Black rectangles denote the bins used for the computation of conditional probabilities. (b) Curves for conditional probability $P[V({{\boldsymbol r}}_{\bar {z}})|S({\boldsymbol r}_\eta )]$ of being inside a vortex, given that there is a scar at the surface at the same $x,y$ -position, sorted by scar area so that each curve represents one bin in panel (a). The vertical dashed line marks the limit of the viscous boundary layer. The dash-dotted black line is the unconditional probability, $P[V({{\boldsymbol r}}_{\bar {z}})]$ . (c) Corresponding cumulative probabilities, by scar area. Inset shows zoom-in on the region where curves cross from the viscous layer to the blockage layer.

By the same binning and averaging that we presented for dimples above, results for conditional and cumulative vortex probabilities beneath scars are computed for different scar sizes and shown in figure 9, along with binning details. Since the scars in our simulations are less numerous and longer lasting (on average) than the dimples, we require more instances in each bin for scars to avoid the danger of single bins containing only a small amount of long-lasting scars (recall that each count denotes an imprint at a single time step). Therefore, each bin contains at least $10\,000$ counts to ensure statistical convergence, except for the one containing the largest scars, which does not reach the minimum number and is excluded from the analysis.

The difference between figures 7 and 9 is striking. Unlike the probability of finding vortices a finite distance beneath a dimple which depends strongly on dimple area, the curves in figure 9(b,c) collapse near perfectly and show no dependence on scar size. The results very clearly support the conclusion that below every scar, there is a vortex situated near the lower edge of the viscous layer, $\bar {z}\approx {L_\nu }$ , regardless of the size of the scar. Towards the formulation of a PCD, we conclude that a true scar is an imprint of a sub-surface vortex which is not connected to the surface. Further statistical insight into these sub-scar vortices comes to light in § 5.6.

5.5. Reynolds number and Weber number effects on sub-feature statistics

Figure 10.

Conditional probability for dimples, $P[V({{\boldsymbol r}}_{\bar {z}})|V({\boldsymbol r}_\eta )]$ (blue, solid line) and scars, $P[V({{\boldsymbol r}}_{\bar {z}})|S({\boldsymbol r}_\eta )]$ (red, solid line), for cases 1–6 in panels (a)–(f), respectively. In each panel, the vertical dashed line marks the limit of the viscous boundary layer and the dash-dotted black line is the unconditional probability, $P[V({{\boldsymbol r}}_{\bar {z}})]$ .

Consider figure 10, depicting the conditional probabilities $f^{\textrm{V}}_{\eta }(\bar {z})$ and $f^{\textrm{S}}_{\eta }(\bar {z})$ for all six cases detailed in table 1. As for case 1, also for cases 2–6, $f^{\textrm{V}}_{\eta }(\bar {z})$ falls off approximately as a Gaussian through the blockage layer, whereas $f^{\textrm{S}}_{\eta }(\bar {z})$ increases rapidly in the vicinity of the surface, with a maximum near $\bar {z} = L_\nu$ .

Comparing the results for different Weber numbers reveals a negligible effect of surface tension on the conditional probabilities (with the caveat that we have considered only three values of $We$ ). The Reynolds number, however, does affect the statistical depth of the sub-surface vortices for both dimples and scars.

A different scaling of the abscissa points in the direction of what drives the Reynolds number effect. Observe in figure 11, where the conditional probabilities for dimples, $f^{\textrm{V}}_{\eta }(\bar {z})$ , and scars, $f^{\textrm{S}}_{\eta }(\bar {z})$ , are plotted as a function of depth scaled by $L_\nu$ (or $\lambda _{\textrm{T}}$ , discussed below) rather than $L_\infty$ , that all curves for $f^{\textrm{V}}_{\eta }(\bar {z}; Re,We) = P[V({{\boldsymbol r}}_{\bar {z}})|V({\boldsymbol r}_\eta ),{Re},{We}]$ collapse onto one curve $f^{\textrm{V}}_{\eta }(\bar {z}/L_\nu ) = P[V({{\boldsymbol r}}_{\bar {z}}/L_\nu )|V({\boldsymbol r}_\eta )]$ in the vicinity of the surface. The same effect is evident for scars, by the peak in $f^{\textrm{S}}_{\eta }(\bar {z})$ around $\bar {z} = {L_\nu }$ observed for all Reynolds and Weber numbers. Hence, $f^{\textrm{S}}_{\eta }(\bar {z}; Re, We)$ collapses onto $f^{\textrm{S}}_{\eta }(\bar {z}/L_\nu )$ , which peaks at $\bar {z}/L_\nu \approx 1$ for all cases. Further from the surface, some discrepancy between the results for $Re= 2500$ and $Re= 1000$ can be observed. This is due to the rapid decay of the turbulence rising from the deeper bulk for the lower Reynolds number cases.

Figure 11.

Conditional probabilities for (a) dimples, $P[V({{\boldsymbol r}}_{\bar {z}})|V({\boldsymbol r}_\eta )]$ and (b) scars, $P[V({{\boldsymbol r}}_{\bar {z}})|S({\boldsymbol r}_\eta )]$ for different Reynolds numbers and Weber numbers. Scaling by Taylor microscale ( $\lambda _{\textrm{T}}$ ) or viscous layer thickness ( $L_\nu$ ).

Recalling that the thickness of the viscous layer is computed by $L_\nu = 2 Re_\infty ^{-1/2} L_\infty$ and inserting the definitions of $L_\infty$ and $Re_\infty$ (see § 3.2) yields:

(5.1) \begin{align} L_\nu & = \frac {\lambda _{\textrm{T}}}{\sqrt {15}}\, \approx 0.26\,\lambda _{\textrm{T}}, \end{align}

i.e. a purely Taylor microscale-dependent number. Hence, our study supports the conclusion that the vertical reach of a surface-attached vortex beneath a dimple and the depth of the vortical structure below scars depends on the Taylor microscale of the turbulence in the bulk only, independent of surface tension. In terms of conditional probability, this may be expressed as

(5.2a) \begin{align} f^{\textrm{V}}_{\eta }(\bar {z}, Re,We) \to & f^{\textrm{V}}_{\eta }(\bar {z}/\lambda _{\textrm{T}}) , \end{align}

(5.2b) \begin{align} f^{\textrm{S}}_{\eta }(\bar {z}, Re,We) \to & f^{\textrm{S}}_{\eta }(\bar {z}/\lambda _{\textrm{T}}) , \end{align}

with max $\{f^{\textrm{S}}_{\eta }(\bar {z}/\lambda _{\textrm{T}})\}$ at $\bar {z} \approx 0.26 \lambda _{\textrm{T}}$ . Note, however, that this conclusion is shown for a small sample of $Re$ to hold within a limited range of Reynolds and Weber numbers. Further exploration is called for to establish the limits of applicability.

5.6. Inclination angles for vortices below scars and dimples

The conditional probabilities for finding vortices directly beneath dimples and scars give insight into the prevalence of vortices, but not their properties. Since turbulence near the surface is highly anisotropic, a key property of a vortex is its orientation. We therefore now consider the inclination angles of vortices below dimples and scars.

Vortex inclination angles $\theta$ are computed from (4.1), now considering only grid points that are covered (from above) by a scar or a dimple, i.e. from components of the vorticity field ${\boldsymbol \omega }(\hat {x},\hat {y},\bar {z},t)$ , where $\hat {x}, \hat {y}$ denote the time-dependent subset of $x,y$ -coordinates for which $V({\boldsymbol r}_\eta ) = 1$ for dimples and $S({\boldsymbol r}_\eta ) = 1$ for scars at time $t$ .

Figures 12 and 13 depict the resulting weighted inclination angle distributions at different depth levels beneath dimples and scars, respectively, for case 1 (similar figures for the other cases can be found in the supplementary material). The most obvious conclusion from the figures is – as predicted above – that the vorticity is predominantly vertical beneath dimples and predominantly horizontal beneath scars.

Figure 12.

Normalised histograms of weighted vortex inclination angles for case 1, at increasing average depths spanning the surface and blockage layers, including only regions below vortex dimples. All depths are given in multiples of the viscous layer thickness.

Figure 13.

The same as figure 12, but now for regions beneath scars. The added red line represents results for scars that do not overlap with dimples.

The trend from the results of the inclination angles is clear when compared with the inclination angles of the vorticity vector in the flow in general (figure 4). Recall that we observed a transition from one preferred vortex inclination angle (vertical; $\theta = 0$ and $\pi$ ) to another (horizontal; $\theta = \pi /2$ ) when starting at the free surface and moving through the viscous layer. More precisely, the unconditional distribution of vorticity inclination shows that inside the viscous layer, the vorticity field is either near-horizontal or near-vertical with only a small fraction of points where $|{\boldsymbol \omega }|$ is significant having intermediate orientations. At the surface, all vorticity is vertical (or, more properly, surface-orthogonal) due to the viscous effect, gradually shifting to being dominated by horizontal orientation for $\bar {z}\sim L_\nu$ . The vorticity field remains anisotropic throughout the blockage layer, with vertical and horizontal orientations being more probable.

As observed in figures 12 and 13, the picture is very different below dimples and scars, respectively. The vorticity at the surface for areas classified as ‘true’ dimples is surface normal by definition, and we observe that the preference for horizontal orientation is entirely absent. The statistical indication is that surface-attached vortices extend directly downwards, typically penetrating deep into the blockage layer before curving away, and deflect away any horizontal vortices in the vicinity from entering the space beneath them.

A consideration of vorticity inclination below scars reveals the opposite trend. There is a strong bias towards horizontal orientation throughout the viscous layer, and horizontal vortex inclination remains the preferred inclination well into the blockage layer. In contrast to the inclination angle results for the flow in general, shown in figure 4, and for the flow below dimples, we notice that the vertical component of the vorticity is very small directly below the surface and all but extinguished below scars as early as at depth $\bar {z} / L_\nu = 0.3$ .

The blue bars in figure 13 display a significant peak around the vertical orientations for the two topmost depths – at and very close beneath the surface. This initially puzzling fact is readily explained from imperfections in our scar detection. Figure 5 shows an example situation where a dimple lies close to a scar. Occasionally, an area identified as a scar also contains a dimple. Moreover, although there is a high peak in the blue bars near $0$ and $\pi$ for the two top depths, the vorticity itself is very small here so the impact of a very small number of dimples among the scars will be high. Possible scenarios are vortex rings and vortex pairs connection to the surface (e.g. Bernal & Kwon Reference Bernal and Kwon1989; Dommermuth Reference Dommermuth1993; Weigand & Gharib Reference Weigand and Gharib1995; Gharib & Weigand Reference Gharib and Weigand1996; Sarpkaya Reference Sarpkaya1996; Willert & Gharib Reference Willert and Gharib1997; Zhang et al. Reference Zhang, Shen and Yue1999; Nagaosa & Handler Reference Nagaosa and Handler2003; Terrington et al. Reference Terrington, Hourigan and Thompson2022). During the connection process, the vertical vorticity is enhanced, which can generate dimples that overlap with scars.

After excluding the approximately $7\,\%$ of scars (identified with the SOC) that contain a true dimple (identified with the $\lambda _2$ criterion), the red curves in figure 13 results. (A more sophisticated implementation of scar identification from the set of SOC would have mitigated this problem without reference to sub-surface flow – we allow ourselves this pragmatic fix to demonstrate the physics involved.) Except at the very surface itself, which is ambiguous, beneath the ‘corrected scars’, the vertical vorticity preference is completely extinguished throughout the blockage layer, until the turbulence eventually becomes isotropic. Hence, we conclude towards the formulation of a PCD that the sub-surface vortex which imprints the surface as a scar is horizontally aligned.

Furthermore, considering the orientation of vortices below scars in light of the results in § 5.4 points to a cross-sectional length of $\mathcal {O}(L_\nu )$ of vortices located at the lower edge of the viscous layer below scars, regardless of the scar size. The longitudinal lengths of these vortices are, however, not clearly correlated to the scars: the surface imprint is very sensitive to the exact distance between vortex core and surface, and vortices winding slightly up and down in a serpentine manner may manifest on the surface as multiple scars when approaching the viscous boundary layer from below in several places along its length. Inversely, scars are often far shorter in length than their concomitant horizontal vortices.