2.1 Numerical methods

The data sets of stratified turbulent wakes and the far-field evolution of the wake-emitted IWs are generated numerically by a highly scalable parallel spectral multidomain penalty incompressible Navier–Stokes solver developed by Diamessis, Domaradzki & Hesthaven (Reference Diamessis, Domaradzki and Hesthaven2005). This flow solver has been used to investigate stratified wakes by Diamessis et al. (Reference Diamessis, Spedding and Domaradzki2011) (hereinafter referred to as DSD) and the wake-emitted IWs in the near field by A&D through implicit large-eddy simulations. Critical to the suitability of this solver in the present study is its spatial discretization: Fourier in the horizontal and Legendre multidomain penalty-based subdomains in the vertical, which allows one to attain high resolution in particular vertically localized flow regions of interest. Such vertical discretization also enables flexibility in choosing boundary conditions in the top and bottom boundaries. Spectral filtering and a penalty scheme provide the numerical stability needed for simulating the larger scales of high- $Re$ flow without having to resolve the full spectrum of turbulent motions.

Figure 1. Computational domain for implicit large-eddy simulation of a temporally evolving, stratified towed-sphere wake (Diamessis et al. Reference Diamessis, Domaradzki and Hesthaven2005, Reference Diamessis, Spedding and Domaradzki2011; Abdilghanie & Diamessis Reference Abdilghanie and Diamessis2013). The wake centreline is at $(y,z)=(0,0)$ . The sphere is towed in the direction of positive $x$ , which is the only statistically homogeneous direction in the computational set-up. $x=0$ is set to the midpoint of the domain in the homogeneous streamwise direction. The effect of the towed sphere is not computed explicitly but rather introduced as a complex two-stage turbulent wake initialization procedure (Diamessis et al. Reference Diamessis, Spedding and Domaradzki2011). The wave sponge layer (Abdilghanie Reference Abdilghanie2010) operates on the bottom and lateral boundaries, and the top boundary is a free-slip rigid lid.

2.2 Problem geometry

The flow configuration (see figure 1) and the implementation of numerical methods are almost identical to those of DSD, to which readers interested in the details are referred. Here the necessary modifications applied to the DSD set-up for the purpose of pursuing the objectives of the present study are summarized:

Table 1. Summary of numerical simulations at various wake Reynolds and Froude numbers. $N_{x}\times N_{y}\times N_{z}$ is the number of grid points; $N_{p}$ is the number of MPI processors used.

1. (i) Domain dimensions. The computational domain of dimensions $L_{x}\times L_{y}\times L_{z}$ is made taller and wider (see dimensions in table 1) than in DSD. This is to provide more space in the wake ambient fluid for the radiated IWs to evolve and reflect off the top boundary, i.e. the model sea surface, over a sufficiently wide spanwise extent. The wake centreline, corresponding to $(y,z)=(0,0)$ , is always located at a distance of $9D$ below the top boundary. As a result, the model sea surface is at a fixed vertical distance above the wake centreline in all simulations. The fixed distance $9D$ is chosen such that it provides the IWs enough space to evolve in the far field before being reflected. Such a domain depth is also feasible to simulate with available computing resources.

2. (ii) Wave sponge layers. A&D set up sponge layers (Abdilghanie Reference Abdilghanie2010) to wrap up completely the top/bottom ( $xy$ ) and lateral ( $xz$ ) boundaries to prevent the IWs from re-entering the domain by reflecting or through the periodic boundaries. In the present study, however, the sponge layer on the top boundary is removed to allow the IWs to reflect from the model sea surface (see figure 2).

3. (iii) Boundary conditions. Similar to DSD, periodic boundary conditions are employed in the horizontal directions. The top surface is free-slip to model the sea surface under low wind, and has zero vertical gradient of density perturbation. The model sea surface is flat and non-deformable in the vertical, following scaling arguments (Moum & Smyth Reference Moum and Smyth2006; Zhou & Diamessis Reference Zhou and Diamessis2013) which predict that the surface deflections due to IWs are much smaller than the horizontal length scales of the waves, and hence the surface slope is negligible.

4. (iv) Vertical resolution. Utilizing the spatial adaptivity of the spatial discretization in the vertical, the subdomain deployment focuses on both the turbulent wake core and the IW-reflecting subsurface zone in order for them to be resolved adequately.

Figure 2. Horizontal divergence ${\it\Delta}_{z}$ (normalized by buoyancy frequency $N$ ) fields sampled at the O $yz$ centreplane at various times (marked on top of the colour maps) at $Re=5\times 10^{3}$ and $Fr=4$ (R5F4). The interface between the fluid and wave-absorbing virtual sponge layer is marked by white dashed lines.

2.3 Initial condition

The initial condition of the simulations is a reasonable approximation of the near wake of a towed sphere in a stratified water tank. The reader may find the full details of the set-up in DSD. Here we provide a brief discussion of the most important aspects of the design and implementation of the initial condition. Once initialized, the flow field starts to evolve in time under the assumption that the turbulence is statistically homogeneous in the streamwise direction. Such an approach to simulate a temporally evolving shear flow, introduced in the pioneering work of Orszag & Pao (Reference Orszag and Pao1975), is commonly used in the stratified sphere-wake literature (Gourlay et al. Reference Gourlay, Arendt, Fritts and Werne2001; Dommermuth et al. Reference Dommermuth, Rottman, Innis and Novikov2002; Brucker & Sarkar Reference Brucker and Sarkar2010; Diamessis et al. Reference Diamessis, Spedding and Domaradzki2011; Redford, Lund & Coleman Reference Redford, Lund and Coleman2015) to avoid explicitly accounting for the sphere in the calculation. Alternatively, computing the flow explicitly around a sphere (Pasquetti Reference Pasquetti2011; Dairay, Obligado & Vassilicos Reference Dairay, Obligado and Vassilicos2015; Orr, Spedding & Constantinescu Reference Orr, Domaradzki, Spedding and Constantinescu2015) and simulating a spatially evolving wake (Pal, de Stadler & Sarkar Reference Pal, de Stadler and Sarkar2013) have also been implemented recently. However, the computational cost of such approaches is prohibitive for the present study, particularly at the higher $Re$ considered here.

Given the temporally evolving wake set-up, it is possible to map the time $t$ associated with an instantaneous flow field in the simulation to the corresponding downstream distance $X$ from the sphere to the computational domain (the reader should not confuse $X$ with the local coordinate $x$ of the computational domain, see the caption of figure 1). Specifically, in analogy with the field of view of the corresponding laboratory experiments (e.g. Spedding et al. Reference Spedding, Browand and Fincham1996), at time $t$ in the simulation, the right end of the computational domain is assumed to be at a distance $X=Ut$ from the centre of the sphere, which is assumed to already have travelled, from left to right, through the domain in the positive streamwise direction (Orszag & Pao Reference Orszag and Pao1975).

In setting up the initial flow field for our simulations, the mean and r.m.s. velocity profiles are prescribed by the model by Meunier, Diamessis & Spedding (Reference Meunier, Diamessis and Spedding2006) at a downstream distance $X/D=2$ where a self-similar flow is assumed. By virtue of (1.1) and the choice of $X/D=2$ as the initial condition, the initial dimensionless time $(Nt)_{0}$ is $1$ , $0.25$ and $0.0625$ for the simulations at $Fr=4$ , $16$ and $64$ , respectively. These initial $Nt$ values are smaller than $Nt\approx 2$ , which is the typical transition point in time from the three-dimensional (3-D) regime to the NEQ regime (Spedding Reference Spedding1997). It is after this transition that buoyancy effects begin to impact the wake significantly and the turbulence-driven IW radiation takes place. Therefore, the simulated turbulence-emitted IW field is expected be a faithful representation of its experimental counterpart, although the numerical model does not explicitly compute the wake at $Nt<(Nt)_{0}$ (or $X/D<2$ ) and a self-similar wake is assumed from the onset of the simulation.

2.4 Summary of numerical simulations

The six simulations performed for the present study have the same wake parameter space investigated in DSD, i.e. $Re\in \{5\times 10^{3},10^{5}\}$ and $Fr\in \{4,16,64\}$ . A summary of these simulations is in table 1. Hereinafter, each simulation will be labelled as R $a$ F $b$ , where $a=Re/10^{3}$ and $b=Fr$ . The computational cost increases dramatically with both $Re$ and $Fr$ . A regridding strategy similar to that of DSD is applied to R100F64 in the lateral direction $y$ whenever the half-width of the wake $L_{H}$ exceeds 10 % of the width of the domain.

4.1 Arrival time of peak wave impact

Since the passage of the wake-generating body, there exists a dimensionless time, $N\hat{t}$ , at which the most energetic waves arrive at the surface. We are to determine this critical $N\hat{t}$ , and then determine how this arrival time of peak wave impact changes over the wake parameter space $(Re,Fr)$ . To this end, a model is developed based on linear IW theory and verified by numerical simulations.

One can apply linear IW theory (Sutherland Reference Sutherland2010) to estimate the time for an IW to propagate from $z_{0}$ to $z$ . This time lag can be written as:

(4.1) $$\begin{eqnarray}t-t_{0}=\frac{z-z_{0}}{c_{g,z}}=\frac{z-z_{0}}{{\displaystyle \frac{N}{2{\rm\pi}/{\it\lambda}_{H}}}\sin {\it\theta}\cos ^{2}{\it\theta}},\end{eqnarray}$$

where $(z_{0},t_{0})$ can be considered as the vertical location and time from and at which the waves are emitted, i.e. the virtual origin of the waves, and $c_{g,z}$ is the vertical component of the wave group velocity. Additionally, ${\it\lambda}_{H}$ is the horizontal wavelength, and ${\it\theta}$ is the angle between the wave’s group velocity and the vertical direction.

The critical time $\hat{t}$ at which the surface undergoes the peak wave impact can be considered as the time at which waves of wavelength $\hat{{\it\lambda}}_{H}$ , the most energetic wavelength, arrive at the surface. Noting that the surface is located at $z=9D$ in all simulations, and assuming that in the far field of the wake, $\hat{t}\gg t_{0}$ , i.e. the waves are emitted at a considerably early $t_{0}$ and it takes a considerable amount of time for a wave to propagate in space to reach the observation plane, one can apply (4.1) for waves of wavelength $\hat{{\it\lambda}}_{H}$ arriving at the surface at $\hat{t}$ and obtain

(4.2) $$\begin{eqnarray}N\hat{t}\approx \frac{2{\rm\pi}}{\sin {\it\theta}\cos ^{2}{\it\theta}}\left(\frac{9D-z_{0}}{D}\right)\left(\frac{\hat{{\it\lambda}}_{H}}{D}\right)^{-1}.\end{eqnarray}$$

It is found in § 4.2.5 that

(4.3) $$\begin{eqnarray}\hat{{\it\lambda}}_{H}/D\propto Fr^{1/3},\end{eqnarray}$$

i.e. the most energetic wavelength scales as $Fr^{1/3}$ , and in § 4.2.3 that the virtual origin $z_{0}$ can be considered as a constant across all simulations. Moreover, the dependence of ${\it\theta}$ in $Fr$ is found to be rather weak (see § 4.3 for more details). As a result, one can deduce from (4.2) that

(4.4) $$\begin{eqnarray}N\hat{t}\propto Fr^{-1/3},\end{eqnarray}$$

which predicts that the peak arrival time in $Nt$ scales as $Fr^{-1/3}$ .

Figure 4. Time series of surface r.m.s. horizontal divergence ${\it\Delta}_{z,rms}$ normalized by the buoyancy frequency $N$ : (a,b) at $Re=5\times 10^{3}$ and (c,d) at $Re=10^{5}$ . The peak wave amplitude occurs at a narrow range of $Fr^{1/3}Nt$ for each $Re$ , as marked by the shaded area in (b,d).

To verify (4.4), time series of the spatial r.m.s. values (sampled over the surface $xy$ plane) of horizontal divergence ${\it\Delta}_{z}$ (a quantity denoted by ${\it\Delta}_{z,rms}$ as a measure of wave-induced strain) are computed for all simulations. The ${\it\Delta}_{z,rms}$ amplitude is normalized by $N$ , which forms a dimensionless measure of the waves’ steepness (A&D) and bears implications for the visibility of these waves (see § 5.2). The results are shown in figure 4, depicting the time evolution of the surface-observed wave strain amplitude. A delay of peak arrival time in $Nt$ with decreasing $Fr$ can be observed at both values of $Re$ . When the horizontal axis $Nt$ is rescaled by a factor of $Fr^{1/3}$ , the peak arrival times collapse into a considerably narrow range of $Fr^{1/3}Nt$ values, i.e. $[120,130]$ at $Re=5\times 10^{3}$ and $[230,330]$ at $Re=10^{5}$ , as highlighted by the shaded areas in figure 4(b,d). This collapse implies that the arrival time of peak wave impact at the surface in terms of $Nt$ scales approximately as $Fr^{-1/3}$ at a given $Re$ . This result is consistent with (4.4), a prediction due to linear wave theory.

One may alternatively consider the downstream distance from the towed sphere $\hat{X}/D$ at which the strongest waves arrive at the surface. Noting that the conversion from $Nt$ to $X/D$ involves a linear factor of $Fr/2$ as in (1.1), one can obtain from (4.4) that

(4.5) $$\begin{eqnarray}\frac{\hat{X}}{D}=\frac{Fr}{2}N\hat{t}\propto Fr^{2/3},\end{eqnarray}$$

which implies that the peak wave strains arrive at a shorter downstream distance normalized by the sphere diameter $X/D$ (which, however, corresponds to a later dimensionless time in buoyancy units, i.e. $Nt$ ) at a lower $Fr$ .

On the other hand, while the time (or distance) to peak wave impact at the surface can be modelled using simple kinematic arguments of IWs, it is inconclusive as to how the peak wave magnitude itself scales with $Fr$ , although an increase of ${\it\Delta}_{z}/N$ with $Fr$ can clearly be observed. The waves emitted by the higher- $Re$ wakes are also recorded to have higher amplitude when arriving at the surface. Similar observations for the amplitudes have also been reported in the wake’s near field (Abdilghanie Reference Abdilghanie2010). Understanding the detailed mechanism by which the waves extract momentum and energy from the turbulence would shed light on how the energy content of the turbulence-emitted waves vary with $(Re,Fr)$ . This mechanism is not known completely, and is deferred to further study focusing on turbulence–wave interaction in the wake’s near field.

In this subsection, the time evolution of the surface-observed wave amplitude has been examined. A key ingredient of the model for the time (or distance) to peak wave impact is the scaling of the most energetic wavelength, i.e. (4.3). In the following subsection, the analyses on the length scales of these waves are presented.

4.2 Wavelength

4.2.2 Joint distribution of wave events in wavelength–time space

The variation of the wavelengths ${\it\lambda}_{H}$ corresponding to the most energetic wave packets is found to be much weaker in the spanwise direction $y$ (as can be seen from sample snapshots in figure 3) than the variation in $Nt$ . One can therefore proceed by also collapsing the data in $y/D$ as well, and focus on the evolution of the length scale ${\it\lambda}_{H}$ in time. The joint PDF of these wave events constructed in the $({\it\lambda}_{H}/D,Nt)$ space is shown in figure 5. In agreement with figure 4, at a given $Re$ , energetic packets with the most frequently occurring wavelength reach the surface at a later $Nt$ as $Fr$ is reduced. At the higher $Re$ , the arrival time of the most frequently occurring wavelength is shifted in time towards a later $Nt$ . The joint PDFs appear to lie along a universal line, with the centres of the PDF contours varying as a function of $(Re,Fr)$ . In § 4.2.3, the universal trend in which the prevalent wavelength ${\it\lambda}_{H}$ evolves with $Nt$ , characterized by the evolution of the conditionally averaged wavelength $\overline{{\it\lambda}}_{H}$ at each time (see the time variation of dominant wavelength in figure 3), will be further discussed. In § 4.2.5, the most energetic wavelength $\hat{{\it\lambda}}_{H}$ , i.e. the wavelength ${\it\lambda}_{H}$ corresponding to the peaks of the joint PDFs, will be parameterized by the wake’s $Re$ and $Fr$ .

Figure 5. Joint probability density contours of most energetic wave packets in the $({\it\lambda}_{H}/D,Nt)$ sample space at (a) $Re=5\times 10^{3}$ and (b) $Re=10^{5}$ , respectively, for $Fr\in \{4,16,64\}$ .

Here we highlight two differences between the analyses performed by A&D (who employed the same continuous wavelet transform techniques) and the results presented here: (1) A&D focused on the wake’s near field, whereas we examine the far field; and (2) A&D considered exclusively the variation of conditional mean $\overline{{\it\lambda}}_{H}$ in time (which we show in figure 6), whereas we also discuss the most energetic wavelength $\hat{{\it\lambda}}_{H}$ , which is identified from the joint probability distribution of wave events in the $({\it\lambda}_{H}/D,Nt)$ space (figure 5).

4.2.3 Evolution of mean observable wavelengths

Figure 6 shows the conditional mean IW horizontal wavelength, $\overline{{\it\lambda}}_{H}$ , as a function of $Nt$ for all simulations. $\overline{{\it\lambda}}_{H}$ is the mean ${\it\lambda}_{H}$ computed based on the wave events extracted in the $({\it\lambda}_{H}/D,y/D)$ sample space at a given $Nt$ (see § 4.2.1(b)). For the range of $Nt$ values at which a significant probability density is present (figure 5), the time evolution of $\overline{{\it\lambda}}_{H}$ in $Nt$ seems to follow a universal line across all simulations covering the wake parameter $(Re,Fr)$ . The decay of $\overline{{\it\lambda}}_{H}$ in $Nt$ follows a power law of $(Nt)^{-1}$ which is consistent with the laboratory experiments and analysis by Bonneton et al. (Reference Bonneton, Chomaz and Hopfinger1993) and Spedding et al. (Reference Spedding, Browand, Bell and Chen2000). This observation motivates the linear wave model that is to be constructed and calibrated below.

Figure 6. Time series of conditionally averaged wavelengths $\overline{{\it\lambda}}_{H}/D$ observed at the sea surface. Best fit to the linear propagation model (4.7) is shown as dashed line.

One may consider the linear propagation of IWs and rewrite (4.1) as

(4.6) $$\begin{eqnarray}\frac{{\it\lambda}_{H}}{D}=\frac{z-z_{0}}{D}\frac{2{\rm\pi}}{\sin {\it\theta}\cos ^{2}{\it\theta}}\frac{1}{N(t-t_{0})}.\end{eqnarray}$$

Again, consider the wake’s far field and invoke the far-field condition $t\gg t_{0}$ :

(4.7) $$\begin{eqnarray}\frac{{\it\lambda}_{H}}{D}\approx \frac{z-z_{0}}{D}\frac{2{\rm\pi}}{\sin {\it\theta}\cos ^{2}{\it\theta}}\frac{1}{Nt}.\end{eqnarray}$$

Hence the power law $(Nt)^{-1}$ is recovered in (4.7), in agreement with the observations shown in figure 6 and experiments (Bonneton et al. Reference Bonneton, Chomaz and Hopfinger1993; Spedding et al. Reference Spedding, Browand, Bell and Chen2000). The effectiveness of linear theory here indicates that these IWs emitted by turbulence act in the far field as linear waves, at least within the parameter range and uniform stratification examined here. The key physical process driving the time evolution of wavelength at the surface is the waves’ dispersion, i.e. at a given ${\it\theta}$ , waves of greater wavelength propagate faster and thus arrive at the observation plane at an early $Nt$ , and vice versa.

Figure 7. Schematic of key concepts associated with the linear wave model (4.7) evaluated at the sea surface ( $z=z_{sfc}$ ). Internal waves (dashed arrows) of various wavelengths ${\it\lambda}_{H}$ are observed at the surface as if they had been emitted from the virtual origin $z_{0}$ at $t_{0}\ll t$ . Equation (4.7) relates the surface-observed mean ${\it\lambda}_{H}$ to the submerged turbulence characterized by the dimensionless time $Nt$ since the passage of the sphere and the corresponding dimensionless distance $X/D$ to the wake-generating body as in (1.1).

The linear propagation model in (4.7) can be applied specifically at the sea surface located at $z=z_{sfc}$ , where $z_{sfc}$ is the depth of the wake centreline below the surface (figure 7). It is desirable to first complete the linear propagation model in (4.7) by calibrating the virtual origin $z_{0}$ of the waves based on the observations at $z_{sfc}/D=9$ in our simulations. Note that the intended applicability of the model and the calibration parameter $z_{0}$ is restricted to the wake’s far field, for which $(z,t)\gg (z_{0},t_{0})$ , and to conditions in which the waves are indeed linear and given enough room to disperse when they arrive at the vertical location $z$ . In the wake’s near field, however, the waves’ interaction with turbulence is possibly nonlinear (Sutherland & Linden Reference Sutherland and Linden1998), and the waves are not fully dispersed in space, being close to the turbulent source. Therefore, the model (4.7) and any fitting parameter are not expected to give a reliable description of the waves in the near field. In fact, A&D reported a power-law decay of $\overline{{\it\lambda}}_{H}$ in $Nt$ with exponents different from $-1$ , based on near-field sampling of these waves.

In order to calibrate $z_{0}$ for the model (4.7), a least squares fit can be applied to the data points in figure 6 contributed by all simulations. The best fit to the data following (4.7) is found to be (with 95 % confidence level):

(4.8) $$\begin{eqnarray}\frac{{\it\lambda}_{H}}{D}=\frac{149\pm 2}{Nt},\end{eqnarray}$$

based on observations at a sea surface located at $z/D=9$ . By substituting ${\it\theta}=45^{\circ }$ , a characteristic prevalent wave emission angle reported by A&D, into (4.7), and utilizing the fit in (4.8), the virtual origin (based on observations in the far field) is found to be

(4.9) $$\begin{eqnarray}\frac{z_{0}}{D}=0.62\pm 0.11.\end{eqnarray}$$

$z_{0}$ is strictly a calibration parameter for the far-field linear wave model (4.7) which enables the prediction of wavelength in the far field for an arbitrary $z$ (see in § 4.2.4). As such, $z_{0}$ is not to bear implications for the local dynamics at $z=z_{0}$ in the wake’s near field, and it should not be confused with the location of the wake edge, itself a function of $Fr$ (figure 7 of DSD).

Interestingly, Brandt & Rottier (Reference Brandt and Rottier2015) built upon the theories of Gilreath & Brandt (Reference Gilreath and Brandt1985), Voisin (Reference Voisin, Redondo and Métais1995) and Dupont & Voisin (Reference Dupont and Voisin1996), and deduced that, at $Fr=4$ , the origin of the turbulence-generated waves would be at $z_{0}/D\approx 1$ and $Nt_{0}\approx 3$ , which is comparable to the numerical results presented here for the given ${\it\theta}=45^{\circ }$ . However, there is a considerable degree of sensitivity of the calibrated $z_{0}/D$ value in (4.9) to the choice of polar angle ${\it\theta}$ ; for example, as ${\it\theta}$ varies from $45^{\circ }$ to $55^{\circ }$ , a range in which the IWs are observed to be the most energetic (see figure 10 in § 4.3), the $z_{0}$ value may vary by approximately $2D$ .

4.2.4 Application of the linear model

Figure 8. Time evolution of ${\it\lambda}_{H}/D$ in $Nt$ for various wake depths, $z_{sfc}/D\in [9,18,36]$ , as predicted by (4.7) using the calibrated virtual origin $z_{0}/D=0.62$ for waves propagating at ${\it\theta}=45^{\circ }$ to the vertical. A given ${\it\lambda}_{H}/D$ observed at the surface (marked by the horizontal grey line) may correspond to different $Nt$ values (where the vertical grey arrows point at, for each $z_{sfc}/D$ considered in the figure) of the submerged turbulence, which may be located at various depths.

As sketched in figure 7, the specific application of linear model (4.7) that is of interest is to connect surface observations of ${\it\lambda}_{H}/D$ at a given point in time with the information associated with the submerged turbulence. This information may include dimensionless time $Nt$ elapsed since the passage of the object, and/or the dimensionless distance $X/D$ to the object. The evolution of wavelength in space and time $(z/D,Nt)$ specified by (4.7) is independent of the wake parameters $(Re,Fr)$ , if one ignores the weak dependence of ${\it\theta}$ in $Re$ (see more in § 4.3). Therefore, knowing these parameters a priori may not be a prerequisite for the connection that one would like to make. Figure 8 illustrates some of the implications of (4.7): when the depth $z_{sfc}$ and the $Nt$ value of the submerged turbulence are both known, one can readily compute the corresponding ${\it\lambda}_{H}/D$ to expect at the surface, with (4.7) and the calibrated $z_{0}$ in (4.9). However, inferring $Nt$ solely from surface observations of ${\it\lambda}_{H}/D$ is not as straightforward, because ${\it\lambda}_{H}/D$ is not uniquely determined by $Nt$ but rather the combination of $z_{sfc}/D$ and $Nt$ . One can deduce from observations of ${\it\lambda}_{H}/D$ the combinations of $z_{sfc}/D$ and $Nt$ that would produce such a wavelength, but neither $z_{sfc}/D$ nor $Nt$ individually (see figure 8).

4.2.5 Selection of peak wavelength by wake turbulence

While the time evolution of the mean observable wavelength at a given $(z/D,Nt)$ is universally determined by the kinematic properties of IWs prescribed by the linear theory, the most probable wavelength $\hat{{\it\lambda}}_{H}$ to be observed, i.e. the most energetic one, does depend on the parameters of the wave-emitting wake, i.e. $(Re,Fr)$ (see figure 5). Figure 9 shows the most frequently occurring wavelength identified from figure 5, corresponding to the peaks of the joint PDFs, as a function of $Fr$ . For each value of $Re$ , an approximate $Fr^{1/3}$ -scaling can be observed for $\hat{{\it\lambda}}_{H}/D$ , which is consistent with previous studies (Spedding et al. Reference Spedding, Browand, Bell and Chen2000; Abdilghanie & Diamessis Reference Abdilghanie and Diamessis2013). A dependence on $Re$ can also be observed: increasing $Re$ from $5\times 10^{3}$ to $10^{5}$ reduces $\hat{{\it\lambda}}_{H}$ by 50 % at a given $Fr$ . More discussion on this specific length-scale selection by wake turbulence is given in § 6.2.

Figure 9. $Fr$ - and $Re$ -dependence of the most energetic wavelength $\hat{{\it\lambda}}_{H}/D$ observed at the sea surface.

4.3 Wave period

Similar to the length-scale computation outlined in § 4.2, a one-dimensional wavelet transform (with the Morlet mother wavelet) has been applied to time series of horizontal divergence ${\it\Delta}_{z}$ signals, sampled at every grid point on the surface, to interrogate the period of these waves. One can then isolate, in a given time series of the transform modulus, the most energetic wave packets (with transform modulus above 35 % of the maximum value at all values of $t$ for all locations) as a function of time of occurrence $Nt$ and wave period $T$ . Again, the $x$ -dependence is to be ignored, i.e. the data are to be collapsed in the statistically homogeneous $x$ axis. Effectively, for each $Re$ and $Fr$ , a large population of energetic wave events are sampled in the three-dimensional sample space $(y/D,Nt,T/(2{\rm\pi}/N))$ . One can then collapse these three-dimensional data in two ways, either along the $Nt$ axis onto the $y/D$ – $T/(2{\rm\pi}/N)$ plane, or along the $T$ axis on the $y/D$ – $Nt$ plane. On each of these planes, joint PDFs may be constructed, as shown in figures 10–12.

Figure 10. Joint probability density contours of the most energetic wave packets in the $(T/(2{\rm\pi}/N),Nt)$ sample space: (a) R5F4 and (b) R100F4. $T/(2{\rm\pi}/N)$ values corresponding to polar propagation angles ${\it\theta}$ of $\{30^{\circ },45^{\circ },55^{\circ },60^{\circ },65^{\circ }\}$ are marked in dashed lines for reference. For R5F4, the PDF peaks at polar angle ${\it\theta}=55^{\circ }$ , consistent with the prediction of Voisin (Reference Voisin, Redondo and Métais1995) for the maximum emission angle. For R100F4, the peak polar angle moves noticeably towards a lower value falling in between $45^{\circ }$ and $55^{\circ }$ .

Figure 11. Distribution of most energetic wave packets in the $(|y|/D,T/(2{\rm\pi}/N))$ sample space for $Re=5\times 10^{3}$ , $Fr=4$ : (a) probability density contours, (b) marginal PDF versus $|y|/D$ and (c) marginal PDF versus $T/(2{\rm\pi}/N)$ . Marginal PDFs are obtained by integrating the joint PDF over the entire range of one of the random variables, e.g. see Devore (Reference Devore2015).

Figure 12. Distribution of most energetic wave packets in the $(|y|/D,T/(2{\rm\pi}/N))$ sample space for $Re=10^{5}$ , $Fr=4$ : (a) probability density contours, (b) marginal PDF versus $|y|/D$ and (c) marginal PDF versus $T/(2{\rm\pi}/N)$ .

Figure 10 shows the distribution of most energetic wave packets at $Fr=4$ in $(y/D,Nt)$ sample space. Results at $Fr\in \{16,64\}$ on the distribution of wave packets with respect to $y/D$ are qualitatively similar, and thus not shown. The structure of the joint PDF indicates that the most energetic wave packets arrive at a later time and over a longer window in time at the higher $Re$ , in agreement with figure 5. Moreover, unlike the strong variation of wavelengths with $Nt$ , the variability of wave periods reaching the surface with time is rather weak, especially after the peak wave amplitude has been reached at $Nt=70$ and $180$ for R5F4 and R100F4, respectively.

However, a strong variation of wave periods can be observed in space; specifically, in the spanwise direction $y$ . Figures 11 and 12 illustrate this spatial variability of wave periods and clearly show the dispersive character of the wake-radiated IWs. Regardless of $Re$ , the wave period reaches the peak likelihood at approximately 1.5 times the buoyancy period $2{\rm\pi}/N$ . Longer periods that are progressively smaller in likelihood (as shown by the marginal PDFs versus $T/(2{\rm\pi}/N)$ in figures 11(c) and 12(c)), arrive at larger spanwise offsets. This is because these waves have a group velocity vector that forms a larger angle ${\it\theta}$ with respect to the vertical, following the dispersion relation of IWs. This distinct correlation between the spanwise location $y$ at which the waves arrive at the surface and the waves’ period $T$ is found to be independent of $Fr$ . At the higher $Re$ shown at $Fr=4$ , i.e. R100F4 in figure 12(a), the $y/D$ – $T/(2{\rm\pi}/N)$ correlation is stronger as compared to the lower- $Re$ counterpart, i.e. R5F4 in figure 11(a). Moreover, at R100F4, the marginal PDF in $y/D$ , as shown in figure 12(b), starts to decay as far as $15D$ away from the centreline, whereas at R5F4 the PDF decays immediately after it peaks at $|y|=7D$ , as shown in figure 11(b). Consequently, the waves at higher wake $Re$ are more likely to impact larger spanwise distances from the wake centreline.

Near-field measurements of the wake-emitted waves by A&D reported one relatively narrow band of prevalent wave frequencies (and thus propagation angles) for each of the two $Re$ values, respectively. It is also reported that a selective viscous decay mechanism, similar to that proposed for waves from a turbulent bottom Ekman layer (Taylor & Sarkar Reference Taylor and Sarkar2007), acts in the wake’s near field to further narrow down the frequency band, at the lower Reynolds number ( $Re=5\times 10^{3}$ ). In the far field, however, the waves are in fact observed to be more broadband in the frequency domain, as compared to laboratory results corresponding to polar angles ${\it\theta}=55^{\circ }$ – $65^{\circ }$ sampled at vertical distances up to $z=4D$ via visual observations of isopycnal displacements (Brandt & Rottier Reference Brandt and Rottier2015). This is because the various frequencies are given enough distance to disperse when they are allowed to propagate over a greater distance into the far field. Dispersion of these waves in the wake’s lateral direction due to the variability in propagation angles ${\it\theta}$ widens the spatial extent in $y$ in which the waves may manifest themselves. (The polar angle ${\it\theta}$ of the wave’s group velocity is not to be confused with the azimuthal angle ${\it\alpha}$ associated with the wave strain field on the surface, the latter angle being the subject of the next subsection.) In light of this, one would naturally expect the waves to disperse more in space (specifically, in $y$ ), should there be a greater vertical offset from the wave source to the observation level. Therefore, the degree of wave dispersion in the spanwise direction $y$ (characterized by the width of wave impact in $y$ ) on the observational plane may give an estimate of the depth of the turbulent source, especially when such information is not available.

4.4 Wave orientation

Information on the location and the orientation at which the IWs are most likely to reach the surface may assist the recognition of these wave patterns from remotely acquired sea-surface images. At any given time and at each point on the sea surface, a two-dimensional eigenvalue/vector decomposition may be applied to the surface rate-of-strain tensor associated with the horizontal velocities $(u,v)$ (the vertical velocity $w$ vanishes at the rigid lid), which provides the magnitudes of the two principal strain rates and the directions of the two corresponding principal axes (PAs) on the $xy$ plane. The dominant PA (with the higher strain-rate magnitude) tends to align perpendicular to the waves’ isophase lines, across which velocity due to the wave varies the most rapidly. Therefore, the direction of the dominant PA gives a convenient measure of the waves’ orientation at a given location on the surface. The joint PDF of the orientation (azimuthal angle ${\it\alpha}$ with respect to the wake axis $x$ ) and the spanwise offset $y$ of the wave strains exceeding a certain threshold (80 % of the peak ${\it\Delta}_{z,rms}$ at all $Nt$ values) may then be computed. These results are shown in figure 13 for $Fr=4$ at two values of $Re$ . The results at $Fr\in \{16,64\}$ are similar, and thus not shown. Two prevalent preferred directions can be identified in figure 13: one at $|{\it\alpha}|=0.15$ –0.25 rad ( $9^{\circ }$ – $14^{\circ }$ ) which occurs closer to the wake centreline, and the other at $|{\it\alpha}|=0.7$ –0.8 rad ( $40^{\circ }$ – $46^{\circ }$ ) at a further distance ( $|y|>5D$ ). At the higher $Re$ , a noticeable protrusion of the probability density into greater distances ( $|y|>15D$ ) at azimuthal angle $1<|{\it\alpha}|<1.5$  rad can be observed.

Figure 13. Joint probability density contours of significant wave events in the $(|y|/D,|{\it\alpha}|)$ sample space: (a) R5F4 and (b) R100F4, where ${\it\alpha}$ is the azimuthal angle of the dominant principal strain rate with respect to the $x$ axis. $|{\it\alpha}|=30^{\circ }$ and $60^{\circ }$ are marked as vertical dashed lines for reference.

5.1 A general discussion

In this section, possible sea-surface IW manifestation mechanisms are studied. In order for the turbulence-emitted waves to self-manifest, the surface hydrodynamic condition must act in favour of one or more of the various surface modulation mechanisms of remote sensing significance (see § 1.2). Two outstanding issues, which need to be addressed from remote sensing perspectives, are contrast and spatial scales. First, to generate enough contrast on remote sensing images to be perceivable, e.g. to create enough Bragg scattering for the features to appear on the radar image (Alpers Reference Alpers1985), or to generate surface slicks of sufficient thickness (Moum, Carlson & Cowles Reference Moum, Carlson and Cowles1990), the characteristic surface strain rate, i.e. ${\it\Delta}_{z}$ , must exceed a certain threshold specified by the modulation mechanism. Second, the spatial scales of these surface anomalies must be small enough to fit the field of view of a certain remote sensor, and large enough compared to the sensor’s resolution. As a first attempt to assess the potential remote visibility of these waves, in § 5.2, the feasibility of generating surface slicks is examined and, in § 5.3, wave-induced surface mass transport is discussed.

5.2 Wave strains and surface slicks

5.2.1 Model formulation

One of the most common mechanisms through which surface velocity strains can modulate the sea surface, and thus create anomalies on remote sensing images, is through the organic-rich surfactant materials within the surface microlayer (Ewing Reference Ewing1950a ,Reference Ewing b ; Alpers Reference Alpers1985; Moum et al. Reference Moum, Carlson and Cowles1990; Ermakov et al. Reference Ermakov, Salashin and Panchenko1992; Soloviev & Lukas Reference Soloviev and Lukas2006; Klemas Reference Klemas2012). The microlayer can increase in thickness due to the accumulation of surfactant under favourable conditions, i.e. flow convergence and low mixing/dispersion due to wind. The surfactant can potentially form surface ‘slicks’ when the microlayer thickness exceeds a certain threshold. These slicks smooth out surface capillary–gravity waves and decrease the surface reflectance to light and microwaves. As a result, the slicks appear as dark bands on optical and radar images, and can possibly reveal the compressive surface strains in the surface flow. Moum et al. (Reference Moum, Carlson and Cowles1990) performed in situ measurements of the surface chemistry and the subsurface hydrodynamics simultaneously at the slicks. They were able to relate the formation of the slicks to the flow convergence, and reported a 5–25 % increase in the surfactant concentration across the slicks relative to the background.

In an attempt to predict the feasibility of slick formation under surface strains due to turbulence-radiated IWs, simulations of the two-dimensional advection equation for the surfactant concentration are performed at the topmost free-slip surface of the computational domain (see figure 1). These auxiliary simulations of the passive scalar field are forced by the surface velocity field generated by the wake simulations. The scalar simulations employ second-order finite difference in space and fourth-order Adams–Bashforth–Moulton (Press et al. Reference Press, Teukolsky, Vetterling and Flannery2007) in time. More realistic models which incorporate the diffusion/relaxation of surfactant gradients (Ermakov et al. Reference Ermakov, Salashin and Panchenko1992) and coupling with surface wave actions (Alpers Reference Alpers1985; Ermakov et al. Reference Ermakov, Salashin and Panchenko1992) have been proposed, and can be readily incorporated should the site-specific parameters involved in these models become available.

Figure 14. Surface snapshot of (a) wave-induced horizontal divergence ${\it\Delta}_{z}$ and (b) strain-driven surfactant concentration fields at $Nt=70$ for R100F64. Concentration of ${\it\Gamma}/{\it\Gamma}_{0}>1.05$ may correspond to surface slicks that can be detectable via remote sensing. Isolines with ${\it\Gamma}/{\it\Gamma}_{0}=1.05$ are marked in grey, delineating possible locations of slicks.

5.2.3 Linear analysis

It is desirable to correlate the ratio ${\it\Gamma}/{\it\Gamma}_{0}$ with the surface velocity strain magnitude. A scaling analysis of the two-dimensional governing equation is performed here. Decompose concentration ${\it\Gamma}$ as

(5.1) $$\begin{eqnarray}{\it\Gamma}(x,y,t)={\it\Gamma}_{0}+{\it\Gamma}^{\prime }(x,y,t),\end{eqnarray}$$

where ${\it\Gamma}^{\prime }$ is the perturbation to the background state ${\it\Gamma}_{0}$ . With this decomposition, the two-dimensional advection equation

(5.2) $$\begin{eqnarray}\frac{\partial {\it\Gamma}}{\partial t}+\frac{\partial ({\it\Gamma}u)}{\partial x}+\frac{\partial ({\it\Gamma}v)}{\partial y}=0,\end{eqnarray}$$

can be rewritten as

(5.3) $$\begin{eqnarray}\frac{\partial {\it\Gamma}^{\prime }}{\partial t}+{\it\Delta}_{z}({\it\Gamma}_{0}+{\it\Gamma}^{\prime })+u\frac{\partial {\it\Gamma}^{\prime }}{\partial x}+v\frac{\partial {\it\Gamma}^{\prime }}{\partial y}=0.\end{eqnarray}$$

With small perturbations

(5.4) $$\begin{eqnarray}{\it\Gamma}^{\prime }\ll {\it\Gamma}_{0},\end{eqnarray}$$

and the particle’s horizontal velocity magnitude much smaller than the waves’ celerity, (5.3) reduces to

(5.5) $$\begin{eqnarray}\frac{\partial {\it\Gamma}^{\prime }}{\partial t}+{\it\Delta}_{z}{\it\Gamma}_{0}=0.\end{eqnarray}$$

Given that

(5.6) $$\begin{eqnarray}\frac{\partial }{\partial t}\sim \frac{1}{T}\sim {\it\omega}\sim N,\end{eqnarray}$$

the perturbation ratio ${\it\Gamma}^{\prime }/{\it\Gamma}_{0}$ scales as

(5.7) $$\begin{eqnarray}\frac{{\it\Gamma}^{\prime }}{{\it\Gamma}_{0}}\sim \frac{{\it\Delta}_{z}}{N}.\end{eqnarray}$$

This simple analytical scaling is further confirmed by simulation data (figure 15), which indeed show that the ratio ${\it\Gamma}^{\prime }/{\it\Gamma}_{0}$ scales with ${\it\Delta}_{z}/N$ , a measure of the IW steepness (A&D). The surface-measured ${\it\Delta}_{z}/N$ , and subsequently ${\it\Gamma}^{\prime }/{\it\Gamma}_{0}$ , show a distinct increase with both $Re$ and $Fr$ (see § 4.1), suggesting that at operational parameter values (see § 6.5), at least for an idealized uniform stratification, these turbulence-emitted waves can produce remotely visible surface patterns.

Figure 15. Peak surfactant perturbation ratio ${\it\Gamma}^{\prime }/{\it\Gamma}_{0}$ versus the peak dimensionless surface strain $|{\it\Delta}_{z}|/N$ observed in all six simulations. Only ${\it\Gamma}^{\prime }>0$ , i.e. surfactant enrichment, and ${\it\Delta}_{z}<0$ , i.e. flow convergence, are considered. The characteristic enrichment ratio ${\it\Gamma}^{\prime }/{\it\Gamma}_{0}=0.05$ (dashed line) across observable slicks according to field measurements (Moum et al. Reference Moum, Carlson and Cowles1990) is exceeded by all cases except for R5F4. A scaling of ${\it\Gamma}^{\prime }/{\it\Gamma}_{0}\sim {\it\Delta}_{z}/N$ (solid line) can be observed.

5.3 Wave-driven Lagrangian flows

Here we study the mass transport driven by the reflecting IWs at the free-slip surface (with vertical velocity $w=0$ ). To investigate it, two-dimensional tracking of Lagrangian particles forced by the surface $(u,v)$ velocity field from the wake simulations has been performed. These Lagrangian tracers are inserted into the surface two-dimensional flow field at the earliest $Nt$ of each simulation when the surface flow velocities are virtually zero. Large number of tracers are inserted at regularly distributed locations $(x_{0}^{+},y_{0}^{+})$ at the surface at approximately every $0.1D$ . The particle-tracking scheme employs fourth-order cubic spline interpolation to obtain the particles’ instantaneous velocities from the simulation data, and marches in time using a locally fourth-order Adams–Bashforth–Moulton method (Press et al. Reference Press, Teukolsky, Vetterling and Flannery2007). In particular, we are interested in the mass transport in the spanwise direction $y$ due to the wave impact. To this end, the lateral displacement

(5.8) $$\begin{eqnarray}{\it\Delta}|y^{+}|\equiv \left\{\begin{array}{@{}ll@{}}y^{+}-y_{0}^{+}\quad & \text{ for }y_{0}^{+}\geqslant 0,\\ y_{0}^{+}-y^{+}\quad & \text{ for }y_{0}^{+}<0,\end{array}\right.\end{eqnarray}$$

can be defined, where $y^{+}(t)$ is the particle’s instantaneous position. Note that ${\it\Delta}|y^{+}|$ is displacement with respect to the initial position $y_{0}^{+}$ , and is constructed to be positive when the particle moves away from the wake centreline at $y=0$ .

Figure 16. PDF of lateral particle displacements away from the wake centreline at (a) R100F4 and (b) R100F16. Dashed lines denote the PDF before the most energetic waves interact with the surface, and solid lines the PDF after the most energetic wave impact.

We focus on the lateral transport due to the nonlinear effects of the waves during reflection (Zhou & Diamessis Reference Zhou and Diamessis2015). This transport is found to be more significant at the higher $Re=10^{5}$ with IWs of higher amplitude interacting with the surface, and the results are shown for R100F4 and R100F16 in figure 16. (Similar effects were observed for R5F64, but no distinct directionality of the mean drift was observed for R5F4 or R5F16, where the wave amplitudes reaching the surface are much weaker in amplitude and shorter in duration.) Before the most energetic waves reach the surface (sampled at $Nt=100$ and $40$ for R100F4 and R100F16, respectively), the PDF of ${\it\Delta}|y^{+}|$ centres around zero with a weak diffusion. After the most intense wave impact (sampled at $Nt=300$ and $250$ for R100F4 and R100F16, respectively), the PDF diffuses out more and drifts significantly towards the positive values of ${\it\Delta}|y^{+}|$ , which implies that the waves advect the Lagrangian tracers away from the wake centreline. This collective net motion of tracers at both sides of the wake centreline creates a local divergence in the lateral mass transport at $y=0$ . The magnitude of the drift increases significantly with $Fr$ , i.e. the maximum displacement approximately quadruples as $Fr$ increases from 4 to 16 at $Re=10^{5}$ . Following the same trend, much larger particle displacements are observed at R100F64. However, the particles at R100F64, especially the ones inserted close to the lateral boundaries of the computational domain, move far enough to reach the lateral sponge layers, where motions are artificially damped out. The particles then get trapped at the sponge layers, which makes the interpretation of this particular Lagrangian data set ambiguous. Future work is needed in this regard to investigate the Lagrangian dynamics of IW-impacted sea surfaces at high $Fr$ of the wave-emitting turbulence, with an even wider domain in the lateral direction to accommodate increasing particle excursion length scales with $Fr$ . Although the Lagrangian effects reported here may be considerably weak, i.e. net drifts on $O(0.01)$ – $O(0.1)$ of the diameter $D$ of the turbulence-generating body, the effect seems to scale strongly with $Fr$ , and may potentially be significant at operational parameter values (see § 6.5).