2 Description of the physical vortex-ring-induced mixing experiments

We attempt to model the mixing produced by a large number of independent vortex-ring-induced mixing experiments. Physical measurements of such a system have been recently reported in Olsthoorn & Dalziel (Reference Olsthoorn and Dalziel2015). In each of those experiments, a tank was initially filled with a stable density stratification consisting of two nearly homogeneous layers with a sharp density transition between them. For the remainder of this paper, we will denote this type of stratification as a ‘continuous two-layer stratification’, with the understanding that the stratification is approximately two layer with a continuous transition region between them. A sequence of vortex rings were then generated within the tank, such that they propagated along the direction of gravity. The maximum distance below the interface that any vortex ring penetrated was small compared with the depth of the lower layer fluid, such that the bottom of the tank did not significantly affect the dynamics of the flow. As each vortex ring translated under its self-induced velocity, it displaced the isosurfaces of the density field. The perturbation to the density field resulted in the production of secondary vorticity through a baroclinic torque. This secondary vorticity was produced directly at the location of the density interface. Lawrie & Dalziel (Reference Lawrie and Dalziel2011) have previously argued that the co-location of vorticity with the peak density gradients, as was the case in these experiments, will lead to a high mixing efficiency. Further, in a recent publication, Olsthoorn & Dalziel (Reference Olsthoorn and Dalziel2017) have demonstrated that the coupling of the secondary vorticity with the impinging vortex-ring results in an instability that rapidly generates turbulence. Thus, to review, each propagating vortex ring displaces the isopycnal surfaces, which produces secondary vorticity that, through an interaction with the vortex ring, is unstable to an instability identified in Olsthoorn & Dalziel (Reference Olsthoorn and Dalziel2017). The subsequent turbulent production further enhances the stirring of the density field, generating density fluctuations down to the Kolmogorov scale. In the experiments of Olsthoorn & Dalziel (Reference Olsthoorn and Dalziel2015), the time interval between the generation of each vortex ring was sufficient to allow the fluid within the tank to become nearly quiescent (except for thermal fluctuations). By measuring the density field between a sequence of these mixing events, Olsthoorn & Dalziel (Reference Olsthoorn and Dalziel2015) quantified the mixing induced by each vortex ring.

Figure 1 presents representative snapshots of a single stratified vortex-ring experiment. Although presented slightly differently, the experiment shown in figure 1 is the same as one of those presented in figure 2 of Olsthoorn & Dalziel (Reference Olsthoorn and Dalziel2017), to which the reader is referred for details on the experimental set-up. Here, figure 1 shows the computed azimuthal vorticity field with overlaid velocity field vectors (left) and the evolution of the density field (right) within a vertical laser sheet. These equally spaced snapshots highlight the propagation of the vortex ring (figure 1 a), the displacement of the density interface (figure 1 b), followed by the production of secondary vorticity (figure 1 c), the instability of the vortex ring (figure 1 d) and the slow transition back to quiescence (figure 1 e). While this figure has been generated from a single vortex-ring experiment, each of these five steps are characteristic of the vortex-ring experiments considered here. Note that the height of the mixing region is comparable to the diameter of the impacting vortex ring.

The above description of the vortex-ring experiment’s mechanics highlights the stirring (and the associated production of strong density gradients) of the density field. This is only one component of mixing, which occurs through a combination of stirring and diffusion. In the vortex-ring experiments, the smallest-scale features of the flow were produced through the turbulent eddies. Thus, we argue, that the majority of the fine-scale stirring, and consequently the mixing, induced by the vortex ring will only occur once the flow becomes unstable to the instability discussed in Olsthoorn & Dalziel (Reference Olsthoorn and Dalziel2017). As we will see below, we construct our model such that the growth rate of the vortex-ring instability will limit the mixing rate.

Both a velocity and length scare are required in order to parameterize the vortex-ring-induced mixing. For the vortex-ring experiments, it is natural to select the vortex-ring propagation velocity $U$ as the characteristic velocity, and the vortex-ring diameter $a$ as the characteristic length scale. This paper focuses on three dimensionless parameters: the Reynolds number ( $Re$ , the ratio of inertia to viscous forces), the Richardson number ( $Ri$ , the ratio of buoyancy to advective forces) and the Schmidt number ( $Sc$ , the ratio of viscous to molecular diffusion). These are defined as

(2.1a-c ) $$\begin{eqnarray}\displaystyle & \displaystyle \mathit{Re}=\frac{Ua}{\unicode[STIX]{x1D708}},\quad \mathit{Ri}=\frac{g\left(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1}\right)}{\unicode[STIX]{x1D70C}_{1}}\frac{a}{U^{2}},\quad \mathit{Sc}=\frac{\unicode[STIX]{x1D708}}{\unicode[STIX]{x1D705}}. & \displaystyle\end{eqnarray}$$

Here, $g$ is the acceleration due to gravity, $\unicode[STIX]{x1D70C}_{1},\unicode[STIX]{x1D70C}_{2}$ are densities associated with the stratification, $\unicode[STIX]{x1D708}$ is the kinematic viscosity (here, $\unicode[STIX]{x1D708}=1\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$ ) and $\unicode[STIX]{x1D705}$ is the coefficient of mass diffusion. In this paper, we will model the vortex-ring-induced mixing produced in the physical experiments and in numerical simulations (presented in § 4). With each of these methodologies, we will ensure consistency by comparing the Reynolds, Richardson and Schmidt numbers.

Another important parameter is the kinetic energy of each vortex ring ( $\text{KE}_{ring}$ ). According to Norbury (Reference Norbury1973), the kinetic energy of a vortex ring is given as

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \text{KE}_{ring}=C_{KE}\left(\frac{1}{2}\unicode[STIX]{x1D70C}_{1}U^{2}\right)\left(\frac{4}{3}\unicode[STIX]{x03C0}\left(\frac{a}{2}\right)^{3}\right). & \displaystyle\end{eqnarray}$$

The constant $C_{KE}$ is a function of the vortex-ring aspect ratio (the ratio of the core width to ring diameter). For the vortex rings used in Olsthoorn & Dalziel (Reference Olsthoorn and Dalziel2015), a value of $C_{KE}=6.5$ has been estimated, and thus we use this value when discussing the model experiments below.

Figure 2 presents a diagram of the energy pathways resulting from the input of $\text{KE}_{ring}$ . The propagating vortex ring produces available potential energy (APE, see Winters et al. Reference Winters, Lombard, Riley and D’Asaro1995) by displacing the isopycnal surfaces from their equilibrium position that, in turn, generates the kinetic energy (KE) associated with the secondary vorticity. In general, KE will produce APE in the system (and vice versa) via a reversible buoyancy flux ( $\mathbb{H}$ ). The instability that results from the coupling of the primary vortex ring and secondary vorticity then produces turbulent kinetic energy (TKE) at a rate $\mathbb{I}$ . This TKE production is coupled to (both generates and is generated by) the APE associated with the small-scale stirring of the density field that, through diffusion, mixes the stratification at some rate $\mathbb{M}$ , increasing the background potential energy (BPE) of the system. KE (predominantly through TKE) also viscously dissipates at a rate $\unicode[STIX]{x1D716}$ , acting as a source for internal energy. Diffusion of the background density profile ( $\mathbb{D}_{\unicode[STIX]{x1D70C}}$ ) will also slowly increase the BPE of the system. There are additional partitions of energy in the system, such as internal waves, that are not specifically labelled within figure 2 as these other energy reservoirs do not significantly contribute to the dominant mixing mechanism. This description of the energy pathways is consistent with the work of Winters et al. (Reference Winters, Lombard, Riley and D’Asaro1995), who considered mixing in an oceanographic context. The total KE and APE, excluding $\text{KE}_{ring}$ , will be denoted as stirring energy. Here, stirring is a result of both the coherent (reproducible) motions and the turbulent fluctuations. We denote the energy associated with the coherent motions as ‘coherent stirring energy’, and likewise, we denote the energy associated with the turbulent fluctuations as ‘turbulent stirring energy’. We argue that the turbulent fluctuations (turbulent stirring energy) are the dominant contributor to the change in the BPE of the system. We will return to this when we discuss the model construction.

In the experiments, once the transient stirring energy in the system had sufficiently dissipated, another vortex ring was generated and the cycle repeated. This process continued until the desired number of vortex rings was generated. We will model the repetitive generation of vortex rings below.

Figure 2. Diagram of the simplified energy pathway for the vortex-ring experiments. The input of kinetic energy from the vortex ring ( $\text{KE}_{ring}$ ) will lead to an increase in the gravitational potential energy (BPE) of the system. The size of each energy reservoir does not correspond to its relative contribution.

Figure 3. Diagram of the energy pathways in the model construction. The kinetic energy of the vortex ring ( $T$ ) (a) breaks down into stirring energy ( $e$ ) (b), which subsequently mixes the density profile ( $\unicode[STIX]{x1D70C}$ ) (c) changing the background potential energy of the system (BPE).

3 Model of the vortex-ring-induced mixing experiments

The purpose of the present model is to predict the mixing produced by isolated vortex-ring-induced mixing events within a stratified fluid. To characterize this system, we model the coherent vortex-ring energy density ( $T$ ), the stirring energy density ( $e$ ) and the background density field ( $\unicode[STIX]{x1D70C}$ ). We consider horizontally averaged quantities such that each variable is only a function of a single spatial (vertical) dimension and time. This model builds upon the conceptualization introduced in figure 2. As the majority of the mixing and dissipation will result from turbulent motions, we model $e$ as a turbulent quantity. That is, the model includes the coherent stirring energy implicitly through the model breakdown parameter $\mathbb{I}_{M}$ . Figure 3 presents a cartoon of the simplified model. Note that the evolution of $e$ and $T$ are dependent on the density field creating a coupled dynamical system for $T,e$ and $\unicode[STIX]{x1D70C}$ . As discussed above, the coherent energy $T$ does not directly mix the density field ( $\unicode[STIX]{x1D70C}$ ), but acts as a propagating source for $e$ .

This model can be written as a system of three coupled differential equations:

(3.1a-c ) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x2202}_{t}T=\mathbb{A}-\mathbb{I}_{M}+\mathbb{S},\quad \unicode[STIX]{x2202}_{t}e=\mathbb{D}_{e}-\unicode[STIX]{x1D716}-g\mathbb{B}+\mathbb{I}_{M},\quad \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D70C}=-\unicode[STIX]{x2202}_{z}\mathbb{B}+\mathbb{D}_{\unicode[STIX]{x1D70C}}. & \displaystyle\end{eqnarray}$$

The non-turbulent vortex-ring energy density ( $T$ ) is produced ( $\mathbb{S}$ ), is advected ( $\mathbb{A}$ ) and will, in the presence of the stratification, feed the stirring energy, $e$ , at a rate $\mathbb{I}_{M}$ . The stirring energy then diffuses ( $\mathbb{D}_{e}$ ), dissipates ( $\unicode[STIX]{x1D716}$ ), and produces BPE via an irreversible buoyancy flux ( $\mathbb{B}$ ) that raises the centre of mass of the density field $\unicode[STIX]{x1D70C}$ . That is, $\mathbb{B}$ is positive semi-definite. Finally, the density field also diffuses ( $\mathbb{D}_{\unicode[STIX]{x1D70C}}$ ). Each of the operators, described above, will vary with the vertical coordinate $z$ and time $t$ .

Balmforth et al. (Reference Balmforth, Llewellyn Smith and Young1998) constructed a turbulence model that coupled the horizontally averaged turbulent kinetic energy and the density (buoyancy) field. That model depends critically on a mixing length scale $l$ , over which the turbulent eddies can mix the surrounding fluid. We follow an approach similar to Balmforth et al. (Reference Balmforth, Llewellyn Smith and Young1998) to model $e$ and $\unicode[STIX]{x1D70C}$ . In this formulation, we write:

(3.2a-c ) $$\begin{eqnarray}\displaystyle & \displaystyle \mathbb{D}_{e}=\unicode[STIX]{x2202}_{z}[(\unicode[STIX]{x1D708}_{e}+\unicode[STIX]{x1D708})\unicode[STIX]{x2202}_{z}e],\quad \unicode[STIX]{x1D716}=\unicode[STIX]{x1D6FD}\frac{e^{3/2}}{l},\quad \mathbb{B}=-\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D708}_{e}\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D70C}, & \displaystyle\end{eqnarray}$$

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle \mathbb{D}_{\unicode[STIX]{x1D70C}}=\unicode[STIX]{x1D705}\unicode[STIX]{x2202}_{z}^{2}\unicode[STIX]{x1D70C}. & \displaystyle\end{eqnarray}$$

Here, we have augmented the previous model with an explicit kinematic viscosity ( $\unicode[STIX]{x1D708}$ ) and molecular diffusivity ( $\unicode[STIX]{x1D705}$ ). Both $e$ and $\unicode[STIX]{x1D70C}$ are primarily driven by eddy diffusion, defined in terms of a turbulent viscosity that, on dimensional grounds, is given as $\unicode[STIX]{x1D708}_{e}=l\sqrt{e}$ . The turbulent dissipation ( $\unicode[STIX]{x1D716}$ ) and buoyancy flux ( $\mathbb{B}$ ) are similarly constructed. The parameters $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$ are model constants and will be discussed below. Finally, the non-dimensional turbulent length scale ( $l$ ) will depend on the local density gradient. For a nearly uniform density field, this length scale will be set by the vortex-ring diameter ( $l=a$ ). However, where there is a strong density gradient, the vertical length scales are constrained. The experimental work of Park, Whitehead & Gnanadeskian (Reference Park, Whitehead and Gnanadeskian1994) suggested that, in a strongly stratified environment, the turbulent length scale will be proportional to $e/\sqrt{g|\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D70C}|}$ . As such, Balmforth et al. (Reference Balmforth, Llewellyn Smith and Young1998) proposed a simple model for the length scale that preserves these limits,

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle l=\frac{a\sqrt{e}}{\sqrt{e-\unicode[STIX]{x1D6FE}g\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D70C}}}, & \displaystyle\end{eqnarray}$$

with free parameter $\unicode[STIX]{x1D6FE}$ .

To model the vortex-ring system, we need to augment this model with the input of energy from the vortex ring, $T$ . We define the advection $\mathbb{A}$ and breakdown $\mathbb{I}_{M}$ terms as:

(3.5a,b ) $$\begin{eqnarray}\displaystyle & \displaystyle \mathbb{A}=U\unicode[STIX]{x2202}_{z}T,\quad \mathbb{I}_{M}=\unicode[STIX]{x1D706}g\left(\frac{\unicode[STIX]{x1D70C}-\unicode[STIX]{x1D70C}_{1}}{\unicode[STIX]{x1D70C}_{0}}\right)\sqrt{T}. & \displaystyle\end{eqnarray}$$

The density $\unicode[STIX]{x1D70C}_{1}=\unicode[STIX]{x1D70C}(z=z_{0})$ is the density at the vortex-ring initialization height $z_{0}$ . The advection term ( $\mathbb{A}$ ) prescribes that $T$ is transported vertically downward at the constant propagation speed $U$ . Based upon the work of Olsthoorn & Dalziel (Reference Olsthoorn and Dalziel2017), we know that the stratified vortex-ring system is unstable, with a growth rate proportional to the bulk Richardson number of the flow. The parameterization of $\mathbb{I}_{M}$ , which is constructed on dimensional grounds, captures this dependence (see below) with constant $\unicode[STIX]{x1D706}$ , a free parameter that we will set to unity.

Finally, we must prescribe the generation rate $\mathbb{S}$ of the vortex rings. In this model, $T$ is forced periodically and instantaneously. That is, after each time interval $\unicode[STIX]{x0394}t=\unicode[STIX]{x1D70F}_{R}$ , a vortex ring is instantaneously introduced into the system. Mathematically, this is written as

(3.6) $$\begin{eqnarray}\displaystyle & \displaystyle \mathbb{S}=\mathop{\sum }_{n=0}^{N}\frac{\text{KE}_{Ring}}{Aa}f(z-z_{0})\unicode[STIX]{x1D6FF}(t-n\unicode[STIX]{x1D70F}_{R}). & \displaystyle\end{eqnarray}$$

Here, $A$ is the plan area of the stratified tank and $\unicode[STIX]{x1D6FF}$ is a Dirac delta function. The system is periodically forced for a specified number of iterations $N$ . The index $n$ is the vortex-ring generation number, which identifies the number of vortex rings that have been input into the system. The functional form of $f(z)$ is defined below.

We non-dimensionalize the physical parameters as

(3.7a,b ) $$\begin{eqnarray}\displaystyle & \displaystyle z^{\prime }=\frac{z}{a},\quad t^{\prime }=\frac{U}{a}t, & \displaystyle\end{eqnarray}$$

(3.8a-c ) $$\begin{eqnarray}\displaystyle & \displaystyle T^{\prime }=\frac{T}{U^{2}},\quad e^{\prime }=\frac{e}{U^{2}},\quad \unicode[STIX]{x1D70C}^{\prime }=\frac{\unicode[STIX]{x1D70C}-\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{0}}, & \displaystyle\end{eqnarray}$$

(3.9a-c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\prime }=\frac{\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x1D70C}_{0}},\quad \unicode[STIX]{x1D70F}=\frac{U\unicode[STIX]{x1D70F}_{R}}{a},\quad K^{\prime }=\frac{\text{KE}_{Ring}}{\unicode[STIX]{x1D70C}_{0}U^{2}Aa}, & \displaystyle\end{eqnarray}$$

where the reference density $\unicode[STIX]{x1D70C}_{0}$ is selected to be the initial minimum density of the system. Similarly, $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\prime }$ is defined as the difference between the initial maximum and minimum density of the system, from which we specify an initial Richardson number $Ri_{0}=g(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{0})/\unicode[STIX]{x1D70C}_{0}(a/U^{2})$ . Finally, $K^{\prime }$ is the non-dimensionalized kinetic energy of the vortex ring. The model then reduces to the following, dropping the primes for convenience:

(3.10) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}T=\unicode[STIX]{x2202}_{z}T-\unicode[STIX]{x1D706}\mathit{Ri}_{0}(\unicode[STIX]{x1D70C}-\unicode[STIX]{x1D70C}_{1})\sqrt{T}+\mathop{\sum }_{n=0}^{N}Kf(z-z_{0})\unicode[STIX]{x1D6FF}\left(t-n\unicode[STIX]{x1D70F}\right), & \displaystyle\end{eqnarray}$$

(3.11) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}e=\unicode[STIX]{x2202}_{z}\left[\left(\unicode[STIX]{x1D708}_{e}+\frac{1}{\mathit{Re}}\right)\unicode[STIX]{x2202}_{z}e\right]-\unicode[STIX]{x1D6FD}\frac{e^{3/2}}{l}+\unicode[STIX]{x1D6FC}\mathit{Ri}_{0}\unicode[STIX]{x1D708}_{e}\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D70C}+\unicode[STIX]{x1D706}\mathit{Ri}_{0}(\unicode[STIX]{x1D70C}-\unicode[STIX]{x1D70C}_{1})\sqrt{T}, & \displaystyle\end{eqnarray}$$

(3.12) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D70C}=\unicode[STIX]{x2202}_{z}\left[\left(\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D708}_{e}+\frac{1}{Re~Sc}\right)\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D70C}\right]. & \displaystyle\end{eqnarray}$$

The functional form of $f$ is then given:

(3.13) $$\begin{eqnarray}\displaystyle & \displaystyle f(z)=\frac{1}{\sqrt{2\unicode[STIX]{x03C0}\unicode[STIX]{x1D70E}^{2}}}\exp \left[-\frac{z^{2}}{2\unicode[STIX]{x1D70E}^{2}}\right]. & \displaystyle\end{eqnarray}$$

Here, $\unicode[STIX]{x1D70E}=1/4$ such that the width of the forcing equals the vortex-ring size.

This model has four free parameters. The work of Tominaga & Stathopoulos (Reference Tominaga and Stathopoulos2007) has shown that the turbulent Schmidt (Prandtl) number $\unicode[STIX]{x1D6FC}$ has a typical value of 0.2–1.3, depending on the flow structure. For the purposes of this model, we set $\unicode[STIX]{x1D6FC}=1$ . As reported in Vassilicos (Reference Vassilicos2015) for decaying turbulence, the dissipation parameter $\unicode[STIX]{x1D6FD}$ , where it is constant, has a value near one, and thus we set $\unicode[STIX]{x1D6FD}=1$ . With reference to Park et al. (Reference Park, Whitehead and Gnanadeskian1994), the parameter $\unicode[STIX]{x1D6FE}$ is of order one, and thus we set this parameter to one. Based upon the work of Olsthoorn & Dalziel (Reference Olsthoorn and Dalziel2017), we suggest that the value of $\unicode[STIX]{x1D706}$ is also $O(1)$ . The vortex breakdown parameter $\unicode[STIX]{x1D706}$ is therefore also set to one. Thus, in this paper we restrict ourselves to the case where the free parameters are all set to unity. We return to this later.

The model was implemented on a uniform grid, using pseudospectral spatial derivatives and a first-order semi-implicit time stepping scheme. The computational domain was defined with 1024 grid points. Varying the number of grid points demonstrated that this resolution was sufficient for the parameter sets presented here. The code was shown to preserve mass to near machine precision. Adaptive time stepping was used to control the total energy conservation, which had a relative energy loss typically within $O(10^{-4})$ . A spectral filter was also used to limit the aliasing of the Fourier modes.

We ran a set of model experiments (runs) in a manner similar to that described for the physical experiments in § 2. Four parameter cases were performed, which prescribe the functional form of the stratification and the vortex-ring parameters. We label these model cases M1–M4. For each of these cases, four different stratification strengths ( $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=\{0.01,0.02,0.04,0.08\}$ ) were set, resulting in a total of 16 runs. As described in § 2, each model run will comprise of sequentially generated vortex rings enumerated as $n=\{0,1,2,\ldots ,N\}$ , for a total of $N=500$ generated vortex rings in each of the 16 runs, as prescribed by (3.13). The parameters associated with each of these runs are presented in table 1. In the first three model cases (M1–M3), similar to the experiments of Olsthoorn & Dalziel (Reference Olsthoorn and Dalziel2015), a continuous two-layer density profile was specified using a $\tanh$ function, with an initial interface height of $H_{0}$ and an interface thickness of $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D70C}}$ . These density profiles are prescribed as

(3.14) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}_{M1-M3}(z,t=0)=\frac{1}{2}\left(1-\tanh \left[\frac{z-H_{0}}{\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D70C}}}\right]\right). & \displaystyle\end{eqnarray}$$

The dimensionalized initial conditions are given in table 2 and were selected to approximate the physical experiments performed with salt water in Olsthoorn & Dalziel (Reference Olsthoorn and Dalziel2015). The fourth case (M4) was initialized with a linear stratification, given

(3.15) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}_{M4}(z,t=0)=(z-L). & \displaystyle\end{eqnarray}$$

Here, $L$ is the height of the domain.

Table 1. Table of the relevant characteristic parameters of the different vortex-ring cases. Data from the physical experiments were taken from Olsthoorn & Dalziel (Reference Olsthoorn and Dalziel2015). Note that the diameter of the Hill’s vortex ring (shown in parentheses) is defined as $2R$ .

Table 2. Table of the dimensional domain parameters for the model and numerical simulation.

Figure 4. Plot of the evolution of the density profiles for one run ( $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=0.01$ ) from each of the model cases. Here, model experiments (a) M1, (b) M2: high $Ri$ , (c) M3: low $Sc$ , and (d) M4: linear stratification, are all presented. Density profiles are plotted just prior to the generation of vortex ring $n=\{0,125,250,375\}$ .

Figure 4 shows the evolution of the density profiles for one run ( $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=0.01$ ) from each of the four model cases. Density profiles were plotted just prior to the generation of vortex ring $n=\{0,125,250,375\}$ . The results show excellent qualitative agreement with the physical experiments. We observe that, as in Olsthoorn & Dalziel (Reference Olsthoorn and Dalziel2015), the evolution of the density profiles is defined by three generic characteristics. First, the vortex rings sharpen the density interface. Second, the vortex rings generate a middle fluid layer that is near homogeneously mixed. Third, the growth of the middle fluid layer is limited by the vortex-ring injection height.

Comparison of figure 4(a) (Model Exp. M1) and figure 4(b)(M2) demonstrates how the density field evolution changes for different vortex-ring parameters ( $\text{M1}:\,\mathit{Re}=1600,\mathit{Ri}=2.5\text{ versus }\text{M2}:\,\mathit{Re}=400,\mathit{Ri}=4.9$ ). The same characteristic evolution of the density profiles is observed for M2, though only a small amount of scouring of the density interface has occurred due to the decrease in kinetic energy input. Figure 4(c) (M3) varies the molecular diffusivity ( $\unicode[STIX]{x1D705}$ ) of the stratification ( $\mathit{Sc}=3$ versus $\mathit{Sc}=1000$ ). Again, the same features of the density profiles are observable, except that the vortex rings are no longer able to effectively sharpen the lower interface as it diffuses. Due to the finite domain size, significant diffusion of the density interface can limit the run time of the model. This will be important when discussing the numerical simulations below. Finally, figure 4(d) (M4) encapsulates the effect of a different initial background stratification (linear profile), demonstrating the same characteristic evolution, although we have no matching physical experiments against which to compare these runs.

Figure 5 shows the partition of energy for the model experiments for a single mixing event ( $n=250$ ) of M1 ( $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=0.01$ ). This shows the integrated vortex-ring energy (VRE $=\int _{0}^{L}T\,\text{d}z^{\prime }$ ) and integrated stirring energy ( $\text{SE}=\int _{0}^{L}e\,\text{d}z^{\prime }$ ) along with the change in background potential energy ( $\unicode[STIX]{x0394}$ BPE) of the system (correcting for the background diffusion $\mathbb{D}_{\unicode[STIX]{x1D70C}}$ ). The integrated dissipation (D $=\int _{t_{n}}^{t}\int _{0}^{L}\unicode[STIX]{x1D716}\,\text{d}z^{\prime }\,\text{d}t^{\prime }$ ) has also been plotted. We observe that the mixing is temporally confined near the peak in $e$ , as emphasized in § 2. However, data on the time-dependent dissipation and mixing rates are not available for the physical experiments and thus comparison is limited to that of the density profiles $\unicode[STIX]{x1D70C}(z,t)$ .

Figure 5. Plot of the partition of energy between the integrated vortex-ring energy (VRE) and stirring energy (SE) for vortex ring $n=250$ . The change in potential energy ( $\unicode[STIX]{x0394}$ BPE) and the integrated dissipation (D) are also plotted. Time has been normalized by the inter vortex-ring spacing $\unicode[STIX]{x1D70F}$ . Data have been plotted for M1 with $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=0.01$ .

4 Simulation of the vortex-ring-induced mixing experiments

We validate the 1-D model results using a 3-D pseudospectral numerical solver (SPINS; see Subich, Lamb & Stastna Reference Subich, Lamb and Stastna2013) to solve the incompressible Navier–Stokes equations under the Boussinesq approximation. These equations can be written as

(4.1) $$\begin{eqnarray}\displaystyle & \displaystyle \left(\unicode[STIX]{x2202}_{t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\right)\boldsymbol{u}=-\unicode[STIX]{x1D735}p-\mathit{Ri}\unicode[STIX]{x1D70C}g\hat{\boldsymbol{z}}+\frac{1}{\mathit{Re}}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}, & \displaystyle\end{eqnarray}$$

(4.2) $$\begin{eqnarray}\displaystyle & \displaystyle \left(\unicode[STIX]{x2202}_{t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\right)\unicode[STIX]{x1D70C}=\frac{1}{Re~Sc}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D70C}, & \displaystyle\end{eqnarray}$$

(4.3) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0. & \displaystyle\end{eqnarray}$$

Here, $\boldsymbol{u}$ and $p$ are the velocity and pressure fields, respectively. Boldface variables denote vector quantities.

Experimental visualizations of the interaction of a vortex ring with a stratified interface (see Olsthoorn & Dalziel Reference Olsthoorn and Dalziel2017) demonstrate that, where the vortex ring propagates parallel to the direction of gravity, the flow field remains predominantly axisymmetric about the vortex-ring axis throughout the majority of the interaction, despite the formation of a three-dimensional instability. In order to facilitate the numerical computations, we take partial advantage of this symmetry by simulating a quarter ring in a triply periodic, free-slip (cosine transform) domain. The parameters associated with the numerical simulations can be found in table 1. Grid resolution studies determined that the resolution (see table 1) was sufficient to estimate the mixing efficiency, although we note that we do not resolve down to the Batchelor scale of the flow. As with the model results, a high molecular diffusivity results in a thick density interface, which will eventually violate the continuous two-layer set-up considered here, and will limit the run time of each experiment. Thus, we desire the lowest diffusivity that is computationally viable. In these simulations, we select $Sc=3$ . Four sets of numerical simulations were performed. The initial density stratification for each case was defined via a tanh profile similar to (3.14). See table 2 for the initial conditions.

Three sets of simulations (Num. Exp. N1–N3) were initialized with a Hill’s spherical vortex as it is a classical vortex-ring solution. The Hill’s vortex can be written as

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle \left.\begin{array}{@{}c@{}}\displaystyle u_{r}=\left\{\begin{array}{@{}ll@{}}\displaystyle \frac{3}{2}U\frac{zr}{R^{2}}\quad & r\leqslant R,\\ \displaystyle \frac{3}{2}U\frac{zr}{R^{2}}\left(\frac{R^{2}}{z^{2}+r^{2}}\right)^{5/2}\quad & r>R,\end{array}\right.\\ \displaystyle u_{z}=\left\{\begin{array}{@{}ll@{}}\displaystyle \frac{3}{2}U\left(\frac{5}{3}-\frac{2r^{2}+z^{2}}{R^{2}}\right)\quad & r\leqslant R\\ \displaystyle U\left[\left(\frac{R^{2}}{r^{2}+z^{2}}\right)^{5/2}\left(\frac{2z^{2}-r^{2}}{2R^{2}}\right)-1\right]\quad & r>R.\end{array}\right.\end{array}\right\} & \displaystyle\end{eqnarray}$$

In this paper, $R$ is the radius of the Hill’s vortex, and $U<0$ is its propagation speed. We note that there is a mismatch between the definition of the Hill’s vortex diameter ( $2R$ ) and the experimentally measured vortex-ring diameter that was defined as the distance between vorticity centroids (see Olsthoorn & Dalziel Reference Olsthoorn and Dalziel2015). A random initial velocity perturbation of $O(10^{-4})$ , relative to the vortex-ring propagation speed, was added to the numerical simulations in order to trigger any instabilities in the interaction between the vortex ring and the stratification.

For the fourth set of simulations (N4), a different initial condition was used to assess the dependence of the results on the ring aspect ratio (core size/ring diameter). Similar to Archer, Thomas & Coleman (Reference Archer, Thomas and Coleman2009), a vorticity distribution was initialized into the numerical solver (we used an azimuthally rotated shielded dipole) that, when time evolved, produced a coherent vortex ring. This resultant non-spherical vortex ring was then used as the initial condition for the numerical mixing simulations. Fitting the vortex core to a Gaussian distribution, we estimate the aspect ratio of this vortex ring to be 0.17. In the physical experiments, the vortex rings had an aspect ratio of ${\approx}0.1$ .

As with the experimental setup described in § 2, and the model set-up described in § 3, the simulations were run by generating vortex rings that interact with the stratification. The flow was then evolved until the velocity field dissipated sufficiently. After a delay ( $\unicode[STIX]{x1D70F}$ ), the velocity field for a new vortex ring was superimposed (by addition) onto the residual velocity field. This cycle was repeated until the desired number of iterations was achieved. We set the end time to be 100 vortex rings. Thirteen different parameter cases (requiring 1300 simulations) were performed for various Reynolds and Richardson numbers, the details of which can be found in table 1. As with the experimental results of Olsthoorn & Dalziel (Reference Olsthoorn and Dalziel2015) and the model results above, there is an initial set-up period, within which the functional form of the stratification varies. After this setup period, the stratification tends to a self-similar form and the mixing rate is nearly constant, and it is this value that is reported.

Figure 6 shows the sorted density profile every 20 vortex-ring generations for one numerical simulation (N3: $Re=3500$ , $Ri=2.75$ ). Here, we again observe the same characteristic features of the background density field evolution. Note the similarity between figure 4(c) and figure 6. As before, the diffusion of the background stratification is significant and must be accounted for when considering the mixing rate of each vortex ring.

Figure 6. Plot of the sorted density profile every twenty vortex-ring generations of one numerical simulation (N3: $Re=3500$ , $Ri=2.75$ ). A solid line is drawn at $z=z_{0}$ .

Additionally, figure 7(a) shows the evolution of the distribution of energy into its various compartments for the first vortex ring of one numerical simulation (N3: $Re=3500$ , $Ri=2.75$ ). This figure is reminiscent of figure 5 from the model results. In particular, the time dependence of the mixing (M) and the total dissipation (D) provided in figure 7(a) are similar to those found previously in figure 5, though their relative values are different. Unlike the model, the numerical simulations explicitly support the generation of APE. Associated with this APE is the generation of internal waves that manifest as oscillations between the APE and the kinetic energy (KE). Figure 7(b) shows the relative energy distribution prior to the generation of a new vortex ring, for all vortex-ring generations. This plot quantifies the incremental change to the mixing rate of each subsequent vortex ring. Both panels (a–b) have been normalized by the initial vortex-ring energy (E) and the interval $\unicode[STIX]{x1D70F}$ between vortex rings for comparison with the model results. Typical net relative energy loss at time $t=\unicode[STIX]{x1D70F}$ is $5\times 10^{-3}$ . Due to the late-time exponential decay of the internal waves generated, there will always be some residual stirring energy ( $\text{RE}=\text{KE}+\text{APE}$ at $t=\unicode[STIX]{x1D70F}$ ) in the system prior to subsequent vortex-ring generations. As there is a practical limitation on the length of each numerical simulation, we terminated the simulations when $\text{RE}/\text{KE}_{ring}=O(10\,\%)$ . The physical experiments also have RE, though it is much less than that of the numerical simulations, as we can wait longer between vortex-ring generations at almost no cost. As this residual energy remains nearly constant with subsequent vortex-ring generations, the RE will have a small, near constant, contribution to the increase in BPE of the system, when compared to the mean mixing rate of each vortex ring. As mentioned above, these simulations demonstrate an initialization period, after which the change in potential energy of the system is near constant.

Figure 7. (a) Plot of the energy partition for the first vortex-ring interaction of Num. Exp. N3 (N3: $Re=3500$ , $Ri=2.75$ ). Note that the KE is initially slightly above 1 as a result of the initial perturbation. (b) Plot of the energy partition at the end of each subsequent $n$ vortex-ring interactions.

5 Discussion and conclusion

For each experiment, we compute the ratio of the change in background potential energy ( $\unicode[STIX]{x0394}$ BPE) between successive vortex rings ( $\unicode[STIX]{x0394}t=\unicode[STIX]{x1D70F}$ ) versus the energy of the input vortex ring ( $\text{KE}_{ring}$ ). We define this ratio as the mixing efficiency ( $\unicode[STIX]{x1D702}$ ), indicating the amount of background potential energy change for a given energy input. This definition of the mixing efficiency is unambiguous where $\text{RE}=0$ . Where $\text{RE}\neq 0$ , provided that the RE is constant between vortex-ring generations, the associated mixing will also be constant and the interpretation of the mixing efficiency remains well defined over the interval between vortex rings. The mixing efficiency is then computed as

(5.1) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}=\frac{\unicode[STIX]{x0394}\text{BPE}-\unicode[STIX]{x0394}\text{PE}_{\unicode[STIX]{x1D705}}}{\text{KE}_{ring}},\quad \text{where }\unicode[STIX]{x0394}\text{BPE}=gA\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}a^{2}\int \left[\unicode[STIX]{x1D70C}_{s}^{(n+1)}-\unicode[STIX]{x1D70C}_{s}^{(n)}\right]z\,\text{d}z. & \displaystyle\end{eqnarray}$$

Here, $\unicode[STIX]{x1D70C}_{s}^{(n)}$ is the sorted density profile after $n$ vortex rings have been produced. The change in BPE is corrected for the diffusive increase in potential energy ( $\unicode[STIX]{x0394}\text{PE}_{\unicode[STIX]{x1D705}}=gA\unicode[STIX]{x1D705}\unicode[STIX]{x1D70F}(\unicode[STIX]{x1D70C}(0)-\unicode[STIX]{x1D70C}(L))$ ) as we are interested only in the contribution due to the vortex ring.

Figure 8. Mixing efficiency as a function of Richardson number for the model, numerical simulations and physical experiments found in table 1. The black solid line corresponds to the mean mixing efficiency of the physical experiments ( $\unicode[STIX]{x1D702}_{0}=0.42$ ). The dashed red line is the estimated mixing efficiency of the numerical simulations ( $\unicode[STIX]{x1D702}_{N}=0.45$ ).

Comparing the model, numerical simulations and physical experiments, figure 8 shows the mixing efficiency determined for all cases identified in table 1. In this plot, the mixing efficiency is near constant with $Ri$ . Error bars are computed as the root mean squared error from the associated mean mixing efficiency value, once the system has completed its initial set-up period. We observe that the mixing efficiency of the numerical simulations ( $\unicode[STIX]{x1D702}_{N}\approx 0.45$ ) is slightly higher than the physical, salt-water experiments ( $\unicode[STIX]{x1D702}_{0}\approx 0.42$ ), as would be expected from the lower $Sc$ . (The computational resources necessary to use the experimental value of $Sc$ were not available.) The mixing efficiencies found in the model are consistently higher still ( $\unicode[STIX]{x1D702}_{M}\approx 0.49$ ) than the numerical simulations or physical experiments, though it is still within the experimental uncertainty of the physical experiments.

As we have noted previously, the model is dependent on four free parameters ( $\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D706}$ ). Figure 9 presents the linear parameter analysis for the four free model parameters. In each of the four plots, one of the four parameters was varied while the other three were held constant. In each case the mixing efficiency was computed for one run of case M1 ( $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=0.01$ ) at $n=100$ . We note that an increase in the dissipation parameter $\unicode[STIX]{x1D6FD}$ by ${\approx}15\,\%$ (see figure 9 b) would account for the difference between the experimental value of the mixing efficiency and the model runs. It is worth noting that that the parameters associated with the dissipation rate are where most of the sensitivity of the model resides. A more precise parameter selection is left for future work.

Figure 9. The dependence of the mixing efficiency on the model parameters (a) $\unicode[STIX]{x1D6FC}$ , (b) $\unicode[STIX]{x1D6FD}$ , (c) $\unicode[STIX]{x1D6FE}$ and (d) $\unicode[STIX]{x1D706}$ . A dashed line has been plotted at $\unicode[STIX]{x1D702}=0.5$ . Mixing efficiency values were generated from case M1 ( $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=0.01$ ) at $n=100$ .

As the Richardson number decreases below $O(1)$ , the mixing efficiency dependence on $Ri$ becomes more ambiguous. Indeed, recent work by Shrinivas & Hunt (Reference Shrinivas and Hunt2015) has indicated that the vertical confinement of turbulent mixing may change the mixing efficiency dependence on $Ri$ . This confinement will be especially pronounced at low $Ri$ due to the deep penetration of the vortex rings into the lower layer. The confinement is entirely omitted in the model due to its 1-D construction. As the effect of confinement will influence the three-dimensional structure of the flow, one might initially model it by modifying the propagation speed ( $U$ ) of the vortex rings near the boundaries. We do not attempt this here.

Finally, we want to highlight that the mixing efficiency is defined here as an aggregate, time-independent quantity. That is, the net fluid mixing that results from a given energy input. This definition suggests that vortex rings are effective mixers as they transport energy directly to the density interface (with minimal dissipation), produce vorticity directly at the location of the peak in the density gradient (which Lawrie & Dalziel (Reference Lawrie and Dalziel2011) argued would result in a high mixing efficiency), and, through a flow instability (see Olsthoorn & Dalziel Reference Olsthoorn and Dalziel2017), generate turbulence. This series of events enables each vortex ring to create a near optimal mixing state such that nearly all the vortex-ring energy is deposited directly at the location of peak mixing. That is, the kinetic energy of the vortex ring produces stirring energy at the location of peak density gradient. In addition, turbulent stirring energy is generated at a rate proportional to the bulk Richardson number of the system, which is essential for the system to establish a self-similar density profile. The model, which is a simplification of the vortex-ring system, emphasizes this picture by only generating stirring energy (at a rate proportional to the bulk Richardson number) where the density field is not constant; where mixing can occur. This is in contrast to grid-generated turbulence, which dissipates significantly before reaching the density interface.

This paper presents a model for isolated vortex-ring-induced stratified mixing experiments. This work has been shown to provide qualitative and quantitative agreement with both physical experiments and numerical simulations. At moderate $Ri$ , the mixing efficiency of the vortex rings has been shown, in all three methodologies, to be near constant after an initialization period with very similar asymptotic values. This constant mixing efficiency regime of the vortex-ring experiments has been previously reported in Olsthoorn & Dalziel (Reference Olsthoorn and Dalziel2015), although the numerical work found in this paper demonstrates this regime at a much lower Schmidt number ( $\mathit{Sc}=3\text{ versus }\mathit{Sc}=700$ ). The 1-D model constructed in the paper encapsulates the essential features of the energy pathways for the vortex-ring-induced mixing experiments. In particular, this work highlights the important contribution of the vortex-ring breakdown being proportional to the Richardson number ( $\mathbb{I}_{M}\propto Ri_{0}$ ). As demonstrated in Olsthoorn & Dalziel (Reference Olsthoorn and Dalziel2017), the dominant vortex-ring instability in the strongly stratified system ( $Ri>O(1)$ ) has a growth rate proportional to the bulk Richardson number of the flow. This model demonstrates that the identified vortex-ring instability plays a key role in establishing the constant mixing efficiency regime.

We study vortex-ring-induced mixing in analogy to large-scale turbulent-eddy mixing events. However, it should be clear that stratified turbulence is characterized by its large range of length scales and complex flow structures. As such, a natural extension of the present model would investigate a convolution of the individual mixing events discussed here. Future work will investigate the application of this model to a mixing box experiment, similar to the one described in Turner (Reference Turner1968).