2.1 Computational set-up

We perform DNS of breaking internal waves on slopes in an idealized, laboratory-scale domain using the parallel Navier–Stokes code of Cui (Reference Cui1999). This code employs the fractional-step method of Zang, Street & Koseff (Reference Zang, Street and Koseff1994), and has been used extensively in the past to study geophysical flows at the laboratory scale (see Venayagamoorthy & Fringer Reference Venayagamoorthy and Fringer2007; Chou & Fringer Reference Chou and Fringer2010; Arthur & Fringer Reference Arthur and Fringer2014, Reference Arthur and Fringer2016). The computational set-up used here is based on that of Arthur & Fringer (Reference Arthur and Fringer2014, Reference Arthur and Fringer2016), and is summarized in figure 1.

Figure 1.

The computational set-up used to study breaking internal waves on slopes. ( $x_{1}$ , $x_{2}$ , $x_{3}$ ) are Cartesian coordinates in the cross-shore, lateral (along shore, into the page) and vertical directions, respectively. (a) A schematic of the domain, where the length $L=4~\text{m}$ , the total depth $H=0.56~\text{m}$ , the upper-layer depth $h_{1}=0.3~\text{m}$ , the lower-layer depth $h_{2}=0.26~\text{m}$ , the upper-layer density $\unicode[STIX]{x1D70C}_{1}=985~\text{kg}~\text{m}^{-3}$ , the lower-layer density $\unicode[STIX]{x1D70C}_{2}=1015~\text{kg}~\text{m}^{-3}$ , the slope $s=0.218$ , the amplitude of the initial half-Gaussian $a_{0}=0.1~\text{m}$ and the length scale of the initial half-Gaussian $L_{0}=0.7~\text{m}$ . The length of the flat region before the slope begins is $L_{s}=1.675~\text{m}$ and the radius of curvature of the rounded bottom at the beginning of the slope is 3 m. The pycnocline thickness $\unicode[STIX]{x1D6FF}$ varies for each case and is shown in table 1. The lateral ( $x_{2}$ , into the page) width of the domain is $W=0.1$  m. (b) The initial density profile for each case, as defined in (2.1). The stratified region around the interface $x_{3}=-h_{1}\pm \unicode[STIX]{x1D6FF}/4$ is highlighted by the intersection of the vertical dotted lines at $\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D70C}_{0}=0.987,1.013$ with each density profile. (c) The maximum buoyancy frequency for each case $N_{max}$ as a function of $\unicode[STIX]{x1D6FF}$ , including a scaling line of $N_{max}\sim \unicode[STIX]{x1D6FF}^{-1/2}$ (Troy & Koseff Reference Troy and Koseff2005).

Computations are made on an orthogonal curvilinear grid of size $n_{1}\times n_{2}\times n_{3}=1152\times 96\times 128\approx 14\times 10^{6}$ points that was generated using the open-source software gridgen (available online at https://code.google.com/p/gridgen-c/). Grid stretching of the form $\unicode[STIX]{x0394}x_{i}^{k}=r\unicode[STIX]{x0394}x_{i}^{k+1}$ , where $r$ is the stretching factor and $k$ is the index of the grid point, is employed in the horizontal ( $x_{1}$ ) and vertical ( $x_{3}$ ) directions to increase resolution in the breaking region. Grid spacing in the lateral ( $x_{2}$ ) direction is constant. In the horizontal direction, stretching is applied from $x_{1}=0$ to $L_{s}$ with $r=1.01$ , concentrating 1024 points (approximately 90 % of $n_{1}$ ) into the sloping region. The grid is also stretched slightly in $x_{1}$ in the sloping region to maintain orthogonality. In the vertical direction, stretching is performed twice. First, grid points are concentrated toward the bottom with $r=1.02$ in order to resolve the flow near the wall. Second, 100 of 128 points (78 %) are concentrated in the bottom 35 % of the domain with $r=1.07$ in order to resolve the pycnocline and lower layer during breaking. The resulting grid resolution at the pycnocline in the breaking region is approximately $\unicode[STIX]{x0394}x_{1}\times \unicode[STIX]{x0394}x_{2}\times \unicode[STIX]{x0394}x_{3}=4~\text{mm}\times 1~\text{mm}\times 2~\text{mm}$ . Near the bottom wall, the vertical coordinate is measured in wall units $x_{3}^{+}=x_{3}/\unicode[STIX]{x1D6FF}_{S}$ , where $\unicode[STIX]{x1D6FF}_{S}=\sqrt{2\unicode[STIX]{x1D708}/\unicode[STIX]{x1D714}}$ (see § 2.2 for the definition of the wave frequency $\unicode[STIX]{x1D714}$ ). The vertical resolution near the wall is therefore $\unicode[STIX]{x0394}x_{3}^{+}=\unicode[STIX]{x0394}x_{3}/\unicode[STIX]{x1D6FF}_{S}<1$ . In the worst case scenario, the grid spacing is found to be within approximately one order of magnitude of the Kolmogorov length scale $\unicode[STIX]{x1D702}_{k}$ , providing sufficient resolution for DNS (see Arthur & Fringer Reference Arthur and Fringer2014). The boundary conditions for all velocity components are no slip on the bottom wall and free slip on the top, left and right walls. The density field has a gradient-free boundary condition on the top, bottom, left and right walls. All variables are periodic in the lateral ( $x_{2}$ ) direction.

A time step of $\unicode[STIX]{x0394}t=0.003$ s is used for all simulations. The number of time steps varies by case and ranges from 18 000 for case 1 to for 27 000 for case 8. Simulations were run on the US Army Research Lab Department of Defense Supercomputing Research Center (ARL DSRC) supercomputer Pershing using 432 processors. With a computation time of approximately 20 s per time step, this resulted in wall-clock simulation times ranging from 100 h (43 000 processor h; case 1) to 150 h (65 000 processor h; case 8).

Table 1.

The breaking wave cases considered in this study in terms of the pycnocline thickness $\unicode[STIX]{x1D6FF}$ , the wave amplitude $a$ , the wave speed $c$ and the wave length scale $L_{w}$ , the normalized interface thickness $k\unicode[STIX]{x1D6FF}$ , the wave steepness $ka$ and the wave Reynolds number  $Re_{w}$ .

2.2 Initial wave parameters

We consider eight breaking wave cases with varying pycnocline thickness $\unicode[STIX]{x1D6FF}$ , as summarized in table 1. The stratification within the domain is initialized as

where the reference density $\unicode[STIX]{x1D70C}_{0}=1000~\text{kg}~\text{m}^{-3}$ , the density difference between the top and bottom layers is $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1}=30~\text{kg}~\text{m}^{-3}$ , the upper-layer depth $h_{1}=0.3$  m, and $\unicode[STIX]{x1D6FC}=0.99$ . The pycnocline thickness $\unicode[STIX]{x1D6FF}$ determines the buoyancy frequency $N=\sqrt{-(g/\unicode[STIX]{x1D70C}_{0})\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}/\unicode[STIX]{x2202}x_{3}}$ of the initial stratification; the maximum buoyancy frequency $N_{max}$ scales with $\unicode[STIX]{x1D6FF}^{-1/2}$ (Troy & Koseff Reference Troy and Koseff2005), as is shown for each case in figure 1(c). Throughout the paper, we refer to time using the variable $\unicode[STIX]{x1D70F}$ in order to avoid confusion with $t$ , which stands for ‘turbulent’ in the notation of Scotti & White (Reference Scotti and White2014).

A wave of depression is created in the domain by initializing a half-Gaussian at the left boundary. The shape of the interface is given by

where $a_{0}=0.1~\text{m}$ and $L_{0}=0.7~\text{m}$ are the initial amplitude and length scale of the half-Gaussian, respectively. An initial perturbation $\unicode[STIX]{x1D701}^{\prime }R$ , where $\unicode[STIX]{x1D701}^{\prime }=1~\text{mm}$ and $R\in \{-1,1\}$ is a uniformly distributed random number, is also added to the interface to trigger lateral instabilities during breaking. Relaxation of the initial half-Gaussian produces an internal wave of depression that propagates toward the slope and breaks. We note that the initial wave is solitary-like, but that it does not necessarily satisfy the Dubreil–Jacotin–Long (DJL) equation for internal solitary waves (see Lamb Reference Lamb2002, and references therein). We therefore refer to the waves in this study more generally as ‘internal waves’. Further discussion of the difference between the waves in the present study and true DJL internal solitary waves can be found in Arthur & Fringer (Reference Arthur and Fringer2016).

Here, case 1 corresponds to case 3 from Arthur & Fringer (Reference Arthur and Fringer2014) and cases 3t1/3t2 from Arthur & Fringer (Reference Arthur and Fringer2016). The additional cases 2–8 differ from this ‘base case’ only in the value of $\unicode[STIX]{x1D6FF}$ ; all other initial parameters are held constant. As a result of changing $\unicode[STIX]{x1D6FF}$ , the other properties of the incoming wave (after it has formed from the initial half-Gaussian) change slightly relative to the base case (see table 1). These properties, including the wave amplitude $a$ , the wave speed $c$ and the wave length scale $L_{w}=1/a\int _{0}^{L_{s}}\unicode[STIX]{x1D702}(x_{1})\,\text{d}x_{1}$ , where $\unicode[STIX]{x1D702}$ is the displacement of the $\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{0}$ isopycnal from $x_{3}=-h_{1}$ (Michallet & Ivey Reference Michallet and Ivey1999), are calculated after the wave has formed but before it begins to interact with the slope.

In order to quantify the effect of changing the pycnocline thickness $\unicode[STIX]{x1D6FF}$ on the incoming wave, we classify each breaking wave case in terms of the relevant non-dimensional parameters, including $k\unicode[STIX]{x1D6FF}$ , $ka$ and $Re_{w}$ (Troy & Koseff Reference Troy and Koseff2005), as shown in table 1: $k\unicode[STIX]{x1D6FF}$ is the normalized pycnocline thickness, where the horizontal wavenumber $k=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D706}$ and $\unicode[STIX]{x1D706}\approx 2L_{w}$ ; $ka$ is the wave steepness; $Re_{w}=a^{2}\unicode[STIX]{x1D714}/\unicode[STIX]{x1D708}$ is the wave Reynolds number, where the wave frequency $\unicode[STIX]{x1D714}=c/\unicode[STIX]{x1D706}$ . We note that the value of $Re_{w}=146$ reported here for case 1 differs from the value of $Re_{w}=208$ reported for the same case in Arthur & Fringer (Reference Arthur and Fringer2014) and Arthur & Fringer (Reference Arthur and Fringer2016). Here, we use the actual wave speed $c$ as opposed to the linear phase speed in a two-layer system that was used previously. As shown in table 1, the normalized pycnocline thickness varies by nearly two orders of magnitude from case 1 to case 8. However, $ka$ and $Re_{w}$ change by no more than 20 % from the base case. The bathymetric slope $s$ is also held constant. We therefore focus our analysis on how the energetics of wave breaking change with the normalized pycnocline thickness $k\unicode[STIX]{x1D6FF}$ . Troy & Koseff (Reference Troy and Koseff2005) considered the combined effect of the pycnocline thickness $k\unicode[STIX]{x1D6FF}$ and the wave steepness $ka$ using the wave Richardson number $Ri_{w}\sim k\unicode[STIX]{x1D6FF}/(ka)^{2}$ . However, we do not consider $Ri_{w}$ here because $ka$ is roughly constant and therefore $Ri_{w}\sim k\unicode[STIX]{x1D6FF}$ .

4.1 Local definitions of turbulent dissipation, mixing and mixing efficiency

The local mixing efficiency of stratified turbulent flows has traditionally been quantified using the flux Richardson number $R_{f}={\mathcal{B}}/({\mathcal{B}}+\unicode[STIX]{x1D716})$ (1.1), where ${\mathcal{B}}=(g/\unicode[STIX]{x1D70C}_{0})\overline{\unicode[STIX]{x1D70C}^{\prime }u_{3}^{\prime }}$ is the turbulent buoyancy flux and $\unicode[STIX]{x1D716}=2\unicode[STIX]{x1D708}\overline{\unicode[STIX]{x1D61A}_{ij}^{\prime }\unicode[STIX]{x1D61A}_{ij}^{\prime }}$ is the turbulent dissipation (e.g. Osborn Reference Osborn1980; Ivey & Imberger Reference Ivey and Imberger1991; Shih et al. Reference Shih, Koseff, Ivey and Ferziger2005). The turbulent rate of strain tensor is $\unicode[STIX]{x1D61A}_{ij}^{\prime }=((\unicode[STIX]{x2202}u_{i}^{\prime }/\unicode[STIX]{x2202}x_{j})+(\unicode[STIX]{x2202}u_{j}^{\prime }/\unicode[STIX]{x2202}x_{i}))/2$ . The overbar denotes a Reynolds average, while the prime denotes a departure from that average. As discussed by Venayagamoorthy & Koseff (Reference Venayagamoorthy and Koseff2016), the turbulent buoyancy flux ${\mathcal{B}}$ can include reversible countergradient fluxes, leading to low or negative values of $R_{f}$ , especially in weakly turbulent conditions. Interpretation of $R_{f}$ using this method is further complicated in the presence of unsteady, inhomogeneous turbulence due to the difficulty in choosing the appropriate Reynolds average (e.g. Hult et al. Reference Hult, Troy and Koseff2011). Thus, even in the current DNS study, it is difficult to interpret the meaning of $R_{f}$ .

To address the problems associated with quantifying turbulent dissipation, mixing and mixing efficiency in regions of unsteady, inhomogeneous turbulence, Scotti & White (Reference Scotti and White2014) developed a new method of analysis based on a locally defined background potential energy. Their method builds upon the background potential energy framework of Winters et al. (Reference Winters, Lombard, Riley and D’Asaro1995) to define local turbulent dissipation and irreversible turbulent mixing quantities. Turbulent dissipation is defined as above,

where we have adopted the notation of Scotti & White (Reference Scotti and White2014) for clarity. Irreversible turbulent mixing is defined as the turbulent dissipation of available potential energy,

where $b(x_{1},x_{2},x_{3},\unicode[STIX]{x1D70F})=g(\unicode[STIX]{x1D70C}-\unicode[STIX]{x1D70C}_{0})/\unicode[STIX]{x1D70C}_{0}$ is the buoyancy field and $N_{\ast }^{2}(x_{3},\unicode[STIX]{x1D70F})$ is the squared buoyancy frequency of the background density field $\unicode[STIX]{x1D70C}_{\ast }(x_{3},\unicode[STIX]{x1D70F})$ . The background density field represents the lowest possible potential energy state of the system if it were to be adiabatically rearranged (Winters et al. Reference Winters, Lombard, Riley and D’Asaro1995), and is obtained numerically by sorting the density field $\unicode[STIX]{x1D70C}$ . Changes in the background potential energy can only occur due to molecular diffusion, and are therefore irreversible. In the present study, the Reynolds average overbar denotes a lateral average (in the $x_{2}$ direction), while the prime denotes a departure from that average. Our calculations of $\unicode[STIX]{x1D716}_{k}^{t}$ and $\unicode[STIX]{x1D716}_{p}^{t}$ are therefore functions of $x_{1}$ , $x_{3}$ and $\unicode[STIX]{x1D70F}$ . In (4.2), $N_{\ast }^{2}(b)$ is found by mapping the given local value of $b$ to a corresponding value of $N_{\ast }^{2}$ in the background buoyancy frequency profile and is therefore a function of $x_{1}$ , $x_{2}$ , $x_{3}$ and $\unicode[STIX]{x1D70F}$ . Alternatively, $N_{\ast }^{2}(\overline{b})$ must be found by interpolating $N_{\ast }^{2}$ to the given value of $\overline{b}$ , and is therefore a function of $x_{1}$ , $x_{3}$ and  $\unicode[STIX]{x1D70F}$ .

Figure 6.

(a) The background density field $\unicode[STIX]{x1D70C}_{\ast }/\unicode[STIX]{x1D70C}_{0}$ and (b) the background buoyancy frequency field $N_{\ast }^{2}$ used in the calculation of $\unicode[STIX]{x1D716}_{p}^{t}$ . Data are shown for case 5 ( $k\unicode[STIX]{x1D6FF}\approx 1$ ) at the time of maximum turbulent mixing $(\unicode[STIX]{x1D716}_{p}^{t})_{max}^{V}$ (5.2), which corresponds to the time shown in figures 2(e), 4 and 5(m). The background buoyancy frequency is shown at $\unicode[STIX]{x1D70F}=0$ and $\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D70F}_{f}$ in (b) as well. The horizontal dotted line denotes the vertical location of the maximum background buoyancy frequency $N_{\ast }^{2}$ at $\unicode[STIX]{x1D70F}/T_{1}=0.36$ , while the horizontal dash-dot line denotes the vertical location of the initial pycnocline $x_{3}=-h_{1}$ .

Irreversible mixing leads to a decrease in the maximum background buoyancy frequency over time. This is shown in figure 6, which presents example calculations of the background density field $\unicode[STIX]{x1D70C}_{\ast }$ and background buoyancy frequency field $N_{\ast }^{2}$ for case 5. The background buoyancy frequency is shown for the initial condition $\unicode[STIX]{x1D70F}=0$ , during the upslope surge at $\unicode[STIX]{x1D70F}/T_{1}=0.36$ (the time shown in figures 2 e, 4 and 5 m), and at the end of the breaking event $\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D70F}_{f}$ . The background density field is only shown for the intermediate time $\unicode[STIX]{x1D70F}/T_{1}=0.36$ because the change from $\unicode[STIX]{x1D70F}=0$ to $\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D70F}_{f}$ is not discernible in the figure. The horizontal dotted line denotes the vertical location of the maximum buoyancy frequency at $\unicode[STIX]{x1D70F}/T_{1}=0.36$ , which is slightly below $x_{3}=-h_{1}$ because there is initially less lower-layer fluid ( $\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{2}$ ) than upper-layer fluid ( $\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{1}$ ) in the domain due to the sloping bottom. Layering is also evident in the background density field (figure 6 a), along with the corresponding wiggles in $N_{\ast }^{2}$ (figure 6 b).

Using the definitions of turbulent dissipation and irreversible turbulent mixing in (4.1) and (4.2) above, the local irreversible flux Richardson number can be calculated with

This measure of the flux Richardson number is not affected by countergradient buoyancy fluxes, and is therefore a preferable measure of the local mixing efficiency to $R_{f}$ (Scotti & White Reference Scotti and White2014; Venayagamoorthy & Koseff Reference Venayagamoorthy and Koseff2016).

Figure 7.

Example local turbulent energetics calculations for case 5 ( $k\unicode[STIX]{x1D6FF}\approx 1$ ) at the time of maximum turbulent mixing $(\unicode[STIX]{x1D716}_{p}^{t})_{max}^{V}$ ((5.2), corresponding to the time shown in figures 2(e), 4 and 5(m)). Shown are (a) the density field $\overline{\unicode[STIX]{x1D70C}}/\unicode[STIX]{x1D70C}_{0}$ , (b) the turbulent dissipation $\unicode[STIX]{x1D716}_{k}^{t}$ , (c) the turbulent mixing $\unicode[STIX]{x1D716}_{p}^{t}$ , and (d) the local irreversible flux Richardson number $R_{f}^{\ast }$ . In (b) and (c), the turbulent quantities are normalized by the maximum turbulent dissipation at the time shown, $\unicode[STIX]{x1D716}_{k,max}^{t}$ . Note that in (d), the $R_{f}^{\ast }$ field is restricted to locations where $Re_{L}>10$ for clarity ( $Re_{L}$ is defined in (4.7)), while (a), (b) and (c) include the full field.

Local turbulent dissipation and irreversible mixing quantities are calculated over time for each breaking wave case. For illustrative purposes, we consider these quantities for case 5 during the upslope surge of dense fluid after breaking, when shear instabilities (billows) are present at the interface (figure 7 a). At this time, turbulent dissipation $\unicode[STIX]{x1D716}_{k}^{t}$ and mixing $\unicode[STIX]{x1D716}_{p}^{t}$ are elevated throughout the upslope surge (figure 7 b,c). However, mixing is more concentrated along the roll-up regions of the billows and the nose of the upslope surge (figure 7 c), where sharp density gradients are present. A similar result was noted by Smyth et al. (Reference Smyth, Moum and Caldwell2001), who used the dissipation rate of scalar variance $\unicode[STIX]{x1D712}=2\unicode[STIX]{x1D705}\overline{|\unicode[STIX]{x1D6FB}\unicode[STIX]{x1D70C}^{\prime }|^{2}}$ to examine turbulent mixing during the breakdown of Kelvin–Helmholtz billows. The quantity $\unicode[STIX]{x1D712}$ is related to $\unicode[STIX]{x1D716}_{p}^{t}$ by the background buoyancy frequency $N_{\ast }^{2}$ , viz.

The local irreversible flux Richardson number $R_{f}^{\ast }$ (figure 7 d) generally follows $\unicode[STIX]{x1D716}_{p}^{t}$ . It approaches 1 in the billow roll-up regions and at the nose of the upslope surge, where $\unicode[STIX]{x1D716}_{p}^{t}$ is elevated, and approaches 0 within the billow cores and throughout the rest of the turbulent region, where $\unicode[STIX]{x1D716}_{p}^{t}$ is smaller relative to $\unicode[STIX]{x1D716}_{k}^{t}$ . Although $R_{f}^{\ast }$ clearly varies with the local characteristics of the flow, the bulk (spatially and time-averaged) $R_{f}^{\ast }$ is found to be roughly 0.3 for each case. We therefore present an analysis of $R_{f}^{\ast }$ using the parameter space of Mater & Venayagamoorthy (Reference Mater and Venayagamoorthy2014) in the next section.

4.2 Local mixing efficiency analysis

Mater & Venayagamoorthy (Reference Mater and Venayagamoorthy2014) defined three non-dimensional parameters that quantify the effects of shear, stratification and viscosity in stratified turbulence. They include a shear strength parameter $ST_{L}$ , a turbulent Froude number $Fr_{k}$ and a turbulent Reynolds number $Re_{L}$ , as presented in § 1. Here, we define these parameters as

where $\overline{S}=|\unicode[STIX]{x2202}\overline{u}_{1}/\unicode[STIX]{x2202}x_{3}|$ is the absolute value of the laterally averaged vertical shear and $k^{t}=\overline{(u_{i}^{\prime }u_{i}^{\prime })}/2$ is the turbulent kinetic energy,

where $N_{\ast }(\overline{b})$ is the buoyancy frequency of the background density field for the local value of $\overline{b}$ , as in (4.2), and

Figure 8.

$R_{f}^{\ast }$ versus $Fr_{k}$ results for breaking wave cases 1–8. Data are restricted to $4<ST_{L}<9$ and $100<Re_{L}<400$ . (a) The overall mean irreversible flux Richardson number $\langle R_{f}^{\ast }\rangle _{all}$ along with $\langle R_{f}^{\ast }\rangle _{c}$ for cases 1–8. (b) The area-weighted fraction of total turbulent mixing as a function of $Fr_{k}$ (i.e. a binned sum of $\unicode[STIX]{x1D716}_{p}^{t}A_{cell}W/(\unicode[STIX]{x1D716}_{p}^{t})_{tot}^{V}$ in each $Fr_{k}$ bin, where $A_{cell}$ is the ( $x_{1}$ , $x_{3}$ ) area of the corresponding grid cell and $W$ is the $x_{2}$ width of the domain), for cases 1, 3, 5 and 8. The quantity $(\unicode[STIX]{x1D716}_{p}^{t})_{tot}^{V}$ used in (b) is the total volume-integrated turbulent mixing (defined in (5.4)).

Because we are interested in the effect of stratification on the local turbulent mixing efficiency, we consider how $R_{f}^{\ast }$ varies with the turbulent Froude number $Fr_{k}$ . For each breaking wave case, we calculate a mean irreversible flux Richardson number $\langle R_{f}^{\ast }\rangle _{c}$ in bins of $Fr_{k}$ (figure 8 a). $\langle R_{f}^{\ast }\rangle _{c}$ reaches a local maximum between 0.2 and 0.45 for $0.5<Fr_{k}<0.9$ . We also define an overall mean irreversible flux Richardson number $\langle R_{f}^{\ast }\rangle _{all}$ as a function of $Fr_{k}$ by calculating the mean $R_{f}^{\ast }$ value among all cases in each $Fr_{k}$ bin. The $\langle R_{f}^{\ast }\rangle _{all}$ versus $Fr_{k}$ curve captures the general trends of the $\langle R_{f}^{\ast }\rangle _{c}$ versus $Fr_{k}$ results (figure 8 a). It reaches a local maximum of approximately 0.36 for $Fr_{k}\approx 0.7$ , where the balance of turbulence and stratification leads to the optimal mixing efficiency. $\langle R_{f}^{\ast }\rangle _{all}$ decreases as $Fr_{k}$ increases above 0.7 because there is less stratification to be mixed. As $Fr_{k}$ decreases below 0.7 until roughly 0.1, $\langle R_{f}^{\ast }\rangle _{all}$ decreases because the increased stratification suppresses turbulence. For $Fr_{k}<0.1$ , $\langle R_{f}^{\ast }\rangle _{all}$ increases as $Fr_{k}$ decreases because $\unicode[STIX]{x1D716}_{k}^{t}$ approaches zero. Despite the relatively large values of $\langle R_{f}^{\ast }\rangle _{c}$ for $Fr_{k}<0.1$ , only a small fraction of the total mixing occurs in these regions. This is evidenced by the result in figure 8(b), which shows the total irreversible mixing $\unicode[STIX]{x1D716}_{p}^{t}$ that occurs within each $Fr_{k}$ bin for cases 1, 3, 5 and 8.

We note that the data in this analysis are restricted to $4<ST_{L}<9$ and $100<Re_{L}<400$ , which is where the most total mixing occurs over the course of each breaking wave event. The $\langle R_{f}^{\ast }\rangle _{all}$ versus $Fr_{k}$ curve in figure 8(b) shows similarities to $R_{f}$ and $R_{f}^{\ast }$ versus $Fr_{k}$ relationships found in previous studies (Shih Reference Shih2003; Mater & Venayagamoorthy Reference Mater and Venayagamoorthy2014), particularly for $Fr_{k}>0.1$ . However, a complete analysis of $R_{f}^{\ast }$ in the ( $ST_{L}$ , $Fr_{k}$ , $Re_{L}$ ) parameter space, which would facilitate direct comparisons to these studies, is left to future work.

5.1 Volume-integrated turbulent dissipation and mixing

By volume-integrating local dissipation and mixing quantities during wave breaking, we examine the changes in turbulent energetics as a function of the initial interface thickness $k\unicode[STIX]{x1D6FF}$ . We define the volume-integrated turbulent mixing as

where the volume of integration $V_{s}$ is defined as the sloping region of the domain ( $L_{s}\leqslant x_{1}\leqslant L$ , see figure 1) in order to better isolate the effects of breaking (Arthur & Fringer Reference Arthur and Fringer2014). Maximum, cumulative and total volume-integrated turbulent mixing quantities are defined as

and

respectively. Similar quantities are defined for turbulent dissipation by replacing $\unicode[STIX]{x1D716}_{p}^{t}$ with $\unicode[STIX]{x1D716}_{k}^{t}$ in (5.1)–(5.4). Turbulent dissipation and mixing both begin to rise after wave breaking occurs, and continue until the end of the breaking event $\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D70F}_{f}$ (figure 9 a,b, points d–i for case 5). Note that the volume-integrated turbulent mixing efficiency measures in figure 9(c,d) are defined and discussed in § 5.3.

Figure 9.

Volume-integrated turbulent energetics results for breaking wave cases 1–8. Shown are (a) cumulative dissipation $(\unicode[STIX]{x1D716}_{k}^{t})_{cum}^{V}$ , (b) cumulative mixing $(\unicode[STIX]{x1D716}_{p}^{t})_{cum}^{V}$ , (c) cumulative mixing efficiency $\unicode[STIX]{x1D702}_{cum}^{t}$ and (d) instantaneous mixing efficiency $\unicode[STIX]{x1D702}^{t}$ over time for cases 1, 3, 5 and 8. In (a) and (b), $(\unicode[STIX]{x1D716}_{k}^{t})_{cum}^{V}$ and $(\unicode[STIX]{x1D716}_{p}^{t})_{cum}^{V}$ are normalized by the total time-integrated values $(\unicode[STIX]{x1D716}_{k}^{t})_{tot}^{V}$ and $(\unicode[STIX]{x1D716}_{p}^{t})_{tot}^{V}$ , respectively. The critical $R_{f}$ value of 0.17 is shown by a horizontal dotted line in (c) and (d). Labels a–i in (a–d) correspond to the times shown in figure 2 for case 5. Also shown are (e) maximum and (f) total (time-integrated) values of dissipation and mixing as a function of $k\unicode[STIX]{x1D6FF}$ for all cases. Values in (e) are normalized by the maximum instantaneous turbulent dissipation for all cases $(\unicode[STIX]{x1D716}_{k}^{t})_{max}^{V,all}$ , while values in (f) are normalized by the maximum total turbulent dissipation for all cases $((\unicode[STIX]{x1D716}_{k}^{t})_{tot}^{V})_{max}^{all}$ . Vertical dotted lines in (e) and (f) represent $k\unicode[STIX]{x1D6FF}$ regimes I and II.

The variation in volume-integrated turbulent dissipation and mixing with $k\unicode[STIX]{x1D6FF}$ can be summarized by examining the maximum and total values of these quantities, as shown in figures 9(e) and 9(f), respectively. Here, we define two regimes of $k\unicode[STIX]{x1D6FF}$ based on the peaked shape of $(\unicode[STIX]{x1D716}_{p}^{t})_{max}^{V}$ (figure 9 e) and the asymptotic shape of $(\unicode[STIX]{x1D716}_{p}^{t})_{tot}^{V}$ (figure 9 f). Regime I is defined for $k\unicode[STIX]{x1D6FF}<1$ , where both $(\unicode[STIX]{x1D716}_{p}^{t})_{max}^{V}$ and $(\unicode[STIX]{x1D716}_{p}^{t})_{tot}^{V}$ increase monotonically. Regime II is defined for $k\unicode[STIX]{x1D6FF}\geqslant 1$ , where $(\unicode[STIX]{x1D716}_{p}^{t})_{max}^{V}$ decreases monotonically and $(\unicode[STIX]{x1D716}_{p}^{t})_{tot}^{V}$ asymptotes. The maximum volume-integrated turbulent dissipation $(\unicode[STIX]{x1D716}_{k}^{t})_{max}^{V}$ does not follow the peaked shape of $(\unicode[STIX]{x1D716}_{p}^{t})_{max}^{V}$ . Rather, $(\unicode[STIX]{x1D716}_{k}^{t})_{max}^{V}$ increases with $k\unicode[STIX]{x1D6FF}$ for low $k\unicode[STIX]{x1D6FF}$ values in regime I, and is relatively scattered as $k\unicode[STIX]{x1D6FF}$ increases into regime II. The total volume-integrated turbulent dissipation $(\unicode[STIX]{x1D716}_{k}^{t})_{tot}^{V}$ , however, increases monotonically with $k\unicode[STIX]{x1D6FF}$ through both regimes I and II.

5.2 Turbulent contribution to dissipation and mixing

Due to the computational restrictions associated with DNS, the Reynolds number $Re_{w}$ of our simulations is limited. In order to contextualize our results, it is therefore important to quantify the turbulent contribution to the full energy budget, which includes both laminar and turbulent effects. In this section, we use volume-integrated dissipation and mixing quantities to consider the turbulent contribution to the full dissipation and mixing that occurs during wave breaking. Scotti & White (Reference Scotti and White2014) define the local dissipation of the full kinetic energy as

where $\unicode[STIX]{x1D61A}_{ij}=((\unicode[STIX]{x2202}u_{i}/\unicode[STIX]{x2202}x_{j})+(\unicode[STIX]{x2202}u_{j}/\unicode[STIX]{x2202}x_{i}))/2$ is the rate of strain tensor, see (1.1) therein. Thus, we define the full volume-integrated dissipation as

where the volume of integration $V_{s}$ is defined as in § 5.1. We define the full volume-integrated mixing as in Winters et al. (Reference Winters, Lombard, Riley and D’Asaro1995),

where $\unicode[STIX]{x1D719}_{d}$ is the volume-integrated rate of change in the background potential energy and $\unicode[STIX]{x1D719}_{i}$ is the irreversible rate of change from internal to potential energy. Here, $\unicode[STIX]{x1D70C}_{top}$ and $\unicode[STIX]{x1D70C}_{bottom}$ are the densities on the top and bottom boundaries of the domain, respectively, and $A$ is the corresponding surface area of the boundary. We note that the full volume-integrated dissipation ${\mathcal{D}}$ and mixing ${\mathcal{M}}$ measures defined here are equal to those used previously by Arthur & Fringer (Reference Arthur and Fringer2014). The dissipation $\unicode[STIX]{x1D716}$ in Arthur & Fringer (Reference Arthur and Fringer2014) is equal to ${\mathcal{D}}$ and the effective irreversible mixing $\unicode[STIX]{x1D719}_{e}$ in Arthur & Fringer (Reference Arthur and Fringer2014) is equal to ${\mathcal{M}}$ . The values ${\mathcal{D}}$ and ${\mathcal{M}}$ include both laminar and turbulent effects, while $(\unicode[STIX]{x1D716}_{k}^{t})^{V}$ and $(\unicode[STIX]{x1D716}_{p}^{t})^{V}$ (5.1) include only turbulent effects.

Figure 10.

The turbulent contribution to dissipation and mixing. Volume-integrated (a) dissipation ( $(\unicode[STIX]{x1D716}_{k}^{t})^{V}$ , turbulent and ${\mathcal{D}}$ , full) and (b) mixing ( $(\unicode[STIX]{x1D716}_{p}^{t})^{V}$ , turbulent and ${\mathcal{M}}$ , full) over time for breaking wave case 5 ( $k\unicode[STIX]{x1D6FF}\approx 1$ ). Values in (a) and (b) are normalized by ${\mathcal{D}}_{max}$ , the maximum full dissipation for case 5. The fraction $f_{turb}$ of both dissipation and mixing (c) over the entire breaking wave event (‘Total’) and (d) during the turbulent portion of the event (‘Post-turbulent’), as a function of $k\unicode[STIX]{x1D6FF}$ for all cases. The vertical dotted lines in (c) and (d) represent $k\unicode[STIX]{x1D6FF}$ regimes I and II defined in § 5.1.

The initial peak in the full dissipation ${\mathcal{D}}$ (shown for case 5 in figure 10 a) is associated with near-bottom downslope flow in the lower layer prior to breaking (figure 2 a). As noted in Arthur & Fringer (Reference Arthur and Fringer2014), this dissipation is laminar because turbulence does not develop in the domain until wave breaking occurs. The turbulent dissipation $(\unicode[STIX]{x1D716}_{k}^{t})^{V}$ , therefore, does not begin to rise until breaking and does not include this pre-breaking laminar dissipation. The shapes of the turbulent mixing $(\unicode[STIX]{x1D716}_{p}^{t})^{V}$ and full mixing ${\mathcal{M}}$ curves (shown for case 5 in figure 10 b) are similar. This is expected because the turbulence generated by overturning during the breaking event is the primary driver of mixing in the domain. Overall, however, $(\unicode[STIX]{x1D716}_{k}^{t})^{V}$ and $(\unicode[STIX]{x1D716}_{p}^{t})^{V}$ are generally lower than ${\mathcal{D}}$ and ${\mathcal{M}}$ , respectively, throughout the breaking event.

In order to quantify the difference between the turbulent and full measures of dissipation and mixing, we define $f_{turb}$ as the fraction of either dissipation or mixing due to turbulence over the course of the breaking event. Thus $f_{turb}=(\unicode[STIX]{x1D716}_{k}^{t})_{tot}^{V}/{\mathcal{D}}_{tot}$ for dissipation and $f_{turb}=(\unicode[STIX]{x1D716}_{p}^{t})_{tot}^{V}/{\mathcal{M}}_{tot}$ for mixing, where ${\mathcal{D}}_{tot}$ and ${\mathcal{M}}_{tot}$ are the time-integrated values of ${\mathcal{D}}$ and ${\mathcal{M}}$ , respectively, defined as in (5.4) for $(\unicode[STIX]{x1D716}_{p}^{t})_{tot}^{V}$ . The turbulent contribution to dissipation (figure 10 c) increases within regime I from approximately 17 % for case 1 to approximately 35 % for case 5 ( $k\unicode[STIX]{x1D6FF}\approx 1$ ). Within regime II, $f_{turb}$ remains relatively constant. The turbulent contribution to mixing (figure 10 c) is approximately 34 % for case 1, but is then roughly constant between 40 %–50 % for larger $k\unicode[STIX]{x1D6FF}$ . The onset of turbulence in the domain (when $(\unicode[STIX]{x1D716}_{k}^{t})^{V}$ and $(\unicode[STIX]{x1D716}_{p}^{t})^{V}$ begin to rise) is indicated by a dot in figure 10(a,b). Another calculation of $f_{turb}$ , this time over only the post-turbulent portion of the event (as opposed to the entire event) is shown in figure 10(d). The turbulent contributions to dissipation and mixing are higher in the post-turbulent stage because the laminar pre-breaking stage is not included, reaching approximately 45 % for dissipation and 45 %–55 % for mixing. However, similar trends are seen with  $k\unicode[STIX]{x1D6FF}$ .

Overall, the turbulent contribution to dissipation and mixing is relatively low ( ${<}60$  %). This is a consequence of the low Reynolds number $Re_{w}$ used in our simulations, which is necessitated by the computational restrictions associated with DNS. Although we reach local $Re_{b}$ values of the order of $10^{2}{-}10^{3}$ (in the energetic regime of Shih et al. Reference Shih, Koseff, Ivey and Ferziger2005), these regions make up only a small fraction of the domain. Arthur & Fringer (Reference Arthur and Fringer2014) found that as $Re_{w}$ increases, the turbulent contribution to dissipation and mixing during breaking internal wave events increases because turbulence develops earlier in the event. Thus, we expect that as the Reynolds number of the cases considered in this study increases, the turbulent contribution to dissipation and mixing will also increase. The potential effect of increasing the Reynolds number on the mixing efficiency is discussed in the next section. However, further analysis at larger Reynolds number is left to future work.

5.3 Volume-integrated mixing efficiency

Using the volume-integrated calculations of turbulent dissipation and mixing above, volume-integrated measures of the turbulent mixing efficiency are made. We define a cumulative volume-integrated turbulent mixing efficiency

as well as a bulk volume-integrated turbulent mixing efficiency

Thus, $\unicode[STIX]{x1D702}_{cum}^{t}(\unicode[STIX]{x1D70F}_{f})=\unicode[STIX]{x1D702}_{B}^{t}$ . As in Arthur & Fringer (Reference Arthur and Fringer2014), we also define a full bulk mixing efficiency, which includes both laminar and turbulent effects, as

The cumulative mixing efficiency $\unicode[STIX]{x1D702}_{cum}^{t}$ (figure 9 c) is initially large due to the offset between the growth of mixing and dissipation in the early stages of turbulence. That is, as sharp density surfaces are created by turbulent strain, mixing begins to increase before turbulent dissipation occurs. This effect is discussed in Smyth et al. (Reference Smyth, Moum and Caldwell2001), who observed a similar phenomenon in stratified shear flow turbulence. After turbulence has developed, $\unicode[STIX]{x1D702}_{cum}^{t}$ decreases, approaching $\unicode[STIX]{x1D702}_{B}^{t}$ . The bulk volume-integrated turbulent mixing efficiency $\unicode[STIX]{x1D702}_{B}^{t}$ varies in the $k\unicode[STIX]{x1D6FF}$ regimes defined in § 5.1 (figure 11). $\unicode[STIX]{x1D702}_{B}^{t}$ is low (roughly 0.25) for small $k\unicode[STIX]{x1D6FF}$ (case 1). As $k\unicode[STIX]{x1D6FF}$ increases within regime I, the mixing efficiency increases, reaching a peak of 0.37 when $k\unicode[STIX]{x1D6FF}\approx 1$ (case 5). Within regime II, $\unicode[STIX]{x1D702}_{B}^{t}$ levels out around 0.35.

Figure 11.

Bulk mixing efficiency measures $\unicode[STIX]{x1D702}_{B}^{t}$ and $\unicode[STIX]{x1D702}_{B}$ as functions of $k\unicode[STIX]{x1D6FF}$ . The vertical dotted line represents $k\unicode[STIX]{x1D6FF}$ regimes I and II defined in § 5.1. Also shown are the bulk mixing efficiency results of Fringer & Street (Reference Fringer and Street2003).

The volume-integrated bulk mixing efficiency values calculated for breaking wave cases 1–8 are consistently higher than the critical value of the flux Richardson number, $R_{f}^{crit}\approx 0.17$ . However, we note that the instantaneous volume-integrated turbulent mixing efficiency, defined as

and shown in figure 9(d), does approach $R_{f}^{crit}\approx 0.17$ towards the end of the breaking event, especially for low values of $k\unicode[STIX]{x1D6FF}$ . This asymptotic behaviour occurs for each case, indicating that the turbulence after breaking reaches an equilibrium state representative of what one might expect in turbulence that is predominantly shear driven. A similar trend was observed by Smyth et al. (Reference Smyth, Moum and Caldwell2001) for the flux coefficient, which is analogous to the mixing efficiency.

Table 2.

Mixing efficiency values from previous studies.

For further comparison to the present results, table 2 shows mixing efficiency values that have been calculated for various flows in previous studies. It is similar to table 1 in Fringer & Street (Reference Fringer and Street2003), but with the addition of several more recent studies. The range of bulk turbulent mixing efficiency values calculated here ( $\unicode[STIX]{x1D702}_{B}^{t}=0.25-0.37$ ) is on the upper end of the values reported in other studies of breaking internal waves on slopes (Michallet & Ivey Reference Michallet and Ivey1999; Arthur & Fringer Reference Arthur and Fringer2014). This may be because we consider only the turbulent contribution to the bulk mixing efficiency, as discussed in § 5.2. For the relatively low Reynolds number in our simulations, laminar dissipation occurs due to strong downslope flow near the bottom boundary prior to wave breaking. This laminar dissipation reduces the overall mixing efficiency of the breaking event $\unicode[STIX]{x1D702}_{B}$ (see figure 11), but is not captured in $\unicode[STIX]{x1D702}_{B}^{t}$ , which includes only turbulent effects. At higher Reynolds number, we would expect turbulence to develop earlier in the simulation (i.e. during the downslope flow of lower-layer fluid before breaking), increasing turbulent dissipation and thus lowering the turbulent mixing efficiency. Therefore, at increased Reynolds number, we would expect $\unicode[STIX]{x1D702}_{B}^{t}$ and $\unicode[STIX]{x1D702}_{B}$ to approach one another because the turbulent contribution to dissipation and mixing would overwhelm the laminar contribution.

The bulk mixing efficiency values calculated here are also larger than those presented by Hult et al. (Reference Hult, Troy and Koseff2011) for breaking internal waves over a ridge. This is due, in part, to the presence of turbulence in the unstratified bottom boundary layer of their domain. Conversely, in the present breaking wave cases, the stratified region intersects the bottom boundary layer on the slope, allowing for near-bottom mixing. Because the method of Scotti & White (Reference Scotti and White2014) defines turbulent dissipation and mixing locally, $\unicode[STIX]{x1D716}_{k}^{t}$ and $\unicode[STIX]{x1D716}_{p}^{t}$ can be integrated over the bottom boundary and interior regions separately. We define the bottom boundary layer as in § 2.1 as $x_{3}<-d+\unicode[STIX]{x1D6FF}_{S}$ , where $d$ is the local depth, and the interior as the remainder of the domain. Although we find that 80 %–90 % of the total dissipation and mixing occur in the interior region, the bulk turbulent mixing efficiencies in both layers are similar to the overall value $\unicode[STIX]{x1D702}_{B}^{t}$ (see Arthur Reference Arthur2015).

The present bulk mixing efficiency values are lower than those calculated by Fringer & Street (Reference Fringer and Street2003) for breaking progressive internal waves in deep water because of the absence of bottom boundary-layer dissipation in their domain. Fringer & Street (Reference Fringer and Street2003) also found a different shape of the bulk mixing efficiency as a function of $k\unicode[STIX]{x1D6FF}$ (see figure 11). In their study, $ka$ is relatively large; this allows $k\unicode[STIX]{x1D6FF}$ to control the mixing efficiency because billows are not restricted in size by $ka$ . Thus, weak convective breaking occurs for large $k\unicode[STIX]{x1D6FF}$ , causing the mixing efficiency to decrease. In our simulations, $ka$ is relatively small. Thus, the mixing efficiency is controlled by $ka$ for large $k\unicode[STIX]{x1D6FF}$ (regime II), and is relatively constant. For further discussion of the billow scale in our breaking wave cases, see § 6 below.

An important consideration when comparing the present mixing efficiency calculations to those of other flows, especially in the laboratory or the field, is the Prandtl number $Pr=\unicode[STIX]{x1D708}/\unicode[STIX]{x1D705}$ . In water, $Pr=7$ for temperature and $Pr=700$ for salt. As $Pr$ increases ( $\unicode[STIX]{x1D705}$ decreases), one would expect the mixing efficiency to decrease because the time scale necessary for mixing to occur increases. However, $Pr=1$ in the current simulations due to the computational restrictions associated with resolving mixing for smaller values of $\unicode[STIX]{x1D705}$ . This provides another potential explanation for the relatively high mixing efficiencies in the present study.

5.4 Numerical effects

Numerical effects also play a role in the energetics, and are a function of the chosen numerical methods and discretization. There are, in fact, three sources of dissipation in our simulations: that due to mean (or laminar) strains, that due to turbulent strains and that due to numerical effects. We estimate the relative magnitude of the numerical dissipation, in a domain- and time-integrated sense, by comparing the change in total energy over the course of each simulation to that captured by resolved quantities. Following Winters et al. (Reference Winters, Lombard, Riley and D’Asaro1995) and Arthur & Fringer (Reference Arthur and Fringer2014), the total energy is defined as

where $V_{d}$ is the full volume of the domain, $E_{K}$ is the kinetic energy and $E_{P}$ is the potential energy. The evolution of $E_{T}$ should satisfy

because there are no boundary fluxes (Winters et al. Reference Winters, Lombard, Riley and D’Asaro1995). Numerical dissipation can therefore be estimated as the difference between the total energy change and the resolved energy change,

As shown in figure 12(a) for case 1, the majority of the numerical dissipation occurs prior to wave breaking ( $\unicode[STIX]{x1D70F}/T_{1}<0.2$ ), as the initial wave evolves on the coarse, offshore portion of the grid.

Figure 12.

(a) Cumulatively integrated $E_{T}$ evolution terms from (5.14) over time for case 1, normalized by $({\mathcal{M}}_{tot}+{\mathcal{D}}_{tot})$ . (b) Domain- and time-integrated estimate of numerical dissipation ${\mathcal{N}}_{tot}$ relative to the resolved full dissipation and mixing $({\mathcal{M}}_{tot}+{\mathcal{D}}_{tot})$ , as a function of $k\unicode[STIX]{x1D6FF}$ .

The total domain- and time-integrated numerical dissipation can be estimated as ${\mathcal{N}}_{tot}={\mathcal{N}}_{cum}(\unicode[STIX]{x1D70F}_{f})$ , allowing us to estimate the effect of numerical dissipation on the full bulk mixing efficiency $\unicode[STIX]{x1D702}_{B}$ (5.10). We define an effective full bulk mixing efficiency as

Because $1/(1-r)=1+r+r^{2}/2+\cdots \,$ , equation (5.17) can be simplified to

For cases 2–8, the quantity $|{\mathcal{N}}_{tot}/({\mathcal{M}}_{tot}+{\mathcal{D}}_{tot})|$ varies between 0.06–0.11 (figure 12 b), indicating that $\unicode[STIX]{x1D702}_{B,eff}$ is within approximately 6 %–11 % of $\unicode[STIX]{x1D702}_{B}$ . For case 1, which has the smallest value of $k\unicode[STIX]{x1D6FF}$ , the value of ${\mathcal{N}}_{tot}/({\mathcal{M}}_{tot}+{\mathcal{D}}_{tot})$ is slightly larger in magnitude (approximately 0.2). However, as noted above in reference to figure 12(a), the majority of the numerical dissipation occurs prior to wave breaking. Although case 1 is shown in figure 12(a) because it has the largest amount of numerical dissipation, cases 2–8 show qualitatively similar results as a function of time. The overall findings of this paper related to the bulk mixing efficiency would not change despite the relatively large effect of numerical dissipation in case 1. In particular, the numerical dissipation would act to further decrease the bulk mixing efficiency for case 1, which already possesses the smallest bulk mixing efficiency among all cases.

The numerical dissipation in our simulations is due to the non-conservative energy behaviour of non-staggered grid methods, which do not allow the discrete form of the energy flux related to $pu_{i}$ (where $p$ is the pressure) to be written in conservative form (Ferziger & Perić Reference Ferziger and Perić2002). It is therefore related to the projection scheme used by Zang et al. (Reference Zang, Street and Koseff1994) and not to numerical diffusion of momentum, which is expected because momentum advection is computed with the dispersive QUICK scheme (Leonard Reference Leonard1979). Based on numerical tests (not shown), $\unicode[STIX]{x1D719}_{i}$ and ${\mathcal{D}}$ are relatively weak functions of the grid spacing and time step size, while ${\mathcal{N}}$ is a relatively strong function of these parameters. This suggests that if we had the computational resources to perform the simulations with higher resolution and a smaller time step size, the effective bulk mixing efficiency $\unicode[STIX]{x1D702}_{B,eff}$ would converge to the reported bulk mixing efficiency $\unicode[STIX]{x1D702}_{B}$ . Such higher-resolution simulations would also allow for a grid-resolution study of the dissipation spectra, a more rigorous method for assessing numerical dissipation.