2.1. Governing equations

We consider an incompressible fluid flow driven by small density differences, such that the Boussinesq approximation can be employed. The fluid motion is described by the three-dimensional Navier–Stokes equations for an incompressible flow in the Boussinesq limit, and by an advection–diffusion equation for the transport of the buoyancy scalar. The governing equations are given by

(2.1)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}=0, \end{gather}$$

(2.2)$$\begin{gather}\frac{\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} )\boldsymbol{u}={-}\frac{1}{\rho_0}\boldsymbol{\nabla} p + \nu{\rm \Delta} \boldsymbol{u} -b\boldsymbol{e}_z, \end{gather}$$

(2.3)$$\begin{gather}\frac{\partial b}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} b = \kappa\nabla^{2}b + g\mathcal{S}, \end{gather}$$

where $\boldsymbol {u}$ is the fluid velocity decomposed in Cartesian coordinates as $\boldsymbol {u} = u\boldsymbol {e}_x + v\boldsymbol {e}_y +w\boldsymbol {e}_z$, $p$ is the pressure, $\rho _0$ is the reference fluid density, $\nu$ is the kinematic viscosity, $b=g(({\rho -\rho _0})/{\rho _0})$ is relative buoyancy, $\kappa$ is the scalar diffusivity and $\mathcal {S}$ is a relative density source term. In geophysical applications of gravity currents, the buoyancy scalar is often the result of salinity, such that $\kappa$ is the diffusivity of salt ions in water. This modelling approach can also be used to investigate particle-driven currents as long as the particles can be considered non-inertial, such that their transport can be modelled in the equilibrium-Eulerian framework as a scalar quantity. Then, the advection velocity of the particles is equal to the sum of the fluid velocity $\boldsymbol{u}$ and the constant Stokes particle settling velocity $\boldsymbol{V}$. This approach has previously been successfully employed to investigate numerically the effects of settling on the dynamics of particle-laden flows (Necker et al. Reference Necker, Härtel, Kleiser and Meiburg2002), as well as their dissipation and mixing properties (Necker et al. Reference Necker, Härtel, Kleiser and Meiburg2005). In the present study, we consider a generic buoyancy term $b$ which can represent continuous quantities such as salinity, or non-settling particles that satisfy the equilibrium-Eulerian assumptions. While the role of settling is not investigated here, it can be simply included in this modelling approach, provided that the equilibrium-Eulerian assumption can be made. The validity of this assumption, which depends on the Stokes number of the particles, might not be valid in the immediate wake of a deep-sea mining collector where considerable turbulence is expected, but has been successfully employed in the study of turbidity currents (Necker et al. Reference Necker, Härtel, Kleiser and Meiburg2005; Nasr-Azadani, Meiburg & Kneller Reference Nasr-Azadani, Meiburg and Kneller2018; Ouillon, Meiburg & Sutherland Reference Ouillon, Meiburg and Sutherland2019).

2.2. Numerical set-up and non-dimensional equations

We consider a spherical source of mass, of diameter $D$, moving with velocity $U_s$. Preliminary simulations (not presented) revealed that, in the reference frame of the moving source, there exists a statistically steady-state regime for sufficiently large source velocities. In order to study this regime for longer times without requiring a very large numerical domain, we conduct the numerical simulations in the reference frame of the source. The centre of the source is located at a stationary position $z=D/2$ in the reference frame introduced in figure 1. An inflow and outflow condition are imposed at the left and right boundaries, respectively, with the inflow condition being set to $u=U_s$. The left inflow boundary is at a distance $x_s$ from the source. At the top and bottom boundary, a wall velocity $U_w=U_s$ is imposed. Periodic boundary conditions are imposed in the $y$-direction. Provided that the flow perturbations introduced by the source are sufficiently far from the lateral boundaries, this set-up is equivalent to that of a moving source in a closed domain with bottom and top no-slip boundary conditions. The source term in the advection–diffusion equation (2.3) is assumed to be $\mathcal {S}=S\xi (\boldsymbol {x})$ where $S$ is the intensity of the source and $\xi$ is the function that defines the sphere in three-dimensional space. The intensity of the source is a function of its volume $V$ and of the influx of mass $\dot {m}$ (in $\mathrm {kg}\ \mathrm {s}^{-1}$), i.e.

(2.4)\begin{equation} S=\frac{\dot{m}}{\rho_0V}=\frac{\dot m}{\rho_0\dfrac{\rm \pi}{6}D^{3}}.\end{equation}

An error function is used to avoid large gradients at the edge of the source, so that

(2.5)\begin{equation} \xi(\boldsymbol{x},t) = \frac{1}{2}\left(1-\textrm{erf}\left[\frac{\|\boldsymbol{x}\|-\dfrac{D}{2}}{\delta}\right]\right), \end{equation}

where $\delta$ is the thickness of the error function, typically chosen as $\delta = 2h$ where $h$ is the grid size.

Figure 1.

Two-dimensional sketch of the three-dimensional simulations of a gravity current from a spherical moving source, in the reference frame of the moving source.

The characteristic velocity of dense currents propagating on a flat bottom generally scales with a buoyancy velocity based on the height and buoyancy of the current head (von Kármán Reference von Kármán1940; Benjamin Reference Benjamin1968; Härtel, Carlsson & Thunblom Reference Härtel, Carlsson and Thunblom2000a; Huppert Reference Huppert2006; Cantero et al. Reference Cantero, Lee, Balachandar and Garcia2007). In the moving source problem, the negative buoyancy accumulated by the source depends not only on advection away from the source by the gravity current itself, but also on how fast the source is moving. We define a characteristic buoyancy velocity $U_b=\sqrt {BD}$, where $D$ is the diameter of the source, and $B=gS\tau$ is the characteristic buoyancy of the gravity current, with $\tau$ the characteristic accumulation time. When the source moves slowly with respect to the buoyancy velocity, we expect the characteristic accumulation time to scale as $\tau \sim {D}/{U_b}$. On the other hand, when the source moves rapidly, we expect the characteristic accumulation time to scale as $\tau \sim {D}/{U_s}$. In order to make the equations non-dimensional, we choose the characteristic buoyancy $B={gSD}/{U_b}$, such that the buoyancy velocity is defined as

(2.6)\begin{equation} U_b=(gSD^{2})^{1/3}. \end{equation}

With $U_b$ the reference velocity, $D$ the reference length, $P=\rho _0U_b^{2}$ the reference pressure and $B={gSD}/{U_b}$ the reference buoyancy, the non-dimensional governing equations are given by

(2.7)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}=0, \end{gather}$$

(2.8)$$\begin{gather}\frac{\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} )\boldsymbol{u} ={-}\boldsymbol{\nabla} p + \frac{1}{Re_b}{\rm \Delta} \boldsymbol{u} - b\boldsymbol{e}_z, \end{gather}$$

(2.9)$$\begin{gather}\frac{\partial b}{\partial t} +\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} b = \frac{1}{Pe_b}{\rm \Delta} b +\xi(\boldsymbol{x}), \end{gather}$$

where $Re_b = {U_bD}/{\nu }$ is the buoyancy Reynolds number and $Pe_b={U_bD}/{\kappa }$ is the Péclet number. The non-dimensional source speed is given by $a={U_s}/{U_b}$. The non-dimensional source speed $a$ controls the buoyancy available for the gravity current to form. When $a$ is small, the non-dimensional buoyancy is expected to be independent of $a$ such that $b\sim 1$. When $a$ is large, and the source speed controls the accumulation time, then the non-dimensional buoyancy depends on $a$ such that $b\sim {1}/{a}$.

We conduct a series of DNS, for source velocities $a=0$, 0.126, 0.252, 0.378, 0.63, 0.945, 1.26, 1.89 and 2.52. The position of the source within the numerical domain is adjusted as the source speed increases in order to make the best use of the numerical domain and maximize the region in which the gravity current can be studied. Similarly, for the two largest source velocities, a narrower but longer domain is employed, as the region occupied by the gravity current narrows as the source speed increases. The domain height is set to $L_z=1.5$ for all simulations. The front velocity of planar lock-release gravity currents in the slumping phase converges asymptotically with increasing Reynolds numbers (Härtel et al. Reference Härtel, Meiburg and Necker2000b; Cantero et al. Reference Cantero, Lee, Balachandar and Garcia2007), and gravity currents have thus been successfully investigated at a reduced scale from the original problem, both in the laboratory and in numerical simulations. Preliminary simulations (not shown) of the moving source problem revealed that when $a\geq 0.63$, the non-dimensional buoyancy starts decreasing as $b\sim {1}/{a}$. In order to keep the effective Reynolds number constant between simulations, we therefore keep $Re_b$ constant for $a\leq 0.63$ and increase its value with $\sqrt {{1}/{a}}$ for $a>0.63$. The role of the Péclet number is not investigated here and it is set to be equal to the Reynolds number. The Reynolds number is sufficiently large that the Péclet number is expected to play a limited role on the velocity of the gravity current front (Bonometti & Balachandar Reference Bonometti and Balachandar2008). The simulation parameters are reported in table 1.

Table 1.

Numerical parameters of the simulation campaign. The Reynolds number $Re_b$ is constant for $a\leq 0.63$, then increases with $\sqrt {{1}/{a}}$.

The equations are solved by our finite difference code, originally developed under the name TURBINS (Nasr-Azadani & Meiburg Reference Nasr-Azadani and Meiburg2011). A small $O(10^{-3})$, spatially varying random perturbation is added to the $u$-component of the velocity field at $t=0$, in order to facilitate the three-dimensional evolution of the flow. The equations are integrated in time until $t=63$, which is sufficient for the gravity current to reach one of the boundaries in cases where $a<0.63$ or for a statistically steady state to be obtained and studied in cases where $a\geq 0.63$. The grid size is set to ${\rm \Delta} x = 0.01$ in order to fully resolve the flow features down to the smallest scale, resulting in a numerical domain typically comprised of ${\sim}6\times 10^{8}$ grid points. The simulations are run for ${\sim}48$ hours on 384 cores on the Knights Landing (known as KNL) nodes of the Stampede2 supercomputer. The simulation parameters are summarized in table 1.

In § 3, we investigate the qualitative properties of the gravity current that forms from the dense scalar released at the source, as a function of the ratio $a$ of source speed $U_s$ to buoyancy velocity $U_b$. In § 4, we discuss the physical mechanisms that lead, above a critical source speed, to the establishment of a statistically steady state in the reference frame of the moving source, in what we refer to as the supercritical regime. In § 5, we demonstrate that for sufficiently large values of $a$ in the supercritical regime, the front propagation normal to the direction of the source motion is equivalent to the front propagation in a classic lock-release gravity current. We further show that existing models for the propagation speed of lock-release gravity currents in the constant velocity phase, as well as predictions of the time at which the front velocity starts to decrease, can be directly applied to the moving source current in this high-$a$ regime. Proof-of-concept experiments of a dense discharge from a moving object in a tow-tank are then presented and the results are compared with the numerical model. Finally, we explore the implications of these findings to sediment plumes released by the mining of polymetallic nodules in the deep ocean, in the Appendix (A).

5.1. Available buoyancy

When $a\gg 1$, in the reference frame of the moving source, (2.3) is dominated by advection by the flow past the source in the $x$-direction. The release of the dense scalar, of which the buoyancy is expected to scale with $1/a$ in the supercritical regime, thus predominantly results in gradients in the $y$-direction, and therefore an increase in momentum in the $y$-direction. If we consider a fixed vertical plane normal to the direction of motion of the source, the passage of the source through that plane when $a\gg 1$ results in a local increase in buoyancy which then forms a gravity current approximately constrained to that plane since gradients in the spanwise $x$-direction become small. (‘Spanwise’ is used to refer to the direction normal to the predominant direction of propagation of the gravity current, consistent with the literature.) With the exception of the shape of the lock, the problem in a fixed $y$–$z$ plane is analogous to the instantaneous lock-release problem of a gravity current in a rectangular channel. In the limit where the time scale associated with the source motion is much smaller than the time scale associated with buoyancy, i.e. ${D}/{U_s}\ll {D}/{U_b}$, the relative mass per unit length released by the passage of the source is simply given by the analytical, non-dimensional expression

(5.1)\begin{equation} m_0 = \frac{\rm \pi}{6}\frac{gSD^{3}}{U_s}\cdot\frac{1}{BD^{2}} =\frac{\rm \pi}{6a},\end{equation}

where we recall that $B={gSD}/{U_b}$ is the reference buoyancy used to non-dimensionalize the equations. Under the same assumptions, the maximum buoyancy generated by the source is found along the centreline of the source trajectory, and is given by

(5.2)\begin{equation} b_0 = \frac{gSD}{U_s}\cdot\frac{1}{BD}=\frac{1}{a}.\end{equation}

The relative mass per unit length can be directly compared with the measured mass per unit length computed from the simulations as

(5.3)\begin{equation} m_x(x) =\iint_{y,z}b(x,y,z)\,\textrm{d} y\,\textrm{d} z, \end{equation}

and plotted in figure 7(a) for $t=63$. The mass per unit length is seen to reach a local maximum at some distance behind the source for all values of $a$ in the subcritical regime. This supports the previous observation that in the subcritical regime, the dense current spreads in all directions and thus the mass per unit length in the direction of motion of the source varies with distance from the source. In the supercritical regime, however, the buoyancy increases at the source and for a short distance behind the source before reaching a constant value independent of $x$. This conservation of mass along the $x$-axis in the supercritical regime, on the other hand, demonstrates that the negative buoyancy released by the source spreads over a constant area.

Figure 7.

(a) Mass per unit length $m_x(x) = \iint _{y,z}b(x,y,z)\,\textrm {d} y\,\textrm {d} z$ function of $x\in [-1,4]$ for different values of $a$ at $t=63$. (b) Maximum mass per unit length for different values of $a$ at $t=63$, plotted against the analytical value in the absence of buoyancy advection or diffusion (see (5.1)).

In figure 7(b), the global maxima of $m_x$ and the analytical solution $m_0$ (see (5.1)) in the limit of no diffusion or advection are plotted as a function of $a$, and show excellent agreement for all values of $a$ in the supercritical regime, except for $a =0.63$, where the analytical formula underestimates the buoyancy, suggesting that other accumulation processes allow the buoyancy to build up at this source velocity.

The excellent agreement between the theoretical mass per unit length in the absence of buoyancy-driven advection or diffusion, expressed in (5.1), and the measured mass per unit length in the simulations when $a\geq 0.945$ shows that for such high values of $a$, the source acts as an instantaneous release of buoyancy in the $y$–$z$ plane. In § 5.2, we show that this release in buoyancy at high values of $a$ drives a flow analogous to small-release gravity currents in the rectilinear configuration.

5.2. Normal front velocity

As shown in § 5.1, the buoyancy generated by the source remains constant within the vertical $y$–$z$ plane, which suggests that in the regime where $a>1$, the source does not act as a constant-flux release, but as a constant-volume release constrained to the vertical plane. In the following, we therefore compare the propagation of a rectilinear lock-release (RLR) gravity current with the propagation of the MS gravity currents in the $y$-direction. For simplicity, we consider an RLR gravity current with a rectangular lock of height $h=1$ and buoyancy $b=1$. The length $l$ of the lock is set so that the lock has the same area as the MS gravity current in the vertical plane. By symmetry, the propagation of the RLR gravity current is only considered in one half of the vertical plane, such that $l=\frac {1}{2}{m_0}/{b_0}={{\rm \pi} }/{12}$. The initial condition of the RLR is sketched in figure 8 alongside a sketch of the vertical slice of buoyancy imparted by the moving source.

Figure 8.

Sketch of the RLR gravity current initial condition alongside a sketch of the vertical slice of buoyancy imparted by the moving source in a fixed vertical plane normal to the direction of motion of the source.

Using the same numerical framework, we simulate the RLR gravity current in a channel of height $1.5$, length $10$ and width $0.5$. The Reynolds number is set to $10\,000$ so as to reproduce the turbulent conditions of the moving source problem. The normal front velocity $V_f$ calculated from the moving source simulations can be directly compared with the results of the RLR gravity current simulation by rescaling $V_f$ with the buoyancy velocity of the equivalent RLR gravity current, i.e. $\sqrt {b_0}=\sqrt {{1}/{a}}$. Considering the moving source problem in the stationary reference frame, the front evolves as a function of time at a given position which relates directly to the distance from the source in the simulations as $t'=x/a$. Rescaled with the modified buoyancy velocity, the equivalent spreading time is $\tilde t = \sqrt {{1}/{a}}t'=\sqrt {{1}/{a^{3}}}x$.

We now compute the normal front velocity $V_f$ of the moving source gravity current using (4.3). Here $V_f$ is plotted as a function of the equivalent lock-release time $\tilde t$ in figure 9. Note that the front position is smoothed over a moving window in order to avoid discontinuities in the spatial derivatives that inevitably arise as a result of the discrete nature of the front position in the simulations. The front velocities for $a<0.63$, shown as dashed lines, are plotted for completeness but are outside of the regime of applicability of (4.3) since a steady state is not observed. The front velocity displays a very similar behaviour to small-release gravity currents in rectangular channels (Rottman & Simpson Reference Rottman and Simpson1983; Cantero et al. Reference Cantero, Lee, Balachandar and Garcia2007). For instance, at $a=1.26$ and $a=1.89$, the front velocity initially increases to a maximum value, then slightly decreases as it enters a phase of constant velocity, and finally starts decreasing after a finite time. At $a=2.52$, the numerical domain is too short for the current to enter this decreasing velocity regime and the constant velocity regime is observed until the maximum resolved time. At $a=0.945$, the front also reaches a maximum velocity followed by a phase of relatively constant velocity. This phase is, however, followed by a sharp increase in front velocity which precedes an abrupt drop. This unexpected change in front speed is most likely due to the presence of shear-induced, large-scale three-dimensional structures at the front, which do not exist in classical planar currents and are not investigated further here. At $a=0.63$, the front velocity quickly increases to a maximum, and promptly enters a phase of velocity decrease. A period of relatively constant velocity is only observed for a short time, for $0.6< t'<3.6$. For all values of $a> 0.63$, and to a lesser extent for $a=0.63$, the rescaled front velocity curves collapse onto a single curve, showing that both the temporal scaling of $\tilde t$ and the velocity scaling of $\sqrt {{1}/{a}}U_b$ adequately capture the behaviour of the normal front. In addition, the front velocity calculated from the simulation of the RLR gravity current in the rectilinear channel, with lock height $1$ and lock length ${\rm \pi} /12$, is plotted in figure 9. A time shift of 0.5 non-dimensional times is applied to the RLR in order to align the beginning of the acceleration phase with the moving source simulations. We note that the initial behaviour of the RLR gravity current might differ from that of the moving source as a result of the small lock length to lock height aspect ratio in the RLR simulation, leading to a column-collapse behaviour, which might explain the need for the time shift. The RLR simulation achieves a similar velocity maximum at the end of the acceleration phase and enters a phase of relative constant velocity, followed by a phase of velocity decrease.

Figure 9.

Normal front velocity $V_f$ as a function of the relative time $\tilde t$ for various values of $a$. The velocity is computed from the time averaged front position and is plotted as a dashed line for $a< 0.63$ since the front position has not converged. The vertical dotted lines mark the estimated start and end times of the constant velocity phase for $a\geq 0.63$. The black full line represents the front velocity calculated from the RLR DNS with lock height $1$ and lock length ${\rm \pi} /12$.

The normal front velocity of the MS gravity current is averaged over the time interval of relatively constant front velocity $\tilde t\in [\tilde t_1,\ \tilde t_2]$. The values of $\tilde t_1$ and $\tilde t_2$ are estimated as $\tilde t_1 \approx 1.64$, 1.79, 1.93, 2.11, 2.28, and $\tilde t_2 \approx 3.6$, 8.0, 9.4, 9.7, 10 for $a=0.63$, 0.945, 1.26, 1.89, 2.52, respectively. These times are reported in figure 9 as a short vertical dotted line. As for the moving source problem, the beginning and end times of the constant velocity phase are visually estimated and added to figure 9. Close agreement is observed between the RLR gravity current front velocity and the MS gravity current normal front velocity for $a>1$. The duration of the constant velocity phase also closely matches that of the simulation at the largest value of $a$ tested. The time averaged normal front velocity $\bar V_f$ of the MS gravity current is plotted as a function of $a$ in figure 10(a), and rescaled with the effective buoyancy $\sqrt {{1}/{a}}$ to be directly compared with the time averaged front velocity of the RLR gravity current. In addition, both the MS and RLR gravity current DNS results are compared with the front velocity predicted by the shallow water theory of Rottman & Simpson (Reference Rottman and Simpson1983) for a lock-release gravity current with the same initial lock height to channel height ratio as the RLR gravity current DNS. Great quantitative agreement is found between the normal front velocities of the MS gravity currents for values of $a\geq 0.945$ and the RLR gravity current DNS, as observed in figure 9. The observed front velocities of the MS and RLR gravity current DNS are slightly smaller than those predicted by the shallow water model of Rottman & Simpson (Reference Rottman and Simpson1983). This small discrepancy could be a result of the model parameter $\beta$, here set to one, or due to the fact that some of the assumptions break down for ratios of initial lock height to channel height larger than $\frac {1}{2}$. We note that the agreement between the MS and RLR DNS is good for $a\geq 0.945$, but significantly poorer at $a=0.63$, which further suggests that while the gravity current enters the supercritical regime at $a\approx 0.63$, the released gravity current only behaves as a lock-release current when $a\approx 1$ or larger.

Figure 10.

(a) Mean front velocity in the constant velocity regime as a function of $a$ for $a\geq 0.63$. The MS gravity current average front velocity is compared with the RLR simulation, and to the solution obtained from the shallow-water equation following Rottman & Simpson (Reference Rottman and Simpson1983) with $\beta =1$. (b) Front position relative to initial lock length at the transition time $\tilde t_2$ when the current front starts slowing down. The MS simulations and RLR simulation are compared with the experimental observations of Rottman & Simpson (Reference Rottman and Simpson1983).

In addition to the averaged front velocity, we compare the front position of the MS and RLR gravity current fronts at the time when the current starts slowing down from the constant velocity phase. Rottman & Simpson (Reference Rottman and Simpson1983) estimated that the front velocity of gravity currents produced by instantaneous releases in a rectangular channel started decreasing when the front position reached three to 10 lock-lengths depending on the initial lock to channel height ratio, increasing almost linearly with this ratio. For the present lock to channel height ratio of $\frac {2}{3}$, the findings of Rottman & Simpson (Reference Rottman and Simpson1983) suggest that the current front will start slowing down when it reaches six to nine times the initial lock length. The front position at time $\tilde t_2$, is shown as a function of $a$ for the MS gravity current simulations in figure 10(b). It is compared with the front position at the transition time $\tilde t_2$ computed from the RLR simulation, as well as with the transition distances observed in the experiments of Rottman & Simpson (Reference Rottman and Simpson1983), between $6x_0$ and $9x_0$ where we recall that $x_0 ={{\rm \pi} }/{12}$ is the lock length of the RLR gravity current simulation. The transition distance of the MS gravity current increases with $a$ and converges for values of $a\geq 1.26$, for which great agreement is found with the RLR gravity current simulation. For $a=0.63$, and to a lesser extent $a=0.945$, the gravity current travels shorter distances before starting to slow down. Both the RLR simulation and the MS simulations for all values of $a$ except $a=0.63$ find a longer transition distance than in the experiments of Rottman & Simpson (Reference Rottman and Simpson1983). The relative increase over the upper bound of $9x_0$ is approximately $30\,\%$ for the RLR simulation. We do not further investigate the cause of this discrepancy.

In summary, the temporal evolution of the MS gravity current front normal to the direction of motion of the source is found to agree well with that of the RLR gravity current with initial buoyancy $\sqrt {{1}/{a}}$, lock height $1$ and lock length ${{\rm \pi} }/{12}$. The lock length is chosen such that the mass of the dense fluid contained within the lock is equal to the mass per unit length produced by the spherical moving source in the limit of large values of $a$. We also find that while the gravity current enters the supercritical regime when $a= 0.63$, and its front velocity in the constant velocity phase is well predicted by the theory of lock-release gravity currents for $a= 0.945$, its behaviour is most similar to that of the equivalent lock-release problem when $a\geq 1.26$ in our simulations. Above this threshold, the normal front velocity of the MS gravity current in the constant velocity phase is, in dimensional units, approximately equal to $Fr\sqrt {{1}/{a}}U_b$ where the Froude number $Fr$ can be determined from the equivalent RLR gravity current problem, as shown here.

5.3. Velocity field and transition to lock-release in the supercritical regime

The good quantitative agreement between the front position and velocity of a lock-release gravity current, and the MS gravity current in the constant velocity phase for large values of $a$, is not observed for $a=0.63$ despite the flow having transitioned to the supercritical regime. The horizontal velocity averaged over a height of $\frac {1}{2}$ from the bottom,

(5.4)\begin{equation} \langle(u,v)\rangle_{{1}/{2}}= 2\int_0^{{1}/{2}} (u,v)\, \textrm{d} z, \end{equation}

is computed and plotted in figure 11 for $a=0.63$, 0.945 and 2.52 for the $y>0$ half of the domain, removing $a$ from the $x$-component of the velocity so as to visualize the velocity in the stationary reference frame. It is superimposed onto a pseudocolour plot of the magnitude of the horizontal velocity. The contour of buoyancy $\langle b \rangle _{{1}/{2}}=0.1$ is shown in red to identify the position of the head of the current. For all values of $a$ in the supercritical regime, the velocity magnitude is maximal in the head of the current. The direction of the flow, on the other hand, is markedly different at $a=0.63$ when compared with larger values of $a$. In the head of the current, the flow is initially predominantly aligned with the $y$ direction close to the source, but the velocity vectors quickly align with the direction normal to the front. The realignment of the velocity vectors with the normal to the front is characteristic of axisymmetric dense fronts spreading under their own buoyancy. Thus, despite the front position being statistically invariant in the reference frame of the moving source at $a=0.63$, the dense fluid is transported normally to the front in the head of the current, and not along the $y$-direction.

Figure 11.

Quiver plot of the horizontal velocity expressed in the stationary reference frame, vertically averaged over $z\in [0, \frac {1}{2}]$, $\langle (u,v)\rangle _{{1}/{2}}$, superimposed to a pseudocolour plot of the velocity magnitude $\|\langle (u,v)\rangle _{{1}/{2}}\|$. The velocity vectors are plotted at every 50 grid points for readability. The $0.1$ contour of buoyancy vertically averaged over $z\in [0, \frac {1}{2}]$, $\langle b\rangle _{{1}/{2}}$ is added in red to visualize the front of the current. Only the $y>0$ half of the domain is shown.

For $a\geq 0.945$, on the other hand, the velocity vectors in the head of the current are not aligned with the normal to the front, and are instead predominantly aligned with the $y$-axis in the entirety of the numerical domain behind the source. In a Lagrangian sense, fluid parcels that make up the head and that contain the dense scalar released by the source are thus transported within the $y$–$z$ plane throughout the propagation of the current. We finally conclude that despite the transition to the supercritical regime in which the source moves more rapidly than the dense front ahead of the source at $a=0.63$, the MS gravity current only undergoes the transition to a unidirectional, instantaneous lock-release current when $a\approx 1$.