2.1 Reynolds equations

The flow with which we are concerned is a round turbulent plume that is swirl free and statistically axisymmetric. The plume is comprised of fluid that is of lower density than its surroundings, which are assumed to be of uniform density ${\it\rho}_{0}$ , and therefore has positive buoyancy

where ${\it\rho}$ is the fluid density. Using the Boussinesq approximation, the equations governing the transport of volume, longitudinal specific momentum and buoyancy are

where $\boldsymbol{u}\equiv (u,v,w)$ is the velocity, with components in the radial, azimuthal and longitudinal directions, respectively, $p$ is the kinematic pressure relative to a hydrostatic balance and ${\it\nu}$ and ${\it\kappa}$ are the kinematic viscosity and buoyancy diffusion coefficient, respectively. Noting the statistical axisymmetry of the flow, an ensemble and azimuthal average of (2.2)–(2.4) in the limit of high Reynolds and Péclet numbers yields

Here $\overline{{\it\chi}}$ denotes the ensemble average of the variable ${\it\chi}$ , so that ${\it\chi}\equiv \overline{{\it\chi}}+{\it\chi}^{\prime }$ , where $\overline{{\it\chi}^{\prime }}=0$ . Multiplication of (2.6) by $2\overline{w}$ yields an equation for the mean longitudinal kinetic energy:

Due to the fact that $\overline{w}\gg \overline{u},\overline{v}$ we will omit ‘longitudinal’ in describing (2.8) hereafter, referring to it instead as the equation for the mean kinetic energy. For further details pertaining to the derivation of (2.5)–(2.8) the reader is referred to Craske & van Reeuwijk (Reference Craske and van Reeuwijk2015a ).

Having been obtained via invertible manipulations, the momentum–energy system (2.6) and (2.8), is equivalent to the volume–momentum system, (2.5) and (2.6). Both satisfy local volume and momentum conservation, but whereas (2.5) is a diagnostic relation or constraint, (2.8) is a prognostic equation. In the following section we will show that the momentum–energy system constitutes the natural choice for developing and understanding integral models of unsteady plumes.

2.2 Integral equations

Integration of (2.6)–(2.8) with respect to $r$ from zero to the radial extent of the plume $r_{d}$ , yields the integral equations

where the dependent variables are the longitudinal volume flux $Q$ , the longitudinal specific momentum flux $M$ (hereafter referred to as the momentum flux for brevity) and the integral buoyancy $B$ :

The Greek letters in (2.9)–(2.11) are dimensionless profile coefficients that allow us to express unknown integrals in terms of the dependent variables $Q$ , $M$ and $B$ , and are defined as

We refer to (2.9)–(2.11) as generalised unsteady plume equations, because we have not yet made an assumption about the way in which the dimensionless profile coefficients might depend on $z$ , $t$ , $Q$ , $M$ or $B$ . Physically, the profile coefficients ${\it\beta}_{g}$ , ${\it\gamma}_{g}$ , ${\it\theta}_{g}$ and ${\it\delta}_{g}$ correspond to gross dimensionless fluxes of momentum, energy and buoyancy, and the gross dimensionless production of turbulence kinetic energy (including the redistribution of energy via pressure), respectively. Note that ${\it\theta}_{g}$ is used in (2.11) because it accounts for the total transport of buoyancy, whereas ${\it\theta}_{m}$ is used in (2.10) because the forcing of the mean flow energy by buoyancy does not include turbulent transport ${\it\theta}_{f}$ . In addition, note that ${\it\beta}_{m}$ in (2.13), which we include for completeness, is unity by definition. In obtaining (2.9)–(2.11) it was assumed that the ensemble and azimuthally averaged velocity $\overline{w}$ and buoyancy $\overline{b}$ are equal to zero at $r=r_{d}$ .

Before proceeding, we note that the integral buoyancy $B$ is not used in steady plume theory, in which the effects of buoyancy are typically expressed in terms of the mean buoyancy flux $F\equiv {\it\theta}_{m}MB/Q$ (see, e.g. Hunt & van den Bremer Reference Hunt and van den Bremer2011, and references therein):

However, in the unsteady integral equations (2.9)–(2.11), the use of $B$ is convenient from both a mathematical and a conceptual perspective. With $B$ as a dependent variable, the temporal derivatives in (2.9)–(2.11) are decoupled and each obeys a classical conservation equation. Moreover, it is natural to view the buoyancy flux $F$ as an unknown quantity requiring assumptions because it depends on the correlation of $\overline{w}$ and $\overline{b}$ . In this regard, we also note that while we do not advocate using the continuity equation (2.5) directly, the behaviour of the volume flux $Q$ nevertheless plays a crucial role in the system (2.9)–(2.11) as a dependent variable.

2.3 Modelling assumptions

The central message of this paper is that each of the dimensionless profile coefficients appearing in (2.9)–(2.11) play an independent and non-trivial role in determining the structure of the governing equations. In a steady state, for which $\partial _{t}=0$ , the situation is different because the profile coefficients no longer play an independent role. Therefore, the form of the classical steady-state power-law solutions, which will be discussed in § 2.5, are essentially independent of the values of the profile coefficients. In a steady state the influence that the profile coefficients have on the behaviour of the system is hidden.

Figure 1. The mean dimensionless energy flux ${\it\gamma}_{m}$ associated with different radial profiles of the longitudinal velocity. The top-hat profile (a) is a limiting case for which the dimensionless energy flux is minimal, while in the case of a Gaussian profile (c)  ${\it\gamma}_{m}=4/3$ .

By definition, the assumed shape of the mean velocity profile does not affect the volume flux or the momentum flux in the plume, i.e. the volume flux and the momentum flux are equal to $Q$ and $M$ , respectively, regardless of the behaviour of $\overline{w}$ . Similarly the shape of the buoyancy profile does not, by definition, affect the integral buoyancy $B$ . In self-similar steady-state models of jets and plumes it is therefore conventional (Turner Reference Turner1973), and entails no loss of generality, to regard $\overline{w}$ as having a uniform distribution of amplitude $w_{m}$ for $r\leqslant r_{m}$ , which is known as a top-hat distribution, as illustrated in figure 1(a). However, the assumed velocity profile does affect the mean energy flux ${\it\gamma}_{m}M^{2}/Q$ in the plume, because ${\it\gamma}_{m}$ depends on the radial dependence of the velocity field. In unsteady jets and plumes, unlike their statistically steady counterparts, the energy flux plays an independent role in the governing equations and therefore, the assumption of a particular velocity profile does entail a loss of generality. A useful result in this regard is that for a Gaussian mean velocity profile $\overline{w}(r)$ , illustrated in figure 1(c), the dimensionless energy flux ${\it\gamma}_{m}$ can be determined exactly (Craske & van Reeuwijk Reference Craske and van Reeuwijk2015a ) as

The mean energy flux associated with top-hat velocity profiles, on the other hand, is minimal (assuming that $\overline{w}>0$ ), such that ${\it\gamma}_{m}=1$ . It is evident from (2.14) that any variation in the radial dependence of the velocity profile will result in ${\it\gamma}_{m}>1$ , as illustrated schematically in figure 1. Although Hewitt & Bonnebaigt (Reference Hewitt and Bonnebaigt2014) remark that integral models of turbulent plumes remove information pertaining to radial dependence, we find that this is not necessarily the case. Use of the mean kinetic energy equation exposes the fact that ${\it\gamma}_{m}$ appears as an independent parameter, determining the response of the plume to source perturbations and the trajectories of its characteristic curves.

2.4 Volume conservation

The motivation for working with the mean kinetic energy equation (2.8) rather than the continuity equation (2.2) is that the radial integral of (2.2),

does not provide an explicit prognostic equation (for definiteness and consistency with the classical approach, here we have assumed that $r_{d}\rightarrow \infty$ in (2.12)). Instead, (2.19), which is ostensibly identical to the steady-state volume flux equation given by Morton et al. (Reference Morton, Taylor and Turner1956), relates the volume flux in the plume to the induced radial velocity in the environment. Indeed, (2.19) can be used in place of the integral energy equation (2.10), but the modeller is then faced with the task of providing a closure for $ru|_{\infty }$ . In the steady state, Morton et al. (Reference Morton, Taylor and Turner1956) replace the right-hand side of (2.19) with $2{\it\alpha}_{0}M^{1/2}$ , thereby defining the classical entrainment coefficient ${\it\alpha}_{0}$ , to obtain a closed system of equations. In the unsteady case, however, an induced radial flow in the environment does not necessarily correspond to entrainment into the plume, because the radius of the plume might change with respect to time (cf. the temporal jet studied by van Reeuwijk & Holzner Reference van Reeuwijk and Holzner2014, for example). Conversely, efforts to obtain a prognostic equation for the area $r_{m}^{2}\equiv Q^{2}/M$ of the plume rely on particular assumptions regarding its velocity profile. In all but the simplest cases (e.g. the rigorously derived top-hat model of Scase et al. (Reference Scase, Caulfield, Dalziel and Hunt2006b )), such models have questionable foundations, because typically, ensemble-averaged quantities in plumes do not have a well-defined edge and depend continuously on the radial coordinate. The use of a generalised plume theory, resulting in (2.9)–(2.11), circumvents these issues altogether and allows one to understand unsteady plume models in a broader context. Consequently, the theory provides a physics-based means of understanding and remedying the problems associated with existing unsteady plume models (for details of these problems see Scase & Hewitt (Reference Scase and Hewitt2012)).

From the integral equations for momentum (2.9) and mean energy (2.10), one can obtain a prognostic equation for the area of the plume without assumption. To appreciate this, note that the continuity equation (2.2) was used to obtain the mean energy equation (2.8) from the mean momentum equation (2.6). Together, (2.6) and (2.8) therefore imply the continuity equation (2.2). Similarly, one can recover an integral volume conservation equation by combining integral equations for momentum (2.9) and mean energy (2.10). A consistent equation for the area of the plume is therefore obtained by expanding $\partial _{t}(Q^{2}/M)$ and substituting for $\partial _{t}Q$ and $\partial _{t}M$ using (2.9) and (2.10) respectively:

where

The advantage of (2.20) and (2.21) in comparison with (2.19) is that not only do they account explicitly for the temporal change of the area of the jet, they account for the process of entrainment in terms of known physical processes such as the production of turbulence kinetic energy and the transport of mean kinetic energy. As described in van Reeuwijk & Craske (Reference van Reeuwijk and Craske2015), (2.21) is an entrainment relation that ensures consistency between equations for volume, momentum and energy conservation. Equation (2.21) is not an entrainment model because we have not assigned values/functions to the profile coefficients ${\it\beta}_{g}$ , ${\it\gamma}_{g}$ , ${\it\theta}_{g}$ and ${\it\delta}_{g}$ by invoking assumptions.

There are two significant features of (2.20). First is the factor $1/{\it\gamma}_{g}$ preceding $\partial _{t}(Q^{2}/M)$ , which, by implying an effective area $Q^{2}/(M{\it\gamma}_{g})$ , provides a definition of the plume edge for continuous velocity profiles. Reassuringly, for top-hat profiles ${\it\gamma}_{g}=1$ , which makes the left-hand side of (2.20) consistent with the top-hat model of Scase et al. (Reference Scase, Caulfield, Dalziel and Hunt2006b ). Second, is the entrainment coefficient ${\it\alpha}$ , which is not necessarily constant. The first contribution on the right-hand side of (2.21) accounts for the production of turbulence kinetic energy and typically dominates the remaining terms. When the profile coefficients are constant, it is useful to regard the second and third terms on the right-hand side of (2.21) as relating to the integral advective acceleration of the plume $\partial _{z}M$ . Note, however, that due to entrainment $\partial _{z}M>0$ does not always imply a local advective acceleration $\partial _{z}w_{m}>0$ , but nevertheless accounts for the way in which the plume is forced in an integral sense. How entrainment is affected by the plume’s integral acceleration depends on the radial dependence of the velocity, turbulence and pressure. For top-hat profiles ${\it\gamma}_{g}={\it\beta}_{g}=1$ , the entrainment coefficient is insensitive to the integral acceleration of the plume and the second and third terms on the right-hand side of (2.21) are equal to zero. For further details regarding the physical interpretation of the first three terms on the right-hand side of (2.21), the reader is referred to Craske & van Reeuwijk (Reference Craske and van Reeuwijk2015a ), as they also appear in the equations describing an unsteady jet. The final contribution in (2.21) is only found in plumes and is proportional to the flux-balance parameter ${\it\Gamma}$ of Morton (Reference Morton1959). If ${\it\theta}_{m}>1$ , the buoyancy provides slightly more forcing in the energy equation (2.10) than one would expect from identically distributed $\overline{b}$ and $\overline{w}$ . Noting that the area $Q^{2}/M$ is inversely proportional to the momentum flux $M$ , the effect of ${\it\theta}_{m}>1$ is to reduce the entrainment coefficient. Conversely, when ${\it\theta}_{m}<1$ the buoyancy provides slightly less forcing in the energy equation than one might expect, and the entrainment coefficient increases.

In appendix A we derive generalised time-dependent similarity solutions of (2.9)–(2.11) involving an algebraic decrease in the plume’s source buoyancy flux. The solutions show that for Gaussian plumes the entrainment coefficient increases, relative to its steady-state value ${\it\alpha}_{0}$ , in such a way that the plume’s radius retains its steady-state dependence on $z$ and is independent of time. Conversely, in plumes with a top-hat velocity profile, the solutions predict that the entrainment coefficient is always equal to ${\it\alpha}_{0}$ and that the spreading rate of the plume decreases relative to that associated with the steady state due to the time dependence of the source conditions, the latter result being first obtained by Scase et al. (Reference Scase, Caulfield, Dalziel and Hunt2006b ).

2.5 The steady state

Using the momentum–energy formulation (2.9)–(2.11), it is interesting to see how the profile coefficients can be incorporated into the classical steady-state plume solutions. For constant source buoyancy flux $F_{s}$ and constant profile coefficients, the latter assumption being consistent with far-field self-similarity, a steady-state solution to (2.9)–(2.11), which we denote $(Q_{0},M_{0},B_{0})$ , is

where the constant steady-state entrainment coefficient, which follows from substitution of (2.22) into (2.10), is

Noting that the steady-state mean buoyancy flux $F$ is related to $F_{s}$ according to $F={\it\theta}_{m}F_{s}/{\it\theta}_{g}$ , and that $B=FQ/({\it\theta}_{m}M)$ , the solution to the steady-state integral buoyancy $B_{0}$ is

The effects of ${\it\theta}_{g}$ and ${\it\beta}_{g}$ are felt only via ${\it\alpha}_{0}$ and an effective buoyancy flux $F_{E}\equiv F_{s}/({\it\beta}_{g}{\it\theta}_{g})$ (van Reeuwijk & Craske Reference van Reeuwijk and Craske2015), and do not affect the classical power-law scaling of the steady-state solutions. Due to the fact that $F_{s}={\it\theta}_{g}F/{\it\theta}_{m}$ , (2.22) can also be expressed in terms of the mean buoyancy flux $F$ . Noting that the use of the total source buoyancy flux $F_{s}$ , rather than the effective source buoyancy flux $F_{E}$ , would lead to an over-estimation of $Q$ and $M$ , the solutions (2.22) might be useful to experimentalists who wish to compare data to theory using a known source buoyancy flux $F_{s}$ . Comparison of the system (2.22) to the classical plume solutions of Morton et al. (Reference Morton, Taylor and Turner1956) shows that the flux-balance parameter of Morton (Reference Morton1959) is

which characterises the relative importance of buoyancy compared with inertia in the flow. When $0<{\it\Gamma}<1$ , the plume is dominated by inertia and referred to as being ‘forced’; when $1<{\it\Gamma}$ the plume is dominated by buoyancy and is referred to as being ‘lazy’. Jets and pure plumes correspond to the special cases for which ${\it\Gamma}=0$ or ${\it\Gamma}=1$ , respectively. Thus, for the far-field similarity solutions (2.22), describing a pure plume, the fluxes of volume, momentum and buoyancy are balanced such that ${\it\Gamma}=1$ . In general ${\it\theta}_{m}\approx 1$ and ${\it\beta}_{g}>1$ , hence one would expect the classical definition of ${\it\Gamma}$ , namely $(5Q^{2}F)/(8{\it\alpha}_{0}M^{5/2})$ , to be slightly greater than unity, based on observations of $Q$ , $M$ and $F$ from a pure plume.

4.1 Steady plumes

Before discussing the simulations of unsteady plumes we discuss those of steady plumes, focusing on the values of the dimensionless profile coefficients that were introduced in § 2. The profile coefficients, defined by (2.13)–(2.16), determine the relative importance of each term in the governing integral equations (2.9)–(2.11). Physically, they can be viewed as dimensionless flux and turbulence production/redistribution terms. As explained in §§ 2.3 and 2.5, while the profile coefficients play a passive role in a steady state, we will demonstrate in § 5 that in an unsteady state they play an active role in determining the structure of the governing equations. In order to make predictions about unsteady plumes, it is therefore useful to establish the values of the profile coefficients to leading order by inspecting the steady-state plume data. Whether the profile coefficients are themselves affected by unsteadiness is a higher-order question that we defer until § 6.1.

Table 2 displays the values of the profile coefficients evaluated from the steady-state data provided by simulations L and H. The values were obtained by averaging over the interval $z/r_{s}\in [20,40]$ , in which the profile coefficients were observed to have reached an approximately constant value. The near field, in which the profile coefficients exhibit appreciable variation, is not investigated in this work. To begin, it is useful to recall that ${\it\beta}_{m}$ , ${\it\gamma}_{m}$ and ${\it\theta}_{m}$ are the leading-order dimensionless mean fluxes of momentum, energy and buoyancy. The dimensionless production ${\it\delta}_{m}$ is also of leading order but, due to the difference between the longitudinal and radial length scale of the plume, is smaller than ${\it\beta}_{m}$ , ${\it\gamma}_{m}$ and ${\it\theta}_{m}$ by a factor $O({\it\alpha})$ . The approximate equality between the profile coefficients in L and H supports the view that the simulations are well resolved and independent of Reynolds number. The validity of the simulations is further supported by the close agreement of the radial profiles of velocity, buoyancy and turbulent transport that we discuss in appendix B. Were the Reynolds number smaller, one might expect to see a difference between the data obtained from L and H, due to a possible dependence on the Reynolds number. We note that the observed entrainment coefficient ${\it\alpha}\approx 0.11$ (inferred by evaluating $\text{d}Q/\text{d}z/(2M^{1/2})$ ) is at the low end of the values that are reported from laboratory experiments, which typically range from 0.12 to 0.17 (Carazzo, Kaminski & Tait Reference Carazzo, Kaminski and Tait2006).

Noteworthy for the present study is the fact that ${\it\gamma}_{m}\approx 1.28$ is close to the value $4/3$ , associated with Gaussian velocity profiles (2.18), rather than 1, which is the value associated with top-hat velocity profiles. We remind readers that ${\it\gamma}_{m}=4/3$ is an exact result for the dimensionless energy flux associated with a Gaussian velocity profile, and is due to the cubic term $\overline{w}^{3}$ appearing in the integrand of (2.18). In contrast, one would expect the dimensionless mean buoyancy flux ${\it\theta}_{m}$ for Gaussian plumes to be equal to unity, which is consistent with what we observe in the DNS data. As expected, the gross dimensionless buoyancy flux exceeds that arising solely from the mean flow, resulting in ${\it\theta}_{g}>{\it\theta}_{m}$ . Indeed, contributions to momentum, energy and buoyancy transport from turbulence (see ${\it\beta}_{f}$ , ${\it\gamma}_{f}$ and ${\it\theta}_{f}$ , respectively) are of the order of 20 % and are therefore not insignificant.

4.2 Unsteady plumes

Following a step change in the buoyancy flux at the source of an otherwise steady plume a disturbance, or wave, propagates in the direction of the mean flow. In referring to this disturbance as a wave, we follow Whitham (Reference Whitham1974) and regard it as a recognisable signal propagating with a certain velocity. Here the notion of a wave is perhaps more appropriate than the notion of a front, which was employed in Craske & van Reeuwijk (Reference Craske and van Reeuwijk2015a ), as it encompasses a wider variety of possible disturbances, which can be comprised of several fronts. More precisely, in § 5.1 we will see that the wave is comprised of characteristic surfaces along which wave fronts travel.

Figure 2. Isoregions of the ensemble and azimuthally averaged buoyancy (red) and threshold of the instantaneous enstrophy (blue) in an unsteady turbulent plume at times $t_{n}\approx 1.95{\it\tau}_{s}n$ , $n=3,\ldots ,9$ , where ${\it\tau}_{s}\equiv r_{s}^{4/3}(F_{s}^{A})^{-1/3}$ . The buoyancy displayed in the figure has been non-dimensionalised using the local characteristic buoyancy scale $b_{m0}\equiv B_{0}M_{0}/Q_{0}^{2}$ of a steady plume with buoyancy flux $F^{B}$ .

The unsteady plume is illustrated in figure 2, which displays the azimuthal and ensemble-averaged normalised buoyancy field at several time instances, in addition to the instantaneous boundary of the plume, corresponding to a single member of the ensemble. The normalised buoyancy field was obtained by dividing $\overline{b}$ by a characteristic steady-state buoyancy $b_{m0}(z)\equiv B_{0}M_{0}/Q_{0}^{2}$ , based on observations from L (i.e. a steady plume with source buoyancy flux $F^{B}$ ). The boundary of the plume is defined by an isoline associated with a relatively small value of enstrophy, and therefore separates turbulent and non-turbulent parts of the flow. In the averaged buoyancy field one observes a smoothly defined cigar-shaped region penetrating progressively further into the domain. Although the size of the region increases with respect to time, its shape is approximately invariant over the range of times displayed in figure 2. Looking at the buoyancy field and the boundary of the plume in figure 2, one does not get the impression that the width or radial extent of the plume is strongly affected by the step change in the buoyancy flux. In this regard, we note that experimental observations of a plume, whose source buoyancy flux was suddenly reduced at the source, suggested that such plumes become narrower in the vicinity of the step change (Scase et al. Reference Scase, Caulfield and Dalziel2008). However, consistent with our interpretation of figure 2, the change in plume width is not readily discernible from the radial extent of the passive scalar field presented in figures 3 and 4 of Scase et al. (Reference Scase, Caulfield and Dalziel2008). Only by quantifying the width of the plume with a top-hat width based on the concentration of the passive scalar do Scase et al. (Reference Scase, Caulfield and Dalziel2008) find that the plume becomes narrower (figure 6 of their study). The interpretation and observed behaviour of the plume radius will be discussed in further detail below.

Figure 3. (a) Dimensionless buoyancy flux, (b) dimensionless momentum flux and (c) dimensionless volume flux, corresponding to individual members of the ensemble (thin lines) and their ensemble average (thick line) at $t/{\it\tau}_{s}=16$ . The width of the line denoting the ensemble average is equal to twice the standard deviation of the sample, divided by $\sqrt{n}$ , where $n=24$ is the number of members of the ensemble. The dashed lines correspond to steady-state data before and after the step change in the source buoyancy flux.

The primary focus of this study is the behaviour of integral quantities such as the volume flux, mean momentum flux and the mean buoyancy flux. Figure 3 displays the longitudinal dependence of these quantities some time after the step change in the buoyancy flux at the source. To begin, we choose to display the mean buoyancy flux $F$ , rather than the integral buoyancy $B$ , in figure 3. Being equal to the constants $F^{A}$ and $F^{B}$ upstream and downstream of the travelling wave respectively, the behaviour of the mean buoyancy flux is easier to interpret than that of the integral buoyancy. Integrals from each member of the ensemble show significant variation in comparison with the relatively smooth profile that is obtained from their ensemble average. In the mean buoyancy flux $F$ the step change has propagated to approximately $z/r_{s}=30$ . Since $F$ appears as a forcing term in the governing differential equation for the momentum flux $M$ , which is an integral of $\overline{w}^{2}$ rather than $\overline{w}$ , one observes that the behaviour of $Q$ is generally smoother than $F$ and $M$ , making the step change in figure 3(c) comparatively difficult to discern.

Figure 4. DNS data and model prediction of the volume flux $Q$ (a,b); the momentum flux $M$ (c,d); and the integral buoyancy $B$ (e,f). (a,c,e) Display the initial conditions, while (b,d,f) display the data/predictions at times $t_{n}\approx 1.95{\it\tau}_{s}n$ . The thick solid line corresponds to the DNS data. GPM refers to the Gaussian plume model described in § 6 of this paper, while TPM refers to the top-hat plume model described by Scase & Hewitt (Reference Scase and Hewitt2012).

In figure 4 we examine the volume flux $Q$ , the momentum flux $M$ and the integral buoyancy $B$ , which are the dependent variables of the system (2.9)–(2.11). The thick black lines correspond to data obtained from the simulation at a given time. Figure 4 also includes model predictions, which will be discussed in § 6. In each variable, one can observe a wave that travels in the direction of positive $z$ for increasing $t$ . As one would expect, knowing that the plume is driven by buoyancy, the position of the wave is approximately the same in each of the dependent variables. At a given time, upstream and downstream of the wave, the integrals $Q$ , $M$ and $B$ satisfy a quasi-steady balance, and therefore the classical power-law solutions $B\sim z^{1/3}$ , $M\sim z^{4/3}$ and $Q\sim z^{5/3}$ are valid (cf. the dashed lines in figure 3), $F_{s}$ in (2.22) being given by either $F_{s}^{A}$ (upstream) or $F_{s}^{B}$ (downstream) accordingly. In figure 4(b) and (f) for $z/r_{s}>42$ one can discern a small increase in $Q$ and $B$ , respectively. The outflow boundary condition we employ results in an increase in $r_{d}$ just beneath the top of the domain, and therefore a corresponding increase in $Q$ and $B$ (note that, being proportional to $\overline{w}^{2}$ , rather than $\overline{w}$ , the momentum flux $M$ is not affected). In the case of $B$ , it is evident that the effects of the outflow are small in comparison with the amplitude of the travelling wave.

The wave is perhaps most clearly seen in the integral buoyancy $B$ , which typically has a relatively weak power-law dependence on $z$ . To within a constant rescaling factor, the longitudinal dependence of the integrals does not appear to alter significantly with respect to time. For example, the $z$ -dependence of $B(z,t_{10})$ is qualitatively similar to that of $B(z,t_{4})$ , if a suitable rescaling is applied. Although the main focus of this study is on the integrals $Q$ , $M$ and $B$ , we note for future reference that the pressure integral ${\it\beta}_{p}M$ warrants further attention. As one might expect, the dimensionless pressure integral ${\it\beta}_{p}$ increases at the leading edge of the travelling wave and is therefore expected to influence the dynamics of the wave. Further understanding of this aspect of the flow would require a detailed analysis of the lateral components of the turbulence kinetic energy and are beyond the scope of this paper.

Figure 5. (a) Normalised buoyancy flux at the times $t_{n}\approx 1.95{\it\tau}_{s}n$ . The black circle denotes $z^{\ast }(t)$ , which corresponds to the location of the wave. Profiles $t_{2},\ldots ,t_{5}$ appear to be influenced by near-field effects and are therefore excluded from the plot in (b), which illustrates the self-similarity of the normalised buoyancy flux.

Based on the behaviour of the dependent variables evident in figure 4, one wishes to understand the scaling associated with the travelling wave. To this end, we consider the buoyancy flux $F$ in figure 5(a) at several time instances. Being constant upstream and downstream of the wave, the buoyancy flux is a convenient variable to analyse in this regard. We now turn our attention to the scaling of the position of the step change $z^{\ast }(t)$ , leading eventually to the collapsed data from several time instances that is plotted with respect to the similarity variable $z/z^{\ast }$ in figure 5(b).

Using the integrals (2.12), it is useful to define the following characteristic length, velocity and buoyancy scales respectively:

Thus, according to classical plume theory (2.22), when ${\it\theta}_{m}={\it\beta}_{g}=1$ , the characteristic velocity in a steady plume, whose physical source is located at $z=0$ , is

where $z_{v}$ is the location of an asymptotic virtual source. One therefore presumes that, sufficiently far from the finite source, the location $z^{\ast }$ of the propagating wave obeys

where ${\it\lambda}^{\ast }$ is a constant of proportionality and $F^{\ast }$ is a characteristic buoyancy flux such that $F^{B}\leqslant F^{\ast }\leqslant F^{A}$ . We define the characteristic buoyancy flux by imposing buoyancy conservation in the wave and solving the corresponding Rankine–Hugoniot jump conditions, as outlined in appendix C:

As one would expect, (4.4) implies that for infinitesimal changes in $F$ , $F^{\ast }\sim F^{A}\sim F^{B}$ .

One can integrate equation (4.3) to find that

where $t_{v}$ is the time for which $z^{\ast }(t_{v})=z_{v}$ .

We determine $z^{\ast }(t)$ from the DNS data in a simple and reliable manner according to

Consequently, if $F$ is monotonic, according to (4.6) $z^{\ast }$ is single-valued function of time. To ensure a unique definition of $z^{\ast }$ in situations in which $F$ is not monotonic, we take the maximum value of $z^{\ast }$ satisfying (4.6). With the exception of a small region in the near field (see $t_{2}$ and $t_{3}$ in figure 5), we find that (4.6) is sufficient in defining a unique value of $z^{\ast }$ . The position $z^{\ast }(t)$ , determined according to (4.6), is displayed in figure 5(a) with respect to the buoyancy flux at several time instances and evidently provides a useful indication of the wave’s position. In figure 5(b) we rescale the longitudinal coordinate $z$ using the observed front position $z^{\ast }(t)$ . Plotting the mean buoyancy flux $F$ at times $t\geqslant t_{6}$ with respect to $z/z^{\ast }$ , we observe an approximate collapse of the data and therefore self-similarity. Here, self-similarity implies that sufficiently far from the source (or, equivalently, at sufficiently large times) the process reaches a state of ‘equilibrium’ in which $z$ and $t$ cease to play independent roles. We also note that if $z^{\ast }\sim t^{3/4}$ then self-similarity implies that the longitudinal extent (i.e. the spreading rate) of the front also scales according to $t^{3/4}$ . For small times (see $t_{2}$ , $t_{3}$ and $t_{4}$ in figure 5(a)) we observe that the disturbance in the buoyancy flux is oscillatory and that the local peak in the buoyancy flux at approximately $z/r_{s}=9$ breaks down before $t_{4}$ and results in a significant reduction in the steepness of the wave. Consequently, between $t_{3}$ and $t_{6}$ the wave front appears to become steeper as it adjusts to a far-field equilibrium.

Figure 6. Dimensionless buoyancy $b_{m}(z,t)/b_{m0}(z)$ in the plume, where $b_{m0}$ is the steady-state buoyancy corresponding to simulation L with buoyancy flux $F^{B}$ , in addition to points denoting the location of the travelling wave. The crosses correspond to the observed position of the front, while the line denotes a best fit to the front position, of the form $z^{\ast }-z_{v}\propto (t-t_{v})^{3/4}$ .

In figure 6 we test the hypothesis that $z^{\ast }\sim t^{3/4}$ and determine the constant of proportionality ${\it\lambda}^{\ast }$ from the observed location of the front $z^{\ast }(t)$ . Evident from the normalised characteristic buoyancy $b_{m}/b_{m0}$ displayed in figure 6 is that in the far field the wave’s propagation adheres to the predicted $z^{\ast }\sim t^{3/4}$ scaling. In addition, we find that the constant of proportionality ${\it\lambda}^{\ast }\approx 1.9$ in (4.5). In other words, with reference to (4.3), we observe that the front propagates at nearly twice the local top-hat velocity $w_{m}$ . This behaviour is consistent with unsteady turbulent jets (Craske & van Reeuwijk Reference Craske and van Reeuwijk2015a ), and will be explained in § 5 in terms of the dimensionless profile coefficients.

Figure 7. DNS data and model prediction of the plume radius $r_{m}$ (a) and the flux-balance parameter ${\it\Gamma}$ (b) at times $t_{n}\approx 1.95{\it\tau}_{s}n$ . The thick solid line corresponds to the DNS data. GPM refers to the Gaussian plume model described in § 6 of this paper and TPM refers to the top-hat plume model described by Scase & Hewitt (Reference Scase and Hewitt2012). Note that the profiles in (b) are separated by a distance of 1 unit, as indicated by the scale in the bottom left corner.

Compared with $Q$ , $M$ and $B$ , it is perhaps the derived quantities $r_{m}\equiv Q/M^{1/2}$ and ${\it\Gamma}\propto QB/M^{3/2}$ that reveal more about the dynamics of the system, because, unlike $Q$ , $M$ and $B$ , the characteristic plume radius $r_{m}$ and flux-balance parameter ${\it\Gamma}$ are independent of $F$ . Without scrutinising the governing equations, it is therefore difficult to predict whether $r_{m}$ and ${\it\Gamma}$ will increase or decrease in the vicinity of the propagating wave. Figure 7 demonstrates that to leading-order, the behaviour of $r_{m}$ and ${\it\Gamma}$ is practically unaffected by the step change in the buoyancy flux. One is inclined to conclude that $r_{m}$ and ${\it\Gamma}$ are slightly reduced in the vicinity of the wave although the observed change, being not more than 15 % of their steady-state values, is relatively small. We will discuss the physics determining the behaviour of $r_{m}$ and ${\it\Gamma}$ in §§ 5.2–5.3 and demonstrate why they are only weakly sensitive to changes in the plume’s buoyancy flux.

5.1 The hyperbolic system

A logical starting point to understand the leading-order physics associated with the system (2.9)–(2.11) is to analyse its characteristic curves. Without assigning a value to the profile coefficients ${\it\beta}_{g}$ , ${\it\gamma}_{g}$ and ${\it\theta}_{g}$ , which correspond to dimensionless fluxes of momentum, energy and buoyancy, respectively, the unsteady plume equations can be expressed as

where $\boldsymbol{R}$ consists of the right-hand side sink or source terms appearing in (2.9)–(2.11). Characteristic curves of this system satisfy the following relation

where ${\it\lambda}^{\ast }$ is a dimensionless velocity defined by

The relation (5.2) implies that

To date, the assumption of top-hat velocity profiles and the omission of turbulence has resulted in the unsteady plume equations being regarded as a parabolic system (Scase et al. Reference Scase, Aspden and Caulfield2009). However, use of the generalised framework described in § 2 reveals that, in general, the unsteady plume equations comprise a hyperbolic system, even when higher-order turbulent transport terms are neglected from (2.6) and (2.8) (resulting in ${\it\beta}_{g}={\it\beta}_{m}$ and ${\it\gamma}_{g}={\it\gamma}_{m}$ ). Evident from (5.4) is that in the ‘top-hat’ limit, ${\it\gamma}_{g},{\it\beta}_{g},{\it\theta}_{g}\rightarrow 1$ , ${\it\lambda}^{\ast }=1$ , the characteristic curves collapse onto a single family propagating at the local characteristic velocity. In that case, the system cannot be decomposed using linearly independent eigenvectors and should indeed be regarded as parabolic (see Whitham Reference Whitham1974). The top-hat unsteady plume formulation is therefore a degenerate case amongst a wide variety of possible alternatives, each possessing quite different dynamical properties and underlying physics. In fact, top-hat models represent a non-physical limit in the sense that a discontinuous mean velocity field is not realisable in a real turbulent plume. However, we would like to point out that in comparing the radial dependence of different velocity profiles in the light of (5.4) one is not interested in whether they possess compact support per se, but in the extent to which they are non-uniform and/or account for turbulence and pressure. As described in § 2.3, a non-uniform radial dependence of longitudinal velocity always results in a higher dimensionless energy flux ${\it\gamma}_{m}$ , and therefore a separation of characteristic curves, and an unsteady equation for the area of the plume (2.19) that is fundamentally different from that which is associated with top-hat velocity profiles.

A leading-order representation of a Gaussian unsteady plume is obtained by neglecting turbulent transport (hence ${\it\beta}_{g}=1$ , ${\it\theta}_{g}=1$ and ${\it\gamma}_{g}={\it\gamma}_{m}$ ), and using the fact that ${\it\gamma}_{m}=4/3$ , as shown in (2.18). Consequently, (5.4) shows that for Gaussian profiles ${\it\lambda}_{1}=2$ , which implies that the fastest characteristic propagates at twice the local characteristic velocity. This prediction is in reasonably close agreement with the observations reported in § 4 which suggest that ${\it\lambda}^{\ast }\approx 1.9$ . We tentatively attribute the fact that the observed propagation velocity is slightly less than one would expect to find in a Gaussian plume ( ${\it\lambda}^{\ast }\approx 2.0$ ) to pressure, which typically being negative inside the plume (see, e.g. ${\it\beta}_{p}$ and ${\it\gamma}_{p}$ in table 2), leads to a reduction in both ${\it\beta}_{g}$ and  ${\it\gamma}_{g}$ .

5.2 The behaviour of the plume radius and flux-balance parameter

As pointed out in § 4, it is difficult to predict the behaviour of the characteristic radius $r_{m}\equiv Q/M^{1/2}$ and the flux-balance parameter ${\it\Gamma}$ , which are derived quantities, in an unsteady plume a priori. It is therefore useful to understand how their behaviour is affected by the value of the profile coefficients.

Substitution of the relation $Q=r_{m}M^{1/2}$ into (2.20), reveals that

If one assumes that at a particular time the plume is straight sided (conical), such that $\partial _{z}r_{m}=6{\it\alpha}_{0}/5$ in the right-hand side of (5.5) then, using (2.23),

Inspection of (5.6) implies that

When the profile coefficients are such that (a) ${\it\gamma}_{g}/{\it\beta}_{g}=4/3$ and (b) ${\it\theta}_{m}=1$ , the plume remains straight sided for all time, provided that it has a constant source area. Physically, these conditions correspond to (a) the gross dimensionless fluxes of energy and momentum in the plume being in the ratio of $4/3$ and (b) the mean dimensionless flux of buoyancy in the plume being equal to unity. The latter condition is satisfied if the mean buoyancy profile has the same shape and width as the mean velocity profile.

The difference between straight sidedness in jets (see Craske & van Reeuwijk Reference Craske and van Reeuwijk2015b ) and plumes is an additional, independent, condition on ${\it\theta}_{m}$ in the case of the latter. Indeed, ${\it\theta}_{m}$ places further control on the behaviour of $\partial _{t}M$ relative to $\partial _{t}Q$ (see (2.9)–(2.10)) and therefore affects the behaviour of $Q^{2}/M$ . The relations (5.7) are useful because they allow one to relate the ‘internal’ properties of a plume, such as ${\it\gamma}_{g}$ , to an ‘external’ observable such as $r_{m}$ at the integral level. This contrasts with steady-state plumes, for which it is not possible to distinguish between a Gaussian or a top-hat distribution of velocity by observing $Q,M$ and  $F$ .

If we assume that the conditions (5.7) for straight sidedness are satisfied, then, using (2.25),

Substitution of (5.8) into (2.11) and using (2.10), in addition to ${\it\beta}_{g}=3{\it\gamma}_{g}/4$ , ${\it\theta}_{m}=1$ for straight sidedness, one obtains

provided that ${\it\Gamma}=1\,\forall z$ is an initial condition and that ${\it\Gamma}=1$ at the source. Therefore, only if the gross dimensionless buoyancy flux ${\it\theta}_{g}$ (i.e. inclusive of turbulent transport) is equal to the gross dimensionless energy flux ${\it\gamma}_{g}$ will the plume remain pure, such that ${\it\Gamma}=1$ , regardless of the temporal dependence of $M(z,t)$ . Noting that ${\it\theta}_{m}=1$ for straight sidedness and that ${\it\gamma}_{g}\approx 4/3$ , it is only possible that ${\it\theta}_{g}={\it\gamma}_{g}$ if approximately $1/4$ of the total buoyancy flux is transported by turbulence. Table 2 suggests that ${\it\theta}_{g}\approx 1.2$ , while ${\it\gamma}_{g}\approx 1.4$ , and therefore we would not expect for the plume to remain precisely pure. However, in situations for which the prediction of local changes in ${\it\Gamma}$ is not crucial, the assumption of equality between ${\it\gamma}_{g}$ and ${\it\theta}_{g}$ will prove to be a useful idealisation in the context of modelling, which we discuss in § 6.

5.3 The structure of waves in a plume

The way in which the properties of travelling waves in the plume depend on the profile coefficients can be established by examining the behaviour of the invariants associated with (5.1). The left eigenvectors corresponding to the eigenvalues (5.4) are

where ${\it\phi}\equiv \sqrt{1-{\it\beta}_{g}/{\it\gamma}_{g}}$ ,

In general, $L$ is invertible, which means that the unsteady plume system (5.1) can be decomposed into a system of ordinary differential equations. Each differential equation corresponds to the derivative of a quasi-invariant quantity $Y_{1},Y_{2}$ or $Y_{3}$ along a characteristic curve. In particular, introducing the integrating factor $Q^{c_{1}}/M^{c_{2}}$ :

one finds that

The remaining invariants are identical to those found in unsteady jets (Craske & van Reeuwijk Reference Craske and van Reeuwijk2015a ):

When the system is forced due to buoyancy and the production of turbulence kinetic energy, $\boldsymbol{R}\neq 0$ in (5.1) and the ‘invariants’ need not be constant along characteristic curves. It is nevertheless instructive to consider the homogeneous problem, for which $\boldsymbol{R}=0$ , such that $Y_{1}$ , $Y_{2}$ and $Y_{3}$ are constant along characteristics and therefore determine the behaviour of the original dependent variables $Q,M$ and $B$ . The value of $Y_{1}$ , $Y_{2}$ and $Y_{3}$ at a given point in the domain can be determined by tracing their respective characteristic curves to the source at $z=0$ . We consider the case for which ${\it\lambda}_{3}<{\it\lambda}_{2}<{\it\lambda}_{1}$ . We denote with $S$ , $S_{1}$ and $S_{2}$ the regions $[{\it\lambda}_{3},{\it\lambda}_{1}]$ , $[{\it\lambda}_{2},{\it\lambda}_{1}]$ and $[{\it\lambda}_{3},{\it\lambda}_{2}]$ respectively, as indicated in figure 8. Hence

and

Figure 8. The three distinct characteristic curves and associated eigenvalues ${\it\lambda}_{1}$ , ${\it\lambda}_{2}$ and ${\it\lambda}_{3}$ in a generalised formulation of the unsteady plume equations. Regions $A$ and $B$ denote points in $(z,t)$ space upstream and downstream of a wave, while regions $S_{1}$ and $S_{2}$ denote regions of the wave. The depicted behaviour of $F$ in the wave is for schematic purposes only.

The values of $M$ and $Q$ in the wave are not affected by the location of the second characteristic associated with ${\it\lambda}_{2}$ . The buoyancy $B$ , however, takes distinct values in $S_{1}$ and $S_{2}$ according to the position of the second characteristic. If it is assumed that both the source Richardson number and the source radius are fixed, we may say that $Q^{A}\propto F^{A\,1/3}$ , $Q^{B}\propto F^{B\,1/3}$ , $M^{A}\propto F^{A\,2/3}$ and $M^{B}\propto F^{B\,2/3}$ . Hence (5.15) and (5.16) imply that

Therefore, following a step change in the source buoyancy flux $F_{s}^{B}\mapsto F_{s}^{A}$ , the parameter ${\it\phi}$ , which depends on ${\it\gamma}_{g}$ and ${\it\beta}_{g}$ , determines the behaviour of $Q^{S}$ and $M^{S}$ . In particular, for a given $F^{A}$ and $F^{B}$ , ${\it\phi}$ will determine the step change in the plume velocity $w_{m}$ in $S$ and therefore whether each characteristic is associated with rarefaction or compression behaviour. A rarefaction wave occurs when the characteristic velocity $w_{m}$ undergoes a positive step change in the positive $z$ direction, whereas a compression wave occurs when $w_{m}$ undergoes a negative step change. For Gaussian plumes or, more generally, when ${\it\gamma}_{g}/{\it\beta}_{g}=4/3$ (which implies that ${\it\phi}=1/2$ ), the behaviour of the wave is particularly simple, because $Q^{S}=Q^{A}$ and $M^{S}=M^{A}$ . In that case the leading characteristic associated with ${\it\lambda}_{1}$ is a rarefaction wave or a compression wave, according to whether $F^{A}<F^{B}$ or $F^{A}>F^{B}$ , respectively. Furthermore, in agreement with the deductions made in the previous section, when ${\it\phi}=1/2$ the plume radius $r_{m}$ does not deviate from its steady-state value in the region $S$ , although, as noted previously, for this to be the case it is necessary that ${\it\theta}_{m}=1$ . Since these properties are inherited from unsteady jets, the reader is referred to Craske & van Reeuwijk (Reference Craske and van Reeuwijk2015b ) for further details. Here we will focus on the additional properties resulting from buoyancy, as manifest in the characteristic curve associated with  ${\it\lambda}_{2}$ .

Using (5.17) one can determine how the flux-balance parameter ${\it\Gamma}$ behaves in regions $S_{1}$ and $S_{2}$ :

where

In ${\it\Gamma}$ , unlike $Q$ and $M$ , the characteristic associated with ${\it\lambda}_{2}$ , which is determined by the dimensionless buoyancy flux ${\it\theta}_{g}$ , can result in a discontinuity. In general, distinct values of ${\it\Gamma}$ are found in regions $S_{1}$ and $S_{2}$ and depend on ${\it\gamma}_{g},{\it\beta}_{g}$ and ${\it\theta}_{g}$ .

When ${\it\gamma}_{g}/{\it\beta}_{g}=4/3$ the plume is straight sided, ${\it\phi}=1/2$ and

Therefore, substituting for $c_{1}$ and $c_{2}$ using (5.11) and assuming ${\it\Gamma}^{A}={\it\Gamma}^{B}=1$ ,

where

such that the behaviour of ${\it\Gamma}$ over the wave depends only on the ratio of the dimensionless buoyancy flux ${\it\theta}_{g}$ to the dimensionless energy flux ${\it\gamma}_{g}$ . While the plume remains pure in $S_{2}$ (note that ${\it\Gamma}^{S_{2}}=1$ ), ${\it\Gamma}^{S_{1}}$ depends on the location of the second characteristic curve. If ${\it\theta}_{g}/{\it\gamma}_{g}<1$ , there is a deficit of buoyancy in $S_{1}$ and the plume becomes ‘forced’ ( ${\it\Gamma}^{S1}<1$ ). Conversely, if ${\it\theta}_{g}/{\it\gamma}_{g}>1$ surplus buoyancy enters $S_{1}$ and the plume becomes ‘lazy’. These results are consistent with the condition established in § 5.2 that ${\it\theta}_{g}/{\it\gamma}_{g}=1$ in order for the plume to remain pure, and are illustrated in figure 9. For details about the behaviour of ${\it\Gamma}$ when the fluxes in the plume are reduced according to a self-similar scaling, the reader is referred to appendix A.

Figure 9. Wave structure in a straight-sided plume ( ${\it\gamma}_{g}/{\it\beta}_{g}=4/3$ , ${\it\theta}_{m}=1$ ). The special case for which ${\it\Gamma}=1$ in the wave corresponds to ${\it\theta}_{g}={\it\gamma}_{g}$ . The variable ${\it\lambda}$ is a scaled $z$ coordinate, which is defined rigorously in (6.21).

5.4 Response to source perturbations

The final property of the generalised unsteady plume equations that we consider is their response to harmonic perturbations. From a modelling perspective it is necessary to understand whether an unsteady plume model is stable with respect to source perturbations. Scase & Hewitt (Reference Scase and Hewitt2012) showed that top-hat unsteady plume models are ill posed because they admit the unbounded (exponential) growth of high-frequency source perturbations, making it impossible to obtain convergence in numerical simulations of the governing integral equations. Using the framework described in this paper, it is possible to adopt a more general approach and analyse the response of the unsteady plume equations without making an assumption about the velocity profile, turbulent transport or pressure. It is consequently possible to establish whether a given plume model is well posed. More generally, one can establish a relation between the dimensionless profile coefficients that ensures well posedness and understand how this relates to the underlying physics of an unsteady plume.

Using the dimensionless variables

with the characteristic steady-state velocity

we consider a harmonic source perturbation of amplitude ${\it\epsilon}$ and seek solutions of the form

At $O({\it\epsilon})$ the unsteady plume equations can be expressed as

where

In (5.27) and (5.28) we have used (2.23), rearranged as

to eliminate ${\it\delta}_{g}$ .

Here we restrict our attention to a plume with velocity and buoyancy profiles of Gaussian form and equal width, and neglect the longitudinal turbulent transport of momentum. Hence, we assume that ${\it\theta}_{m}={\it\beta}_{g}=1$ and ${\it\gamma}_{g}={\it\gamma}_{m}=4/3$ and focus our attention on the way in which ${\it\theta}_{g}$ affects the downstream development of source perturbations, complementing the predictions pertaining to ${\it\Gamma}$ that were obtained in § 5.3. There are several reasons for focusing on straight-sided unsteady plumes at this point. First, theoretical and observational evidence suggests that if one has to choose a value of ${\it\gamma}_{m}$ , then $4/3$ is the most appropriate choice. Second, straight sidedness is a distinguished case and therefore worthy of investigation in its own right. In the absence of compelling evidence to the contrary, we feel that it would be unhelpful and less general to restrict the parameter space in an alternative way. Third, source perturbation analysis for jets, for which straight sidedness was not assumed, was conducted in Craske & van Reeuwijk (Reference Craske and van Reeuwijk2015b ). Therefore, by focusing on the effects of buoyancy transport in the present section, we augment existing work on the subject.

The assumed values of ${\it\beta}_{g}$ , ${\it\theta}_{m}$ and ${\it\gamma}_{g}$ ensure that the plume remains straight sided (see § 5.2) such that $\mathscr{M}_{1}=2\mathscr{Q}_{1}$ for the linearized perturbation, and thus (5.27) simplifies to

We assume an oscillatory solution of the form

Substitution of (5.31) into (5.30) and elimination of $\hat{\mathscr{B}}_{1}$ results in

Transformation of the coordinate system according to ${\it\zeta}_{\ast }\equiv \text{i}({\it\theta}_{g}-2){\it\zeta}/(2{\it\theta}_{g})$ , and defining the dependent variable $\hat{\mathscr{M}}_{\ast }\equiv \hat{\mathscr{M}}_{1}\exp (\text{i}{\it\zeta}/2)$ , results in

where

Equation (5.33) is a confluent hypergeometric equation (Abramowitz & Stegun Reference Abramowitz and Stegun1970, § 13, p. 504), whose solution can be expressed as

where the independent solutions $\text{M}$ and $\text{U}$ are Kummer functions. At the source, (5.35) implies that $\mathscr{B}_{1}(0)=\mathscr{M}_{1}(0)~\forall c_{1},c_{2}$ , which means that to leading order ${\it\Gamma}(0)=1$ . To ensure finite $\hat{\mathscr{M}}_{1}(0)$ , we set $c_{2}=0$ and assume that $c_{1}=1$ , without loss of generality. The asymptotic expansion of (5.35), for ${\it\zeta}\rightarrow \infty$ and $0\leqslant {\it\theta}_{g}<2$ , is

where $\unicode[STIX]{x1D758}$ is the Gamma function. The exponent $a({\it\theta}_{g})$ therefore determines the growth rate of the perturbations, as illustrated in figure 10. When $1\leqslant {\it\theta}_{g}<4/3$ (forced waves, according to § 5.3), the growth rate is negative, but when $4/3<{\it\theta}_{g}<2$ (lazy waves), the growth rate is positive. Between growth and decay of source perturbations lies the special behaviour associated with ${\it\theta}_{g}=4/3$ , for which the plume is pure and exhibits a neutral response to source perturbations. We note that ${\it\theta}_{g}\equiv {\it\theta}_{m}+{\it\theta}_{f}$ is not likely to exceed $4/3$ in practice (see e.g. table 2) because ${\it\theta}_{m}\approx 1$ and the dimensionless turbulent buoyancy flux ${\it\theta}_{f}\approx 0.2$ . Therefore, this analysis demonstrates that the generalised unsteady plume equations are well posed for physically realistic values of the dimensionless profile coefficients.

Figure 10. Stability of the system when ${\it\gamma}_{m}=4/3$ in response to source perturbations with respect to the dimensionless longitudinal coordinate ${\it\zeta}\propto z^{4/3}{\it\sigma}$ for (a)  ${\it\theta}_{g}=1$ ; (b)  ${\it\theta}_{g}=5/4$ ; (c)  ${\it\theta}_{g}=4/3$ and (d)  ${\it\theta}_{g}=5/3$ . The thin solid lines correspond to the exact solution of the linearized problem (5.35) and the thin dashed lines to the asymptotic solution (5.36). The thick lines denote the envelope of the asymptotic solution. The dimensionless buoyancy flux ${\it\theta}_{g}=4/3$ is the special case for which perturbations exhibit neutral growth. For models employing realistic values ${\it\theta}_{g}\leqslant 4/3$ , the system is well posed in the sense that source perturbations are bounded and it is possible to obtain convergent numerical approximations.

6.1 Shear-flow dispersion

In the previous section we established the response of a plume when the profile coefficients are constant. In that case the plume integrals are, in general, discontinuous along characteristic curves, as illustrated schematically in figure 8. However, it is readily apparent from figure 4 that the plume integrals are continuous, varying smoothly from their value upstream of the disturbance to their value downstream of the disturbance. The longitudinal gradients that are produced by unsteady changes in the plume are capable of inducing a local departure from self-similarity. Noting that constant profile coefficients are a necessary condition for self-similarity, one would expect the profile coefficients to deviate from their steady-state value in the vicinity of a sudden change in the plume’s integral quantities. Taylor (Reference Taylor1953) analysed an equivalent situation in pipe flow, showing that in the vicinity of a step increase in the radially averaged scalar concentration there exists a deformation in the otherwise radially uniform concentration field. In correlation with a non-uniform velocity field, the deformation results in a local increase in the scalar flux, which, to leading order, satisfies a dispersive analogue of Fick’s law. Similarly, in Craske & van Reeuwijk (Reference Craske and van Reeuwijk2015b ), a dispersion closure was developed for energy transport in jets, which we will apply here to model the higher-order transport of buoyancy and energy in unsteady plumes.

In order to focus on longitudinal changes that constitute a departure from a steady state, it is convenient to examine dimensionless quantities. For a given quantity $f_{m}(z,t)$ , we define $\mathscr{F}\equiv f_{m}/f_{m0}$ , where $f_{m0}$ is the steady-state value of $f_{m}$ . To quantify the departure from the steady state, we can compute the longitudinal gradient of  $\mathscr{F}$ :

Noting that $r_{m}=6{\it\alpha}_{0}z/5$ in the steady state, it is convenient to define the steady-state function $f_{m0}$ locally:

where $f_{m0}\propto z^{n}$ , which means that the second term on the right-hand side of (6.1) is equal to $-6{\it\alpha}_{0}n/5$ . This is a sensible choice in general, because it means that departures from a steady state are defined relative to the local (in  $z$ ) size of the plume rather than the size it should have with respect to a fixed virtual source. Of course, if the plume is straight sided, such that $r_{m}\propto z~\forall t$ , the definitions coincide. In particular, we define the dimensionless velocity and buoyancy:

where $w_{m0}\propto z^{-1/3}$ and $b_{m0}\propto z^{-5/3}$ are the steady-state power-law solutions for the velocity and buoyancy respectively. Using (6.1), with $\mathscr{F}$ replaced with either $\mathscr{U}$ or $\mathscr{B}$ , we follow Craske & van Reeuwijk (Reference Craske and van Reeuwijk2015b ) and propose a closure for the dimensionless energy flux ${\it\gamma}_{m}$ and the dimensionless buoyancy flux ${\it\theta}_{m}$ :

Note that the fractions $1/3$ and $5/3$ in (6.4) and (6.5), respectively, result from the fact that $w_{m0}\sim z^{-1/3}$ and $b_{m0}\sim z^{-5/3}$ , respectively. Although the terms inside the parenthesis in the definition of ${\it\gamma}_{m}$ appear to differ from those presented in Craske & van Reeuwijk (Reference Craske and van Reeuwijk2015b ), the difference is superficial, as ${\it\gamma}_{1}$ can be redefined to render the expressions equivalent. A notable virtue of the dispersion closure is that the incorporation of (6.5) does not affect the steady-state solutions. In contrast, the diffusive mixing term for momentum proposed by Scase & Hewitt (Reference Scase and Hewitt2012) was not extended to buoyancy transport, because it would have modified the exponents associated with the steady-state solutions.

6.2 Unsteady plume model

For the purposes of establishing a relatively simple representation of an unsteady plume, we neglect turbulent transport and assume ${\it\gamma}_{g}={\it\gamma}_{m}$ , ${\it\theta}_{g}={\it\theta}_{m}$ and ${\it\beta}_{g}={\it\beta}_{m}$ , resulting in the following model

In the interest of brevity, ${\it\theta}_{m}$ on the right-hand side of (6.7) has not been expanded and for consistency should therefore be replaced with (6.5).

We solve the system (6.6)–(6.8) using fourth-order accurate central differences and a fourth-order Runge–Kutta method for integration with respect to  $t$ . To adequately describe both near- and far-field scales in the plume we employ a stretched computational grid in space, for which the ratio of the first and last computational cells is approximately equal to 0.1. At the first and last computational cells the values of the dependent variables are held constant (i.e. independent of time), and by letting the height of the domain be equal to $64r_{s}$ we ensure that the outflow condition has a negligible influence in those parts of the domain with which we are primarily concerned. As indicated in figure 4(a,c,e), smooth initial conditions for $Q,M$ and $B$ were obtained by fitting an error function to the observed buoyancy flux $F$ at $t=t_{1}$ and by substituting $F(z,t_{1})$ into the steady-state solutions (2.22). Although this method of determining initial conditions is approximate, at $t=t_{1}$ the unsteady part of the flow is very close to the source and a precise description of the rapid change in $F$ from $F^{B}$ to $F^{A}$ makes little difference to predictions at later times. For $t>t_{1}$ , $Q,M$ and $B$ are held constant at the source.

We compare our model’s predictions, which we refer to as the Gaussian plume model (GPM), to those of the adjusted top-hat plume model (TPM) proposed by Scase & Hewitt (Reference Scase and Hewitt2012). For the parameters of GPM we assume Gaussian profiles and neglect turbulence transport, setting ${\it\gamma}_{0}=4/3$ and ${\it\theta}_{0}=1$ . For the higher-order mixing terms we use ${\it\gamma}_{1}={\it\theta}_{1}=0.0017$ . For the parameter ${\it\varepsilon}$ that appears in the regularisation proposed by Scase & Hewitt (Reference Scase and Hewitt2012) we set ${\it\varepsilon}=5{\it\gamma}_{1}$ , which ensures that the diffusive flux of energy in TPM is equal to the dispersive flux of energy in GPM (see Craske & van Reeuwijk Reference Craske and van Reeuwijk2015b ). To obtain predictions from TPM it was necessary to employ a flux-limiting scheme, and to this end we followed Kurganov & Tadmor (Reference Kurganov and Tadmor2000). For TPM it was also necessary to use 4000 computational points to obtain a converged prediction, which was verified by halving the number of points to 2000. In contrast, a converged prediction was obtained from GPM using only 400 computational points, and therefore took significantly less time to compute than TPM.

Figures 4 and 7 compare the predictions of (6.6)–(6.8) to those of the regularised top-hat model described by Scase & Hewitt (Reference Scase and Hewitt2012). The predictions of $Q$ and $M$ using GPM displayed in figure 4 are in reasonably good agreement with the DNS data. In particular, GPM correctly represents the position and the spreading rate of the wave front in $M$ . In the volume flux $Q$ it appears that GPM predicts that the wave front travels slightly faster than the DNS observations suggest. Nevertheless, despite using an order of magnitude fewer computational points, GPM provides a better overall agreement with the DNS observations than TPM, which predicts a large overshoot in $Q,M$ and  $B$ .

In figure 7 the predictions from GPM exhibit a good agreement with both the radius of the plume and the flux-balance parameter. TPM predicts a local increase in both $r_{m}$ and ${\it\Gamma}$ , which we do not observe in the DNS simulations. Consistent with the analysis of the hyperbolic problem in § 5.1, for ${\it\theta}_{g}<4/3$ , GPM predicts a small reduction in ${\it\Gamma}$ in the vicinity of the wave. Due to statistical uncertainty, it is not possible to say whether this prediction of a very small reduction in ${\it\Gamma}$ agrees with the DNS observations, which to leading order suggest that ${\it\Gamma}$ is approximately constant.

Figure 11. Wave structure for different values of ${\it\theta}_{g}$ . The dashed lines are normalised integral quantities $\mathscr{Q}\equiv Q/Q_{0}$ , $\mathscr{M}\equiv M/M_{0}$ and the solid line is $\mathscr{B}\equiv B/B_{0}$ . The shaded region denotes the extent to which ${\it\Gamma}(z,t)$ is different from unity. For ${\it\theta}_{g}\approx 4/3$ , ${\it\Gamma}=1$ and the wave is pure.

There is a noticeable difference between the GPM predictions and the DNS observations of the integral buoyancy $B$ in the vicinity of the wave front (see figure 4(f)). The observations suggest that the propagating front in $B$ is steeper and faster than that which is predicted by GPM. While the purpose of this paper is not to investigate the fine tuning of the various profile coefficients to give an optimised agreement with the DNS data, for physical insight and to relate our observations to the theory discussed in § 5.3, it is useful to understand the way in which unsteady plume models depend on the assumed amount of turbulent buoyancy flux. Figure 11 displays the quantities $Q(z,t)$ , $M(z,t)$ and $B(z,t)$ normalised by the scalings associated with the steady state, namely $Q_{0}(z)\sim z^{5/3}$ , $M_{0}(z)\sim z^{4/3}$ and $B_{0}(z)\sim z^{1/3}$ , at a fixed time. In addition, the shaded region in figure 11 indicates the extent to which the flux-balance parameter ${\it\Gamma}$ (2.25) deviates from unity. Each of the predictions displayed in figure 11 corresponds to a different choice of the dimensionless parameter ${\it\theta}_{g}\sim {\it\theta}_{0}+{\it\theta}_{f}$ , and therefore, since ${\it\theta}_{0}=0$ is fixed, to a different amount of turbulent buoyancy transport ${\it\theta}_{f}$ in the plume. More precisely, we add a constant ${\it\theta}_{f}\in [0.5,2.0]$ to ${\it\theta}_{0}$ in (6.8) to account for the turbulent transport of buoyancy. Values of ${\it\theta}_{g}$ much less or greater than unity, although physically unrealistic, have been included in figure 11 to provide a more complete picture of the system’s response.

The wave in the dimensionless buoyancy integral $\mathscr{B}$ in figure 11 appears to propagate faster for large values of ${\it\theta}_{g}$ than for small values of ${\it\theta}_{g}$ , due to the distribution of buoyancy within the wave. In agreement with the analysis of § 5.3, it is also clear that the wave is lazy when ${\it\theta}_{g}>4/3$ and forced when ${\it\theta}_{g}<4/3$ . As ${\it\theta}_{g}$ increases, one observes that the front in the integral buoyancy $B$ becomes steeper, in spite of the fact that the longitudinal mixing (dispersion) parameter ${\it\theta}_{1}$ is identical in each case. Hence, the reason for the disparity between the GPM prediction of $B$ and the DNS results in figure 4 is not the higher-order term on the right-hand side of (6.8), but that by setting ${\it\theta}_{0}=1$ we neglected the turbulent buoyancy flux. As indicated in figure 11, when the turbulent transport of buoyancy is included (e.g. ${\it\theta}_{g}\sim {\it\theta}_{0}=1.2$ ), the wave front in $B$ predicted by GPM becomes steeper and exhibits a better agreement with the features of $B$ observed in figure 4.

6.3 Pure straight-sided unsteady plume model

Motivated by the observations reported in § 4 and the theory developed in § 5, the model we derive in this section is based on two main assumptions; (1) the plume is straight sided $(r_{m}=6{\it\alpha}_{0}z/5~\forall t)$ , hence ${\it\gamma}_{g}/{\it\beta}_{g}=4/3$ and ${\it\theta}_{m}=1$ (see § 5.2); (2) the plume remains pure $({\it\Gamma}=1~\forall t)$ , hence ${\it\theta}_{g}={\it\gamma}_{g}$ (see § 5.2). For definiteness, we will neglect the turbulent transport of momentum and set ${\it\beta}_{g}=1$ , which implies that ${\it\gamma}_{g}=4/3$ . Under these conditions the structure of the wave described in § 5.3 is particularly simple because it consists of a single front propagating along the fastest characteristic curve. Behind the front the plume behaves in accordance with a steady state. Noting from (5.4) that when ${\it\gamma}_{g}=4/3$ and ${\it\beta}_{g}=1$ the dimensionless velocity of the fastest characteristic is ${\it\lambda}^{\ast }={\it\lambda}_{1}=3{\it\gamma}_{g}/2$ , we find that ${\it\theta}_{g}={\it\gamma}_{g}=2{\it\lambda}^{\ast }/3$ for consistency.

In the light of our observations from steady plumes (see table 2), the validity of assuming that ${\it\theta}_{g}={\it\gamma}_{g}$ is questionable. Thus, for situations in which the prediction of changes in ${\it\Gamma}$ in an unsteady plume are crucial, we recommend use of the full Gaussian model described in § 6.1. However, the model we describe below is valuable in its own right as a simple analytical means of obtaining a first approximation to the behaviour of an unsteady plume and for providing physical intuition. It is also of theoretical value as a distinguished case, illustrating that under the aforementioned assumptions the behaviour of an unsteady plume is similar to that of an unsteady jet (cf. Craske & van Reeuwijk Reference Craske and van Reeuwijk2015b ).

We start by considering the transport equation for buoyancy:

Since we are assuming that the plume remains pure and straight sided,

hence

and thus

To include longitudinal spreading of the wave due to dispersion, we employ (6.5), which in this case reduces to

Note that ${\it\theta}_{1}\ll {\it\theta}_{0}$ , implying that, for smoothly varying $B$ , the effect of ${\it\theta}_{1}$ on ${\it\theta}_{g}$ is extremely small and therefore does not significantly violate our assumption that ${\it\theta}_{g}=2{\it\lambda}^{\ast }/3$ , which is required for pure-plume behaviour. The transport equation for the buoyancy integral becomes

For a given solution of (6.14) one can derive the volume flux and momentum flux by inverting (6.10):

In order to use (6.14) in practice, one needs to specify the dispersion parameter ${\it\theta}_{1}$ , in addition to the steady-state entrainment coefficient ${\it\alpha}_{0}$ . The remaining parameter ${\it\lambda}^{\ast }=3{\it\gamma}_{g}/2=2$ is determined from the assumption that ${\it\gamma}_{g}=4/3$ .

6.4 Linearized similarity solution

Guided by the observations of § 4.2 (in particular, see figure 5), one expects unsteady disturbances in a plume to evolve in self-similar fashion sufficiently far from the source; hence we seek a similarity solution to (6.14). To simplify the analysis we will assume that the imposed change in the buoyancy flux is relatively small and consider a perturbation expansion for $B$ , where the small parameter ${\it\epsilon}$ depends on the magnitude of the imposed step change in the source buoyancy flux:

where

and $F^{\ast }$ was defined in (4.4). To first order, the conditions of straight sidedness and constant ${\it\Gamma}$ can be expressed as

At $O({\it\epsilon})$ we find that

Recalling that

along the leading characteristic curves of the original hyperbolic system, it is natural to define a similarity variable

Figure 12. Predictions using a linearized, straight-sided, constant- ${\it\Gamma}$ model, for which ${\it\lambda}^{\ast }=2$ . The circles denote the prediction that is obtained by solving (6.22) numerically and the solid black line denotes the solution for $\mathit{Pe}\gg 1$ . The grey lines correspond to observed values of the normalised buoyancy integral from the DNS at times in the interval $[t_{6},t_{12}]$ (cf. figure 5).

such that ${\it\lambda}(t,z^{\ast })={\it\lambda}^{\ast }$ . Assuming that the process is indeed self-similar, (6.19) can be significantly simplified. Using (6.21),

where ${\it\lambda}^{\ast }$ is the dimensionless velocity of the front and $\mathit{Pe}\equiv {\it\lambda}^{\ast }/(4{\it\theta}_{1})=300$ for ${\it\theta}_{1}=0.0017$ and ${\it\lambda}^{\ast }=2$ . Equation (6.22) has an analytical solution that can be expressed in terms of a hypergeometric function. In practice however ${\it\theta}_{1}\ll 1$ , therefore $\mathit{Pe}\gg 1$ , and it is appropriate to use the asymptotic form of the solution to (6.22):

where G is the hypergeometric function and $\mathscr{B}_{1}^{B}$ and $\mathscr{B}_{1}^{A}$ are the values of $\mathscr{B}_{1}$ before and after the step change, respectively. Figure 12 displays the solution to (6.22) in addition to the asymptotic solution (6.23). Evident from figure 12 is that the similarity solution (6.23) exhibits a reasonably good agreement with both the observed behaviour of $\mathscr{B}_{1}$ and the full differential equation (6.22). The actual front propagates slightly slower than our predictions using ${\it\lambda}^{\ast }=2$ would suggest, which is consistent with the observation that ${\it\lambda}^{\ast }\approx 1.9$ , as reported in § 4.2. Note that the observed behaviour of the integral buoyancy is approximately self-similar, albeit over the limited sample that we were able to provide in the far field. The similarity scaling employed to derive (6.22), and the observation of self-similarity in figure 12, support the view that the length scales in the plume vary according to $t^{3/4}$ . This scaling applies to both the position of a disturbance and its longitudinal extent and therefore provides some resolution of the open question debated in Scase et al. (Reference Scase, Aspden and Caulfield2009), as to the longitudinal scaling of a propagating pulse structure in a plume. However, we note that our observations are based on results from a domain of relatively limited longitudinal extent.