2 An integral model for unsteady turbulent buoyant plumes

We model an unsteady turbulent buoyant plume formed due to the release of a fluid from a point source into an otherwise quiescent ambient fluid of a different density. A cylindrical coordinate system is adopted, with $\hat{\boldsymbol{r}}$ and $\hat{\boldsymbol{z}}$ denoting unit vectors in the radial and vertical directions respectively. The plume and ambient are composed of incompressible fluids, and the plume is assumed to have a circular (time-averaged) cross-section. The velocity field, $\boldsymbol{u}$ , is assumed to be axisymmetric, with $\boldsymbol{u}=u\hat{\boldsymbol{r}}+w\hat{\boldsymbol{z}}$ . We further assume that the plume is slender (as observed in experiments, e.g. Rouse et al. Reference Rouse, Yih and Humphreys1952), such that $R/H\ll 1$ , where $R$ and $H$ denote typical length scales in the radial and vertical directions respectively, that the Reynolds number of the emitted fluid at the source is sufficiently high, such that the turbulent motion is fully developed throughout the flow domain, and that the plume fluid transports a scalar species (such as heat or salt) with the plume density linearly related to the concentration of the species.

Turbulence in the plume is responsible for the entrainment of ambient fluid and mixing of the plume and ambient fluids, and its role is central to the ensuing dynamics. We therefore adopt the Reynolds-averaged Navier–Stokes equations to describe the fluid motion, with a Reynolds-averaged advection–diffusion equation to describe the transport of the scalar species that results in the density difference between the plume and the ambient fluids, taking $u=\overline{u}+u^{\prime }$ , $w=\overline{w}+w^{\prime }$ and $g_{r}=\overline{g_{r}}+g_{r}$ , where $\overline{u}$ and $u^{\prime }$ denote the ensemble average and the fluctuation about the average of $u$ respectively, and similarly for $w$ and $g_{r}$ , and where $g_{r}=g({\it\rho}_{a}-{\it\rho})/{\it\rho}_{0}$ is the reduced gravity, with $g$ denoting the gravitational acceleration, ${\it\rho}$ and ${\it\rho}_{a}$ denoting the density of the plume and ambient fluids respectively, and ${\it\rho}_{0}$ being a characteristic density scale. The equations for the conservation of (the ensemble averaged) mass, axial momentum and reduced gravity are then

It should be noted that in (2.1) the slenderness of the plume has been used to simplify the Reynolds-averaged Navier–Stokes equations, as described for steady plumes by Linden (Reference Linden, Batchelor, Moffatt and Worster2000). In particular, the radial momentum balance (not given) shows that, to leading order in the plume slenderness, there is no radial pressure gradient. Therefore, the pressure in the plume can be taken to be equal to the hydrostatic pressure of the ambient. We have also invoked the Boussinesq approximation whereby differences in density are sufficiently small that they may be neglected except where they are multiplied by the gravitational acceleration. Furthermore, we have neglected molecular diffusion terms under the assumption that the Reynolds number and Péclet number are sufficiently large that mixing is dominantly due to turbulence.

Integral equations are obtained by integrating each of equations (2.1) over a cross-section of the plume. We define a surface $r=b(z,t)$ representing the boundary of the plume over which entrainment of ambient fluid into the plume occurs, and impose the boundary condition on this surface,

where the entrainment velocity, $u_{e}$ , is the radial velocity of the ambient fluid across $r=b(z,t)$ (note that $u_{e}<0$ when ambient fluid is entrained into the plume). This boundary condition represents the advection of the surface with the local velocity of the plume, with entrainment represented as a sink of mass (equivalently volume under the Boussinesq approximation), with the volume flux per unit area of the surface across $r=b(z,t)$ given by $u_{e}$ . We note that this approach differs from that taken by Craske & van Reeuwijk (Reference Craske and van Reeuwijk2015a ,Reference Craske and van Reeuwijk b ), who adopt the characteristic length and velocity scales defined through the cross-sectionally averaged fluxes of mass and axial momentum which themselves are defined in terms of an additional length scale at which the mean axial velocity can be considered negligible in contrast to the mean axial velocity at the centreline. As turbulent plumes typically exhibit non-uniform radial profiles of axial velocity and concentration of species, with decaying velocity and concentration with the radial distance from the plume axis, in laboratory experiments or direct numerical simulations there is a choice as to how to define the plume edge, taking, for examples, the radial distance at which the axial velocity or concentration reaches a specified threshold, or the radial distance at which the mass flux included in a cross-sectional integral captures a specified proportion of the total mass flux induced by the plume. We discuss below how these choices are represented in our model.

Integration of the point-wise conservation equation (2.1) over a plume cross-section, using Leibniz’s theorem for interchanging differentiation and integration (Flanders Reference Flanders1973) together with the boundary condition (2.2), gives

representing the evolution of mass, momentum and buoyancy respectively.

We define the integral volume flux as $Q=b^{2}W$ (note that, under the Boussinesq assumption, $Q$ can also be referred to as the (specific) mass flux), the (specific) momentum flux as $M=b^{2}W^{2}$ and the buoyancy flux as $F=b^{2}WG^{\prime }$ . Here, $W$ and $G^{\prime }$ denote the cross-sectionally averaged mean axial velocity and the reduced gravity of the plume respectively, and are given by

We then have

where $S_{m}$ and $S_{f}$ are momentum and buoyancy ‘shape factors’ respectively, defined as

and

Thus, the shape factors quantify the effect of non-uniform radial profiles of the axial velocity and density on the rates of transport of momentum and scalar species. Crucially, we note that $S_{m}\geqslant 1$ , which is representative of the more rapid transport of momentum than mass in the plume, unless the mean axial velocity is radially invariant within the plume (in which case $S_{m}\equiv 1$ and $S_{f}\equiv 1$ ). It should be noted that this latter condition occurs only for ‘top-hat’ profiles as described below.

The plume edge, given by the surface $r=b(z,t)$ , is chosen to be at a radial distance such that the boundary terms in (2.3) are negligible in comparison to the integral terms (e.g. $\overline{w}(b)\ll W$ , $\overline{u^{\prime }w^{\prime }}(b)\ll W^{2}$ , etc.). Furthermore, the entrainment assumption of Morton et al. (Reference Morton, Taylor and Turner1956) allows the entrainment velocity to be written as $u_{e}=-{\it\alpha}W$ , where ${\it\alpha}$ is the entrainment coefficient. This simplification of the turbulence mixing dynamics assumes a similarity of turbulent structures at each height in the plume. The entrainment coefficient must be determined empirically and, for fully developed plumes far from the source where the radial profiles of the axial velocity have a self-similar Gaussian form, a constant value ${\it\alpha}=0.1\pm 0.01$ is appropriate (Morton et al. Reference Morton, Taylor and Turner1956; List Reference List1982; Woods Reference Woods2010). However, near to the source (or for non-ideal initial conditions), the entrainment coefficient may vary substantially as the flow develops (Wang & Law Reference Wang and Law2002; Kaminski et al. Reference Kaminski, Tait and Carazzo2005; Ezzamel et al. Reference Ezzamel, Salizzoni and Hunt2015). In this study we take a constant entrainment coefficient, with the exception of appendix A, where we also examine the steady solutions of a model with an entrainment coefficient that varies as the flow develops.

We note from (2.6) that the value of the momentum shape factor is tied to the choice of the plume width, $b(z,t)$ . If we assume that the mean axial velocity has a self-similar Gaussian profile in the radial direction, as inferred from laboratory experiments on steady plumes (Papanicolaou & List Reference Papanicolaou and List1988; Shabbir & George Reference Shabbir and George1994), so that the mean axial velocity of the plume can be written as

where $\overline{W}(z,t)$ is the mean axial velocity on the plume axis and $R(z,t)$ is a characteristic length scale for radial variation determined from observations, then the shape factor is given by

The plume width can be related to the length scale of the Gaussian radial profile $R(z,t)$ through a threshold on the mean axial velocity. For example, Morton et al. (Reference Morton, Taylor and Turner1956) and Papanicolaou & List (Reference Papanicolaou and List1988) define the characteristic length scale of a Gaussian plume as the radial position at which the axial velocity is a factor of $1/e$ of the centreline value (i.e. $b=R$ ), and this results in $S_{m}=1.08$ .

The entrainment hypothesis and the use of shape factors for the integral fluxes of momentum and buoyancy allow the integral conservation equation (2.3) to be written in terms of the integral fluxes of volume, momentum and buoyancy as

While the system (2.10) applies to a spatially and temporally varying ambient density field, for the remainder of this study we consider plumes in an unstratified ambient, with ${\it\rho}_{a}$ constant. The shape factors $S_{m}$ and $S_{f}$ could be spatially and temporally varying, as the radial profiles of the axial plume velocity and reduced gravity evolve (Carazzo et al. Reference Carazzo, Kaminski and Tait2006; Craske & van Reeuwijk Reference Craske and van Reeuwijk2015a ,Reference Craske and van Reeuwijk b ; Ezzamel et al. Reference Ezzamel, Salizzoni and Hunt2015). However, for fully developed plumes, laboratory experiments suggest that the axial plume velocity and reduced gravity attain self-similar forms (Morton et al. Reference Morton, Taylor and Turner1956; Papanicolaou & List Reference Papanicolaou and List1988; Shabbir & George Reference Shabbir and George1994; Ezzamel et al. Reference Ezzamel, Salizzoni and Hunt2015), such that the shape factors can be taken to be constants. While temporal changes in the source conditions could result in local variations in the shape factors, we expect the contribution of that change to be small. Therefore, here we consider the simplest case of $S_{m}=S$ with $S$ constant and $S_{f}\equiv 1$ to investigate the structure of the model, but note that we would expect that $S_{f}$ differs from unity. For example, taking a Gaussian profile for the reduced gravity, $\overline{g_{r}}(r,z,t)=\overline{G_{r}}(z,t)\text{e}^{-r^{2}/{\it\lambda}^{2}R^{2}}$ , where $\overline{G_{r}}(z,t)$ is the centreline value and ${\it\lambda}\approx 1.2$ (Papanicolaou & List Reference Papanicolaou and List1988), we find $S_{f}=1.06$ when taking $b=R$ . However, it will be shown below that it is the divergence of the momentum shape factor, $S_{m}$ , from unity that makes the significant difference to the model system (2.10), changing it from parabolic for $S_{m}=1$ to hyperbolic for $S_{m}>1$ , whereas the buoyancy shape factor plays a less significant role, only leading to loss of hyperbolicity if $S_{f}=S_{m}\pm \sqrt{S_{m}^{2}-S_{m}}$ (see § 4).

In the integral equations for conservation of momentum (2.10b ) and conservation of buoyancy (2.10c ) we have retained integral terms that represent turbulent axial diffusion of momentum and buoyancy respectively. Typically, in steady plume models these diffusive terms have been neglected (we refer to the system of integral equations (2.10) without the turbulent diffusive terms as the ‘non-diffusive’ system). Indeed, Morton (Reference Morton1971) argues that

and therefore the contributions of the turbulent diffusion terms to the plume dynamics are of a similar magnitude to boundary terms which are neglected. We note that some studies on momentum-driven jets from maintained sources retain the axial derivatives of quadratic fluctuation terms to model the contribution of a non-hydrostatic axial pressure gradient (see, e.g., Shabbir & George Reference Shabbir and George1994; Wang & Law Reference Wang and Law2002; Yannapoulos Reference Yannapoulos2006). When these turbulent diffusion terms are neglected, the momentum shape factor is taken to be constant ( $S_{m}=S$ with $S$ constant), the buoyancy shape factor $S_{f}\equiv 1$ and the ambient fluid is unstratified, the governing equations (2.10) can be written as

The steady integral equations with pure-plume boundary conditions ( $Q=0$ , $M=0$ , $F=F_{0}$ at $z=0$ ) have well-known power-law solutions (Morton et al. Reference Morton, Taylor and Turner1956), and these solutions are little altered by the momentum shape factor; for $S\geqslant 1$ the steady solutions of the non-diffusive system of equation (2.12) are

where

In the top-hat limit $S\rightarrow 1$ the solution of Morton et al. (Reference Morton, Taylor and Turner1956) is recovered. Furthermore, the effective radius of the plume, $b_{0}(z)=Q_{0}/\sqrt{M_{0}}=6{\it\alpha}z/5$ , is independent of the shape factor. Therefore, the momentum shape factor cannot be determined from measurement of the plume radius alone, in contrast to the entrainment coefficient which can be determined using the spreading rate of the steady plume.

Scase & Hewitt (Reference Scase and Hewitt2012) advocate modelling of the turbulent diffusion terms (particularly the diffusion of momentum) in (2.10) in the unsteady plume model to obtain a well-posed system of equations. Recently, Craske & van Reeuwijk (Reference Craske and van Reeuwijk2015a ,Reference Craske and van Reeuwijk b ) have demonstrated the relatively weak role of turbulent diffusion for momentum-driven jets, but that diffusive effects in the jet occur due to the departure of flow variables from self-similar forms (type II dispersion in the nomenclature of Craske & van Reeuwijk Reference Craske and van Reeuwijk2015a ,Reference Craske and van Reeuwijk b ). However, we suggest that the dominant dynamics for turbulent plumes can be described by the non-diffusive system of equations. Indeed, we show below that the inclusion of diffusive terms modelling turbulent diffusion can lead to difficulties in the integral model. We therefore analyse the non-diffusive system and show that a momentum shape factor that differs from unity is sufficient to obtain a well-posed model.

3 Difficulties associated with the turbulent diffusive terms

In an analysis of a non-diffusive unsteady plume model with shape factors $S\equiv 1$ and $S_{f}\equiv 1$ , corresponding to top-hat radial profiles for the axial velocity and reduced gravity (i.e. $\overline{w}$ and $\overline{g_{r}}$ are radially invariant for $r\leqslant b$ and equal to zero for $r>b$ ), Scase & Hewitt (Reference Scase and Hewitt2012) assess the well-posedness of the system of equations by introducing small harmonic perturbations to the source buoyancy flux, and examine the growth of the perturbations downstream of the source. A linear analysis shows that the perturbations grow with distance from the source, indicating instability, and, more importantly, the growth rate of the perturbations increases without bound as the frequency of the harmonic oscillation of the source buoyancy flux increases (Scase & Hewitt Reference Scase and Hewitt2012). The non-diffusive system of equations with $S\equiv 1$ and $S_{f}\equiv 1$ is therefore ill-posed, as there is a loss of continuous dependence of the solution on the boundary conditions (Joseph & Saut Reference Joseph and Saut1990; Scase & Hewitt Reference Scase and Hewitt2012). The ill-posedness is manifested in numerical solutions of the system of equations by an inability to compute solutions that are independent of the truncation implicit in the numerical scheme (e.g. grid-scale dependence for finite-difference methods).

The ill-posedness identified by Scase & Hewitt (Reference Scase and Hewitt2012) in the non-diffusive system when top-hat profiles are assumed (i.e. when $S\equiv 1$ ) may be due to a missing physical process that provides a mechanism to curtail the unbounded growth of the arbitrarily short-wavelength/high-frequency modes (Joseph & Saut Reference Joseph and Saut1990; Scase & Hewitt Reference Scase and Hewitt2012). In an attempt to regularize the ill-posed system, Scase & Hewitt (Reference Scase and Hewitt2012) introduce a phenomenological model for the diffusive term in the momentum balance (2.10b ), appealing to Prandtl’s mixing-length theory for turbulent eddy diffusion,

where ${\it\kappa}>0$ is a dimensionless parameter characterizing the diffusion of momentum through the action of turbulent eddies whose length scale is set by the radius of the plume. The equation expressing the balance of axial momentum proposed by Scase & Hewitt (Reference Scase and Hewitt2012) is then given by

Scase & Hewitt (Reference Scase and Hewitt2012) provide numerical evidence that the diffusive term leads to a well-posed system of equations. It should be noted that, while Scase & Hewitt (Reference Scase and Hewitt2012) take $S\equiv 1$ , in the analysis presented below we analyse the more general problem of $S\geqslant 1$ .

Scase & Hewitt (Reference Scase and Hewitt2012) show that steady power-law solutions for pure-plume boundary conditions exist for the diffusive system (with ${\it\kappa}>0$ anticipated to be small), and the structural change to the system of equations leads to a small modification of the solutions of Morton et al. (Reference Morton, Taylor and Turner1956) (given by (2.13) with $S=1$ ). The steady solutions of the ‘regularized’ system of equations with turbulent diffusion of momentum are given by (Scase & Hewitt Reference Scase and Hewitt2012)

However, we find that these power-law solutions are spatially unstable. Indeed, taking perturbations of the form

(noting that the buoyancy flux $F(z)\equiv F_{0}$ for a plume in an unstratified ambient), where ${\it\epsilon}>0$ is an ordering parameter, and linearizing the governing ordinary differential equations (for ${\it\epsilon}\ll 1$ ), we obtain

Seeking solutions of the form $Q_{1}^{(SH)}=q_{1}z^{{\it\sigma}}$ and $M_{1}^{(SH)}=m_{1}z^{{\it\sigma}}$ , we find that non-trivial solutions are possible if ${\it\sigma}={\it\sigma}_{1}=-1$ or

At least one of ${\it\sigma}_{\pm }>0$ for ${\it\kappa}>0$ for any $S\geqslant 1$ (figure 1 a) and therefore the steady solutions of the Scase & Hewitt (Reference Scase and Hewitt2012) ‘regularized’ model are potentially spatially unstable. The growing modes can be suppressed by the application of a boundary condition on a truncated domain (Scase & Hewitt Reference Scase and Hewitt2012), but we note that the system of equations for an unstratified ambient is formally applied on a semi-infinite domain without a far-field boundary condition. We note that, in the limit ${\it\kappa}=0$ , we find that ${\it\sigma}_{+}$ no longer appears in the analysis, while ${\it\sigma}_{-}\rightarrow -10/3$ , and therefore the steady solutions become spatially stable (as expected, since the ‘regularized’ system reduces to the classical Morton et al. (Reference Morton, Taylor and Turner1956) plume model). Spatial instability of the steady solutions is also found in a stratified ambient (not shown here), and in this case leads to severe difficulties as steady solutions cannot, in general, be found analytically, and the spatial instability precludes the numerical computation of steady solutions. Therefore, while the diffusive term introduced by Scase & Hewitt (Reference Scase and Hewitt2012) into the momentum balance resolves the ill-posedness in the unsteady plume model, the modification of the governing equations renders the steady solutions spatially unstable. While our present analysis considers pure-plume boundary conditions (i.e. $Q(0)=0$ , $M(0)=0$ and $F(0)=F_{0}$ ), in appendix B we show that spatially stable steady solutions with arbitrary source conditions are not possible.

Figure 1. Spatial growth rates $\text{Re}({\it\sigma})$ as a function of the diffusion coefficient for (a) the Scase & Hewitt (Reference Scase and Hewitt2012) system with momentum diffusion (3.2), (b) the alternative form of the momentum diffusion equation (3.8) and (c) diffusion of buoyancy (3.11). In (a,b) different values of the momentum shape factor are shown, with $S=1$ (solid line), $S=1.5$ (dashed line) and $S=2$ (dash-dot line). In (c) the three branches of the growth rates are distinguished, with one branch corresponding to a real mode (solid line) and two branches that form a complex conjugate pair (dashed and dash-dot lines).

An alternative (although similar) form of the diffusion term can be formed by re-ordering the cross-sectional integration following the mixing-length parameterization of the fluctuation vertical momentum flux (with boundary terms introduced that are assumed small and subsequently neglected), so that we write

in place of (3.1). We then obtain the following expression for conservation of momentum:

The corresponding steady solutions are

and a spatial stability analysis gives growth rates ${\it\sigma}={\it\sigma}_{1}=-1$ or

Thus, as with the Scase & Hewitt (Reference Scase and Hewitt2012) momentum diffusion term, there is a growing mode (figure 1 b), and the steady solutions are spatially unstable.

Finally, we consider the change to the classical steady solutions of Morton et al. (Reference Morton, Taylor and Turner1956) that occurs when a diffusive term is included in the equation for conservation of buoyancy (2.10c ) rather than in the equation for conservation of momentum (i.e. we now take ${\it\kappa}=0$ ). We again appeal to the Prandtl mixing-length theory to model the fluctuation of the flux of buoyancy, and obtain the following diffusive equation for the conservation of buoyancy:

where ${\it\kappa}_{1}$ is a dimensionless parameter that describes the strength of the diffusion of buoyancy. The classical steady solutions (2.13) for pure-plume boundary conditions remain as solutions of the system (and other possible solutions cannot satisfy the pure-plume boundary condition at $z=0$ ). However, a spatial stability analysis of the steady solutions (here including perturbations to the buoyancy flux in addition to the mass and momentum fluxes) shows that small-amplitude perturbations to the steady solutions are amplified unless the magnitude of the turbulence-induced diffusion is relatively large ( ${\it\kappa}_{1}>4/3$ , figure 1 c).

The consequences of these analyses are that, while axial diffusivity effects may capture some important unsteady processes in the mixing of momentum or buoyancy, the form of the parameterization significantly alters the steady states attainable by the system. Indeed, the states equivalent to those established by Morton et al. (Reference Morton, Taylor and Turner1956) have become spatially unstable, and therefore the unsteady model with axial diffusivity is not able to describe the steady states without an ad hoc far-field boundary condition.

4 Well-posedness of the hyperbolic unsteady plume model

As the turbulent diffusive terms, modelled phenomenologically using a Prandtl mixing-length approach, lead to a pathology in the integral plume model whereby the well-established steady states of Morton et al. (Reference Morton, Taylor and Turner1956) cannot be obtained, the regularization through eddy diffusion, while potentially physically appealing, is mathematically problematic. We therefore seek an alternative regularization of the system of unsteady integral equations; specifically, we examine the non-diffusive system of equations (i.e. we neglect the turbulent diffusive terms in (2.10), as suggested by Morton Reference Morton1971) but allow the momentum shape factor to differ from unity to describe the different rates of transport of mass and momentum that result from non-uniform radial profiles. Thus, from here on in we study the non-diffusive system of equations

The system (4.1) is strictly hyperbolic when $S>1$ , with three real characteristic wave speeds given by

We note, however, that in the limit $S\rightarrow 1$ there is a loss of hyperbolicity since the eigenvectors of the characteristic equation for the system (4.1) do not span $\mathbb{R}^{3}$ , although the wave speeds remain real and equal to $c_{0}$ ; the system is formally parabolic for $S=1$ (Scase et al. Reference Scase, Caulfield, Dalziel and Hunt2006b ). As we demonstrate below, the change in the characteristic structure, and so change in type, of the governing equations that occurs for $S=1$ fundamentally alters solutions of the system of equations. (It should be noted that if a shape factor for the buoyancy flux that differs from unity is included in (4.1c ), then $c_{0}=S_{f}M/Q$ while $c_{\pm }$ remain unchanged, and a loss of hyperbolicity occurs if $S_{f}=S\pm \sqrt{S(S-1)}$ .)

In order to establish well-posedness of the system of equations, we analyse the evolution of small perturbations to the steady solutions. It is convenient for subsequent analysis and numerical computations to factor out the steady solution and to introduce a new spatial coordinate, $x=(q_{0}/m_{0})z^{4/3}$ . We consider the stability of the steady solution to small perturbations. Therefore, we introduce an ordering parameter ${\it\epsilon}\ll 1$ and perturbations to the steady solution of the form

and linearize the governing equation (4.1) to obtain

where $\boldsymbol{q}=(q,m,f)^{\text{T}}$ and the matrices $\unicode[STIX]{x1D63C}$ , $\unicode[STIX]{x1D63D}$ and $\unicode[STIX]{x1D63E}$ are given by

When $S\equiv 1$ the system (4.1) was shown to be ill-posed by Scase & Hewitt (Reference Scase and Hewitt2012), who examined the response of the linearized system of equations to small-amplitude harmonic variation of frequency ${\it\omega}$ in the source buoyancy flux, and found a closed form solution for the perturbations in terms of special functions. It was shown that the amplitude of the perturbations grows with the distance from the source as ${\it\omega}^{-7/8}x^{-7/8}\exp (\sqrt{15{\it\omega}x}/4)$ for $x\gg 1$ (Scase & Hewitt Reference Scase and Hewitt2012). Thus, the steady solutions are linearly unstable, as the amplitudes of the perturbations grow as they propagate away from the source, and, furthermore, the system of equations (4.1) with $S=1$ is ill-posed since high-frequency fluctuations increase in amplitude most rapidly with no cutoff (Scase & Hewitt Reference Scase and Hewitt2012).

For $S>1$ we are unable to find closed form solutions of the linear system (4.4). We therefore examine the evolution of perturbations in the far field through an asymptotic analysis of the linear system (4.4) in the neighbourhood of the irregular singular point $x\rightarrow \infty$ (Bender & Orszag Reference Bender and Orszag1978). We consider two stability problems: an initial-value problem where the evolution of an arbitrary initial perturbation with compact support

is investigated, and a boundary-value problem where a fluctuation at the source $x=0$ is imposed and the response of the system downstream of the source is examined. Both of these problems can be analysed conveniently through the use of integral transforms of the linear system (4.4); a Laplace transform in time for the initial-value problem and a Fourier transform in time for the boundary-value problem. Denoting the Laplace transform of $\boldsymbol{q}(x,t)$ as $\hat{\boldsymbol{q}}(x,p)$ (with the Fourier transform obtained by taking $p=\text{i}{\it\omega}$ , where ${\it\omega}$ is the frequency of the harmonic source fluctuation imposed in the boundary-value problem), we obtain

The far-field behaviour is obtained conveniently by letting $\hat{\boldsymbol{q}}=\boldsymbol{f}(x)\text{e}^{g(x)}$ (Bender & Orszag Reference Bender and Orszag1978), with $g(x)$ a singular function and $\boldsymbol{f}(x)$ regular as $x\rightarrow \infty$ , so $\boldsymbol{f}(x)=\boldsymbol{f}_{0}+\boldsymbol{f}_{1}x^{-1}+\boldsymbol{f}_{2}x^{-2}+\cdots \,$ for $x\gg 1$ . The linear system (4.7) can then be written as

A leading-order balance requires $g(x)\sim p{\it\lambda}x$ as $x\rightarrow \infty$ . Then ${\it\lambda}$ and $\boldsymbol{f}_{0}$ satisfy the generalized eigenproblem $\unicode[STIX]{x1D63C}\boldsymbol{f}_{0}=-4{\it\lambda}\unicode[STIX]{x1D63D}\boldsymbol{f}_{0}/3$ , and therefore

where $\boldsymbol{f}_{0}^{L}$ denotes the left eigenvector satisfying $\boldsymbol{f}_{0}^{L}\unicode[STIX]{x1D63C}=-4{\it\lambda}\boldsymbol{f}_{0}^{L}\unicode[STIX]{x1D63D}/3$ .

To proceed further we let $g(x)\sim p{\it\lambda}x+{\it\mu}\log x$ for $x\gg 1$ . Substitution into (4.8) and balancing coefficients of $x$ gives, at order $O(1/x)$ ,

and so, by multiplying on the left by $\boldsymbol{f}_{0}^{L}$ , we find

and therefore

The coefficient vectors ( $\boldsymbol{f}_{1}$ , $\boldsymbol{f}_{2}$ , etc.) can be determined from higher-order terms in the expansion but are not given here.

The leading-order behaviour in the far field is therefore given by

where ${\it\delta}(x)$ denotes the Dirac delta function, corresponding to the propagation of discontinuities whose strength varies algebraically with distance from the source. The amplitude of the perturbations grows if $1<S<25/24$ , whereas the amplitude decays algebraically if $S>25/24$ , as shown in figure 2. Importantly, the algebraic growth rates ${\it\mu}$ do not depend on the transform variable $p$ . Therefore, for the initial-value problem the evolution of perturbations in the far field is independent of the spatial length scale of the initial perturbation, whereas when considering a boundary-value problem the growth rate of the perturbations does not depend on the frequency of the harmonic boundary oscillations. Therefore, the far-field asymptotic analysis shows that the system (4.1) with $S>1$ is well-posed as there is a continuous dependence of the solutions on the initial or boundary data, and small scales (either spatial scales in the case of an initial-value problem or temporal scales for a boundary-value problem) are not amplified more rapidly than longer scales. This is in contrast to the far-field limit of the solution of a boundary-value problem that is determined by Scase & Hewitt (Reference Scase and Hewitt2012), where exponential growth of the amplitude of perturbations is found with a growth rate that increases exponentially with $\sqrt{{\it\omega}}$ , and therefore when $S=1$ the system of equations is ill-posed. The distinguished behaviour for $S=1$ is apparent in the eigenvectors $\boldsymbol{f}_{0+}$ and $\boldsymbol{f}_{0-}$ , which are singular in the limit $S\rightarrow 1$ . Numerical solutions demonstrating the evolution of a localized initial perturbation are presented in § 5.

Figure 2. The exponents of the far-field algebraic growth rates of the amplitude of the volume flux perturbation, ${\it\mu}_{+}$ (dashed line) and ${\it\mu}_{-}$ (solid line), as functions of the shape factor $S$ . The perturbation to the volume flux grows with increasing distance from the source when ${\it\mu}_{-}>0$ , which occurs for $S<25/24$ .

5 Solutions of the well-posed unsteady plume model

We consider now solutions of the nonlinear system of equations (4.1) with $S>1$ . As the system of equations (4.1) is hyperbolic in this parameter regime, solutions may exhibit discontinuities, across which we enforce the following jump conditions that respectively impose the conservation of mass, momentum and buoyancy fluxes over a discontinuity at $z=z_{s}(t)$ , moving with velocity $c=\,\text{d}z_{s}/\,\text{d}t$ ,

These jump conditions admit two different types of discontinuous solutions, which are more easily analysed in terms of the primitive variables $b$ , $W$ and $G^{\prime }$ . Further denoting the values of the dependent variables either side of the discontinuity with superscripts $^{+}$ and $^{-}$ , corresponding to $z=z_{s}^{+}$ and $z=z_{s}^{-}$ respectively, we find that provided the velocity is discontinuous $(W^{+}\neq W^{-})$ then the reduced gravity $G^{\prime }$ is continuous (i.e. $G^{\prime +}=G^{\prime -}$ ). Furthermore, in this case, by eliminating $b^{+}$ and $b^{-}$ , we find that

When the shape factor is equal to unity, the only solution is $c=W^{+}=W^{-}$ , and this contradicts the assumption that there is a discontinuity. Thus, for $S=1$ it is not possible to construct discontinuous solutions. However, when $S>1$ , we find

where the negative root has been chosen so that $W^{-}>c$ , a condition required for causality. The other type of discontinuity occurs when $W^{+}=W^{-}=c$ . Then the radius of the plume is also continuous, $b^{+}=b^{-}$ , but the reduced gravity may be discontinuous and potentially even unbounded. We refer to this latter form as a ‘contact’ discontinuity, which occurs due to a linearly degenerate field in the system of equations (Lax Reference Lax1973).

To calculate numerically solutions of the nonlinear hyperbolic system of equations, we employ the non-oscillatory central scheme of Kurganov & Tadmor (Reference Kurganov and Tadmor2000) with a third-order total variation diminishing Runge–Kutta time stepping scheme (Gottlieb & Shu Reference Gottlieb and Shu1998), using an adaptive time step that ensures that the Courant–Friedrichs–Lewy (CFL) number remains fixed at $1/8$ to maintain numerical stability. The high-resolution central scheme is a shock-capturing numerical method for conservation laws (Kurganov & Tadmor Reference Kurganov and Tadmor2000) and has been used extensively to compute solutions of nonlinear hyperbolic systems. For numerical convenience, we factor out the steady solutions from the dependent variables and compute solutions in the transformed spatial coordinate $x$ as, although the system is not autonomous in this representation, the characteristic wave speeds remain bounded near to the source (i.e. as $x\rightarrow 0$ ) when the steady solution is realized, whereas the wave speeds of the system (4.1) diverge as $z\rightarrow 0$ . Therefore, in the mapped variables, we can take larger time steps while maintaining numerical stability.

We consider first the evolution of a perturbation to the steady solution (2.13) to confirm the far-field asymptotic analysis in § 4. We take as an initial condition a Gaussian perturbation (in the transformed coordinate $x$ ) to the steady momentum flux of the form

while the volume flux and buoyancy flux are not perturbed from the steady values, and take the shape factor to be $S=1.1$ ( ${>}25/24$ ). We compute numerically the evolution of the perturbation using the Kurganov & Tadmor (Reference Kurganov and Tadmor2000) central scheme for the nonlinear system rather than using the linearized equations, since the linearization introduces a subtle structural change in the governing equations: for the linearized system, all fields are (trivially) linearly degenerate and so discontinuities are of the contact discontinuity type (Lax Reference Lax1973), whereas two fields in the nonlinear system are genuinely nonlinear. We find that the Kurganov & Tadmor (Reference Kurganov and Tadmor2000) scheme resolves the shocks and rarefactions associated with the genuinely nonlinear fields accurately, but has less accuracy when contact discontinuities occur for linearly degenerate fields (see, e.g., Kurganov & Petrova Reference Kurganov and Petrova2000). In order to track the evolution of the perturbations to long times (and so large distance from the source), we implement a moving numerical domain, using the characteristic wave speeds of the linearized system to advect the lower and upper grid points, with boundary conditions in the moving domain given by the respective steady states. The numerical domain is then advected downstream, and the number of grid points increases as the spatial extent of the perturbation grows in order to maintain a specified numerical resolution.

In figure 3 we illustrate the evolution of the initial perturbation for dimensionless times $t\leqslant 40$ (noting that, for pure-plume source conditions, the time scale of motion, ${\it\tau}$ , is related to the length scale, $Z$ , and the source buoyancy flux, $F_{0}$ , through ${\it\tau}=Z^{4/3}F_{0}^{-1/3}$ ), showing the perturbations to the plume radius, mass flux and buoyancy flux as functions of the scaled spatial coordinate $x$ . The initial perturbation develops into a pulse whose spatial extent grows as it propagates. The amplitude of the perturbation slowly decreases, demonstrating linear stability of the system of equations for $S=1.1$ . The three discontinuities are apparent in the evolving perturbation, with the contact discontinuity that propagates at the speed of the intermediate characteristic most clearly seen in the perturbation to the reduced gravity, $G^{\prime }$ .

Figure 3. The perturbations to the fluxes of volume, $Q/Q_{0}-1$ (a), momentum, $M/M_{0}-1$ (b), and buoyancy, $F/F_{0}-1$ (c), and to the plume radius, $b/b_{0}-1$ (d), the plume velocity, $W/W_{0}-1$ (e), and the reduced gravity, $G^{\prime }/G_{0}^{\prime }-1$ (f), as functions of the scaled distance from the source, $x$ , for dimensionless times $t=0$ (orange), $t=10$ (red), $t=20$ (green), $t=30$ (blue), and $t=40$ (black). The momentum shape factor is fixed at a value $S=1.1>25/24$ . The numerical solution scheme adopts a moving and expanding spatial grid with a fixed grid spacing of ${\rm\Delta}x=0.01$ .

We introduce integral measures of the amplitude of the perturbations, with

for the volume flux, and similarly for the fluxes of momentum and buoyancy. Figure 4(a) shows the time evolution of $\mathscr{I}(Q)$ , $\mathscr{I}(M)$ and $\mathscr{I}(F)$ when the shape factor $S=1.1$ . Following a short transient reorganization of the initial condition for $t<5$ (not shown in the figure), the perturbations to the steady solution decay as they propagate through the domain. The perturbations to the steady fluxes of volume and momentum closely follow the prediction of the far-field analysis for $t>100$ , at which time the signal of the perturbation that travels with the fastest characteristic has reached $x\approx 197$ . The decay in the perturbation to the buoyancy flux is less rapid, and the rate of decay diminishes as time progresses. This is expected from the far-field asymptotic analysis of the linearized system, where a component of the buoyancy flux is found to be advected at the speed of the plume without change in amplitude (see (4.12a )).

Figure 4. Integral measures of the size of perturbations from the steady solutions for the volume flux, $\mathscr{I}(Q)$ (solid line), momentum flux, $\mathscr{I}(M)$ (dashed line), and buoyancy flux, $\mathscr{I}(F)$ (dotted line), as functions of the dimensionless time $t$ , for a shape factor (a) $S=1.1>25/24$ and (b) $S=1.02<25/24$ . The thick grey line illustrates the predicted growth rate of perturbations of $t^{2{\it\mu}_{+}}$ with (a) ${\it\mu}_{+}\approx -0.49$ for $S=1.1$ and (b) ${\it\mu}_{+}\approx 0.65$ for $S=1.02$ obtained from the far-field asymptotic analysis (4.12b ).

For $1<S<25/24$ our far-field asymptotic analysis shows that the steady solutions are linearly unstable. Numerical solutions of the governing equations support the asymptotic analysis, as illustrated in figure 4(b), which shows the time evolution of the integral measures of the size of the perturbations for a shape factor $S=1.02$ ( ${<}25/24$ ). For early times ( $t<10$ ) the perturbations to the steady solutions grow algebraically with a growth rate that is close to the growth rate predicted from the far-field analysis. We note that at early times the perturbation has not propagated far downstream from the localized initial condition and therefore the far-field asymptotic analysis would not be expected to describe fully the evolution. At later times ( $t>50$ ) the growth of the perturbations deviates substantially from the predictions of the linear far-field analysis as nonlinear effects begin to influence the evolution.

Numerical solutions of the governing equations also support the far-field asymptotic analysis of the linearized system when a harmonic oscillation of the source boundary conditions is imposed. Figure 5 illustrates the spatial structure of the perturbation to the steady volume flux at time $t=99$ for a shape factor of $S=1.05>25/24$ (figure 5 a) and $S=1.03<25/24$ (figure 5 b). The far-field analysis predicts linear stability when $S=1.05$ and instability when $S=1.03$ , and this is observed in the numerical solutions, with the amplitude of perturbations decaying with increasing distance downstream when $S=1.05$ (figure 5 a), in contrast to the growth in the amplitude of the perturbations that is seen when $S=1.03$ (figure 5 b).

Figure 5. The perturbation to the steady volume flux due to harmonic oscillation of the source buoyancy flux with period $2{\rm\pi}$ as a function of the scaled spatial coordinate $x$ at dimensionless time $t=99$ for a shape factor of (a) $S=1.05>25/24$ and (b) $S=1.03<25/24$ . The red dashed line illustrates an algebraic growth of perturbations of the form $cx^{{\it\mu}}$ , where (a) ${\it\mu}=-0.12$ for $S=1.05$ and (b) ${\it\mu}=0.26$ for $S=1.03$ . The prefactor $c$ is selected to provide an envelope of the numerical solution.

6 Similarity solutions for the unsteady evolution of plumes following an abrupt change in the source buoyancy flux

We examine now the nonlinear evolution of plumes following an abrupt change in the source conditions. Specifically, we consider solutions of the system of nonlinear evolution equations (4.1) with $S>25/24$ (so that the steady states are linearly stable) for an initial boundary-value problem in which the magnitude of the buoyancy flux at the source is instantaneously changed from $F_{0}$ to $F_{1}$ at $t=0$ . For $t<0$ the plume is in a steady state given by (2.13) and (2.14). We calculate the unsteady evolution as the plume adjusts to the new steady state in which the buoyancy flux, $F_{0}$ , in (2.13) and (2.14) is replaced by $F_{1}$ . Thus, the initial conditions are given by (2.13) and (2.14), while for $t>0$ the new source conditions for the plume are given by $F(0,t)=F_{1}$ and $Q(0,t)=M(0,t)=0$ .

Figure 6. The scaled width, $B(z)=z\widehat{Q}/\widehat{M}^{1/2}$ , the scaled velocity, $W(z)=z^{-1/3}\widehat{M}/\widehat{Q}$ , and the buoyancy flux, $\widehat{F}$ , as functions of the distance from source, $z$ at dimensionless times $t=2.4$ (orange), $t=4.9$ (red), $t=7.4$ (green) and $t=9.9$ (blue) for an increase in buoyancy flux $(\mathscr{F}=0.05)$ at $t=0$ and with shape factor $S=1.1$ . In (b) the scaled velocity $W=z^{-1/3}$ , corresponding to the steady velocity field of the original buoyancy flux at the source, and the scaled velocity $W=z^{-1/3}20^{1/3}$ , corresponding to the steady velocity field associated with the new buoyancy flux, are plotted with dashed lines. In (c) the original buoyancy flux, $F=F_{0}=0.05$ , and the new steady-state buoyancy flux, $F=F_{1}=1$ , are plotted with dashed lines.

The numerical solutions demonstrate that the adjustment occurs by the advection of an unsteady ‘pulse’ through the environment (see figures 6, 8 and 10 for examples of computations), which is modelled by the nonlinear evolution equations (4.1). In this section we draw out the self-similar adjustment that occurs in the dependent variables.

The similarity variable is established through the following scaling analyses, which require all of the terms in the governing partial differential equations (4.1) to be of the same order. To this end, to balance all of the derivative terms we require $M\sim Qz/t$ . Turning then to the ‘source’ terms, from (4.1a ) we have $Q/z\sim M^{1/2}$ , while from (4.1b ) $M/z\sim QF_{1}/M$ , where we have used $F\sim F_{1}$ as the scale of the buoyancy flux. Together these yield $z^{4}\sim F_{1}t^{3}$ , and we therefore seek similarity solutions to the system in terms of this gearing between the spatial and temporal scales. We note that the existence of this similarity grouping was identified by Scase et al. (Reference Scase, Caulfield, Dalziel and Hunt2006b ) (see also Scase et al. Reference Scase, Caulfield and Dalziel2008), although they did not construct the complete similarity solution; indeed, for their system with $S=1$ it was not possible to do so because the ill-posedness of the system is manifested as irreconcilable divergences in the solution (see § 4). Instead, Scase et al. (Reference Scase, Caulfield, Dalziel and Hunt2006b ) found a separable solution that did not satisfy the boundary conditions but did capture some of the features found in their numerical computation. We discuss in appendix C the counterpart to their separable solution in the now regularized system of governing equations. Before analysing the form of the similarity solutions, we note that there is another reason for anticipating the similarity scaling deduced above: the system is hyperbolic for $S>1$ , and all of the characteristics are proportional to $W=M/Q\propto F_{1}^{1/3}z^{-1/3}$ . The characteristic curves are of the form $\,\text{d}z/\,\text{d}t\propto W$ , thus also admitting the similarity scaling $z^{4}\sim F_{1}t^{3}$ , and it is through the evolution of the characteristics that the system adjusts to its new state.

We now seek similarity solutions for the fluxes of volume, momentum and buoyancy, which we write in the following form:

where the similarity variable is given by

and the similarity functions $\widehat{Q}({\it\eta})$ , $\widehat{M}({\it\eta})$ and $\widehat{F}({\it\eta})$ are to be determined. In this form the steady state given by (2.13) corresponds to constant values of $\widehat{Q}$ , $\widehat{M}$ and $\widehat{F}$ . It is convenient to further substitute $\widehat{Q}={\it\eta}\widehat{q}$ , $\widehat{M}={\it\eta}^{2}\widehat{m}$ and $\widehat{F}={\it\eta}^{3}\widehat{f}$ because this simplifies the governing equations, which are now given by

Symbolically, we write this coupled differential system as $\unicode[STIX]{x1D63F}{\it\eta}\,\text{d}\widehat{\boldsymbol{q}}/\,\text{d}{\it\eta}=\boldsymbol{b}$ , noting that both the matrix $\unicode[STIX]{x1D63F}$ and the vector $\boldsymbol{b}$ are only functions of the dependent variables $\widehat{\boldsymbol{q}}=(\widehat{q},\widehat{m},\widehat{f})$ . We note that the ‘separable’ similarity solutions derived by Scase et al. (Reference Scase, Caulfield, Dalziel and Hunt2006b ) correspond to $\widehat{\boldsymbol{q}}=(\widehat{q},\widehat{m},\widehat{f})=\text{constant}$ , and these are discussed in appendix C. For the motion driven by an abrupt change in the magnitude of the buoyancy flux at source, the vital parameter is the ratio of the initial to final buoyancy fluxes at the source given by $\mathscr{F}=F_{0}/F_{1}$ . We note that the initial conditions correspond to

In terms of the similarity functions, these are enforced in the far field ( ${\it\eta}\rightarrow \infty$ ) and correspond to the region that is unaffected by the change of the buoyancy flux at the source. The source conditions demand that

Constructing the similarity solutions then corresponds to integrating the governing system of ordinary differential equations (6.3), subject to these conditions.

The matrix $\unicode[STIX]{x1D63F}$ is singular when

This corresponds to locations where $\widehat{q}/\widehat{m}=4/3$ and $\widehat{q}/\widehat{m}=4(S\pm \sqrt{S^{2}-S})/3$ . These are singular points of the governing system of equations and are of significance because the dependent variables may have discontinuous gradients at these locations. We show in appendix D how to derive the local behaviour close to the singular points; this is required to initiate numerical integration from these locations.

The similarity solutions may also feature discontinuous solutions, in which case we define the shock position, $z_{s}(t)$ , scaled according to the similarity variables to be given by

where ${\it\eta}_{s}$ is a constant. In this case of discontinuous solutions, the jump conditions relate the dependent variables at ${\it\eta}={\it\eta}_{s}^{+}$ to those at ${\it\eta}={\it\eta}_{s}^{-}$ and are given by

Key locations in the similarity solutions are the points at which the dependent variables transition from the new steady state to an unsteady pulse and then from this unsteady pulse to the original steady state. The family of slowest moving characteristics associated with the new source is given by

The transition between the steady and unsteady portions of the solution must occur at a singular point of $\unicode[STIX]{x1D63F}$ to permit the gradient of the solution to change discontinuously. Thus, in this case $\widehat{q}/\widehat{m}=4(S-\sqrt{S(S-1)})/3$ , and so we find that

This implies that the transition point occurs at a constant value of the similarity variables, ${\it\eta}={\it\eta}_{c1}$ , given by

using the condition (6.5). Likewise, the family of fastest moving characteristics is given by

and so the boundary between the unsteady pulse and the steady solution associated with the original buoyancy flux occurs at a constant value of the similarity variable, ${\it\eta}={\it\eta}_{c2}$ , given by

using the condition (6.4), provided that the motion is due to a decrease in the source strength $(\mathscr{F}>1)$ . If the source strength increases then the characteristics associated with the new release overtake those due to the original source and, as we show below, the motion forms ‘shocks’.

7 Discussion and conclusion

Steady plume models are applicable on time scales that are long compared with the eddy turnover time that characterizes transient turbulent features, and on time scales shorter than the characteristic time for source variations (Scase et al. Reference Scase, Caulfield, Dalziel and Hunt2006b ). In many applications, the effect of variations of the source fluxes of mass, momentum or buoyancy on the plume dynamics is of fundamental importance. Here, unsteady models are essential. For an example, during volcanic eruptions the source conditions can fluctuate in time due to unsteadiness in the physical processes occurring in the volcanic conduit and at the vent. End members of the source unsteadiness are the initiation phase (when an eruption first produces material that becomes buoyant and ascends through the atmosphere) and the waning phase (when the eruption comes to an end through either a gradual reduction in the rate at which material is erupted or an abrupt cessation of the activity) of the eruption.

The unsteady model of turbulent plumes analysed by Scase et al. (Reference Scase, Caulfield, Dalziel and Hunt2006b ) and Scase (Reference Scase2009) identified a key feature of the unsteady response of plumes to abrupt changes in the source buoyancy flux: the plume adjusts to the new source conditions through a ‘pulse’ that propagates through the plume. However, the model of Scase et al. (Reference Scase, Caulfield, Dalziel and Hunt2006b ) adopts top-hat profiles for the mean axial velocity and was shown to be ill-posed in the analysis of Scase & Hewitt (Reference Scase and Hewitt2012), with the numerical solutions of Scase et al. (Reference Scase, Caulfield, Dalziel and Hunt2006b ) and Scase (Reference Scase2009) attained only due to significant numerical viscosity.

The diffusive regularization of the unsteady model proposed by Scase & Hewitt (Reference Scase and Hewitt2012) has an intuitive physical basis: turbulent eddies have a vertical length scale over which different levels of the plume are connected. Numerical solutions (Scase & Hewitt Reference Scase and Hewitt2012) suggest that the phenomenological model of the turbulent diffusion of momentum advocated by Scase & Hewitt (Reference Scase and Hewitt2012) results in a well-posed unsteady model, but the analysis presented here shows that this modification of the governing system of equations renders the classical power-law solution for steady plumes unstable. As the power-law solutions have strong empirical support (e.g. Morton et al. Reference Morton, Taylor and Turner1956; Papanicolaou & List Reference Papanicolaou and List1988; Shabbir & George Reference Shabbir and George1994), we conclude that the diffusive regularization of Scase & Hewitt (Reference Scase and Hewitt2012) is not appropriate.

Our analysis identifies a different physical process, a difference in the rates of transport of the cross-sectionally averaged mass and momentum of the plume (referred to as type I dispersion by Craske & van Reeuwijk Reference Craske and van Reeuwijk2015a ,Reference Craske and van Reeuwijk b ) which occurs due to non-uniform radial profiles of the mean axial velocity in the plume, and we have shown that including this process through a momentum shape factor leads to a well-posed system of equations. However, this requires the top-hat description of the radial profiles for the axial velocity in the plume, which has been applied extensively (for convenience) in models of steady plumes, to be replaced with a description of the radial variation. The resulting integral model for unsteady plumes introduces only one additional parameter, the momentum shape factor, over the classical steady plume model (Morton et al. Reference Morton, Taylor and Turner1956). Furthermore, the (ensemble averaged) radius of the plume remains dependent on the entrainment coefficient alone, so the unsteady model does not preclude calibration of the entrainment coefficient from laboratory experiments that measure the steady plume width (or spreading angle). To determine the momentum shape factor, the radial profile of the axial velocity can be measured directly (e.g. Papanicolaou & List Reference Papanicolaou and List1988; Shabbir & George Reference Shabbir and George1994), and the shape factor computed from (2.6).

Explicitly including a momentum shape factor that differs from unity results in a strictly hyperbolic system of equations that governs the plume dynamics. This is appealing, as laboratory experiments (Scase et al. Reference Scase, Caulfield and Dalziel2008) and numerical modelling (Scase et al. Reference Scase, Aspden and Caulfield2009) show that the adjustment of the plume to changes in the source conditions occurs through the propagation of an unsteady pulse, as predicted by our non-diffusive hyperbolic model. The development of the unsteady pulse following an abrupt change in the source buoyancy flux is described by a similarity solution of the hyperbolic system of equations. The construction of the similarity solution allows us to identify three qualitatively different regimes of the unsteady evolution. Following an increase in the source buoyancy flux, the pulse takes the form of a localized increase in the plume width with a leading discontinuity. If the source buoyancy flux is reduced then the plume width narrows. For a relatively strong reduction in the source buoyancy flux an internal shock occurs in the similarity solution, whereas no internal shock is found if the source buoyancy flux is reduced by a smaller amount. We expect that diffusive processes that are not represented in our formulation will act to locally smooth the sharp gradients that appear in solutions of the hyperbolic model. However, hyperbolic models have been shown to capture the dominant flow dynamics in many settings (Whitham Reference Whitham1974).

Craske & van Reeuwijk (Reference Craske and van Reeuwijk2015a ) identify, from direct numerical simulations, the dispersion of momentum as a fundamental feature of turbulent jets, and construct an integral model to describe unsteady jets which includes a description of the non-uniform radial profile of the vertical velocity (Craske & van Reeuwijk Reference Craske and van Reeuwijk2015b ) through shape factors. In the model of Craske & van Reeuwijk (Reference Craske and van Reeuwijk2015b ) the integral conservation equations are derived from the point-wise momentum and energy conservation equations, following the approach of Priestley & Ball (Reference Priestley and Ball1955), and an energy shape factor rather than a momentum shape factor is introduced. The resulting system of integral conservation equations shares some features with the unsteady integral model proposed here (albeit for a jet rather than a plume), in particular the hyperbolic structure of the system of equations.

For jets the momentum of the flow as it is ejected from a source drives the motion, and there is significant evolution of the radial profile of the axial velocity (referred to as type II dispersion by Craske & van Reeuwijk Reference Craske and van Reeuwijk2015a ,Reference Craske and van Reeuwijk b ) as the flow develops away from the source. The evolution of the shape of the radial profiles of axial velocity and buoyancy has been incorporated into a non-constant entrainment coefficient by Kaminski et al. (Reference Kaminski, Tait and Carazzo2005) and Carazzo et al. (Reference Carazzo, Kaminski and Tait2006). Furthermore, in the unsteady integral model of Craske & van Reeuwijk (Reference Craske and van Reeuwijk2015b ) the deviation of the radial profiles from self-similar forms gives rise to diffusive terms in the system of integral conservation equations. Our analysis shows that a description of momentum dispersion, through a shape factor that differs from unity, is sufficient to obtain a well-posed model of unsteady plumes. Therefore, while for jets it is necessary to account for the evolution of the radial profile of axial velocity to self-similar form, for unsteady pure plumes the classical entrainment closure of Morton et al. (Reference Morton, Taylor and Turner1956) remains applicable.

The mathematical model we present allows predictions of unsteady plume dynamics to be made and, while our study has focused on changes to the source buoyancy flux, the effect of general temporal variations in the source conditions can be examined. Furthermore, our framework can be applied to describe the unsteady dynamics of plumes in industrial and environmental settings, where additional physical processes such as an ambient flow, thermodynamics and particle transport have a strong influence on the evolution of the plume.