Vortex structures under dimples and scars in turbulent free-surface flows

Abstract

Turbulence beneath a free surface leaves characteristic long-lived signatures on the surface, such as upwelling ‘boils’, near-circular ‘dimples’ and elongated ‘scars’, easily identifiable by eye, e.g. in riverine flows. In this paper, we analyse data from direct numerical simulations to explore the connection between these surface signatures and the underlying vortical structures. We investigate dimples, known to be imprints of surface-attached vortices, and scars, which have yet to be extensively studied, by analysing the conditional probabilities that a point beneath a signature is within a vortex core as well as the inclination angles of sub-signature vorticity. The analysis shows that the probability of vortex presence beneath a dimple decreases from the surface down through the viscous and blockage layers. This vertical variation in probability is approximately a Gaussian function of depth and depends on the dimple’s size and the bulk turbulence properties. Conversely, the probability of finding a vortex beneath a scar increases sharply from the surface to a peak at the edge of the viscous layer, regardless of scar size. The probability distributions of the angle between the vorticity vector and the vertical axis also show a clear pattern about vortex orientation: a strong preference for vertical alignment below dimples and an equally strong preference for horizontal alignment below scars. Our findings corroborate previous studies that tie dimples to surface-attached vertical vortices. Moreover, they suggest that scars can be defined as imprints of horizontal vortices that are located approximately a quarter of the Taylor microscale beneath the free surface.

1. Introduction

Turbulence near air–water surfaces plays a key role in a multitude of environmental, biological and engineering processes, e.g. in the fluxes of heat and gas across the air–water interface in rivers and oceans. Due to the difficulty of bulk flow measurements in the field, relating features on the water surface to the sub-surface flow structures that cause them is of great practical importance (Muraro et al. 2021). The region directly beneath the air–water interface, which affects the gas transfer of the low diffusivity gasses that make up the majority of the atmosphere (Jähne & Haußecker 1998; Herlina & Wissink 2014), is of particular interest. A more detailed coupling between identifiable patterns on the surfaces of lakes, rivers and oceans to the near-surface flow dynamics facilitates the use of data collected from inexpensive remote sensing technology as input, e.g. for estimating greenhouse gas discharge from rivers, similar in magnitude to the total gas flux across the ocean surface and to the total forest carbon uptake rate (see, e.g. Brinkerhoff et al. 2022).

In the absence of wind, an air–water interface is well described as a free surface, i.e. a deformable surface with constant pressure and zero shear stress imposed from the air side. The free surface of a turbulent flow exhibits a multitude of deformations that manifest as identifiable surface patterns. Among them, we typically find upwellings, dimples, waves, scars, bubbles, break up of the surface and droplet formations; a useful taxonomy was given by Brocchini & Peregrine (2001), and reviewed and applied by Muraro et al. (2021). The present study focuses on continuous free surfaces without overturning or breakup, where gravity is dominant and surface tension is weak, the most common state in terrestrial bodies of water with flow (see Brocchini & Peregrine 2001 § 7). Here, the turbulence has enough energy to disrupt the free surface, without breaking it up, and the persistent surface deformations are mainly caused by intermittent upwelling and downwelling events, that is, by rushes of fluid towards and away from the surface (see figure 1). At the surface, the upwelling and downwelling events give rise to regions where the surface divergence ( $\beta = \partial _x u + \partial _y v = -\partial _z w$ , when the surface is nearly flat and normal to $z$ ) is strongly positive and strongly negative, respectively. Areas of upwelling take the shape of large, roughly circular ‘boils’, whereas areas of strong downwelling appear around the curves circumscribing the upwellings. Efforts to properly characterise upwelling boils go back a long time, notably by Matthes (1947) (an introduction is Nezu & Nakagawa 1993). Experiments (Rashidi 1997; Kumar et al. 1998; Mandel et al. 2019; Gakhar et al. 2022) and simulations (Nagaosa 1999; Nagaosa & Handler 2003; Khakpour et al. 2012) show how these boils can be caused by turbulent structures ejected from a wall or canopy shear layer, which then evolve and rise towards the surface. Moreover, upwelling and downwelling events are not limited to flow configurations involving a shear layer that generates coherent turbulent structures. Such events have also been observed when isotropic turbulence approaches a free surface from below (Herlina & Jirka 2008; Guo & Shen 2010; Herlina & Wissink 2014, 2019). These observations suggest that, despite differences in turbulence sources, the interactions between near-surface turbulent coherent structures and a free surface can share common dynamics. Therefore, although isotropic turbulence is an idealised turbulence state, it is fundamental to understanding free-surface turbulence dynamics, such as mixing and gas transfer (see, e.g. Thompson & Turner 1975; Hopfinger & Toly 1976; Brumley & Jirka 1987; Herlina & Jirka 2008; Variano & Cowen 2013) and free-surface deformations (see, e.g. Guo & Shen 2010; Babiker et al. 2023). These studies are particularly relevant to the near-surface turbulence boundary layer where there is zero mean shear, such as the water under calm wind, and where the turbulence is generated far from the surface. Moreover, the isotropic turbulence may be considered a common basis for various types of turbulence sources in nature, as turbulence tends to be reasonably isotropic at small scales.

Figure 1. Snapshot of the river Nidelva in Trondheim, exhibiting a multitude of surface deformations. A small selection of dimples (white ovals), scars (red ovals), upwelling boils (orange boxes) and waves (yellow boxes) are marked. The largest dimples (green boxes) are von Kármán vortices shed from a nearby bridge pillar. Photo by Klervie le Bris.

Dimples – persisting near-circular indentations at the free surface – have long been observed at the edges of upwelling boils (see Banerjee 1994; Longuet-Higgins 1996). These dimples emerge when sub-surface vortices break up and attach to the free surface, as observed in studies on vortex ring connection to free surfaces, both experimentally (e.g. Bernal & Kwon 1989; Weigand & Gharib 1995; Gharib & Weigand 1996; Willert & Gharib 1997) and in simulations (Zhang et al. 1999; Terrington et al. 2022). Also appearing on the edges of upwelling regions are the scars, a type of surface pattern that has attracted less attention in the literature. These elongated, valley-like patterns on the free surface mark regions where the fluid is accelerated from a region of strong positive divergence (upwelling) to a region of strong negative divergence (downwelling), causing a local pressure deficit and a concomitant depression of the free surface. Laminar studies of a rising pair of large horizontal vortices have shown how scars develop over the vortex cores with an elevated area in between (Sarpkaya & Suthon 1991; Dommermuth 1993). Herlina & Wissink (2019) observed that in turbulent free-surface flows, slender horizontal vortices are often observed beneath the surface adjacent to ‘downwelling curves’, shifted towards the side of the nearest strong upwelling; since this is the area where scars are expected, it is an indication that scars and sub-surface vortices are connected, a notion we investigate and quantify here.

Dimples and scars share several practically important traits: they are both local, persistent surface depressions imprinted by strong coherent turbulent structures beneath, they are readily detectable on the free surface by eye or computer vision, and both are closely connected to upwelling and downwelling events which are the main drivers of surface renewal facilitating surface fluxes (Kermani & Shen 2009). This motivates the main focus of the present paper: to increase the understanding of dimple and scar patterns by quantifying their sub-surface patterns in turbulent free-surface flows. We consider the ‘inverse’ problem of inferring probable subsurface properties from surface observations; uniting our findings with those from the ‘forward’ approach of inferring the surface imprint of a vortex (say) with given properties is an avenue worth pursuing in a separate study.

The interest in quantifying the turbulence near the free surface for the estimation of heat and gas transfer has inspired new techniques for estimating turbulence properties close to the free surface. In situ measurements that penetrate the surface are routinely employed, and can be used to quantify the free surface and to make statistical inferences about surface properties or the turbulence below. However, these can only provide measurement data at single points or along trajectories, and upscaling to cover many sites over large areas is forbiddingly costly. An attractive alternative is non-contact measurements from above which use only information from the free surface itself; see the the recent overview by Dolcetti et al. (2022). It requires knowledge of what information about the fluid flow the observed surface features can provide. A step in this direction is taken by Babiker et al. (2023), who detect and track the dimples at the free surface by applying wavelet analysis to surface elevation data (from direct numerical simulations). The detected dimples show a strong correlation to variations in the surface divergence, $\beta$ , and may therefore act as a proxy for measuring the effects of upwelling events on surface renewal (Babiker et al. 2023) – a beneficial connection as the surface divergence is linked to the mass transfer velocities of low solubility gases and is often used to model cross-surface mass flux (see, e.g. Banerjee et al. 2004; Turney & Banerjee 2013; Kermani et al. 2011). A further step of statistical inferences from surface observations is to use machine learning. Xuan & Shen (2023) explored how neural networks may be used to reconstruct the entire flow field below a free surface. Such data-driven methods have the potential to greatly improve remote sensing campaigns by reconstructing the turbulence or turbulent statistics with limited measurements, but are still in their infancy for this application.

Of greatest practical interest is the region henceforth denoted the surface-influenced layer, the region below the free surface where the flow is significantly influenced by the free surface, i.e. where turbulent mixing and surface renewal take place. Hunt & Graham (1978) pioneered a highly fruitful approach to conceptualise the surface-influenced layer. They applied a rapid-distortion model to analyse a turbulent boundary layer atop a wall moving with the mean flow velocity and found two regions in the surface-influenced layer dominated by different mechanisms: an inner, viscous region and an outer, blockage (or source) region. Although there are important differences in boundary conditions between a rigid surface and a deformable, free-slip interface, Hunt & Graham’s subdivision of the boundary layer conceptually applies also in free-surface turbulence. Hunt & Graham’s model has been expanded upon in later studies (Hunt 1984; Teixeira & Belcher 2000; Magnaudet 2003) and widely discussed in the literature on the surface-influenced layer (e.g. Brumley & Jirka 1987; Shen et al. 1999, 2000; Nagaosa & Handler 2003; Flores et al. 2017). In short, the two sub-regions of the surface-influenced layer have been tied to the two boundary conditions at the surface: the inner region, the viscous layer, to the dynamic boundary condition, and the outer region, the blockage layer, to the kinematic boundary condition. We return to this in § 3.2. For now, we are content with stating that it is particularly advantageous to link surface observations to the different sub-regions of the surface-influenced layer, as this yields a measure of the relevant length scales at which the effects act and what physical mechanisms can be tied to the observations.

The present study aims to elucidate how and to what extent dimples and scars are coupled to the sub-surface dynamics of the flow. When exploring the nature of scars, we immediately run into the fundamental question of what a scar really is. Although a number of authors offer qualitative descriptions, at present, we know of no precise set of criteria. To place the analysis on more solid ground and ensure that any circular arguments are avoided, we propose a new conceptual framework, detailed in § 2. Dimples can be readily placed in our framework on the basis of known properties, and a main goal of this study is to establish a formal understanding of scars using the same concepts, thus making the two features directly comparable.

1.1. Procedure and outline

We explore how dimples and scars are linked to sub-surface structures by considering the conditional probability that the points directly beneath dimples and scars lie inside a vortex core. To properly analyse the relation between surface and sub-surface structures, we introduce a conceptual framework coupling a ‘pseudo-causal definition’ (PCD) of a surface feature in terms of the flow structure which ‘creates’ it and a ‘surface-only criterion’ (SOC) which a surface feature must satisfy to be classified a ‘dimple’ or a ‘scar’. Next, we consider the statistics of how the vorticity beneath these surface features is inclined, tying vortex tubes below dimples and scars to two different preferred directions for the vortex inclination.

The dimples in our analysis are identified by building upon the knowledge of their PCD–SOC pair. As the equivalent PCD–SOC pair is unknown for scars, we proceed in an explorative manner in scar detection, by suggesting what amounts to a candidate SOC for scars and then studying statistically the flow field beneath the resulting structures. Based on the statistical results, we are able to identify and argue for a successful formal classification of scars in the form of a PCD–SOC pair. The data analysed stem from direct numerical simulation (DNS) for a selection of Reynolds and Weber numbers.

The paper is organised as follows: in § 2, we introduce the conceptual framework for surface-feature identification; in § 3, we provide details on the mathematical formulation of our flow and outline key aspects of the DNS; and in § 4, the DNS results are presented and analysed. We present our main analysis of vortices and scars in the flow in § 5 and draw conclusions in § 6.

2. Concepts for surface-feature identification

We propose two concepts pertaining to a class of surface features: the pseudo-causal definition and a set of surface-only criteria. When a well-matched pair of these has been identified, it prepares the ground for matching surface imprint observations to sub-surface flow phenomena in a quantitative way. The classes of surface features considered here are dimples and scars. For dimples, both a pseudo-causal definition and a set of surface-only criteria can be formulated based on current knowledge, while it is one of our goals herein to propose these for the case of scars.

A PCD of a class of surface features states that a ‘true’ member of said class is the surface imprint of a specified sub-surface structure. It should be formulated such that an observed feature is unambiguously sorted into at most one class. For the case of dimples, the appropriate PCD is that a ‘true’ dimple is the imprint of a surface-attached vortex and appears where said vortex is attached. For practical use, for instance, in DNS data, a quantitative formulation is required, for example, by identifying a vortex using the $\lambda _2$ criterion which we will discuss in detail later. Strictly speaking, the causal narrative that one typically uses is incorrect. Saying, for example, that a dimple is ‘caused by’ the presence of a surface-attached vortex somewhat misrepresents the fact that the two appear together, simultaneously, and are two expressions of one and the same flow event (note that, e.g. propagating surface waves are different in this respect). Hence, the label ‘pseudo-causal’. The PCD is of limited practical use in a field setting because detailed knowledge of the flow is necessary for its quantification. However, it expresses a physical understanding of the flow picture. In practice, one can hardly achieve the ideal level of universality and user independence: for instance, the $\lambda _2$ criterion for vortex identification requires a user-defined threshold value, and a given surface feature is the result of the full nearby flow field, not merely the single most dominant structure. The very concept of dividing (material or Eulerian) sub-volumes of a continuous flow into ‘vortex’ and ‘not vortex’ is inherently ambiguous and non-unique. A level of pragmatism is inescapable also in formulating a PCD, and such is the nature of the problem at hand.

A set of SOC of a class of surface features refers only to the free surface using no information about sub-surface flow. It is heuristic, pragmatic and might be context-dependent, yet should be formulated so that it can also be quantitatively applied to decide whether or not a point on the free surface is part of a surface feature of a particular type, in our case, a ‘dimple’ or a ‘scar’. When speaking of dimples and scars in the following, we refer to the SOC, i.e. only based on surface appearance, in contrast to ‘true dimples’ and ‘true scars’ as decided by the PCD. A suitable SOC for dimples was demonstrated by Babiker et al. (2023) (without calling it such): a dimple is a surface indentation which persists for a minimum time and is sufficiently circular. The SOC should be usable in practical field settings where only the free surface is observed, but is purely descriptive and expresses no understanding of the underlying flow.

Once a successful PCD–SOC pair has been identified, it paves the way for inferring quantitative information about sub-surface flow from surface observations in a structured way, for instance, via statistical investigations, such as we present herein. A successful PCD–SOC pair must have a high level of agreement. For dimples, Babiker et al. (2023) showed that the mentioned PCD and SOC when applied to DNS data produced the same identification for practical purposes. Since the dynamics of surface-attached vortices (which manifest as dimples on the free surface) is comparatively well understood, the success of this PCD–SOC pair is not particularly surprising.

3. Flow problem and DNS

The flow problem analysed is isotropic turbulence interacting with a free surface. While the simulation set-up follows Guo & Shen (2010), the present study considers different Reynolds numbers and Weber numbers to gain insights into the Reynolds number scaling and the effects of surface tension. Details on the flow problem, simulation parameters and grid resolution criteria based on grid independence studies can be found from Guo & Shen (2009, 2010) and Xuan & Shen (2019, 2022). Only a synopsis is given here.

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 1975; Hopfinger & Toly 1976; Brumley & Jirka 1987; McKenna & McGillis 2004; Herlina & Jirka 2008; Chiapponi et al. 2012; Lacassagne et al. 2017), and more recent experiments, where the turbulence is generated by jets with zero net flow (e.g. Variano & Cowen 2013; Jamin et al. 2024). In our simulations, turbulence is generated by a random linear forcing method (Lundgren 2003; Rosales & Meneveau 2005), 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 2005), 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 (2009). 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 2010; Xuan & Shen 2022). 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 2019). 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 2019).

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 1972, § 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 1972, 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 1972, 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 (2005), 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. 1999). 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 (1987), the viscous layer thickness is conventionally taken to be $L_\nu = 2 \textrm{Re}_\infty ^{-1/2} L_\infty$ (see discussion by Calmet & Magnaudet 2003). No single definition for $L_\beta$ has become standard, yet ${L_\beta }=L_\infty$ is often used for convenience (see e.g. Magnaudet 2003; Guo & Shen 2010; Herlina & Wissink 2014). 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. (1999).

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. Flow in the near-surface region

Before discussing surface manifestations, we discuss some general properties of turbulent flow beneath a free surface, particularly the appearance of the blockage and viscous layers. While in a large part known, the concepts and general properties in § 4 provide the backdrop against which further analysis can be compared.

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 2010).

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. 1999). 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 1978; Perot & Moin 1995; Shen et al. 1999), 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 2003).

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. (1999) 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 (1985), 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. (1999) 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. Dimples, scars and the vortex dynamics in the surface-influenced layer

The results in § 4 give a general overview of sub-surface flow properties for the turbulent flow. Our principal goal is to understand the dynamics of the flow beneath dimples and scars, including (but not limited to) the identification of a SOC–PCD pair for scars. To this end, we investigate statistical properties of the flow conditional on being beneath scars as identified from a heuristically proposed SOC. We begin, however, by performing the same procedure for the far better understood case of dimples, whence interesting insights also emerge.

In general discussions of scars and dimples in § 5.3 and § 5.4, we consider only case 1, for which trends are clearly seen. In § 5.5, we compare these observations to the other cases in table 1.

5.1. Identification of dimples and vortices

It was shown by Babiker et al. (2023) 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 (2001), ‘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 (1991) demonstrated this experimentally, generating the counter-rotating vortex pair through a submerged slit (see also Dommermuth 1993; Sarpkaya 1996). Ohring & Lugt (1991) 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 (1991), 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. 2023; Bullee et al. 2024).

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. (2023), 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. (2023) 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 1989; Dommermuth 1993; Weigand & Gharib 1995; Gharib & Weigand 1996; Sarpkaya 1996; Willert & Gharib 1997; Zhang et al. 1999; Nagaosa & Handler 2003; Terrington et al. 2022). 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.

6. Conclusions

We have performed a statistical investigation of the vortex flow beneath two types of surface features: dimples and scars. The former has been much studied and is fairly well understood while the latter has received little attention in the past. We make use of data from six direct numerical simulations, comprising different Reynolds numbers and three Weber numbers. As is well known, the flow beneath the surface has distinct behaviour inside the viscous layer of width $L_\nu$ equal to approximately a quarter of the Taylor microscale $\lambda _{\textrm{T}}$ , and the blockage layer underneath which extends down to a depth of approximately one integral scale of the bulk turbulence.

Our goal is to establish what statistical inferences can be drawn about sub-surface flow from observing the surface only, a question of significant practical importance. Towards this end, we establish a pair of concepts which, when suitably formulated, prepares the ground for interpretation of observed surface features: surface-only criteria (SOC) which allow a type of imprint (here: dimples and scars) to be identified only from information about the free-surface elevation, and a pseudo-causal definition (PCD) which defines an observable surface feature of a given type as the imprint of an identifiable sub-surface flow structure.

For dimples, an SOC–PCD pair can be readily formulated from previous work: a set of SOC is that a dimple as a depression in the surface which is approximately circular and persists for a long time (compared with other relevant scales). Conversely, a ‘true dimple’ can be defined as the imprint of a vertical vortex tube attached to the surface, appearing where the tube terminates (PCD). This SOC–PCD pair has been previously shown to be in excellent agreement (i.e. they identify the same areas of the surface as dimples).

For scars, we here propose a corresponding set. SOC: a scar is a depression in the surface that is strongly elongated and lives for a long time; PCD: a ‘true scar‘ can be defined as being the imprint of a horizontal vortex beneath it, which is centred approximately at the lower edge of the surface viscous layer, a depth equal to approximately a quarter of the Taylor microscale.

We present statistical evidence to demonstrate that this understanding is a successful description. First, we considered the conditional probability that a point in the flow belongs to a vortex core (according to the $\lambda _2$ criterion), given that the point lies directly beneath a dimple or scar. Then, the corresponding statistics of the orientation of the vorticity field beneath the same features were found. We found thus that nearly all vortices extending vertically down from dimples reach throughout the viscous sublayer and well into the blockage layer, and the largest vortices can reach all the way to the bottom of said layer where the flow is approximately isotropic. The vortex probability beneath dimples decreases with depth $\bar {z}$ as a Gaussian with variance increasing approximately linearly with the area of the dimple. For varying Reynolds and Weber numbers, the dependence of the conditional probability on $\bar {z}$ is a function of $\bar {z}/\lambda _{\textrm{T}}$ and independent of surface tension (for the range of values considered here). Moreover, we found that the horizontal vortex beneath scars is invariably centred at the lower edge of the surface viscous layer, independently of the size of the scar. Going downwards beneath the scar, the cumulative probability of having encountered a vortex increases from negligible to $\gtrsim 95\,\%$ (averaged over our different cases) in the range $0\lt \bar {z}\lt 1.5{L_\nu }$ . After removing a small number of regions identified as scars which also contain a surface-attached vortex (when the two are so close that our simple detection method mistakes it for a part of the scar), we observed that there is practically no vertically oriented vorticity beneath scars within the viscous and upper half of the blockage layer, all having been deflected by the horizontal vortex situated there. A preference for horizontal vorticity persists beneath scars throughout the blockage layer.

The statistical relations we have found between surface shape and sub-surface vortex structures not only support the suggested SOC–PCD pairs for dimples and scars, but also serve to elucidate further the nature of the coherent vortex structures in the surface-influenced layer of a turbulent free-surface flow, as well as their connection to imprints on the free surface. While it is well known that the turbulence in a free-surface flow with no mean flow changes gradually through the blockage layer, from being essentially isotropic outside of it to becoming highly anisotropic as the viscous layer is approached, our statistical study of vorticity inclination angles – and by inference, the orientation of vortex tubes – offers a more detailed and quantitative description. Surface-attached vortices, already well known to appear beneath dimples, are found to penetrate deep into the blockage layer, the largest reaching all the way to $\bar {z}\sim L_\infty$ ; this is seen by how the vorticity, which is necessarily normal to the surface at $\bar {z}=0$ , retains a strong vertical preference throughout the blockage layer, with virtually no horizontal vorticity present in its upper half. Correspondingly, horizontally oriented vortex tubes associated with scars cluster at the edge of the viscous boundary layer, around $\bar {z}\approx L_\nu$ , with remarkable consistency, rather than penetrating it. As such, the viscous boundary layer acts as a ‘soft boundary’, its lower edge demarking the depth where vortex tubes are forced to either penetrate the boundary layer and attach to the surface, or otherwise arrange parallel to it.

We are thus able to conclude that scars on the surface can be used to identify the position of strong, horizontally oriented vortex tubes, similar to how dimples show the location of vertical ones. One might say that dimples directly disclose surface-attached vortical structures and scars indirectly disclose sub-surface vortical structures, in the sense that there can be no surface-parallel vorticity at the surface itself.