2.1 Non-particulate plume model

We assume an isolated plume emerges from a buoyant source, with initial total buoyancy $B_{0}$ rising in an unconfined, stationary, stratified ambient fluid of buoyancy frequency $N$ . For simplicity, we focus on the case where fluxes of mass and momentum at the source are zero, so that we have a ‘pure’ plume rising from a point source. In many practical cases with a source of finite size, the results will be very similar, provided that one integrates back to a ‘virtual origin’ and the vertical distances are offset appropriately. This type of flow has been studied by many authors since the seminal paper of Morton et al. (Reference Morton, Taylor and Turner1956), with the problem closed by assuming that the horizontal entrainment velocity of ambient fluid into the plume is proportional to the local vertical plume speed, $W$ , their ratio being called the ‘entrainment coefficient’, $\unicode[STIX]{x1D6FC}$ (Turner Reference Turner1986). Models typically assume either a ‘top-hat’ formulation, where fluid properties are assumed to jump from the ambient values outside the local plume radius, $R$ , to constant local plume values inside, or Gaussian profiles, where velocity, density deficit and (in principle) particle concentration have an $\exp (-r^{2}/R^{2})$ dependence and $W$ and $g^{\prime }$ signify values on the plume centreline. In this section, in order to extract information from past literature, we allow for either scenario, but later sections will use exclusively top-hat profiles. Writing $Q(=\unicode[STIX]{x03C0}R^{2}W)$ for the volume flux, $M(=QW/\unicode[STIX]{x1D6FE})$ for the specific momentum flux and $B(=Qg^{\prime }/\unicode[STIX]{x1D6FE})$ for the specific buoyancy flux, with $\unicode[STIX]{x1D6FE}=1$ for top-hat and $\unicode[STIX]{x1D6FE}=2$ for Gaussian plumes, we have the classic equations:

(2.1) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}{\displaystyle \frac{\text{d}Q}{\text{d}Z}}=E,\\ {\displaystyle \frac{\text{d}M}{\text{d}Z}}=g^{\prime }A,\\ {\displaystyle \frac{\text{d}B}{\text{d}Z}}=-N^{2}Q,\end{array}\right\}\end{eqnarray}$$

where entrainment rate ( $E$ ), cross-sectional area ( $A$ ) and reduced gravity ( $g^{\prime }$ ) are given by

(2.2a-d ) $$\begin{eqnarray}E=2\unicode[STIX]{x03C0}R\cdot \unicode[STIX]{x1D6FC}W,\quad A=\unicode[STIX]{x03C0}R^{2},\quad g^{\prime }=\frac{\unicode[STIX]{x1D70C}_{a}-\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D70C}_{0}}g,\quad N=\left(-\frac{g}{\unicode[STIX]{x1D70C}_{0}}\frac{\text{d}\unicode[STIX]{x1D70C}_{a}}{\text{d}Z}\right)^{1/2}.\end{eqnarray}$$

Note that $g^{\prime }$ is the local reduced gravity, referenced to ambient density $\unicode[STIX]{x1D70C}_{a}$ at that height, and not the global reference density $\unicode[STIX]{x1D70C}_{0}$ . In this section we assume that $N$ is constant.

Reflecting differences in the literature, we non-dimensionalise in two stages: with and without entrainment coefficient $\unicode[STIX]{x1D6FC}$ in the scaling. First we use dimensional scales $B_{0}$ and $N$ to give

(2.3) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}{\displaystyle \frac{\text{d}\tilde{Q}}{\text{d}\tilde{Z}}}=2\unicode[STIX]{x1D6FC}\sqrt{\unicode[STIX]{x1D6FE}\unicode[STIX]{x03C0}\tilde{M}},\\ {\displaystyle \frac{\text{d}\tilde{M}}{\text{d}\tilde{Z}}}={\displaystyle \frac{\tilde{B}\tilde{Q}}{\tilde{M}}},\\ {\displaystyle \frac{\text{d}\tilde{B}}{\text{d}\tilde{Z}}}=-\tilde{Q},\end{array}\right\}\end{eqnarray}$$

where

(2.4a-d ) $$\begin{eqnarray}Q=Q_{s}\tilde{Q},\quad M=M_{s}\tilde{M},\quad B=B_{0}\tilde{B},\quad \boldsymbol{X}=L_{s}\tilde{\boldsymbol{X}},\end{eqnarray}$$

(2.5a-c ) $$\begin{eqnarray}Q_{s}=(B_{0}^{3}/N^{5})^{1/4},\quad M_{s}=B_{0}/N,\quad L_{s}=(B_{0}/N^{3})^{1/4}.\end{eqnarray}$$

The non-dimensional entrainment constant $\unicode[STIX]{x1D6FC}$ appears in the first of (2.3) only. An $\unicode[STIX]{x1D6FC}$ -independent problem may be derived by a second, $\unicode[STIX]{x1D6FC}$ -dependent, scaling as

(2.6a-d ) $$\begin{eqnarray}Q=q_{s}q,\quad M=m_{s}m,\quad B=B_{0}b,\quad Z=l_{s}z,\end{eqnarray}$$

(2.7a-c ) $$\begin{eqnarray}q_{s}=\sqrt{\unicode[STIX]{x1D6FC}}Q_{s},\quad m_{s}=M_{s},\quad l_{s}=\frac{1}{\sqrt{\unicode[STIX]{x1D6FC}}}L_{s}=\unicode[STIX]{x1D6FC}^{-1/2}(B_{0}/N^{3})^{1/4},\end{eqnarray}$$

to give:

(2.8) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}{\displaystyle \frac{\text{d}q}{\text{d}z}}=2\sqrt{\unicode[STIX]{x1D6FE}\unicode[STIX]{x03C0}m},\\ {\displaystyle \frac{\text{d}m}{\text{d}z}}={\displaystyle \frac{bq}{m}},\\ {\displaystyle \frac{\text{d}b}{\text{d}z}}=-q.\end{array}\right\}\end{eqnarray}$$

Note, however, that for this second non-dimensionalisation the horizontal scale is $\unicode[STIX]{x1D6FC}$ times the vertical scale, since absolute radius $R$ and height $Z$ are given by

(2.9a,b ) $$\begin{eqnarray}R=\sqrt{\unicode[STIX]{x1D6FC}}L_{s}\frac{q}{\sqrt{\unicode[STIX]{x1D6FE}\unicode[STIX]{x03C0}m}},\quad Z=\frac{1}{\sqrt{\unicode[STIX]{x1D6FC}}}L_{s}z.\end{eqnarray}$$

Note also that entrainment coefficients reflect the choice of reference velocity for the plume profile. From (2.3),

(2.10) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}_{top\text{-}hat}=\sqrt{2}\times \unicode[STIX]{x1D6FC}_{Gaussian}.\end{eqnarray}$$

Since the turbulent entrainment mechanism is different, $\unicode[STIX]{x1D6FC}$ also depends on the type of injection. Typical values of $\unicode[STIX]{x1D6FC}$ for a Gaussian plume are 0.054 (jets) and 0.083 (plumes) – see Turner (Reference Turner1986). The corresponding values for top-hat profiles are 0.076 and 0.12, respectively.

Figure 2 shows the result of integrating (2.8) (with top-hat profiles: $\unicode[STIX]{x1D6FE}=1$ ) from an origin where the volume and momentum fluxes are zero. The plume has the same density as the ambient fluid at a height of $1.038l_{s}$ , although the plume will overshoot this, as velocity does not become zero until $1.366l_{s}$ . Corresponding values for a Gaussian plume, $\unicode[STIX]{x1D6FE}=2.0$ , are $0.87l_{s}$ and $1.15l_{s}$ , as in Morton et al. (Reference Morton, Taylor and Turner1956) – the different (by a factor of $2^{1/4}$ ) rise values in non-dimensional variables reflecting the different $\unicode[STIX]{x1D6FC}$ values for the different assumed plume profiles and reference velocities. The volume fluxes at the start and end of the overshoot region are $2.49q_{s}$ and $3.52q_{s}$ (for top-hat plumes) whilst the momentum fluxes at corresponding positions are precisely $B_{0}/N$ and 0 and the buoyancy fluxes are 0 and $-B_{0}$ for both types of plume profile (see § 2.2 below).

Figure 2. Top-hat integral model for a pure plume; (a) volume, momentum and buoyancy fluxes; (b) radius, velocity and reduced gravity. Variables scaled according to (2.6).

The plume model is not strictly applicable in the overshoot region, as the plume fluid falls back around the rising plume as a ‘curtain’, raising the density of fluid entrained into the rising column above that of ambient fluid. However, for pure plumes, the effect appears to be small. The laboratory experiments reported in Morton et al. (Reference Morton, Taylor and Turner1956) suggest, when extrapolated to a pure plume, a maximum rise of $Z_{max}=3.79(B_{0}/N^{3})^{1/4}$ and an entrainment coefficient (assuming Gaussian-plume profiles) of 0.093. Translated to the notation of the present paper this implies a maximum height of $1.16l_{s}$ . This is extremely close to that cited above for the Gaussian-plume zero-momentum level, suggesting that, despite the deficiencies of the entrainment assumption in the overshoot region, the maximum height of rise, $Z_{max}$ , is close to that by direct integration: $1.15l_{s}$ for Gaussian, or $1.37l_{s}$ for top-hat, plume profiles, respectively.

The curtain of plume material, being negatively buoyant, falls back around the rising core, entraining some ambient fluid and feeding an intrusion layer where the density matches the surrounds. A model for this has been proposed by Carrazo, Kaminski & Tait (Reference Carrazo, Kaminski and Tait2010), but in the context of negatively buoyant turbulent fountains. However, the available experimental evidence suggests that, for pure plumes, the amount of ambient fluid entrained in the overshoot region is much smaller than that in the rise phase, probably because, unlike forced plumes, relative velocities in this region are smaller and the interface between plume and ambient fluid is stably stratified. For the height of the intrusion layer, $Z_{T}$ , Richards et al. (Reference Richards, Aubourg and Sutherland2014) cite $Z_{T}=2.7(B_{0}/N^{3})^{1/4}$ for a pure plume. Accepting the maximum rise height $Z_{max}$ in the previous paragraph, this implies $Z_{T}/Z_{max}=0.71$ . Mirajkar et al. (Reference Mirajkar, Tirodkar and Balasubramanian2015), with forced plumes, cite similar ratios: 0.68 for single phase, or 0.75 with particle loading, the larger being primarily a consequence of particle re-entrainment depressing the maximum height of rise. The ratio of neutral-buoyancy height to maximum rise height for our naively integrated plume equations is 0.76, not substantially different from these ratios.

We conclude then that, despite its formal inapplicability in the overshoot (fountain) region, for pure plumes a simple integral plume model is consistent with an initial maximum rise predicted from the zero-momentum level ( $W_{T}=0$ ) and an intrusion at the zero-buoyancy level ( $B_{T}=0$ ) in the plume. This is unlikely to be true for buoyant jets, since the enhanced momentum is expected to lead to more substantial entrainment in the overshoot region.

2.2 Invariants

Before addressing particle-laden plumes, we note two invariants in special cases.

In an unstratified environment, the last of (2.1) gives constant buoyancy flux:

(2.11) $$\begin{eqnarray}B=\text{const.}\end{eqnarray}$$

This is widely used, but not relevant to the stratified environment here.

In a uniformly stratified environment ( $N=\text{const.}$ ) the last two equations of (2.8) can be combined and rearranged to give, for the $\unicode[STIX]{x1D6FC}$ -scaled variables:

(2.12) $$\begin{eqnarray}\frac{\text{d}}{\text{d}z}(b^{2}+m^{2})=0.\end{eqnarray}$$

Hence, in $\unicode[STIX]{x1D6FC}$ -scaled variables,

(2.13) $$\begin{eqnarray}b^{2}+m^{2}=\text{const.}\end{eqnarray}$$

The corresponding result in dimensional variables is a second (and surprisingly, but as far as we can tell, not hitherto highlighted) invariant

(2.14) $$\begin{eqnarray}B^{2}+N^{2}M^{2}=\text{const.}\end{eqnarray}$$

As a corollary, for a pure plume ( $m_{0}=0$ , $b_{0}=1$ ), then at the neutral-buoyancy level ( $b=0$ ) one must have

(2.15) $$\begin{eqnarray}m=1,\end{eqnarray}$$

whilst at maximum rise ( $m=0$ ),

(2.16) $$\begin{eqnarray}b=-1.\end{eqnarray}$$

Thus, if the plume equations were to be integrated to the point of zero vertical velocity then the buoyancy here would be equal and opposite to that at the source (figure 2 a).

2.3 Intrusion layer

The concentration of particles within the spreading intrusion is assumed to be axisymmetric, and the local volumetric flux of particles through the bottom of the intrusion is assumed to be given by fall speed $W_{s}$ multiplied by particle concentration (figure 1). Assuming, first of all, that there is negligible entrainment of ambient fluid in the overshoot region or the spreading intrusion, the volumetric concentration of particles in the intrusion, $C(r)$ , can be found (Sparks et al. Reference Sparks, Carey and Sigurdsson1991):

(2.17) $$\begin{eqnarray}\frac{\text{d}(Q_{T}C)}{\text{d}r}=-2\unicode[STIX]{x03C0}rCW_{s},\end{eqnarray}$$

where $r$ is radius and $Q_{T}$ is the volumetric flux in the plume at the intrusion height $Z_{T}$ (hence also the radial volume flux in the intrusion). The solution matching the plume radius $R_{T}$ at the intrusion height is

(2.18) $$\begin{eqnarray}C=C_{T}\text{e}^{-(\unicode[STIX]{x03C0}W_{s}/Q_{T})(r^{2}-R_{T}^{2})}.\end{eqnarray}$$

This profile has been used in the past for intrusions forming a free-surface gravity current, as in laboratory studies of plumes injected into unstratified fluid (Sparks et al. Reference Sparks, Carey and Sigurdsson1991; Veitch & Woods Reference Veitch and Woods2000). However, it assumes that the volume flux entering the intrusion layer is the same as that in the plume just below this layer. For the present, stratification determined, intrusions we can allow for additional particle free, and hence less dense, fluid entrained in the overshoot region. If we suppose an additional volume flux $\unicode[STIX]{x1D700}Q_{T}$ to be entrained here, then the total volume flux entering the intrusion layer is $(1+\unicode[STIX]{x1D700})Q_{T}$ and the initial concentration is $C_{T}/(1+\unicode[STIX]{x1D700})$ . Equation (2.18) is then amended to read:

(2.19) $$\begin{eqnarray}C=\frac{C_{T}}{1+\unicode[STIX]{x1D700}}\text{e}^{-(\unicode[STIX]{x03C0}W_{s}/(1+\unicode[STIX]{x1D700})Q_{T})(r^{2}-R_{T}^{2})}.\end{eqnarray}$$

The fraction $\unicode[STIX]{x1D700}$ links the buoyancy at the top of the plume, $B_{T}$ , and the distance between the maximum height in the fountain and the intrusion level, both of which can, in principle, be parameterised experimentally. If we estimate the averaged entrained density as that in the ambient fluid half-way between these levels and equate the resulting density of plume fluid entering the intrusion to the ambient fluid density at intrusion height then we have, from mass flux divided by volume flux entering the intrusion:

(2.20) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D70C}_{T}Q_{T}+\left(\unicode[STIX]{x1D70C}_{a}(Z_{T})+{\displaystyle \frac{1}{2}}(Z_{max}-Z_{T}){\displaystyle \frac{\text{d}\unicode[STIX]{x1D70C}_{a}}{\text{d}Z}}\right)\unicode[STIX]{x1D700}Q_{T}}{(1+\unicode[STIX]{x1D700})Q_{T}}=\unicode[STIX]{x1D70C}_{a}(Z_{T}),\end{eqnarray}$$

or, after rearrangement and rewriting in terms of buoyancy flux and buoyancy frequency:

(2.21) $$\begin{eqnarray}B_{T}=-{\textstyle \frac{1}{2}}(Z_{max}-Z_{T})(\unicode[STIX]{x1D700}Q_{T})N^{2}.\end{eqnarray}$$

Thus, the height of the intrusion layer $Z_{T}$ (with corresponding buoyancy $B_{T}=B(Z_{T})$ ) and maximum initial height, $Z_{max}$ , are intimately related to the additional volume flux entrained, $\unicode[STIX]{x1D700}Q_{T}$ . If $Z_{max}$ and $Z_{T}$ can be determined experimentally, then the ratio, $\unicode[STIX]{x1D700}$ , of ambient fluid entrained in the overshoot region to that below the intrusion can be obtained. Conversely, for any given $\unicode[STIX]{x1D700}$ , $Z_{T}$ can be determined numerically from (2.21) and the integral plume solution by noting that the quantity

(2.22) $$\begin{eqnarray}B+{\textstyle \frac{1}{2}}(Z_{max}-Z)(\unicode[STIX]{x1D700}Q)N^{2}\end{eqnarray}$$

varies monotonically from positive at the neutral-buoyancy level, to negative at $Z_{max}$ .

In the remainder of the paper we will retain $\unicode[STIX]{x1D700}$ as a free parameter in order to investigate the effect of the overshoot region. However, we also note the evidence presented at the end of § 2.1, which suggests that, for pure plumes at least, $B_{T}$ , and hence $\unicode[STIX]{x1D700}$ , is nearly 0. This is not the case for buoyant jets, because the experimental evidence cited earlier indicates that the additional momentum increases the initial rise, $Z_{max}$ , and is likely to promote greater mixing in the overshoot region. To highlight the effect that entrainment in the overshoot region could have, we plot graphs of non-dimensional intrusion height ( $Z_{T}/l_{s}$ ), intrusion-level plume buoyancy ( $B_{T}/B_{0}$ ) and critical source buoyancy ratio for particles ( $P_{crit}$ – see § 2.4 below) against $\unicode[STIX]{x1D700}$ in figure 3. The graphs show that entrainment in the overshoot region will, by entraining slightly lighter fluid, lead to a modest increase in intrusion height $Z_{T}$ and small negative values of $B_{T}$ . But, since the entrained fluid in this region is completely particle free, the effect on $P_{crit}$ is more substantial.

Figure 3. Effect of $\unicode[STIX]{x1D700}$ (ratio of volume flux entrained in the overshoot region to that at the top of the plume) on (a) non-dimensional intrusion height, $Z_{T}/l_{s}$ , with maximum rise $Z_{max}/l_{s}$ shown by the dashed line; (b) non-dimensional plume buoyancy, $B_{T}/B_{0}$ , at intrusion height; (c) critical source buoyancy ratio $P_{crit}$ (see § 2.4).

2.4 Collapse of particle-laden plumes

In cases where the plume source contains particles, these particles will affect the overall total buoyancy. We assume that the total buoyancy is given by $B=B_{F}+B_{P}$ , where $B_{F}$ is the buoyancy flux due to the plume fluid and $B_{P}$ is the buoyancy flux (assumed negative) due to the particles. We write

(2.23) $$\begin{eqnarray}P=|B_{0,P}/B_{0,F}|\end{eqnarray}$$

for the magnitude of the source buoyancy flux ratio (with $0\leqslant P<1$ ). In many previous studies particles are assumed to have negligible effect on the flow ( $P\ll 1$ ), but we make no such assumption here. In this section (but not for our more detailed model in § 3) we assume that the rising plume is sufficiently strong that particles remain within the plume until it begins to spread as an intrusion.

Particles fall at their settling velocity, $W_{s}$ , through the ambient fluid. As they fall they are drawn towards the rising plume and those that fall from the intrusion at less than a critical radius $R_{c}$ are re-entrained rather than falling to the floor. To calculate $R_{c}$ it is necessary to integrate the effect of the horizontal entrainment of the particles toward the plume centre line as the particles descend. Such an analytical approximation to the critical radius has been described by Sparks et al. (Reference Sparks, Carey and Sigurdsson1991), Lane-Serff (Reference Lane-Serff1995) and Mingotti & Woods (Reference Mingotti and Woods2015) for a linearly stratified ambient medium. By contrast, our formulation here, which was derived independently but turns out to be a rephrased version of that in Veitch & Woods (Reference Veitch and Woods2000), avoids making assumptions about the ambient stratification or the vertical profile of plume radius or entrainment velocity.

Figure 4. Re-entrainment region and particle fluxes.

The ‘entrainment boundary’ distinguishing falling particles which may be re-entrained into the plume from those which escape to be deposited on the ground (figure 4) may be computed from a particle velocity given in radial/vertical coordinates by

(2.24) $$\begin{eqnarray}\boldsymbol{U}=\left(\begin{array}{@{}c@{}}-{\displaystyle \frac{V_{e}R}{r}}\\ -W_{s}\end{array}\right),\end{eqnarray}$$

where $V_{e}$ is the entrainment velocity, $r$ is the local radial coordinate and $R$ is the plume radius at that height. Then

(2.25) $$\begin{eqnarray}\frac{\text{d}r}{\text{d}Z}=\frac{V_{e}R/r}{W_{s}},\end{eqnarray}$$

or

(2.26) $$\begin{eqnarray}2r\frac{\text{d}r}{\text{d}Z}=\frac{2V_{e}R}{W_{s}}=\frac{\text{entrainment rate}}{\unicode[STIX]{x03C0}W_{s}};\end{eqnarray}$$

i.e.

(2.27) $$\begin{eqnarray}\frac{\text{d}r^{2}}{\text{d}Z}=\frac{1}{\unicode[STIX]{x03C0}W_{s}}\frac{\text{d}Q}{\text{d}Z}.\end{eqnarray}$$

Hence, for a pure plume with $r=Q=0$ at source:

(2.28) $$\begin{eqnarray}r^{2}=\frac{Q}{\unicode[STIX]{x03C0}W_{s}}.\end{eqnarray}$$

The critical radius $R_{c}$ at intrusion height $Z_{T}$ is given by

(2.29) $$\begin{eqnarray}R_{c}^{2}=\frac{Q_{T}}{\unicode[STIX]{x03C0}W_{s}}=R_{T}^{2}\frac{W_{T}}{W_{s}},\end{eqnarray}$$

where $W_{T}$ is the upward velocity in the plume at this level. Note that this critical radius depends only on the flow field generated by the plume below the intrusion and not on fluid entrained in the overshoot (except, indirectly, in so far as that may affect the intrusion height). This equation is equally valid in free-surface-limited and stratification-limited cases.

In the intrusion layer, then, equation (2.19) can be written

(2.30) $$\begin{eqnarray}\frac{C}{C_{T}/(1+\unicode[STIX]{x1D700})}=\text{e}^{(W_{s}/W_{T})/(1+\unicode[STIX]{x1D700})}\text{e}^{-(\unicode[STIX]{x03C0}W_{s}r^{2}/(1+\unicode[STIX]{x1D700})Q_{T})}\end{eqnarray}$$

and the proportion of particulate re-entrained is, with F the advective flux of particulate:

(2.31) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}=1-\frac{F(R_{c})}{F_{T}}=1-\exp ^{(1/(1+\unicode[STIX]{x1D700}))(W_{s}/W_{T})}\exp \left(-\frac{1}{1+\unicode[STIX]{x1D700}}\right).\end{eqnarray}$$

In particular, if the particle settling velocity is small ( $W_{s}/W_{T}\ll 1$ , or $R_{T}/R_{c}\ll 1$ ) then the proportion re-entrained is

(2.32) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}=1-\exp \left(-\frac{1}{1+\unicode[STIX]{x1D700}}\right).\end{eqnarray}$$

If there is no entrainment from the fountain, $\unicode[STIX]{x1D700}=0$ , and this gives the proportion 0.632 derived by Veitch & Woods (Reference Veitch and Woods2000). A key feature is that this is independent of plume buoyancy, background stratification or particle settling velocity (provided the last is small).

Equation (2.32) plays a critical role in establishing conditions for plume collapse. We show below that the critical source buoyancy ratio, $P_{crit}$ is related to the residual buoyancy flux $B_{T}$ in the plume at the height of the intrusion, which, as shown in § 2.3, is related in turn to the amount of entrainment of particle-free fluid in the overshoot region.

Consider first the particle flux $F$ in plume or intrusion (figure 4). For heavy particles, $F$ is proportional to particle buoyancy flux $B_{P}$ . Following re-entrainment, a particle flux $F_{0}$ at source gives rise to a different particle flux $F_{T}$ at the top of the plume (i.e. intrusion level). Flux $\unicode[STIX]{x1D6F7}F_{T}$ is recycled into the plume by re-entrainment, whilst $(1-\unicode[STIX]{x1D6F7})F_{T}$ escapes. In steady state the latter balances the source flux, so

(2.33) $$\begin{eqnarray}F_{0}=(1-\unicode[STIX]{x1D6F7})F_{T},\end{eqnarray}$$

or, in terms of particle buoyancy,

(2.34) $$\begin{eqnarray}B_{T,P}=\frac{B_{0,P}}{1-\unicode[STIX]{x1D6F7}},\end{eqnarray}$$

where subscripts $T$ and 0 refer to top and source, subscripts $P$ and $F$ to particle and fluid.

Now consider the total plume buoyancy (fluid $+$ particles) at the intrusion level:

(2.35) $$\begin{eqnarray}B_{T}=B_{T,F}+B_{T,P};\end{eqnarray}$$

i.e.

(2.36) $$\begin{eqnarray}B_{T}=(B_{T,F}-B_{0,F})+B_{0,F}+\frac{B_{0,P}}{1-\unicode[STIX]{x1D6F7}}.\end{eqnarray}$$

At source, by definition of the source buoyancy ratio,

(2.37a,b ) $$\begin{eqnarray}B_{0,F}=\frac{B_{0}}{1-P},\quad B_{0,P}=-P\frac{B_{0}}{1-P}.\end{eqnarray}$$

Substituting (2.37) in (2.36) gives

(2.38) $$\begin{eqnarray}B_{T}=(B_{T,F}-B_{0,F})+\frac{B_{0}}{1-P}\left(1-\frac{P}{1-\unicode[STIX]{x1D6F7}}\right).\end{eqnarray}$$

The first bracketed term on the right-hand side is negative (since the plume entrains denser fluid), but becomes zero when the plume collapses (since then the ‘top’ $T$ is at the ground, 0). The critical value of $P$ is then given, after rearrangement, by

(2.39) $$\begin{eqnarray}P_{crit}=\frac{\unicode[STIX]{x1D6F7}^{\prime }(1-B_{T}/B_{0})}{1-\unicode[STIX]{x1D6F7}^{\prime }B_{T}/B_{0}},\end{eqnarray}$$

where $B_{T}\;({\leqslant}0)$ is the plume buoyancy at the height of the intrusion layer and $\unicode[STIX]{x1D6F7}^{\prime }=1-\unicode[STIX]{x1D6F7}$ is the non-entrained fraction. These are, in turn, related via the volume fraction $\unicode[STIX]{x1D700}$ entrained in the overshoot region, equations (2.21) and (2.32). $P_{crit}$ is plotted as a function of $\unicode[STIX]{x1D700}$ in figure 3(c).

The minimum height that we might conceivably assign to the spreading layer is the neutral-buoyancy height, $B_{T}=0$ , occurring when $\unicode[STIX]{x1D700}=0$ . Then we have a minimum value of $P_{crit}$ for plume collapse:

(2.40) $$\begin{eqnarray}P_{crit}=\unicode[STIX]{x1D6F7}^{\prime }|_{\unicode[STIX]{x1D700}=0}=0.368.\end{eqnarray}$$

We conclude that the plume will not collapse for any source buoyancy ratio smaller than this. Moreover, as explained at the end of § 2.3, this is likely to be the critical buoyancy ratio for pure plumes. We will show in § 3 that this bound is consistent with a more detailed computational model.

For buoyant jets, however, the situation is more complex, as the greater momentum throughout is expected to lead to both greater initial rise (Mirajkar et al. Reference Mirajkar, Tirodkar and Balasubramanian2015) and more mixing with particle-free fluid in the overshoot region. Equations (2.28)–(2.32) require correction for a finite source, whilst $\unicode[STIX]{x1D6F7}^{\prime }$ and $B_{T}$ must be modified. If the source is not large and the change to $Z_{max}$ is not substantial, then figure 3(c) shows the anticipated change in $P_{crit}$ with the entrained fraction from the overshoot region, $\unicode[STIX]{x1D700}$ . This paper will primarily focus on pure plumes, but will investigate the effect of overshoot-region entrainment in § 3.5.

2.5 Effect of particle sedimentation on the intrusion layer

In response to a reviewer’s comment, we address the effect on the intrusion layer of its particle load being continually diminished by sedimentation. As particles drop out, the overall density of the intrusion would be reduced, potentially leading to a layer that was not horizontal but increased in height with radius.

To our knowledge this has never been observed experimentally. Any buoyant tendency to rise would promote entrainment that quickly restored the neutral buoyancy of the layer, albeit thickening it. Even if there were a resultant rise then, in the region below the base of the intrusion, but above the original intrusion height, there would be no horizontal entrainment into the rising plume, and so particles would drop vertically until they reached the original intrusion height. Thus, apart from a small extra fall time, there is no significant effect on the re-entrainment process, and the model of § 2.4 (where we assumed a horizontal intrusion) will give the correct results for the recycled fraction, at least for a quasi-steady flow.

Even if no entrainment into the intrusion took place and it rose with radius, we can place an upper bound on any rise at the critical radius (the maximum distance at which it can have any effect on the re-entraining part of the flow). Any entrainment in the overshoot region inevitably dilutes the particle concentration. Hence, the maximum effect of particle dropout on intrusion density will occur if the intrusion height $Z_{T}$ is where the plume first reaches overall neutral buoyancy, $B_{T}=0$ , In this case we have $B_{T}^{(F)}=-B_{T}^{(P)}$ . At the critical radius, the particle buoyancy in the intrusion is the same as that at the source, $B_{0}^{(P)}$ . Thus, the buoyancy in the intrusion at the critical radius is (see (2.34) and (2.37)):

(2.41) $$\begin{eqnarray}B=B_{T,F}+B_{0,P}=-B_{T,P}+B_{0,P}=-\frac{\unicode[STIX]{x1D6F7}}{1-\unicode[STIX]{x1D6F7}}B_{0,P}=\frac{\unicode[STIX]{x1D6F7}}{1-\unicode[STIX]{x1D6F7}}\frac{P}{1-P}B_{0}.\end{eqnarray}$$

At $P=P_{crit}$ this is just $B_{0}$ . The density contrast is proportional to the buoyancy divided by the volume flow rate. The height the intrusion needs to rise to match the ambient density in this case is then (in $\unicode[STIX]{x1D6FC}$ -scaled units) $z_{rise}=1/q_{T}$ , where $q_{T}=2.49$ is the non-dimensional volume flux where the plume reaches neutral density, and thus $z_{rise}=0.402$ (compared with an initial non-dimensional intrusion height $z_{T}=1.038$ ). Initial particle loads below the critical value will lead to smaller changes in intrusion height.

We reiterate that we do not expect the rise to happen, and that, even if it did, it would have no effect on the criterion for plume collapse and a very minor effect on the time to collapse.

2.6 Time to collapse

The entrainment of particles into the plume reduces the overall buoyancy of the rising plume and depresses the intrusion height from its initial value. The intrusion height will tend to a steady non-zero value if $P$ is smaller than the critical value, while the intrusion height will reach zero in a finite time for values of $P$ larger than the critical value. We can estimate the change in intrusion height by estimating the overall change in buoyancy flux (since the heights scale with the quarter power of the buoyancy flux) while the time scale for the recycling process depends mostly on the fall time (which is much longer than the rise time in the plume or spreading time in the intrusion). The fall time is given by the intrusion height divided by the fall speed and diminishes as the intrusion height diminishes.

In a simple model we imagine successive discrete particle-recycling events, occurring at the time intervals over which a particle falls with velocity $W_{s}$ a distance equal to the varying intrusion height $Z_{T}$ . The particulate contribution to buoyancy at the top of the plume after $n$ recycling events is (cf. (2.34)):

(2.42) $$\begin{eqnarray}B_{T.P}^{(n)}=B_{0,P}(1+\unicode[STIX]{x1D6F7}+\unicode[STIX]{x1D6F7}^{2}+\cdots +\unicode[STIX]{x1D6F7}^{n})=B_{0,P}\left(\frac{1-\unicode[STIX]{x1D6F7}^{n+1}}{1-\unicode[STIX]{x1D6F7}}\right).\end{eqnarray}$$

Since the particle flux is added to the plume throughout its height, it is not straightforward to use the flux at source in an estimate of the changing intrusion height. However, the following is consistent with the requirements that: (i) the intrusion height $Z_{T}^{(n)}$ scales (in some sense) as (buoyancy) $^{1/4}$ ; and (ii) $Z_{T}^{(n)}\rightarrow 0$ as $n\rightarrow \infty$ for plumes that just collapse:

(2.43) $$\begin{eqnarray}Z_{T}^{(n)}=Z_{T}^{(0)}\left(\frac{B_{0,F}+B_{T,P}^{(n)}-B_{T}}{B_{0,F}+B_{0,P}-B_{T}}\right)^{1/4}.\end{eqnarray}$$

If we take the case of $B_{T}=0$ in (2.43) this gives intrusion heights that are the same as those predicted by injecting all the buoyancy flux from the entrained particles at the source, corresponding to the quickest collapse or lowest final intrusion height.

Using (2.37) for $B_{0,F}$ and $B_{0,P}$ , equation (2.39) for $B_{T}$ , and (2.42) for $B_{T.P}^{(n)}$ :

(2.44) $$\begin{eqnarray}Z_{T}^{(n)}=Z_{T}^{(0)}\left(\frac{\unicode[STIX]{x1D6F7}^{n}-\unicode[STIX]{x1D707}}{1-\unicode[STIX]{x1D707}}\right)^{1/4},\end{eqnarray}$$

where, with $P_{crit}$ the critical source buoyancy ratio for collapse given by (2.39),

(2.45) $$\begin{eqnarray}\unicode[STIX]{x1D707}=\frac{1-{\displaystyle \frac{P_{crit}}{P}}}{1-P_{crit}}.\end{eqnarray}$$

Consider first the marginal case where the plume just collapses: $P=P_{crit}$ or $\unicode[STIX]{x1D707}=0$ , and

(2.46) $$\begin{eqnarray}Z_{T}^{(n)}=Z_{T}^{(0)}\unicode[STIX]{x1D6F7}^{n/4}.\end{eqnarray}$$

This appears to approach $Z_{T}^{(n)}=0$ asymptotically. However, as the fall time is reduced at each recycling, complete collapse is actually obtained in finite time. Summing the first $n$ recyclings, with fall time during the $n$ th recycling, $\unicode[STIX]{x0394}t^{(n)}=Z_{T}^{(n-1)}/W_{s}$ :

(2.47) $$\begin{eqnarray}t^{(n)}=\frac{Z_{T}^{(0)}}{W_{s}}\left(\frac{1-\unicode[STIX]{x1D6F7}^{n/4}}{1-\unicode[STIX]{x1D6F7}^{1/4}}\right).\end{eqnarray}$$

Rearranging for $\unicode[STIX]{x1D6F7}^{n/4}$ and substituting in (2.46) gives for the intrusion height $Z_{T}$ :

(2.48) $$\begin{eqnarray}Z_{T}=Z_{T}^{(0)}-(1-\unicode[STIX]{x1D6F7}^{1/4})W_{s}t.\end{eqnarray}$$

With $\unicode[STIX]{x1D6F7}=0.632$ this gives an estimate of time to collapse for the marginal case as

(2.49) $$\begin{eqnarray}t_{collapse}=9.23\frac{Z_{T}^{(0)}}{W_{s}}.\end{eqnarray}$$

Equations (2.44) and (2.45) confirm that if $P<P_{crit}\;(\unicode[STIX]{x1D707}<0)$ there is no collapse, whereas if $P>P_{crit}\;(\unicode[STIX]{x1D707}>0)$ plume collapse is faster than the critical case.

In the non-collapsing case ( $\unicode[STIX]{x1D707}<0$ ) the model produces an estimate of the ratio of final to initial intrusion height:

(2.50) $$\begin{eqnarray}\frac{Z_{T}}{Z_{T}^{(0)}}=\left(\frac{|\unicode[STIX]{x1D707}|}{1+|\unicode[STIX]{x1D707}|}\right)^{1/4}.\end{eqnarray}$$

In the collapsing case, with $P>P_{crit}$ , the time to collapse can be estimated by summing (numerically) fall times $\unicode[STIX]{x0394}t^{(n)}=Z_{T}^{(n-1)}/W_{s}$ . This depends on $P_{crit}$ , which is determined from the volume fraction $\unicode[STIX]{x1D700}$ entrained in the overshoot region via the re-entrained particulate fraction $\unicode[STIX]{x1D6F7}$ and residual buoyancy ratio $B_{T}/B_{0}$ ((2.21), (2.32), (2.39); figure 3). Typical behaviour is illustrated in figure 5. Collapse time increases markedly near the critical value of $P$ , but does not exceed the limiting values from (2.48) for each value of  $\unicode[STIX]{x1D6F7}$ .

Figure 5. Time to collapse versus source buoyancy ratio $P$ for different degrees of entrainment in the overshoot region.

Figure 6. Schematic of PLUME and FIELD regions.

This idealised model makes the assumption that individual recycling events occur at discrete time intervals $Z_{T}^{(n)}/W_{s}$ , whereas, in reality, particles are re-entrained continuously in time over the full depth of the plume, the greatest particle re-entrainment occurring (at the shortest time interval) near the top of the plume. An advanced model seeks to provide a more accurate picture of events in § 3 below.

3.1 Overview

A computational model has been developed to supplement the simple analytical modelling in § 2. This advanced model offers several major advantages over the simple model of § 2; it enables the temporal evolution of plume and particle cloud to be followed, including the time to collapse; it incorporates time-dependent changes to plume density due to particulate re-entrainment; for heavier particles it can accommodate fallout from the plume during the rise phase. (Although not tested here, it would also admit an arbitrary ambient density profile.)

The model has two components (figure 6): a PLUME component consisting of an integral plume model and a horizontal intrusion, where plume rise is limited by ambient stratification; and a FIELD component where a time-dependent advection equation is solved by a two-dimensional, axisymmetric finite-volume approach for the concentration of particulate material falling out of the horizontal layer and being partially re-entrained into the plume.

These two components are coupled. The PLUME model supplies plume radius as a function of height, as well as the height of the spreading horizontal layer, the entrainment velocity at the plume boundary (which determines the horizontal velocity field outside) and the concentration distribution in the horizontal layer (2.19). The time-dependent FIELD model supplies the ambient concentration at the plume boundary, which is used to determine entrainment of particulate into the plume.

Assuming plume rise velocities to be much larger than particle settling velocity, it is sufficient that within each time step there be a single update of the quasi-static integral PLUME equations with current ambient conditions of density and concentration, followed by an update of the FIELD component.

3.2 PLUME model

In this section, the following integral plume equations are solved for volume flux $Q$ , specific momentum flux $M=WQ$ , mass flux ${\dot{m}}=\unicode[STIX]{x1D70C}Q$ and particle volume flux $F=CQ$ :

(3.1) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}{\displaystyle \frac{\text{d}Q}{\text{d}Z}}=E,\\ {\displaystyle \frac{\text{d}M}{\text{d}Z}}=g^{\prime }A,\\ {\displaystyle \frac{\text{d}{\dot{m}}}{\text{d}Z}}=\unicode[STIX]{x1D70C}_{a}E,\\ {\displaystyle \frac{\text{d}F}{\text{d}Z}}=C_{b}E^{\prime }.\end{array}\right\}\end{eqnarray}$$

Here, $E$ is the entrainment rate, given by

(3.2) $$\begin{eqnarray}E=2\unicode[STIX]{x03C0}R\cdot \unicode[STIX]{x1D6FC}W,\end{eqnarray}$$

and other quantities in the first two equations are as in § 2.

In the third equation we solve for the mass flux ${\dot{m}}=\unicode[STIX]{x1D70C}Q$ , rather than the buoyancy flux $B$ used in § 2. These are related by

(3.3) $$\begin{eqnarray}B=\frac{\unicode[STIX]{x1D70C}_{a}-\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D70C}_{0}}gQ=\frac{1}{\unicode[STIX]{x1D70C}_{0}}(\unicode[STIX]{x1D70C}_{a}gQ-\dot{m}g).\end{eqnarray}$$

We solve for mass flux rather than buoyancy flux in this section because it is computationally more tractable and permits the surrounding particulate concentration to change the density of entrained material, $\unicode[STIX]{x1D70C}_{a}$ .

In the last of (3.1) the net entrainment of heavy particles is, for consistency with the treatment at the boundary of the FIELD region, adjusted for the tendency to fall out of the plume:

(3.4) $$\begin{eqnarray}E^{\prime }=E-W_{s}\frac{\text{d}A}{\text{d}Z},\end{eqnarray}$$

and the boundary value $C_{b}$ is given by

(3.5) $$\begin{eqnarray}C_{b}=\left\{\begin{array}{@{}ll@{}}C_{a}\quad & (\text{ambient concentration})\text{ if }E^{\prime }>0,\\ C\quad & (\text{plume concentration})\text{ if }E^{\prime }<0.\end{array}\right.\end{eqnarray}$$

Finally, we note that, since $A=Q^{2}/M$ , it is computationally advantageous to avoid division by zero by solving the second of (3.1) for $M^{2}$ rather than $M$ :

(3.6) $$\begin{eqnarray}\frac{\text{d}M^{2}}{\text{d}Z}=2g^{\prime }Q^{2}.\end{eqnarray}$$

With these qualifications, the system of (3.1) is written in vector form as

(3.7) $$\begin{eqnarray}\frac{\text{d}\boldsymbol{Y}}{\text{d}Z}=\boldsymbol{F},\end{eqnarray}$$

and solved by the second-order implicit Euler method:

(3.8) $$\begin{eqnarray}\boldsymbol{Y}(Z+\unicode[STIX]{x0394}Z)=\boldsymbol{Y}(Z)+\frac{\unicode[STIX]{x0394}Z}{2}[\boldsymbol{F}(Z,\boldsymbol{Y}(Z))+\boldsymbol{F}(Z+\unicode[STIX]{x0394}Z,\boldsymbol{Y}(Z+\unicode[STIX]{x0394}Z))].\end{eqnarray}$$

The height of the horizontal intrusion is determined as in § 2.3, following (2.21). It is formally a function of the ratio, $\unicode[STIX]{x1D700}$ , of the volume of fluid entrained in the overshoot region to that below the intrusion; however, the majority of calculations to be reported here use $\unicode[STIX]{x1D700}=0$ (corresponding to $B_{T}=0$ ), with an examination of the effect of above-intrusion entrainment in § 3.5.

Although quasi-static, the plume equations have to be solved afresh at each time step as the ambient particulate concentration changes and affects the overall density of material entrained into the plume:

(3.9) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{a}=(1-C_{a})\unicode[STIX]{x1D70C}_{a,F}+C_{a}\unicode[STIX]{x1D70C}_{P},\end{eqnarray}$$

where $\unicode[STIX]{x1D70C}_{a,F}$ is the ambient fluid density, $\unicode[STIX]{x1D70C}_{P}$ is particle density and $C_{a}$ is the ambient concentration (by volume). The last is supplied at the plume boundary by the FIELD model.

For the integral plume calculation we used a step size $\unicode[STIX]{x0394}z/L_{s}=0.0005$ , which was found to be adequate for grid-independent solutions for an isolated plume. Integration was carried out to the height of zero rise velocity, beyond that necessary for the intrusion height.

3.3 FIELD model

In the FIELD model a time-dependent axisymmetric finite-volume approach, with curvilinear grid, is used to calculate concentration in the region external to the plume:

(3.10) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}(rU_{r}C)+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}(rU_{Z}C)=0.\end{eqnarray}$$

The mesh (shown schematically in figure 6) is delimited by the plume boundary, the intrusion layer (whose height may change with time), an outer boundary (fixed, by trial, at the greater of $10L_{s}$ and the critical radius) and the ground. Dimensions, determined by grid-independence checks in collapsing and non-collapsing cases, was $800\times 80$ (horizontal $\times$  vertical), with uniform distribution along each horizontal and near-vertical line.

The relative velocity field of particulate material is determined by the entrainment velocity at the plume boundary $R$ (but falling off as $1/r$ ) and the settling velocity of particles:

(3.11) $$\begin{eqnarray}\left(\begin{array}{@{}c@{}}U_{r}\\ U_{Z}\end{array}\right)=\left(\begin{array}{@{}c@{}}-V_{e}(z){\displaystyle \frac{R(z)}{r}}\\ -W_{s}\end{array}\right).\end{eqnarray}$$

The small downward velocity of cell faces as the intrusion layer descends could be included in a ‘moving-grid’ approach. However, this is considered sufficiently small here to be neglected. It does not affect the end state.

Dirichlet boundary conditions are supplied at inflow boundaries:

(3.12) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}C=C_{top}(r)\quad \text{at upper boundary (2.19)};\\ C=0\quad \text{at outer boundary.}\end{array}\right\}\end{eqnarray}$$

To maintain conservation and for consistency with the plume model (3.4) and (3.5), where there is net positive particle inflow $E^{\prime }$ to the plume zero-gradient boundary conditions are applied to the FIELD region at the plume boundary and supply the ‘ambient’ concentration $C_{a}$ ; conversely, where, due to significant fallout, there is net efflux from the plume this switches locally to Dirichlet conditions, with boundary concentrations those in the plume. At the ground, Neumann boundary condition conditions give the deposition flux.

Equation (3.10) is solved for cell-centre concentrations by a standard iterative approach (alternating direction implicit, with tridiagonal solver along coordinate lines), using first-order upwinding for advection and backward differencing for the time derivatives. With these schemes the system of algebraic equations is very diagonally dominant and converges within a few iterations at each time step. Typical computation times (up to a constant final-layer height or complete collapse) were a few minutes on a desktop PC, longest times being required for source buoyancy ratios near the critical value.

3.5 Effect of entrainment in the overshoot region

For stratification-limited plumes, the final intrusion height is affected by the amount of entrainment of particle-free fluid in the initial overshoot region, since this determines the density of fluid which forms the intrusion. In § 2.3 we formulated a relationship between the intrusion height, $Z_{T}$ , (and hence buoyancy $B_{T}$ ) and the ratio, $\unicode[STIX]{x1D700}$ , of volume entrained in the overshoot region to that in the sub-intrusion rise. We noted that (a) for pure plumes, experimental evidence is that $\unicode[STIX]{x1D700}$ is near zero, but that (b) critical source buoyancy $P_{crit}$ would vary fairly rapidly with $\unicode[STIX]{x1D700}$ , because the entrainment in the overshoot region is, by contrast with that below it, particle free. Thus, for any future investigation of buoyant jets, where, due to momentum, the entrainment in the overshoot region is likely to be more substantial, it is important to examine the potential effect of $\unicode[STIX]{x1D700}$ on collapse behaviour.

Using constant values $W_{s}/U_{s}=0.01$ and $\unicode[STIX]{x1D6FC}=0.12$ , figure 10 shows the variation of final height with source buoyancy ratio $P$ for various overshoot-entrainment ratios $\unicode[STIX]{x1D700}$ . Note that the time-varying maximum height $Z_{max}$ is required to compute the intrusion height (2.21), and this is still determined as the zero-momentum level in the integral plume model. In each case, the plume collapses abruptly as a critical value of $P$ is approached. These critical values are consistent with those deduced earlier, in figure 3(c). The graph confirms that entrainment in the overshoot region makes only a small relative difference to final intrusion height, but a more substantial increase to  $P_{crit}$ .

Figure 10. Variation of intrusion-layer height with source buoyancy ratio for different degrees of entrainment in the overshoot region.