2.1. Rescaling

Our focus is on the dynamics of the intrusion and thus we present results rescaled for buoyancy-driven propagation (see also Ceddia et al. Reference Ceddia, Kerr, Rudman, Settle and Slim2026). Specifically, we use the velocity scale $\sqrt {NQ}$ , length scale $\sqrt {Q/N}$ and time scale $1/N$ to set

(2.5) \begin{equation} \boldsymbol{v} = (u,v) =\hat {\boldsymbol{v}}/\sqrt {NQ}, \quad \boldsymbol{x}= (x,y) = \hat {\boldsymbol{x}}/\sqrt {Q/N} \quad \text{and} \quad t = N\hat {t}, \end{equation}

where unhatted variables are dimensionless. We rescale the density by $\Delta \rho = \rho _0 \sqrt {N^3Q}/g$ , the initial density difference across the typical length scale, and set the density relative to the intrusion in dimensionless form as

(2.6) \begin{align} \rho _{\textit{rel}}(x,y,t) = \frac {\hat {\rho }-\rho _0}{\Delta \rho }. \end{align}

The initial density difference profile is thus $\rho _{\textit{rel}}(x,y,0) = -y$ . Finally, we set the dimensionless pressure relative to the hydrostatic contribution $-\rho _0 g y$ as

(2.7) \begin{align} p_{\textit{rel}} = \frac {\hat {p}+ \rho _0 gy}{\rho _0 NQ}. \end{align}

The governing equations in dimensionless form are now

(2.8a) \begin{align} \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{v} &= 0, \end{align}

(2.8b) \begin{align} \frac {\partial \boldsymbol{v}}{\partial t} + (\boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{v} \ &= -\boldsymbol{\nabla } p_{\textit{rel}} + \frac {1}{\textit{Re}} {\nabla} ^2 \boldsymbol{v} - \rho _{\textit{rel}} \,\boldsymbol{e}_y , \end{align}

(2.8c) \begin{align} \frac {\partial \rho _{\textit{rel}}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }\left (\rho _{\textit{rel}} \boldsymbol{v} \right ) &= \frac {1}{\textit{Re} \, {\textit {Sc}}} {\nabla} ^2 \rho _{\textit{rel}} \qquad \text{and} \end{align}

(2.8d) \begin{align} \frac {\partial c}{\partial {t}}+{\boldsymbol{\nabla }} \boldsymbol{\cdot }\left (c{\boldsymbol{v}}\right ) &= \frac {1}{\textit{Re} \, {\textit {Sc}}}{{\nabla} }^2 c. \end{align}

The rescaling leaves five dimensionless parameters: the Reynolds number ${\textit{Re}}=Q/\nu$ , the Schmidt number ${\textit {Sc}} = \nu /D$ and the three geometric quantities

(2.9) \begin{align} \mathcal{W} = W \sqrt {N/Q}, \quad {\mathcal{H}_a} = H_a \sqrt {N/Q} \quad \text{and} \quad \mathcal{L}_a = L_a \sqrt {N/Q}. \end{align}

We consider a range of Reynolds numbers ${\textit{Re}}\geqslant 200$ and set ${\textit {Sc}}=10$ to ensure viscous diffusion dominates density diffusion. The first of the geometric quantities $\mathcal{W}$ is the dimensionless width of the inlet or alternatively the ratio of the speed of the intrusion $\sqrt {NQ}$ to the speed of the inlet flow $Q/W$ . If $\mathcal{W}$ is relatively small, then the source is jet-like, much narrower and faster than the natural buoyancy scale. Conversely, if $\mathcal{W}$ is relatively large, then the source can be considered ‘lazy’, wider and slower than the natural buoyancy scale. We consider this latter case.

The second geometric quantity $\mathcal{H}_a$ is the dimensionless height of the ambient or alternatively the ratio of the speed of internal waves $NH_a$ to the speed of the intrusion $\sqrt {NQ}$ . If $\mathcal{H}_a$ is relatively small, then the intrusion propagates faster than the fastest internal waves and we would expect the intrusion to be propagating into a minimally disturbed region. Conversely, if $\mathcal{H}_a$ is relatively large, then internal waves can propagate ahead of the intrusion. Note that the speed of mode- $m$ columnar waves in this framework is ${\mathcal{H}_a}/m\pi$ . We consider only symmetric configurations and so only modes with even values of $m$ are excited. We take a range of values ${\mathcal{H}_a}\geqslant \mathcal{W}$ .

The final geometric quantity $\mathcal{L}_a$ is the dimensionless length of the ambient. We restrict our attention to times before the intrusion feels the influence of the far boundary and thus treat this parameter as effectively infinite.

2.2. Numerical method

The equations above were solved using Semtex, a spectral element simulation code of Blackburn et al. (Reference Blackburn, Lee, Albrecht and Singh2019). The spectral element method, as implemented in Semtex, is founded on the discretisation of a two-dimensional planar domain into a series of contiguous, edge-aligned quadrilateral elements, using nodal Lagrange interpolants within each element to form a polynomial function space with $C_0$ continuity across the boundaries of the elements (Deville, Fischer & Mund Reference Deville, Fischer and Mund2002; Karniadakis & Sherwin Reference Karniadakis and Sherwin2005; Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang2007). The original Semtex code was modified to account for the addition of the advection–diffusion relation, (2.4), for the colour function $c$ .

Figure 2.

(a) A sample spectral element mesh used in two-dimensional planar simulations, consisting of 144 elements. The red lines indicate the element boundaries. (b) Legendre polynomial cardinal functions and the corresponding Gauss–Lobatto nodal point distribution for a polynomial degree $p=6$ , illustrating the interpolation basis used in each direction of an element.

Figure 2 demonstrates a sample spectral element mesh designed for our simulations. Each element has shape functions with tensor products of nth-order Lagrange interpolants through the Gauss–Lobatto–Legendre points ( $n=6$ in the figure). The basis polynomials share convergence properties with orthogonal Legendre polynomials (Blackburn et al. Reference Blackburn, Lee, Albrecht and Singh2019), with the primary aim of mesh refinement being to attain exponential convergence of solutions as the polynomial order increases. The polynomial order $p$ was set to be $p = N_p-1$ , where $N_p$ is the number of mesh points along the edge of each element. Note that the actual simulations were run with higher resolutions than shown, with up to $19\,000$ elements and $N_p=10$ .

The simulations were run on the full domain shown in figure 1 with the boundary conditions described above. An outflow condition was implemented on the right boundary, using an approach initially discussed by Dong (Reference Dong2015). In the domain $\mathcal{L}_a-10 \leq x \leq \mathcal{L}_a$ , we also implemented a sponge region that applies a forcing term to steer the flow toward a constant horizontal velocity, $ \boldsymbol{u}_0 = (1/\mathcal{H}_a, 0 )$ , and effectively mitigate potential reflections back into the computational domain.

We used a second-order backward differencing scheme outlined by Karniadakis, Israeli & Orszag (Reference Karniadakis, Israeli and Orszag1991) to integrate the Navier–Stokes equations in time. The time step was set so the Courant–Friedrichs–Lewy number, which is used to determine the numerical stability and physical validity of the simulations, stayed well below unity.

Further details on general aspects of the code, including the computational environment, mesh resolutions and implementation, are given by Blackburn et al. (Reference Blackburn, Lee, Albrecht and Singh2019). For the simulations outlined in § 3, the high-resolution direct numerical simulation runs required approximately 7 days of computational time while running on a single core. The specific processor used was an Intel Xeon Gold 5220R, with multi-threaded execution disabled. Details of code convergence and verification are given in Appendix A. The code was further verified against a version of the two-dimensional conservative finite-volume code of Rudman (Reference Rudman1998). The only difference from Rudman’s original code was the inclusion of a continuously varying density field instead of a discontinuous density interface. The profiles of the intrusions using the two methods show excellent agreement. We also performed selected three-dimensional simulations to confirm that the two-dimensional simulations are representative of the full dynamics. Details of these simulations are given in Appendix B.

2.3. Key diagnostic variables

To describe the evolution of the intrusion, we define $h$ to be the dimensionless half-thickness of the intrusion according to

(2.10) \begin{equation} h(x,t) = \frac {1}{2} \int ^{{\mathcal{H}_a}/2}_{-{\mathcal{H}_a}/2}c(x,y,t)\,\text{d}y. \end{equation}

We introduce variables averaged across the intrusion according to

(2.11) \begin{equation} \bar {f}(x,t) = \frac {1}{2h} \int ^{{\mathcal{H}_a}/2}_{-{\mathcal{H}_a}/2}c(x,y,t)f(x,y,t)\,\text{d}y. \end{equation}

In particular, the vertically averaged horizontal velocity within the intrusion is

(2.12) \begin{equation} \bar u(x,t)=\frac 1{2h}\int ^{{\mathcal{H}_a}/2}_{-{\mathcal{H}_a}/2}cu\,\text{d}y. \end{equation}

We also consider the vertical displacement $\eta (x,y,t)$ of an isopycnal passing through $(x,y)$ at time $t$ from its initial location. This is obtained from the density difference field as

(2.13) \begin{align} \eta (x,y,t)=\rho _{\textit{rel}}(x,y,t)+y. \end{align}

Wave features are evident on the thickness profiles. To analyse these, local maxima of $h$ were identified by first smoothing using a Savitzky–Golay filter to suppress any spurious peaks and then comparing each grid point value $h_i$ against its immediate neighbours $h_{i-1}$ and $h_{i+1}$ .

Finally, the tip location $x_f$ at time $t$ was found by identifying the largest grid point, $\max (x_i)$ , such that $h_i \gt \varepsilon$ , with $\varepsilon = 0.001$ . Reducing $\varepsilon$ did not perceptibly change the tip location identified.

3.1. A typical intrusion

Figures 3–5 show the intrusion thickness, velocity profiles and density profiles respectively for a simulation with ${\textit{Re}}=1000$ , $\mathcal{W}=3$ and ${\mathcal{H}_a}=10$ . Corresponding movies are provided with the Supplementary Material. Figure 6 shows the evolution in time of the density profile along a horizontal slice immediately above the intrusion as well as the position of the intrusion tip and the positions of wave crests on the upper surface of the intrusion in time. The parameters for this simulation approximately match those of experiment 8 of Ceddia et al. (Reference Ceddia, Kerr, Rudman, Settle and Slim2026), with one dimensionless spatial unit corresponding to $4.5\,\text{cm}$ and one dimensionless temporal unit corresponding to $1.0\,\text{s}$ .

Figure 3.

Intrusion profiles for ${\textit{Re}}=1000$ , $\mathcal{W}=3$ , ${\mathcal{H}_a}=10$ and $\mathcal{L}_a=150$ . (a–c) Colour function $c$ for the intruded fluid at the times indicated. A corresponding movie is provided in the Supplementary Material (movie 1). In (a), the (equal aspect ratio) inset shows a detailed view of the highlighted region. (d) Front location against time (solid black curve) and fitted line $x_f = 0.72 t$ (dashed red line). (e) Half-thicknesses at increasing times from $t=10$ to $t=100$ in intervals of $5$ .

Figure 4.

Velocity profiles in $y\geqslant 0$ for the simulation in figure 3. (a) Streamlines (solid black curves) and horizontal velocity (background heatmap) at $t=50$ . The intrusion envelope is shown in bold and extends to $x\approx 38$ . A corresponding movie is provided in the Supplementary Material (movie 2). (b–d) Horizontal velocity profiles at $t=50$ in different locations: (b) near the inlet, (c) elsewhere in the intrusion and (d) ahead of the intrusion. The intrusion envelope is again shown as a bold curve. (e) Horizontal velocity on $y=0$ and (f) vertically averaged horizontal velocity at increasing times from $t=4$ to $t=100$ in intervals of $2$ . The dotted/dashed red lines indicate the speed of mode-4/mode-6 internal waves. In (a–d), the profiles in $y\lt 0$ are the symmetric extensions.

Figure 5.

Density profiles in $y\geqslant 0$ for the simulation in figure 3. (a–c) Contours of the vertical displacement of isopycnals $\eta$ at the times indicated. A corresponding movie is provided in the Supplementary Material (movie 3). (d) Vertical density gradient $\partial \rho _{\textit{rel}}/\partial y$ at time $t=100$ . (e–g) Densities relative to the intrusion on the vertical line segments indicated in (c) at times from $t=5$ to $t=100$ in intervals of $5$ . Dotted and solid coloured curves indicate the displacement before and after the intrusion arrive respectively. The dashed black lines are the unperturbed density profiles $\rho _{\textit{rel}} = -y$ .

Figure 6.

Wave propagation for the simulation in figure 3. The heatmap shows the isopycnal displacement $\eta$ at $y = 1.2$ (indicated by a red line in figure 5 c). White points mark the location of all local crests on the intrusion (see § 2.3 for the calculation). The location of the advancing intrusion tip is indicated by the solid black curve. The solid red line and dashed red line represent the propagation speeds of mode-4 and mode-2 internal waves, respectively.

At a broad scale, the intrusion evolves into a simple, symmetric tapered-wedge shape, decreasing in $x$ from a half-thickness of approximately $1.1$ to zero nearly linearly except for a blunted nose (figure 3 panels a–c, e). The tip propagates at a near-constant velocity of approximately $0.7$ (figure 3 d).

At a finer scale, the half-thickness at the source rapidly adjusts from the inlet width of $\mathcal{W}/2$ to this profile through a prominent Kelvin–Helmholtz billow in $x\lesssim 2$ (figure 3 a–c). On the interface immediately beyond, there are a number of progressively smaller billows. Entrainment nevertheless is limited and the interface remains remarkably sharp (see Appendix C for further, quantitative details on entrainment). Beyond this region of billows, up to approximately the mid-length of the intrusion, the interface appears nearly perfectly linear (figures 3 a–c, e and 6). In the front half of the intrusion, a prominent wave train appears. The waves propagate at a constant speed (figure 6). They are faster than the intrusion tip and, as they reach the tip, they squeeze it and cause an oscillation in its speed (figures 3 d and 6). This interaction with the tip also appears to generate reflected waves. These waves propagate forwards but are markedly slower than the tip and weakly interfere with the incoming wave train (visible as gaps in the crests in figure 6).

Within the intrusion, the velocity is predominantly horizontal, with a minimal vertical component (figure 4 a). The flow rapidly adjusts from the prescribed parabolic inflow profile to a plug-like profile, although in the immediate vicinity of the inlet there is some inversion of the flow profile such that maximum speeds occur near the interface of the intrusion with the ambient (figure 4 b). Near the intrusion tip, a noticeable shear reappears (figure 4 c). As the intrusion thins in $x$ , the flow speeds up nearly linearly, such that the horizontal speed on $y=0$ at a given $x$ -location increases from approximately $0.5$ near the inlet to approximately $0.7$ at the tip at any given instant of time (figure 4 e).

In the ambient fluid, internal waves are clearly apparent (figures 4 and 5). Ahead of the intrusion, these are mode-2 with isopycnals raised in $y\gt 0$ (figure 5 b; note that they have the appearance of mode-1 as only the upper half of the domain is shown). These internal waves propagate with the speed of mode-2 columnar disturbances ${\mathcal{H}_a}/(2\pi )$ (dashed red line in figure 6). Near the tip of the intrusion, mode-4 waves dominate (figure 5 b, c). The isopycnals are raised immediately above the intrusion and lowered near the upper boundary of the domain. Behind this, back to approximately the mid-length of the intrusion, there continues to be internal wave-like motion in the ambient, but the dominant feature becomes a significant uplift of the isopycnals immediately above the intrusion (figure 5 b, c). The internal waves propagate in sync with the waves on the surface of the intrusion described above (solid red line in figure 6). In the near-source half of the intrusion, internal waves are not immediately apparent (neither are surface waves) but there remains a significant, near-constant uplift of the isopycnals reminiscent of evanescent waves (figure 5 b, c). The net effect of this displacement of the isopycnals is that the magnitude of the density gradient is increased by approximately 50 % from $\partial \rho _{\textit{rel}}/\partial y = -1$ to $\partial \rho _{\textit{rel}}/\partial y \approx -1.5$ above the intrusion (figure 5 d–g). Immediately ahead of the intrusion, the density gradient is reduced by approximately 55 % to $-0.45$ (figure 5 d–g).

The internal waves are long and thus the velocity in the ambient is also predominantly horizontal in much of the domain. This horizontal velocity manifests as a sequence of shear layers or alternating jets in $y$ (figure 4 a–d). Ahead of the intrusion there is a single forward and backward propagating pair. Above the intrusion tip, there are two forward propagating regions and an intervening reversed region. Nearer the source there are multiple forward and reversed regions. This behaviour appears qualitatively similar to that described by Manins (Reference Manins1976) and Ceddia et al. (Reference Ceddia, Kerr, Rudman, Settle and Slim2026). Far ahead of the intrusion, the region of forward propagating fluid gradually broadens with $x$ (figure 4 a) and the speed decays commensurately (figure 4 e) until the full layer is moving forward with the uniform velocity $1/{\mathcal{H}_a}$ required to balance the influx at the source.

3.2. Self-similar collapse

Motivated by the observations of Ceddia et al. (Reference Ceddia, Kerr, Rudman, Settle and Slim2026) and the near-constant velocity of the front (figure 3 d), we plot the half-thickness $h$ and vertically averaged horizontal velocity $\bar {u}$ in the intrusion as well as the density profile in the ambient against the similarity variable $\xi = x/t = \hat {x}/\displaystyle \sqrt {NQ}\,\hat {t}$ in figure 7. At a broad scale, the data collapse extremely well. The surface oscillations described above are also clearly evident and appear as secondary, non-self-similar features.

Figure 7.

Thickness, velocity and density profiles in self-similar coordinates. Self-similar collapse of (a) the half-thickness of the intrusion and (b) the vertically averaged horizontal velocity in the intrusion as a function of the similarity variable $\xi = x/t$ . (c–e) Contours of the isopycnal displacement $\eta$ at the times indicated. (f–g) Densities relative to the intrusion on the vertical line segments indicated in (e). Dotted and solid coloured curves indicate the displacement before and after the intrusion arrives, respectively. The times shown range from $t=5$ to $t=100$ in intervals of $1$ for (a) and (f–h). Parameters are as for figure 3.

3.3. Varying parameters

The broad-scale features described in §§ 3.1 and 3.2 do not change as parameters are varied, however, some of the details can change. We now explore the behaviour with varying Reynolds numbers ${\textit{Re}}$ , inlet widths $\mathcal{W}$ and ambient layer thicknesses $\mathcal{H}_a$ .

3.3.1. Varying Reynolds numbers

Figure 8 shows a comparison of the tip locations and intrusion profiles for various Reynolds numbers but otherwise identical parameters to those given above. As noted, the broad-scale features of the flow remain the same. For Reynolds numbers smaller than our baseline case, the intrusions are slightly stouter; for ${\textit{Re}}=200$ they are approximately $10$ % to $15\,\%$ thicker and commensurately shorter. For Reynolds numbers larger than our baseline case, the overall thickness profiles are nearly indistinguishable (figure 8 b). The most significant change is more intense and extensive Kelvin–Helmholtz billows (figure 8 e, g). These billows extend into the front half of the intrusion for larger ${\textit{Re}}$ , forming on top of the surface wave trains (figure 8 g) and generating a thin mixed layer.

For larger ${\textit{Re}}$ , the internal waves in the ambient do not catch the front as quickly and so tip squeezing and the associated pulsing propagation occur slightly later (figure 8 a). The overall structure of the density perturbations in the ambient remains the same, but there is more fine-scale motion, particularly nearer the source, as a result of Kelvin–Helmholtz forcing. Note that the apparent isopycnal displacement along the top boundary of the domain for ${\textit{Re}}=200$ in figure 8(d) results from the initial linear density profile adjusting to one satisfying zero flux. The diffusivity of the density is inversely proportional to the Reynolds numbers and so this effect is more pronounced at lower ${\textit{Re}}$ .

Figure 8.

Varying Reynolds numbers. Comparison of (a) the tip position against time and (b) an instantaneous profile at time $t=50$ for Reynolds numbers ${\textit{Re}}=200$ , $500$ , $1000$ (our baseline case), $2000$ and $5000$ . (c, e, g) Intrusion profiles and (d, f, h) isopycnal displacements at $t=100$ for Reynolds numbers of (c, d) $200$ , (e, f) $2000$ and (g, h) $5000$ (corresponding plots for ${\textit{Re}}=1000$ are figures 3 c and 5 c). Other parameters as in figure 3.

3.3.2. Varying inlet widths $\mathcal{W}$

Figure 9 shows a comparison of the tip locations and intrusion profiles for various inlet widths $\mathcal{W}\geqslant 2.2$ , but otherwise identical parameters to those in § 3.1. The extents and profiles are unchanged, except for a small portion of the profile near the source where larger inlet widths have an increasing undershoot. This results from the anomalously dense source fluid introduced at the top of the inlet falling, and building up momentum as it does. This vertical momentum is then transferred to horizontal as the fluid approaches neutral buoyancy. For larger $\mathcal{W}$ this generates more horizontal momentum and thus generates a faster, narrower near-source region. However, for the range of $\mathcal{W}$ we considered, this momentum was not sufficient to modify the profile away from the source region.

For $\mathcal{W}\lt 2.2$ , the profile becomes markedly different from those shown, with a narrow, steady, jet-like source region that thickens slightly away from the source before developing large-amplitude undulations reminiscent of an undular bore. For narrower inlets still, the motion is unstable to a sinuous instability and rapidly becomes turbulent. We do not consider intrusions with $\mathcal{W} \lt 2.2$ further in this study.

Figure 9.

Varying inlet widths. Comparison of (a) the tip position against time and (b) an instantaneous profile at time $t=50$ for inlet widths $\mathcal{W}= 2.2$ , $2.5$ , $3$ (our baseline case), $4$ and $5$ . (c) Intrusion profiles and (d) isopycnal displacements at $t=100$ for $\mathcal{W} = 5$ . Other parameters as in figure 3.

Figure 10.

Varying ambient layer thicknesses. (a, b) Instantaneous profiles at time (a) $t=50$ and (b) $t=100$ for ambient layer thicknesses ${\mathcal{H}_a}=3$ , $4.5$ , $6$ , $10$ and $20$ . (c) The tip location (bold solid curves) and the location of all local crests on the upper surface of the intrusion (dotted curves). Other parameters as in figure 3. The critical ambient layer thickness is ${\mathcal{H}_a}\approx 4.5$ (purple curves), with internal waves propagating ahead of the intrusion only for ambient layer thicknesses above this.

3.3.3. Varying ambient layer thicknesses $\mathcal{H}_a$

Given the prominent internal wave motion in the ambient, it might be expected that the intrusion profile changes appreciably with the thickness of the ambient, as this dictates the speed of columnar modes. However, this is not the case. Figure 10 shows a comparison of the tip locations and intrusion profiles for various ambient layer thicknesses, but otherwise identical parameters to those in § 3.1. The profiles are essentially unchanged and the overall extent of the intrusions are only slightly affected by the value of $\mathcal{H}_a$ .

Figure 11.

Behaviour in the ambient for varying ambient layer thicknesses. (a,c,e,g,i) Streamlines (solid black curves) and horizontal velocity (background heatmap) at $t=50$ . (b,d,f,h,j) Contours of the isopycnal displacement $\eta$ at time $t=50$ . (k–m) Densities relative to the intrusion on the vertical line segments indicated in (j) at time $t=50$ for varying ambient layer thicknesses. The dashed black lines are the unperturbed density profiles $\rho _{\textit{rel}} = -y$ . Results are shown for (a, b) ${\mathcal{H}_a}=3$ , (c, d) $4.5$ , (e, f) $6$ , (g, h) $10$ and (i, j) $20$ . Other parameters as in figure 3.

Figure 11 shows the velocity field and isopycnal displacements for the same set of ambient layer thicknesses. In contrast to the intrusion, the ambient dynamics changes markedly as $\mathcal{H}_a$ increases. For ${\mathcal{H}_a}=3$ , the speed of the intrusion tip is significantly higher than the speed of the fastest internal waves ( ${\mathcal{H}_a}/2\pi \approx 0.48$ ). In this case, there are no oscillatory modes, density perturbations cannot propagate ahead of the intrusion and the tip propagates at an apparently perfectly constant speed of $0.76$ . For ${\mathcal{H}_a}=4.5$ , the intrusion tip speed almost matches the fastest internal wave speed ${\mathcal{H}_a}/2\pi \approx 0.72$ . The frontal region has a thin elongated nose that propagates faster than in any of the other simulations, although this appears to be purely a result of lengthening of this nose rather than pointing to other structural changes of the intrusion. For ${\mathcal{H}_a}=6$ , mode-2 internal waves are clearly faster than the intrusion whereas mode 4 are slower and barely apparent. For ${\mathcal{H}_a}=10$ , mode-2 and mode-4 waves are faster than the intrusion tip, while for ${\mathcal{H}_a}=20$ , waves up to and including mode-8 all exceed the tip speed.

The wave train on the surface of the intrusion is absent for ${\mathcal{H}_a}=3$ , although the tip has the form of a single appreciable bolus. A wave train is apparent for larger ambient thicknesses, and across all $\mathcal{H}_a$ values that we considered the waves had comparable amplitude, wavelength and the crests moved at very similar speeds (see figure 10 c). This suggests that the speed matching the mode-4 internal waves in our baseline case with ${\mathcal{H}_a}=10$ was coincidental.

4.1. The Mei model: an unperturbed ambient density

Up to this point, all of the assumptions appear justified. However, the system is not yet complete as $\rho _{\textit{rel}}(x,y,t)$ is not known. In Mei’s original shallow-water description, the ambient density is assumed unperturbed from its initial state, or equivalently isopycnals are assumed to simply slide to the right ahead of the advancing intrusion, so $\rho _{\textit{rel}}(x,y,t)=-y$ and $\eta (x,y,t)=0$ for all time. This is not consistent with the observed profiles (see figures 5 and 13 a). Nevertheless adopting this form, from (4.5) we find the pressure in the intrusion is

(4.10) \begin{equation} p_{\textit{int}}(x,t) = p_0 + \frac {1}{2} h^2. \end{equation}

This pressure profile is shown in figure 13(b), calculated from the half-thickness profile obtained from the simulations. Significantly, the pressure gradient is notably too steep (red dashed curve versus black solid). Nevertheless continuing to use this form, (4.8) becomes

(4.11) \begin{equation} \frac {\partial }{\partial t}(\bar {u}h) + \frac {\partial }{\partial x}\left (\bar {u}^2h\right ) + h^2\frac {\partial h}{\partial x} = 0. \end{equation}

Figure 13.

Model fits and predictions for various shallow-water descriptions (coloured curves) together with the simulated data (black/grey curves). Predictions using the Mei shallow-water model are depicted with red dashed curves, those with a piecewise-linear density with blue curves and the linear velocity fit with purple curves. (a) Density profiles together with the simulated data at various $x$ locations within the intrusion at $t=50$ . (b) Pressure profiles at $t=50$ using (4.5) and taking the half-thickness profile from simulations for $h$ in (4.6). (c) The half-thickness profile in similarity form and (d) the vertically averaged horizontal velocity profile in similarity form. Note that the velocity fit (purple curves) can only provide a density prediction for $y$ values at which the intrusion is present, hence the profiles are only available for $|y|\lt h_0$ . The simulated data are from the solution presented in § 3.1.

We now look for solutions of the governing (4.1) and (4.11) subject to the flux boundary condition $\bar {u}h=1/2$ at the source $x=0$ and the tip boundary condition $h=0$ at $x=x_f(t)$ . Based on the simulation results, the additional condition $h=\mathcal{W}/2$ at $x=0$ is not relevant for $\mathcal{W}\gtrsim 2.2$ . Indeed it can be shown that only one characteristic of the hyperbolic system of equations propagates away from the source if it is sufficiently wide and so only one piece of information can be imposed at $x=0$ (see also Gratton & Vigo Reference Gratton and Vigo1994; Ceddia et al. Reference Ceddia, Kerr, Rudman, Settle and Slim2026).

Making the similarity ansatz

(4.12) \begin{equation} h=\mathcal{H}(\xi ) \quad \text{and} \quad \bar {u}=\mathcal{U}(\xi ), \end{equation}

where the similarity variable $\xi = x/t$ , translates the pair of equations (4.1) and (4.11) into

(4.13) \begin{align} (\mathcal{U}-\xi ) \mathcal{H}' + \mathcal{H}\mathcal{U}' = 0 \quad \text{and} \quad (\mathcal{U}-\xi ) \mathcal{U}' + \mathcal{H}\mathcal{H}' = 0. \end{align}

These equations have two types of solutions: a family of trivial, constant solutions and a non-trivial, rarefaction-like solution when the two equations become linearly dependent. This latter solution is the only one that can satisfy the tip condition. Imposing the source-flux condition on this solution gives

(4.14) \begin{equation} u = \frac {1}{\sqrt {2}} + \frac {1}{2}\xi \quad \text{and} \quad h = \frac {1}{\sqrt {2}} - \frac {1}{2}\xi . \end{equation}

These half-thickness and velocity profiles are shown in red in figure 13(c, d), respectively, together with the simulated results. The thickness profile has a similar structural form to the simulations, with a slender, wedge-like shape. However, there is a significant quantitative disagreement, with the similarity solution longer than the simulations by a factor of nearly two and thinner by a factor of approximately $0.7$ . This results directly from the pressure gradient being too steep, generating a stronger flow and hence a longer, thinner intrusion.

Modifying the ambient density profile provides a better approximation to the pressure and the intrusion profiles. We now explore two different approaches: the first is a particularly simple approximation informed by the data while the second uses a solution of the DJL equation to approximate the isopycnal uplift.

4.2. A uniformly perturbed ambient density profile

The simulation data suggest that the density perturbations are nearly steady once established and only weakly dependent on $x$ except for an oscillatory contribution resulting from the internal waves (see the nearly horizontal contours of isopycnal uplift in figure 5(b, c) and the similar structure of the density profiles at different times in figure 5(d–f) once the intrusion has arrived at a given location). Thus we assume the density perturbations are independent of $x$ and $t$ , $\rho _{\textit{rel}}(x,y,t)=\rho _{\textit{rel}}(y)$ and $\eta (x,y,t)=\eta (y)$ . Then (4.9) reduces to

(4.15) \begin{equation} h\overline {\frac {\partial p_{\textit{rel}}}{\partial x}} = - \rho _{\textit{rel}}(h)h\frac {\partial h}{\partial x} = h^2 \frac {\partial h}{\partial x} - \eta (h)h\frac {\partial h}{\partial x}, \end{equation}

and the momentum relation (4.8) becomes

(4.16) \begin{equation} \frac {\partial }{\partial t}(\bar {u}h) + \frac {\partial }{\partial x}\left (\bar {u}^2h\right ) - \rho _{\textit{rel}}(h)h\frac {\partial h}{\partial x} = 0, \end{equation}

or equivalently

(4.17) \begin{equation} \frac {\partial }{\partial t}(\bar {u}h) + \frac {\partial }{\partial x}\left (\bar {u}^2h\right ) + h^2\frac {\partial h}{\partial x} - \eta (h)h\frac {\partial h}{\partial x} = 0 . \end{equation}

When compared with the corresponding Mei version, this last equation illustrates the decelerating role of isopycnal uplift.

Making the same similarity ansatz as in (4.12) translates (4.1) and (4.16) into

(4.18) \begin{equation} (\mathcal{U}-\xi ) \mathcal{H}' + \mathcal{H}\mathcal{U}' = 0 \quad \text{and} \quad (\mathcal{U}-\xi ) \mathcal{U}' - \rho _{\textit{rel}}(\mathcal{H})\mathcal{H}' = 0. \end{equation}

The same boundary conditions apply as above and the solution satisfying the tip condition is again the rarefaction-like one, although in general it must now be found numerically. This is achieved by eliminating $\mathcal{U}$ from one of the equations by imposing a linear dependency condition and solving the resulting ordinary differential equation using Matlab’s bvp4c solver.

We present two different data-inspired approaches for closing the model: using a piecewise-linear fit to the density data and using a linear fit to the velocity profile. The former is a useful direct approach. We present the latter because the velocity profile is particularly clean and simple and so we can infer the density perturbation that provides the best fit to the observed thickness and velocity profiles. Thus we can explore whether a uniform uplift approximation is able to give a good description with a reasonable density profile.

4.2.1. A piecewise-linear fit to the density profile

We begin by taking a piecewise-linear fit to the density profile presented in figure 5 such that

(4.19) \begin{align} \rho _{\textit{rel}}(y) = \begin{cases} - y & \text{for } y\gt h_2,\\ -r_2 y + (r_2-r_1)h_1 & \text{for } h_1 \lt y \lt h_2 = (r_2-r_1)h_1/(r_2-1), \\ -r_1y & \text{for } 0 \lt y \lt h_1. \end{cases} \end{align}

Here, an unperturbed region in $y \gt h_2$ overlies a region with an increased density gradient $r_2\gt 1$ that in turn overlies a region with a reduced density gradient $r_1\lt 1$ in $y\lt h_1$ . The second region corresponds to uplifted denser fluid, while the last corresponds to upstream blocking.

The pressure in the intrusion is now

(4.20) \begin{equation} p_{\textit{int}} = \begin{cases} \dfrac {1}{2}h_2(2h-h_2) + \dfrac {1}{2}r_2(h_2-h)^2 + p_0 & \text{for } h\gt h_1, \text{ }\\[8pt] \dfrac {1}{2}h_2(2h_1-h_2) + \dfrac {1}{2}r_2(h_2-h_1)^2 - \dfrac {1}{2}r_1\left (h_1^2-h^2\right ) + p_0 & \text{for } h \lt h_1. \end{cases} \end{equation}

A density fit and the resulting pressure, thickness and velocity profiles for this description are shown in figure 13. The piecewise-linear density is a reasonably good fit. The resulting pressure profile (from (4.5), where we have used the half-thickness $h$ obtained from the simulations in (4.6)) is significantly closer to the simulated one than that obtained assuming an unperturbed ambient density, although the near-source region $x\lesssim 15$ still has slightly too steep a gradient. The similarity solution thickness and velocity profiles are also significantly closer to the simulated profiles, although they are still longer and thinner. Here, we have set $r_1=1/6$ , $r_2=1.25$ and $h_1=0.6$ , giving $h_2=2.6$ . These values were chosen based on visual agreement with the density gradients observed in the simulation data (figures 5 d and 13 a). We explored optimising the values to find the best fit for the thickness profile. However, this resulted in a significant overestimate of the density and it does not appear possible to match both the density and thickness profiles closely.

We note that a useful alternative approximation for the density perturbation in $y\gt h_1$ is an exponential fit inspired by the evanescent wave-like appearance of the density profile. Such a fit provides comparable pressure, thickness and velocity predictions.

4.2.2. A linear fit to the velocity profile

Alternatively, we can invert the procedure above, fitting the simulated velocity profile and inferring the thickness, and density. This approach is used because the velocity profile is particularly simple and almost linear. The approximation $\bar {u} = u_0 + \xi /3$ with $u_0=0.46$ provides an excellent fit to the data (see the purple line in figure 13 d). Here, we have explicitly set $1/3$ as the slope of the line as it provides a particularly compact form in the subsequent analysis; the slope of a line of best fit gives a minimally different value. Solving the exact (4.1) then gives the self-similar thickness profile $h = h_0\sqrt {1 - 2 \xi /(3u_0)}$ , with $h_0=1/(2u_0)\approx 1.09$ from the source-flux condition. This is also in good agreement with the simulated profiles (see the purple curve in figure 13 c).

We then substitute these forms into (4.16) and solve for the density at the surface of the intrusion $\rho _{\textit{rel}}(h)$ . This yields the remarkably simple form

(4.21) \begin{align} \rho _{\textit{rel}}(h) = -\alpha h^4, \end{align}

where $\alpha = u_0^2/h_0^4=16u_0^6$ . Thus the ambient density in $h\lt y\lt h_0$ can be inferred as

(4.22) \begin{equation} \rho _{\textit{rel}}(y) = -\alpha y^4. \end{equation}

The model does not predict the density in the ambient for $y\gt h_0$ . Figure 13(a) shows this density profile; it has the correct trend although is noticeably higher than the simulated values.

The solutions we have presented in the last two subsections show that, although not a perfect match, a simple description for the ambient density profile can significantly improve the agreement between the simulations and predictions. We believe that the remaining discrepancies primarily result from approximating the isopycnal uplift as independent of $x$ and $t$ .

4.3. An approximation using the Dubreil-Jacotin–Long equation

An alternative, physically motivated description for the density perturbation is possible using a solution of the DJL equation for isopycnal displacements. We emphasise that the DJL equation is only strictly valid for steady flows in which internal waves do not propagate upstream. The wave assumption is valid for our flows if ${\mathcal{H}_a} \lesssim 4.5$ . The steadiness assumption is not valid for our flows, although the tip of the intrusion is propagating at constant velocity and the nearly steady ambient density profile is what we wish to describe and so the approach gives a useful approximation.

In the limit of long-and-thin hydrostatic flows, the DJL equation in dimensionless form is (see, for example, Lamb & Wan Reference Lamb and Wan1998)

(4.23) \begin{equation} \frac {{d}^2{\eta }}{\text{d}{y}^2} + \frac {1}{C^2}\eta = 0. \end{equation}

Here, $C=u_{\!f}-u_{{bg}}$ , where $u_{\!f}$ is the speed of the intrusion tip and $u_{{bg}}=1/{\mathcal{H}_a}$ . The latter is the uniform flow speed far ahead of the intrusion and is non-zero to balance the source flux. It is absent in standard derivations but may readily be included by performing the analysis in a reference frame moving with the background flow speed.

To impose boundary conditions on (4.23), we follow previous authors and allow slow variations in $x$ . We additionally allow slow variations in $t$ . The conditions are then $\eta =0$ on $y={\mathcal{H}_a}/2$ and $\eta =Y(x,t)$ on $y=h(x,t)$ , where $Y(x,t)$ is an as-yet-undetermined isopycnal displacement on the surface of the intrusion. This second condition differs from other work using the DJL equation for topography (Baines Reference Baines2022) and intrusions (Ungarish Reference Ungarish2006), where $\eta (x)=h(x)$ is imposed on $y=h(x)$ corresponding to $\rho _{\textit{rel}}=0$ on $y=h$ . This condition does not match expectations in our geometry, where the initial intrusion surface at $t=0$ contacts ambient isopycnals of different densities, and also does not agree with simulations.

Solutions are thus

(4.24) \begin{equation} \eta (x,y,t) = Y(x,t) \, \frac {\sin \left (\displaystyle \frac {{\mathcal{H}_a}/2-y}{C}\right )}{\sin \left (\displaystyle \frac {{\mathcal{H}_a}/2-h(x,t)}{C}\right )}\,. \end{equation}

The corresponding horizontal velocity profile in the ambient is (Lamb & Wan Reference Lamb and Wan1998)

(4.25) \begin{equation} u(x,y,t) = u_{{bg}} + C \, \frac {\text{d}{\eta }}{\text{d}{y}} \quad \text{for } y\gt h. \end{equation}

Integrating the horizontal velocity in $y$ , assuming a sharp interface and symmetry across $y=0$ , imposing overall conservation of mass and using the definition (2.12) for the vertically averaged horizontal velocity in the intrusion, then gives

(4.26) \begin{equation} \frac {1}{2} = \int _0^{{\mathcal{H}_a}/2} u(x,y',t)\,\text{d}{y'} = \bar {u}h + u_{{bg}} \left ({\mathcal{H}_a}/2-h\right ) -C\,Y ,\end{equation}

and hence, using the value for $u_{{bg}}$ given above,

(4.27) \begin{equation} Y(x,t) = h(x,t) \, \frac {\bar {u}-u_{{bg}}}{C}. \end{equation}

Now imposing $u(x_f,0,t)=u_{\!f}$ in (4.25) and using the definition for $C$ implies $\text{d} \eta /\text{d} y = 1$ at $y=h(x_f,t)=0$ . Differentiating (4.24), substituting for $Y$ and analysing the limit $h\to 0$ gives $C={\mathcal{H}_a}/(2\pi )$ and the tip speed

(4.28) \begin{equation} u_{\!f} =\frac {{\mathcal{H}_a}}{2\pi } + \frac {1}{{\mathcal{H}_a}}. \end{equation}

Substituting (4.27) and (4.28) into (4.24) gives

(4.29) \begin{equation} \eta (x,y,t) = 2\pi \frac {h}{{\mathcal{H}_a}}(\bar {u}-u_{{bg}})\frac {\sin \left (2\pi y/{\mathcal{H}_a}\right )}{\sin \left (2\pi h/{\mathcal{H}_a}\right )} \quad \text{for } y\gt h. \end{equation}

Finally, substituting this form into (4.7) gives the pressure in the intrusion

(4.30) \begin{equation} p_{\textit{int}}(x,t) = p_0 + \frac {1}{2}h^2 + h(\bar {u}-u_{{bg}})\cot \!\left ( \frac {\pi h}{{\mathcal{H}_a}}\right )\!. \end{equation}

This now gives a complete model together with (4.1), (4.8) and (4.9).

Figure 14.

Predictions for the shallow-water model with the DJL description for isopycnal uplift (coloured curves) together with the simulated data. (a) Vertical isopycnal displacements $\eta ^{*}$ at selected positions from $\xi =0.3$ (intrusion body) up to (if valid solutions exist) $\xi = 0.72$ (intrusion tip). To aid comparison across cases, profiles are shifted horizontally: $\eta ^{*} = \eta$ , $\eta + 0.5$ , $\eta + 1.0$ and $\eta + 1.5$ for ${\mathcal{H}_a} = 3$ , $4.5$ , $6$ and $10$ , respectively. Solid curves show $\eta (y\gt h)$ evaluated using the DJL model solutions, while dashed curves show the corresponding simulated $\eta$ . For each $\mathcal{H}_a$ , colour intensity (light to dark) encodes progression in $\xi$ , as indicated by the grey-scale colour bar. (b) Pressure comparisons between self-similar DJL solutions (solid coloured curves) and simulated results (dashed coloured curves) along $y=0$ at $t=50$ . (c) The half-thickness profile in similarity form and (d) the vertically averaged horizontal velocity profile in similarity form. In the last two panels, the simulated data are given by black/grey curves.

A similarity solution of the form (4.12) can again be found numerically using the approach described in § 4.2. A comparison of various expressions given above and the self-similar solutions for $h$ and $\bar {u}$ with simulation data is shown in figure 14. Overall, the agreement is remarkably good. For ${\mathcal{H}_a}=3$ , the isopycnal displacement and pressure profiles are in close agreement. For larger $\mathcal{H}_a$ , the maximum isopycnal displacement is correctly captured, although the profile shape is increasingly mismatched. As a result of this shape mismatch, the pressures have a constant offset, which increases for larger $\mathcal{H}_a$ , although the pressure gradients are more important and these remain in close agreement. For all $\mathcal{H}_a$ , the thickness and velocity profiles are close to the simulations, with those for ${\mathcal{H}_a}=6$ in almost perfect agreement. However, the fact that the solution varies with $\mathcal{H}_a$ is in contrast with the simulation results. We also note that the tip speed obtained in the similarity solution does not match the expression (4.28) used to derive $Y(x,t)$ . The discrepancies between the model and the simulated solutions presumably result from the solutions violating the steadiness assumption underlying the DJL model and that waves can propagate ahead of the intrusion for ${\mathcal{H}_a}\gtrsim 4.5$ .