2.1 Self-similar spreading

For an axisymmetric two-layer current, a constant influx rate corresponds to the critical case $\unicode[STIX]{x1D6FC}=1$ of Lister & Kerr (Reference Lister and Kerr1989), for which there is a similarity solution in which the characteristic thickness of the spreading upper layer is constant (rather than growing or decaying in time). For such a solution it is natural to define a dimensionless total height by $H=h/h_{\infty }$ and a dimensionless upper-layer thickness by $D=d/h_{\infty }$ . A scaling analysis of (2.6a ), governing the evolution of the total height, suggests that the solution may be written as a function of the similarity variable

(2.10) $$\begin{eqnarray}\unicode[STIX]{x1D709}=rt^{-1/2}\left(\frac{3\unicode[STIX]{x1D707}_{2}}{\unicode[STIX]{x1D70C}_{1}gh_{\infty }^{3}}\right)^{1/2}.\end{eqnarray}$$

We write $\unicode[STIX]{x1D709}_{\ast }=\unicode[STIX]{x1D709}(r_{\ast },t)$ for the non-dimensional radial extent of the current. In order to describe the evolution towards self-similarity, we define a rescaled temporal variable $\unicode[STIX]{x1D70F}=\log (t/t_{\ast })$ , where $t_{\ast }$ is an arbitrary reference time scale. The local mass-conservation equations (2.6) then lead to a nonlinear dimensionless system of partial differential equations which consists of a fourth-order system up to the nose coupled to a second-order system beyond the nose, which is given by

(2.11a ) $$\begin{eqnarray}\displaystyle & \displaystyle H_{\unicode[STIX]{x1D70F}}-\frac{1}{2}\unicode[STIX]{x1D709}H_{\unicode[STIX]{x1D709}}+\frac{1}{\unicode[STIX]{x1D709}}(\unicode[STIX]{x1D709}q_{H})_{\unicode[STIX]{x1D709}}=0, & \displaystyle\end{eqnarray}$$

(2.11b ) $$\begin{eqnarray}\displaystyle & \displaystyle D_{\unicode[STIX]{x1D70F}}-\frac{1}{2}\unicode[STIX]{x1D709}D_{\unicode[STIX]{x1D709}}+\frac{1}{\unicode[STIX]{x1D709}}(\unicode[STIX]{x1D709}q_{D})_{\unicode[STIX]{x1D709}}=0, & \displaystyle\end{eqnarray}$$

$q_{H}(\unicode[STIX]{x1D709},\unicode[STIX]{x1D70F})$

$q_{D}(\unicode[STIX]{x1D709},\unicode[STIX]{x1D70F})$

$\unicode[STIX]{x1D709}<\unicode[STIX]{x1D709}_{\ast }$

(2.12a ) $$\begin{eqnarray}\displaystyle & \displaystyle q_{H}=-\left[\left\{H^{3}+(m-1)D^{3}\right\}H_{\unicode[STIX]{x1D709}}+{\textstyle \frac{1}{6}}\unicode[STIX]{x1D700}\left\{(H-D)^{3}\right\}_{\unicode[STIX]{x1D709}}(2H+D)\right], & \displaystyle\end{eqnarray}$$

(2.12b ) $$\begin{eqnarray}\displaystyle & \displaystyle q_{D}=-{\textstyle \frac{1}{2}}\left[\left\{3H^{2}+(2m-3)D^{2}\right\}DH_{\unicode[STIX]{x1D709}}+\unicode[STIX]{x1D700}\left\{(H-D)^{3}\right\}_{\unicode[STIX]{x1D709}}D\right], & \displaystyle\end{eqnarray}$$

and for $\unicode[STIX]{x1D709}>\unicode[STIX]{x1D709}_{\ast }$ by

(2.12c ) $$\begin{eqnarray}\displaystyle q_{H}=-(1+\unicode[STIX]{x1D700})H^{3}H_{\unicode[STIX]{x1D709}}. & & \displaystyle\end{eqnarray}$$

The boundary conditions on (2.11) and (2.12) for $\unicode[STIX]{x1D700}>0$ can be expressed in the dimensionless variables as

(2.13a ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{prescribed influx of upper fluid:}\quad \unicode[STIX]{x1D709}q_{D}\rightarrow {\mathcal{Q}}\quad \text{as }\unicode[STIX]{x1D709}\rightarrow 0, & \displaystyle\end{eqnarray}$$

(2.13b ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{zero influx of lower fluid:}\quad q_{H}-q_{D}\rightarrow 0\quad \text{as }\unicode[STIX]{x1D709}\rightarrow 0, & \displaystyle\end{eqnarray}$$

(2.13c ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{fixed initial height at infinity:}\quad H\rightarrow 1\quad \text{as }\unicode[STIX]{x1D709}\rightarrow \infty , & \displaystyle\end{eqnarray}$$

(2.13d ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{vanishing upper height at nose:}\quad D=0\quad \text{at }\unicode[STIX]{x1D709}=\unicode[STIX]{x1D709}_{\ast }, & \displaystyle\end{eqnarray}$$

(2.13e ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{mass conservation through nose:}\quad (\unicode[STIX]{x1D709}_{\ast })_{\unicode[STIX]{x1D70F}}=\frac{q_{D}}{D}-\frac{1}{2}\unicode[STIX]{x1D709}\quad \text{at }\unicode[STIX]{x1D709}=\unicode[STIX]{x1D709}_{\ast }, & \displaystyle\end{eqnarray}$$

(2.13f ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{continuity of total height:}\quad [H]_{-}^{+}=0\quad \text{at }\unicode[STIX]{x1D709}=\unicode[STIX]{x1D709}_{\ast }, & \displaystyle\end{eqnarray}$$

(2.13g ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{continuity of total flux:}\quad [q_{H}]_{-}^{+}=0\quad \text{at }\unicode[STIX]{x1D709}=\unicode[STIX]{x1D709}_{\ast }, & \displaystyle\end{eqnarray}$$

$(\unicode[STIX]{x1D709}_{\ast })_{\unicode[STIX]{x1D70F}}$

(2.14) $$\begin{eqnarray}{\mathcal{Q}}=\frac{3\unicode[STIX]{x1D707}_{2}Q}{2\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}_{1}gh_{\infty }^{4}}\end{eqnarray}$$

is the non-dimensionalised volumetric influx.

The dimensionless influx ${\mathcal{Q}}$ can be related to the ratio of the height scale for a single-layer current $h_{sl}=[\unicode[STIX]{x1D707}_{1}Q/(\unicode[STIX]{x1D70C}_{1}g)]^{1/4}$ (Huppert Reference Huppert1982b ) and the far-field height $h_{\infty }$ by noting that ${\mathcal{Q}}\propto m(h_{sl}/h_{\infty })^{4}$ .

At late times the evolution of the system becomes self-similar and independent of $\unicode[STIX]{x1D70F}$ . In this limit (2.11) and (2.12) reduce to coupled ordinary differential equations, where the boundary conditions are the same as (2.13) only with a zero nose speed, $(\unicode[STIX]{x1D709}_{\ast })_{\unicode[STIX]{x1D70F}}=0$ .

3.1 The flux function ${\mathcal{F}}$

As we are interested in shock formation, we consider ${\mathcal{F}}$ in the hyperbolic limit $\unicode[STIX]{x1D700}=0$ . For $\unicode[STIX]{x1D700}=0$ the flux function ${\mathcal{F}}$ is equal to the ratio of upper-layer fluid flux to total fluid flux and simplifies to

(3.3) $$\begin{eqnarray}{\mathcal{F}}=\frac{3\unicode[STIX]{x1D706}+(2m-3)\unicode[STIX]{x1D706}^{3}}{2+2(m-1)\unicode[STIX]{x1D706}^{3}}=\frac{q_{D}}{q_{H}}.\end{eqnarray}$$

For the case $m\leqslant 3/2$ we find that the flux function is concave over the entire range $0\leqslant \unicode[STIX]{x1D706}\leqslant 1$ and hence ${\mathcal{F}}^{\prime }$ decreases monotonically as $\unicode[STIX]{x1D706}$ increases, see figure 2. This implies that the thicker the upper layer as a proportion of the current, the slower its local horizontal characteristic velocity $U_{c}$ and vice versa. Therefore, a current with upper-layer thickness decreasing towards the nose is self-spreading and no shocks form in this case.

However, for $m>3/2$ there is a convex region near $\unicode[STIX]{x1D706}=0$ where ${\mathcal{F}}^{\prime }$ increases with $\unicode[STIX]{x1D706}$ . This implies that a part of the current with a thicker upper layer will now have a faster horizontal characteristic velocity than a part of the current with a thinner upper layer nearer the front of the current. Therefore, in some region of the current, the horizontal characteristic velocity will decrease towards the nose, which leads to a convergence of characteristics and hence a steepening of the interface. Ultimately, this steepening mechanism leads to the formation of shocks.

Figure 2.

(a) The flux function ${\mathcal{F}}(\unicode[STIX]{x1D706})$ for viscosity ratios $m=0.1$ , $m=1.5$ and $m=10$ (dashed, dotted and solid respectively); the long-dashed straight line is tangent to ${\mathcal{F}}$ for $m=10$ , and the point of tangency defines the critical value $\unicode[STIX]{x1D706}_{crit}$ . (b) The derivative ${\mathcal{F}}^{\prime }$ of the flux function for the same viscosity ratios $m$ .

To understand the physical origin of the convexity of the flux function ${\mathcal{F}}$ for $m>3/2$ , we consider the variation of ${\mathcal{F}}$ as $\unicode[STIX]{x1D706}\rightarrow 1$ (near the origin) and as $\unicode[STIX]{x1D706}\rightarrow 0$ (near the front of the current). First, near the origin, as $\unicode[STIX]{x1D706}\rightarrow 1$ in (3.3) we obtain the expansion

(3.4) $$\begin{eqnarray}{\mathcal{F}}\sim 1-\frac{3}{2m}(1-\unicode[STIX]{x1D706})^{2}+O(1-\unicode[STIX]{x1D706})^{3}.\end{eqnarray}$$

The leading-order term, ${\mathcal{F}}\sim 1$ , corresponds to the entire current consisting of upper-layer fluid and hence the total fluid flux $q_{H}$ being equal to the upper-fluid flux $q_{D}$ . The next-order correction is quadratic in $(1-\unicode[STIX]{x1D706})$ due to the no-slip condition at the rigid bottom boundary, and comes from integrating the leading-order linear velocity profile in the thin film of lower-layer fluid to obtain a quadratic expression for the lower-layer flux. This correction decreases as $m=\unicode[STIX]{x1D707}_{2}/\unicode[STIX]{x1D707}_{1}$ increases, since a more viscous lower-layer fluid reduces the velocity in the lower-layer fluid. The sign of the correction is such that the flux function is always concave near $\unicode[STIX]{x1D706}=1$ .

In contrast, near the nose, expanding ${\mathcal{F}}$ as $\unicode[STIX]{x1D706}\rightarrow 0$ gives

(3.5) $$\begin{eqnarray}{\mathcal{F}}\sim {\textstyle \frac{3}{2}}\unicode[STIX]{x1D706}+\left(m-{\textstyle \frac{3}{2}}\right)\unicode[STIX]{x1D706}^{3}+O(\unicode[STIX]{x1D706}^{4}).\end{eqnarray}$$

This expansion shows that if the upper layer is only a thin film, then, to leading order, it moves at a constant velocity of $3/2$ , corresponding to the maximum velocity of the underlying parabolic velocity profile. The next-order correction to this is cubic due to the stress-free boundary condition on the velocity at the top surface. The sign of this correction reflects the occurrence of convexity in the flux function for $m>3/2$ . There are two physical processes competing here: first, a larger $\unicode[STIX]{x1D706}$ corresponds to a thinner lower layer and hence a greater influence of the no-slip rigid bottom boundary, which results in a decrease in the average upper-layer velocity; but, secondly, if $m>1$ , the lower-viscosity upper-layer fluid can support a greater shear, which increases its flux in comparison to the lower-layer flux. These effects balance at $m=3/2$ .

In conclusion, we find that for $m\leqslant 3/2$ , no shocks can form, whereas for $m>3/2$ , i.e. a sufficiently lower-viscosity upper-layer fluid, there is a steepening mechanism that leads to the formation of shocks.

The critical viscosity ratio $m=3/2$ can be compared to the known result for the onset of shock formation in miscible two-fluid flows in a Hele-Shaw cell as mentioned in Yang & Yortsos (Reference Yang and Yortsos1997) and Lajeunesse et al. (Reference Lajeunesse, Martin, Rakotomalala and Salin1997). When considering the Hele-Shaw flow, the same flux function ${\mathcal{F}}$ is found. This is due to the fact that the profile of the free-surface flow considered here can be thought of as half the profile of the flow between two rigid parallel plates, which has zero stress at the mid-plane between the two plates due to its symmetry. However, in contrast to the situation in a Hele-Shaw cell, the two-layer gravity current we are considering admits a simple diffusive regularisation via small non-zero density differences. This regularisation allows an analytical derivation of shock conditions, which is not available in the Hele-Shaw flow, for which the regularisation proposed by Yang & Yortsos (Reference Yang and Yortsos1997) relies on solutions to the full Stokes equations around the nose of the intruding finger. As we shall see, this has a significant effect on the location of the shock.

3.2 Shock conditions

We have established a steepening mechanism through which shocks can form for $m>3/2$ . It remains to discuss how these shocks evolve. The general theory of shocks in hyperbolic differential equations, see for example LeFloch (Reference LeFloch2002), characterises shocks according to the behaviour of the associated flux functions and derives conditions on the formation and evolution of these shocks from entropy considerations at the shock and the corresponding regularisations.

Applying this general theory to the equal-density limit of a two-layer viscous gravity current, we find two constraints for any shock: a condition giving the shock speed and a separate entropy condition governing which shocks are allowed according to the directions of the characteristics. Suppose, for example, we have a shock at position $\unicode[STIX]{x1D709}_{s}$ extending from $\unicode[STIX]{x1D706}_{r}$ to $\unicode[STIX]{x1D706}_{l}$ , where $\unicode[STIX]{x1D706}_{r}$ is the limit of $\unicode[STIX]{x1D706}$ as $\unicode[STIX]{x1D709}\rightarrow \unicode[STIX]{x1D709}_{s}$ from the right and $\unicode[STIX]{x1D706}_{l}$ the same limit but from the left.

The first shock condition, often referred to as a Rankine–Hugoniot condition, is given by mass conservation across the shock. Equating the fluxes across the shock in the frame of the shock, gives the speed of the shock as

(3.6) $$\begin{eqnarray}U_{s}\biggm\vert_{\unicode[STIX]{x1D706}_{r},\unicode[STIX]{x1D706}_{l}}=\frac{q_{H}}{H}\frac{{\mathcal{F}}_{r}-{\mathcal{F}}_{l}}{\unicode[STIX]{x1D706}_{r}-\unicode[STIX]{x1D706}_{l}}-\frac{1}{2}\unicode[STIX]{x1D709},\end{eqnarray}$$

where ${\mathcal{F}}_{r}={\mathcal{F}}(\unicode[STIX]{x1D706}_{r})$ and ${\mathcal{F}}_{l}={\mathcal{F}}(\unicode[STIX]{x1D706}_{l})$ .

The second condition is an entropy condition which, in simple hyperbolic systems, is the Lax entropy condition

(3.7) $$\begin{eqnarray}U_{c}|_{\unicode[STIX]{x1D706}_{r}}\leqslant U_{s}\leqslant U_{c}|_{\unicode[STIX]{x1D706}_{l}}.\end{eqnarray}$$

Using (3.6) this can be simplified to

(3.8) $$\begin{eqnarray}{\mathcal{F}}_{r}^{\prime }\leqslant \frac{{\mathcal{F}}_{r}-{\mathcal{F}}_{l}}{\unicode[STIX]{x1D706}_{r}-\unicode[STIX]{x1D706}_{l}}\leqslant {\mathcal{F}}_{l}^{\prime }.\end{eqnarray}$$

This condition states that characteristics at the edge of the shock cannot move away from the shock. For a decreasing jump discontinuity, i.e.  $\unicode[STIX]{x1D706}_{r}<\unicode[STIX]{x1D706}_{l}$ , it implies that a strictly convex flux function leads to a shock, whilst a strictly concave flux function leads to a rarefaction wave, and vice versa for a increasing jump discontinuity $\unicode[STIX]{x1D706}_{r}>\unicode[STIX]{x1D706}_{l}$ .

However, if the flux function has both regions of convexity and concavity then the Lax entropy condition no longer ensures a unique solution, as it only regards properties at the edge of the shock. A shock may split into multiple shocks connected by rarefaction waves depending on the local convexity or concavity of the flux function. Therefore, we instead consider the more general Oleinik entropy condition, which not only restricts the speed of the characteristics at the edge of the shock, but also ensures that the shock is internally consistent. Internal consistency here means that, if we were to consider the shock splitting into multiple smaller shocks, these would remain together, not allowing the formation of rarefaction waves in between them and hence remaining as a single shock. The Oleinik entropy condition can be written as

(3.9) $$\begin{eqnarray}U_{s}|_{\unicode[STIX]{x1D706}_{r},\unicode[STIX]{x1D706}_{i}}\leqslant U_{s}|_{\unicode[STIX]{x1D706}_{r},\unicode[STIX]{x1D706}_{l}}\leqslant U_{s}|_{\unicode[STIX]{x1D706}_{i},\unicode[STIX]{x1D706}_{l}},\quad \text{for all }\unicode[STIX]{x1D706}_{i}\text{ with }\unicode[STIX]{x1D706}_{r}\leqslant \unicode[STIX]{x1D706}_{i}\leqslant \unicode[STIX]{x1D706}_{l},\end{eqnarray}$$

which again can be simplified to

(3.10) $$\begin{eqnarray}\frac{{\mathcal{F}}_{r}-{\mathcal{F}}_{i}}{\unicode[STIX]{x1D706}_{r}-\unicode[STIX]{x1D706}_{i}}\leqslant \frac{{\mathcal{F}}_{r}-{\mathcal{F}}_{l}}{\unicode[STIX]{x1D706}_{r}-\unicode[STIX]{x1D706}_{l}}\leqslant \frac{{\mathcal{F}}_{i}-{\mathcal{F}}_{l}}{\unicode[STIX]{x1D706}_{i}-\unicode[STIX]{x1D706}_{l}},\quad \text{for all }\unicode[STIX]{x1D706}_{i}\text{ with }\unicode[STIX]{x1D706}_{r}\leqslant \unicode[STIX]{x1D706}_{i}\leqslant \unicode[STIX]{x1D706}_{l}.\end{eqnarray}$$

Note that the Oleinik entropy condition (3.9) implies the Lax entropy condition (3.7) by taking the limits $\unicode[STIX]{x1D706}_{i}\rightarrow \unicode[STIX]{x1D706}_{r}$ and $\unicode[STIX]{x1D706}_{i}\rightarrow \unicode[STIX]{x1D706}_{l}$ .

To understand the implication of the Oleinik entropy condition for the heights $\unicode[STIX]{x1D706}$ in which shocks can occur, we can think of a shock as a chord from $(\unicode[STIX]{x1D706}_{r},{\mathcal{F}}_{r})$ to $(\unicode[STIX]{x1D706}_{l},{\mathcal{F}}_{l})$ on the flux function curve, see figure 3(a). This representation as a chord is warranted by the fact that the shock speed $U_{s}$ is related to the slope of these chords in the same way that the speed of characteristics $U_{c}$ is related to the tangent slope ${\mathcal{F}}^{\prime }$ of the flux function. The Oleinik entropy condition implies that chords with $\unicode[STIX]{x1D706}_{r}<\unicode[STIX]{x1D706}_{l}$ must lie entirely above the flux function, whilst chords with $\unicode[STIX]{x1D706}_{r}>\unicode[STIX]{x1D706}_{l}$ must lie entirely below it.

We observe that the upper-layer fluid velocity (2.5a ) increases monotonically from the bottom to the top of the upper layer at every position $\unicode[STIX]{x1D709}$ and hence so must the distance travelled by the upper-layer fluid. This implies that the interface is monotonically sloped upwards and it follows that we must have $\unicode[STIX]{x1D706}_{r}<\unicode[STIX]{x1D706}_{l}$ . (Note that this condition is for a current propagating to the right away from the source; it would need adaptation for a current propagating to the left.)

Hence, as we must have $\unicode[STIX]{x1D706}_{r}<\unicode[STIX]{x1D706}_{l}$ , any chord corresponding to a shock must lie above the flux function and therefore the largest admissible shock is given by the chord that goes through the origin and is tangent to the flux function ${\mathcal{F}}$ . This tangent construction defines a critical height $\unicode[STIX]{x1D706}_{crit}>0$ by

(3.11) $$\begin{eqnarray}{\mathcal{F}}^{\prime }(\unicode[STIX]{x1D706}_{crit})=\frac{{\mathcal{F}}(\unicode[STIX]{x1D706}_{crit})}{\unicode[STIX]{x1D706}_{crit}},\end{eqnarray}$$

which can be solved analytically to obtain

(3.12) $$\begin{eqnarray}\unicode[STIX]{x1D706}_{crit}=2\left(\frac{2m}{3}-1\right)^{-1/2}\sinh \left[\frac{1}{3}\sinh ^{-1}\left\{(m-1)^{-1}\left(\frac{2m}{3}-1\right)^{3/2}\right\}\right].\end{eqnarray}$$

Hence any admissible shock must satisfy $0\leqslant \unicode[STIX]{x1D706}_{r}<\unicode[STIX]{x1D706}_{l}\leqslant \unicode[STIX]{x1D706}_{crit}$ .

Figure 3(a) shows three chords on the flux function which all satisfy the Oleinik entropy condition (3.10), but only two of which correspond to admissible shocks with $\unicode[STIX]{x1D706}_{r}<\unicode[STIX]{x1D706}_{l}$ . The two admissible shocks shown can be associated with distinct types: the chord $\mathbf{QR}$ , which has $0<\unicode[STIX]{x1D706}_{r}<\unicode[STIX]{x1D706}_{l}<\unicode[STIX]{x1D706}_{crit}$ , corresponds to an internal shock; the chord $\mathbf{OP}$ , with $0=\unicode[STIX]{x1D706}_{r}<\unicode[STIX]{x1D706}_{l}<\unicode[STIX]{x1D706}_{crit}$ , corresponds to a frontal shock at the top surface. The third chord $\mathbf{ST}$ does not correspond to an admissible shock, since to satisfy the Oleinik entropy condition (3.10) we must have $\unicode[STIX]{x1D706}_{r}>\unicode[STIX]{x1D706}_{l}$ for this chord, contradicting the requirement that the interface is monotonically sloped upwards.

Figure 3.

(a) The flux function ${\mathcal{F}}$ for $m=100$ with dashed chords representing hypothetical shocks: $\mathbf{OP}$ corresponds to a frontal shock with heights $(\unicode[STIX]{x1D706}_{r},\unicode[STIX]{x1D706}_{l})=(0,0.08)$ , $\mathbf{QR}$ to an internal shock with $(\unicode[STIX]{x1D706}_{r},\unicode[STIX]{x1D706}_{l})=(0.08,0.18)$ and $\mathbf{ST}$ to a shock violating $\unicode[STIX]{x1D706}_{r}<\unicode[STIX]{x1D706}_{l}$ with $(\unicode[STIX]{x1D706}_{r},\unicode[STIX]{x1D706}_{l})=(0.3,0.2)$ . Again, the long-dashed straight line is tangent to ${\mathcal{F}}$ , and defines $\unicode[STIX]{x1D706}_{crit}$ . (b) The corresponding local travelling-wave solutions $\unicode[STIX]{x1D706}(\unicode[STIX]{x1D702})$ and their limits $\unicode[STIX]{x1D706}_{r}$ and $\unicode[STIX]{x1D706}_{l}$ . Note that we show $\unicode[STIX]{x1D706}$ increasing downwards to aid interpretation of the curves as the physical shapes of the interface relative to a top surface at $\unicode[STIX]{x1D706}=0$ .

3.3 Regularised shock structures

To derive the Oleinik entropy condition from physical principles and to understand the local regularised structures of any shock that might occur, we recall that small density differences become significant near regions with a large slope of the interface, i.e. in particular near any shock. These regions can be analysed by defining a local coordinate $\unicode[STIX]{x1D702}$ in a frame travelling with the shock by

(3.13) $$\begin{eqnarray}\unicode[STIX]{x1D709}=\unicode[STIX]{x1D709}_{s}+\unicode[STIX]{x1D700}\unicode[STIX]{x1D702}\left.\frac{H^{4}}{q_{H}}\right|_{\unicode[STIX]{x1D709}=\unicode[STIX]{x1D709}_{s}},\end{eqnarray}$$

where the factor $\unicode[STIX]{x1D700}H^{4}/q_{H}|_{\unicode[STIX]{x1D709}=\unicode[STIX]{x1D709}_{s}}$ is the horizontal scale of the local travelling wave. We express (3.1) governing the relative fluid fraction $\unicode[STIX]{x1D706}$ in terms of this local coordinate, and expand in orders of $\unicode[STIX]{x1D700}$ . Noting that $(\unicode[STIX]{x1D709}_{s})_{\unicode[STIX]{x1D70F}}=U_{s}$ , we obtain at leading order the equation

(3.14) $$\begin{eqnarray}-U_{s}\frac{\text{d}\unicode[STIX]{x1D706}}{\text{d}\unicode[STIX]{x1D702}}+U_{c}\frac{\text{d}\unicode[STIX]{x1D706}}{\text{d}\unicode[STIX]{x1D702}}=\left.\frac{q_{H}}{H}\right|_{\unicode[STIX]{x1D709}=\unicode[STIX]{x1D709}_{s}}\frac{\text{d}}{\text{d}\unicode[STIX]{x1D702}}\left({\mathcal{D}}\frac{\text{d}\unicode[STIX]{x1D706}}{\text{d}\unicode[STIX]{x1D702}}\right).\end{eqnarray}$$

The boundary conditions for this ordinary differential equation are $\unicode[STIX]{x1D706}\rightarrow \unicode[STIX]{x1D706}_{r,l}$ as $\unicode[STIX]{x1D702}\rightarrow \pm \infty$ respectively.

We note that the vertical characteristic velocity $W_{c}$ does not influence the local dynamics near vertical shocks to leading order. This can be explained by noting that $W_{c}$ arises from the divergence in the total flux $q_{H}$ and is independent of the interfacial slope $\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D709}}$ . Hence it remains $O(1)$ even close to a near-vertical interface.

We can simplify (3.14) and integrate once to obtain

(3.15) $$\begin{eqnarray}{\mathcal{D}}\frac{\text{d}\unicode[STIX]{x1D706}}{\text{d}\unicode[STIX]{x1D702}}=-\frac{{\mathcal{F}}_{r}-{\mathcal{F}}_{l}}{\unicode[STIX]{x1D706}_{r}-\unicode[STIX]{x1D706}_{l}}(\unicode[STIX]{x1D706}-\unicode[STIX]{x1D706}_{r})+[{\mathcal{F}}(\unicode[STIX]{x1D706})-{\mathcal{F}}_{r}].\end{eqnarray}$$

This equation has fixed points where the right-hand side is zero. Hence, in particular, $\unicode[STIX]{x1D706}=\unicode[STIX]{x1D706}_{r}$ and $\unicode[STIX]{x1D706}=\unicode[STIX]{x1D706}_{l}$ are both fixed points. To ensure that a heteroclinic connection exists between these fixed points, there must not be any other fixed point $\unicode[STIX]{x1D706}_{i}$ with $\unicode[STIX]{x1D706}_{r}<\unicode[STIX]{x1D706}_{i}<\unicode[STIX]{x1D706}_{l}$ . This is equivalent to the Oleinik entropy condition (3.10), which means that the only admissible solutions under the diffusive regularisation $\unicode[STIX]{x1D700}>0$ , are exactly those with shocks satisfying the Oleinik entropy condition.

Figure 3(b) shows the solutions to (3.15) corresponding to the shocks depicted as chords $\mathbf{OP}$ , $\mathbf{QR}$ and $\mathbf{ST}$ in figure 3(a). Without loss of generality, the origin $\unicode[STIX]{x1D702}=0$ was chosen to be at the average of $\unicode[STIX]{x1D706}_{r}$ and $\unicode[STIX]{x1D706}_{l}$ . In particular, the local travelling-wave solution to $\mathbf{ST}$ highlights the result that for this case $\unicode[STIX]{x1D706}_{r}>\unicode[STIX]{x1D706}_{l}$ and hence the interface slopes downwards rather than upwards.

We note at this point that the regularised shock structures corresponding to a frontal shock with $\unicode[STIX]{x1D706}_{r}=0$ , e.g. the chord $\mathbf{OP}$ , have a square-root singularity at the front. This is consistent with the square-root frontal singularity obtained by Lister & Kerr (Reference Lister and Kerr1989) for two-layer gravity currents with non-zero density differences and is the result corresponding to the cube-root frontal singularity for a single-layer gravity current on a rigid surface as seen in Huppert (Reference Huppert1982b ).

4.1 A simplified model

To illustrate the evolution towards self-similarity and, in particular, the process of shock formation, we consider a simplified system, where we prescribe both the total height $H(\unicode[STIX]{x1D709})$ and the total flux $q_{H}(\unicode[STIX]{x1D709})$ . This corresponds to a fixed evolution of the top surface. To ensure that the total fluid volume is still conserved, we choose $H$ and $q_{H}$ such that they still satisfy

(4.1) $$\begin{eqnarray}{\textstyle \frac{1}{2}}\unicode[STIX]{x1D709}^{2}H_{\unicode[STIX]{x1D709}}=(\unicode[STIX]{x1D709}q_{H})_{\unicode[STIX]{x1D709}}.\end{eqnarray}$$

Prescribing $H$ and $q_{H}$ decouples (3.1) and has the further physical consequence that the pressure is no longer hydrostatic and given by $H$ . Instead the pressure is given by the profile necessary to drive a mass-conserving flow in accordance with the prescribed evolution of $H$ and $q_{H}$ . This means we are effectively solving (3.1) with $\unicode[STIX]{x1D700}=0$ and under the simplifying assumption that both $H$ and $q_{H}$ are known, which we argue captures the essential behaviour of the advective evolution of the interface.

Suppose $H$ and $q_{H}$ are chosen to match the similarity solutions to the full coupled system, assuming that we have calculated these before. Then the interfacial position $\unicode[STIX]{x1D706}(\unicode[STIX]{x1D709},\unicode[STIX]{x1D70F})$ evolves under (3.1) from any initial condition to approach the similarity solution of the coupled system as $\unicode[STIX]{x1D70F}\rightarrow \infty$ . Furthermore, the pressure gradient approaches the hydrostatic pressure gradient in the same limit.

The time-dependent solution for $\unicode[STIX]{x1D706}$ can be found numerically via the method of characteristics for which we introduce characteristic curves $(\unicode[STIX]{x1D709}_{c}(\unicode[STIX]{x1D70F};\unicode[STIX]{x1D6FC}),\unicode[STIX]{x1D706}_{c}(\unicode[STIX]{x1D70F};\unicode[STIX]{x1D6FC}))$ , where $\unicode[STIX]{x1D6FC}$ is a Lagrangian label for the characteristics. Any of these characteristic curves evolve from some initial condition according to

(4.2a,b ) $$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D709}_{c}}=U_{c}(\unicode[STIX]{x1D709}_{c},\unicode[STIX]{x1D706}_{c}),\quad \dot{\unicode[STIX]{x1D706}_{c}}=W_{c}(\unicode[STIX]{x1D709}_{c},\unicode[STIX]{x1D706}_{c}), & & \displaystyle\end{eqnarray}$$

where $U_{c}$ and $W_{c}$ as functions of $\unicode[STIX]{x1D709}$ and $\unicode[STIX]{x1D706}$ are given by (3.1) and (3.2).

As soon as shocks are encountered, as indicated by the intersection of characteristics, the evolution of the shock has to be treated separately to (3.1) taking into consideration the shock conditions (3.6) and (3.10), which determine the shock speed $U_{s}$ and limit the maximal shock height by $\unicode[STIX]{x1D706}_{crit}$ .

4.2 An illustrative example calculation

Figure 4.

A time series of interfacial shapes $\unicode[STIX]{x1D706}(\unicode[STIX]{x1D709},\unicode[STIX]{x1D70F})$ for $m=5$ in the simplified fixed-lid system, with points of interest $\mathbf{A}$ , $\mathbf{B}$ , $\mathbf{C}$ and $\mathbf{N}$ marked (see text). The shapes are shown at equal time intervals $\unicode[STIX]{x0394}\unicode[STIX]{x1D70F}=0.16$ . The boundaries of the shock and the characteristic emanating from $\mathbf{C}$ are highlighted. Note that as in figure 3(b) we have reversed the $\unicode[STIX]{x1D706}$ -axis to aid physical interpretation.

Figure 5.

The characteristics corresponding to the solution shapes in figure 4, focusing on the region near the shock by using coordinates $(\unicode[STIX]{x1D709}-\unicode[STIX]{x1D709}_{\ast },\unicode[STIX]{x1D70F}^{1/2})$ , where $\unicode[STIX]{x1D709}_{\ast }(\unicode[STIX]{x1D70F})$ is the shock position once it has formed and before that the radial position of the characteristic leading up to the shock. The same points of interest $\mathbf{A}$ to $\mathbf{C}$ are marked (see text).

A time series of interfacial shapes $\unicode[STIX]{x1D706}(\unicode[STIX]{x1D709},\unicode[STIX]{x1D70F})$ for the viscosity ratio $m=5$ is shown in figure 4. For illustration, the initial condition was chosen as $\unicode[STIX]{x1D706}(\unicode[STIX]{x1D709},0)=(1-3\unicode[STIX]{x1D709}/2)^{2}$ . The interfacial shape clearly evolves towards a final steady state, which is the anticipated similarity solution. The corresponding characteristics are shown in figure 5, focussing on the region that forms a shock. In both figures, several points of interest have been marked, whose significance is explained below.

Initially the profile is shock free. However, since $m>3/2$ the flux function ${\mathcal{F}}$ has a convex region where ${\mathcal{F}}^{\prime }$ increases with $\unicode[STIX]{x1D706}$ . Therefore, $U_{c}$ also increases with $\unicode[STIX]{x1D706}$ in this region and the steepening mechanism, associated with shock formation, applies in the upper part of the current (see figure 4). The interface steepens until the point $\mathbf{A}$ , where the interface becomes vertical, which constitutes the formation of a shock and corresponds to the first point where the characteristics intersect in figure 5. For the initial condition chosen for this illustration, the shock forms in the interior of the fluid, but other initial conditions can lead to a shock forming at the surface.

Once formed, the shock moves forward with velocity $U_{s}$ . The parts just above and below the shock continue to steepen, eventually also becoming vertical and causing the height of the shock to grow. The edges of the shock $\unicode[STIX]{x1D706}_{l,r}$ are given by the condition of a continuous interface, i.e.  $\unicode[STIX]{x1D706}(\unicode[STIX]{x1D709})\rightarrow \unicode[STIX]{x1D706}_{l,r}$ as $\unicode[STIX]{x1D709}\rightarrow \unicode[STIX]{x1D709}_{s}$ from below or above respectively. Hence, their evolution can be calculated from

(4.3) $$\begin{eqnarray}\dot{\unicode[STIX]{x1D706}}_{l,r}=\{(U_{s}-U_{c})\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D709}}+W_{c}\}|_{\unicode[STIX]{x1D709}=\unicode[STIX]{x1D709}_{\ast },\unicode[STIX]{x1D706}=\unicode[STIX]{x1D706}_{l,r}}.\end{eqnarray}$$

At point $\mathbf{B}$ (figures 4 and 5) the growth of the shock reaches the top surface, i.e.  $\unicode[STIX]{x1D706}_{r}=0$ . This point corresponds to where the shock catches up with the characteristic on the top surface. From $\mathbf{B}$ onwards the upper boundary of the shock is fixed and can no longer grow. Correspondingly there are no more characteristics intersecting the shock from the right in figure 5.

Now only the lower boundary of the shock, $\unicode[STIX]{x1D706}_{l}$ continues to grow until it reaches the critical height $\unicode[STIX]{x1D706}_{crit}$ at point $\mathbf{C}$ . The shock has now reached maximal size, but the shock speed $U_{s}$ at $\mathbf{C}$ is non-zero and the interface continues to evolve.

The evolution from $\mathbf{C}$ onwards contains a new feature, absent from classical shocks: a region associated with a one-sided characteristic shock. The characteristic going through the point $\mathbf{C}$ moves away from the shock and leaves a region between itself and the shock, which needs to be computed differently. This region is filled by characteristics emanating from the shock itself. These initially have the same horizontal speed $U_{c}$ as the shock speed $U_{s}$ , but a non-zero vertical speed $W_{c}$ , which carries them into the region that needs filling. The characteristics emanating from the shock carry away the information $\unicode[STIX]{x1D706}_{l}=\unicode[STIX]{x1D706}_{crit}$ , i.e. that the shock cannot grow beyond $\unicode[STIX]{x1D706}_{crit}$ . Ultimately, the entire shape of the interface is affected by this information that the shock height is limited. We note that if the shock were initially bigger than this critical height $\unicode[STIX]{x1D706}_{crit}$ , then the vertical part below the critical height would fall behind the shock as its characteristic speed is less than the shock speed. In this situation, the interface would thus flatten and again establish the critical shock height.

We recall that the shock speed $U_{s}$ also depends on the position $\unicode[STIX]{x1D709}$ and hence the shock speed continues to change after point $\mathbf{C}$ , even though the shock height is fixed. The speed slows down as $\unicode[STIX]{x1D70F}\rightarrow \infty$ , allowing an approach to a final steady-state similarity solution at point $\mathbf{N}$ .

4.3 Altered boundary conditions for coupled steady similarity solutions in the limit of equal densities

The illustrative example described above shows that for $m>3/2$ the steady-state similarity solution has a frontal shock of height $\unicode[STIX]{x1D706}_{crit}$ . This changes the boundary condition (2.13d ) of zero upper-layer height at the nose. Mathematically, we note that the horizontal velocity of the characteristics at a shock in steady state must be equal to the speed of the shock, since otherwise there would be steepening or flattening of the interface resulting in an unsteady profile. Furthermore, in steady state these velocities have to be zero in the similarity frame and hence $U_{s}=U_{c}=0$ at any shock. The only shock that satisfies this is a frontal shock with height $\unicode[STIX]{x1D706}_{crit}$ . Therefore, we relax the boundary condition (2.13d ) that $\unicode[STIX]{x1D706}=0$ at $\unicode[STIX]{x1D709}=\unicode[STIX]{x1D709}_{\ast }$ and instead replace it by

(4.4) $$\begin{eqnarray}\unicode[STIX]{x1D706}=\left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x1D706}_{crit}(m), & m>3/2\\ 0, & m\leqslant 3/2\end{array}\right.\quad \text{at }\unicode[STIX]{x1D709}=\unicode[STIX]{x1D709}_{\ast }.\end{eqnarray}$$

We also note that the point $(\unicode[STIX]{x1D709},\unicode[STIX]{x1D706})=(0,1)$ is the only stable fixed point of the hyperbolic equation (3.1) with $\unicode[STIX]{x1D700}=0$ and hence all characteristics will tend to that point. This implies that the boundary condition (2.13b ), requiring $\unicode[STIX]{x1D706}=1$ at $\unicode[STIX]{x1D709}=0$ , is satisfied automatically and can therefore be ignored when solving the equations numerically. Having one fewer boundary condition matches the reduction in order of the differential equations in the limit $\unicode[STIX]{x1D700}=0$ . We note that the redundancy of the boundary condition (2.13b ) in the limit of zero density difference is not apparent in the steady similarity equations, but only becomes apparent in the time-dependent system.

6.1 Experimental method

Figure 7.

Schematic diagram of the experimental set-up. Fluid was injected through a hole in the bottom of a tank such that it spreads radially outwards over the top of a lower layer of fluid. The injected fluid was provided at a constant flow rate by means of a peristaltic pump and was dyed green to allow visualisation of the moving nose of the upper spreading layer.

We conducted a suite of experiments in which one fluid was injected into a pre-existing lower layer of another fluid of equal density, but differing viscosity. A schematic diagram of the experimental apparatus is shown in figure 7. A lower-layer fluid was introduced into a tank of square cross-section 48.5 cm $\times$ 48.5 cm, with side walls of height 5 cm, and its free surface allowed to become horizontal under gravity. Once the lower layer had attained a uniform height, a density-matched second fluid was injected through a 1 cm diameter hole in the bottom of the tank at a constant flow rate such that the resulting current spread radially outwards over the top of the lower layer. The flow rate of the injected fluid was set by a peristaltic pump and the flux recorded by a mass balance.

Aqueous solutions of glycerol and salt were used for both fluids. The viscosities of the fluids were varied by altering the concentration of glycerol between 0.1 and 2.3 Pa s. A Bohlin rheometer was used to measure the viscosities as functions of temperature and the temperature in the laboratory was carefully monitored throughout the experiments. The densities of the fluids were matched by adding salt, using a hydrometer to measure the densities. In order to enable the detection of the contact between the two fluid layers the injected fluid was dyed green. The evolution of the nose position of the propagating current could then be determined from images.

The axisymmetric profiles of the total thickness of the current, $h$ , were determined by digitally imaging the deflection of a line, of approximately 1 mm width, projected onto the top surface of the fluids (titanium dioxide was added to the fluids as an opacifier). The line was imaged using a Nikon D300 camera, with a spatial resolution of 0.21 mm, positioned at an angle $\unicode[STIX]{x1D703}=29^{\circ }$ above the horizontal and perpendicular to the projected line. The deflection of the line was measured with respect to a reference image of the projected line on the rigid lower surface with no lower layer, and subpixel accuracy was achieved by fitting a Gaussian profile across the line (see Lister et al. (Reference Lister, Peng and Neufeld2013)). Image analysis was then used to determine the evolution of the total-height profiles of the propagating current. The height and uniformity of the lower layer prior to injection was also checked to within 1 %, using images of the projected line.

To check if lubrication theory is indeed applicable to model these experiments we estimate a modified Reynolds number $\mathit{Re}=Q\unicode[STIX]{x1D70C}h_{\infty }/(\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}_{1}r_{hole}^{2})$ , where the velocity scale is assumed to be the vertical outflow velocity through the hole of radius $r_{hole}=0.5~\text{cm}$ , and the length scale is assumed to be of the order of the initial height $h_{\infty }$ . Using this definition and the values of the physical parameters given in table 1, we find that the Reynolds numbers vary from $0.003$ to $0.02$ between the experiments. This indicates that the experiments are indeed low-Reynolds-number flows and should be well described by lubrication theory.

Furthermore, we calculate Péclet numbers for the experiments as $\mathit{Pe}=\mathit{Sc}\,\mathit{Re}$ , where the Schmidt number is given by $\mathit{Sc}=\unicode[STIX]{x1D707}_{1}/(\unicode[STIX]{x1D70C}D_{glycerol})$ . The diffusivity of water in glycerol is estimated to be $D_{glycerol}\lesssim 10^{-9}~\text{m}^{2}~\text{s}^{-1}$  (D’Errico et al. Reference D’Errico, Ortona, Capuano and Vitagliano2004), and therefore $\mathit{Sc}\gtrsim O(10^{6})$ . Together with the previous estimates for the Reynolds numbers we conclude that the Péclet numbers vary from $5\times 10^{3}$ to $3\times 10^{4}$ and hence it is appropriate to neglect the effects of molecular diffusion of glycerol, and hence the viscosity acts as a passive tracer.

In all of the experiments presented here, the viscosity of the lower layer is greater than or equal to the viscosity of the injected fluid $\unicode[STIX]{x1D707}_{2}\geqslant \unicode[STIX]{x1D707}_{1}$ , such that $m\leqslant 1$ , thereby ensuring that the flow remains axisymmetric. Observations of currents with a less viscous intruding fluid, i.e. $m>1$ , show that the interface can become unstable in a manner analogous to the symmetry-breaking viscous fingering instability in Hele-Shaw cells, and the axisymmetric model derived above is no longer applicable for these cases. Experimental measurements of the growth of this instability and a linear stability analysis is beyond the scope of this study and will instead be the focus of a future communication.

Table 1.

The values in each experiment of the initial height $h_{\infty }$ , influx $Q$ , viscosities $\unicode[STIX]{x1D707}_{i}$ and densities $\unicode[STIX]{x1D70C}_{1}=\unicode[STIX]{x1D70C}_{2}$ , with uncertainties: $\pm 1\,\%$ for $h_{\infty }$ ; $\pm 0.1~\text{cm}^{3}~\text{s}^{-1}$ for $Q$ ; $\pm 5\,\%$ for $\unicode[STIX]{x1D707}_{i}$ ; and $\pm 0.01~\text{g}~\text{cm}^{-3}$ for $\unicode[STIX]{x1D70C}_{i}$ . Also, derived quantities such as the viscosity ratio $m$ , the non-dimensional influx ${\mathcal{Q}}$ and the experimentally measured and numerically computed nose positions, $\unicode[STIX]{x1D709}_{\ast }^{exp}$ and $\unicode[STIX]{x1D709}_{\ast }^{num}$ respectively, and their relative deviation $\unicode[STIX]{x0394}\unicode[STIX]{x1D709}_{\ast }$ . The values $\unicode[STIX]{x1D709}_{\ast }^{exp}$ were obtained from performing a least-squares fit to the time series of measured values of $r_{\ast }(t)$ .

6.2 Comparison of theory and experiments

We first check whether the measured current spreads in a self-similar manner. The experimental measurements show that, with the similarity variable given by (2.10), the profiles of the top surface collapse well onto a self-similar shape, thus supporting the application of the theory developed above. A representative example of this collapse is shown in figure 8.

Figure 9 compares both the self-similar total-height profiles $H(\unicode[STIX]{x1D709})$ obtained from the numerical calculations with the measurements from the experiments. The theory captures the trend of the top-surface profiles very well over a wide range of influxes ${\mathcal{Q}}$ and viscosity ratios $m<1$ and reproduces most features of the flow. In particular, the experiments demonstrate the same trend as observed in the numerical calculations: with increasing influx, ${\mathcal{Q}}$ , the pre-existing fluid layer plays an smaller role such that the injected fluid forms a more compact current, with a sharp, distinct nose. As also predicted by the model, the viscosity ratio $m$ has a significantly smaller effect on the shape of the current than the effect of changes in the influx ${\mathcal{Q}}$ . For several values of ${\mathcal{Q}}$ and $m$ , a small bump is observed between the nose and origin in the experiments, but is absent from the numerical calculations. The cause of the bump is not clear, but it seems to have little effect on the overall profile.

In figure 10 and table 1 we compare the nose positions $\unicode[STIX]{x1D709}_{\ast }$ as measured in the experiments with the values calculated using the model above, which reveals a fairly good fit with relative deviations in the range of 5 to 20 %. The trend is again captured very well, especially considering the simplifying assumptions made, the difficulty resolving the exact nose position from the colour change, and the uncertainties in the physical parameters. Again, both the experimental data and the numerical model show that the nose position depends mainly on the influx ${\mathcal{Q}}$ with only a comparatively small variation with the viscosity ratio $m$ . In particular, both the numerical calculations and the experimental data seem to be approaching the correct asymptotics $\unicode[STIX]{x1D709}_{\ast }\sim {\mathcal{Q}}^{1/2}$ as ${\mathcal{Q}}\rightarrow 0$ and $\unicode[STIX]{x1D709}_{\ast }\sim {\mathcal{Q}}^{3/8}$ as ${\mathcal{Q}}\rightarrow \infty$ .

Overall, the agreement of the numerical and the experimental results suggests that the theory captures the relevant physical processes involved and gives a good description of the experimental results.

Figure 8.

(a) The measured profiles plotted every 10 s for experiment No.  $1$ . (b) The same profiles as functions of the similarity variable $\unicode[STIX]{x1D709}$ . A similar collapse is observed in the other experiments.

Figure 9.

Measured (solid) and computed (dotted) heights $H$ and $D$ as functions of $\unicode[STIX]{x1D709}/\unicode[STIX]{x1D709}_{\ast }^{num}$ , so the similarity variable is rescaled by the numerical nose position. The values of ${\mathcal{Q}}$ and $m$ are as in table 1.

Figure 10.

Log–log plot of the measured (symbols) and computed (curved lines) nose positions, $\unicode[STIX]{x1D709}_{\ast }^{exp}$ and $\unicode[STIX]{x1D709}_{\ast }^{num}$ , respectively, as functions of ${\mathcal{Q}}$ and rescaled by ${\mathcal{Q}}^{1/2}$ with $m$ fixed. The predicted asymptotic dependence on ${\mathcal{Q}}$ for large and small ${\mathcal{Q}}$ is $\unicode[STIX]{x1D709}_{\ast }\sim {\mathcal{Q}}^{3/8}$ and $\unicode[STIX]{x1D709}_{\ast }\sim {\mathcal{Q}}^{1/2}$ respectively (thin straight lines). The experimental data can be found in table 1.