3.1. Centreline statistics

As has been established over the course of many decades, once the plume is sufficiently far from the source, we may deduce asymptotic scaling relationships of centreline statistics of density deficiency, $1-\rho _c/\rho _a=\xi _c \sim z^{-5/3}$ and vertical velocity, $u_c\sim z^{-1/3}$ (Morton, Taylor & Turner Reference Morton, Taylor and Turner1956; Hunt & van den Bremer Reference Hunt and van den Bremer2010). The mechanical pressure is expected to scale with turbulent velocity fluctuations once the flow is self-similar, but the scaling coefficient is expected to be small due to the slenderness of the plume (Shabbir & George Reference Shabbir and George1994), leading to a near-zero value of $\textrm{d}{p}_{c}/\textrm{d}{z}$ far enough downstream. Ensuring the density deficiency and vertical velocity scalings are observed for the present simulations provides further validation of the algorithms used in PeleLMeX for buoyancy-driven flows (Wimer et al. Reference Wimer, Day, Lapointe, Meehan, Makowiecki, Glusman, Daily, Rieker and Hamlington2020; Meehan et al. Reference Meehan, Wimer and Hamlington2022).

Figure 2.

Centreline averages of the (a) density deficiency, $\xi _c$ ; (b) vertical velocity, $u_c$ ; and (c) vertical derivative of the mechanical pressure, $\textrm{d}{p_c}/\textrm{d}{z}$ . In (a) and (b), we include self-similar scaling laws expected for far-field behaviour, and in (c), we include a darkened line at $\textrm{d}{p}/\textrm{d}{z} = 0$ .

We plot the centreline density deficiency, $\xi _c$ , vertical velocity, $u_c$ , and vertical derivative of the mechanical pressure, $\textrm{d}{p_c}/\textrm{d}{z}$ in figure 2. After the helium enters the domain, there is substantial mixing with the ambient fluid, and near $z/D \approx 7$ , the flow transitions to the expected self-similar power-law scalings. Very little research has been done to estimate the parameters dictating this transition, but this downstream location is typical of experiments (Papanicolaou & List Reference Papanicolaou and List1988; Shabbir & George Reference Shabbir and George1994) and simulations (van Reeuwijk et al. Reference van Reeuwijk, Salizzoni, Hunt and Craske2016). Additionally, there is very little difference between the statistics of S1 and S2, but S2 has slightly increased values of $u_c$ because the increased entrainment at the base is pushing more air into the plume.

Upstream of the self-similar power-law scalings, the statistics indicate two additional regimes, which we will classify as the near- and mid-field. The concept of the near- and mid-field was originally proposed by Meehan & Hamlington (Reference Meehan and Hamlington2023), but the authors did not have a sufficiently large domain to provide complete evidence of the three regimes. In the near field, helium is rapidly accelerated in the form of the puffing frequency (see figure 1) and mixes with the air, and as a result, (i) density deficiency decreases, (ii) vertical velocity increases and (iii) a non-constant pressure gradient is observed . Near  $z/D \approx 1$ , the coherent puffing instability transitions into turbulent mixing, but does not immediately transition to self-similar power-law scalings. For $z/D \approx 1-7$ , the scalings of $u_c$ and $\xi _c$ are shallower than far-field scalings, particularly $\xi _c$ . We also see the plume becoming more slender, leading to $\textrm{d}{p_c}/\textrm{d}{z}$ decreasing in this region. We consider this region the mid-field, and its existence may have implications for integral models. It is likely that the existence of the mid-field is largely due to the highly buoyant inflow conditions since it is not generally observed in buoyant jets (Khan & Rao Reference Khan and Rao2023) but may be seen under more lazy conditions (van Reeuwijk et al. Reference van Reeuwijk, Salizzoni, Hunt and Craske2016). Future research is needed to establish the parameters that dictate the mid-field characteristics.

Additional single point statistics, including higher-order fluctuating moments and different source conditions, are available in Meehan & Hamlington (Reference Meehan and Hamlington2023). These are not discussed here since the focus is on statistics pertaining to integral-scale modelling.

3.2. Vertical fluxes

Vertical fluxes are often used as the basis for deriving the plume equations because they are more representative of the global characteristics of the flow. For a non-Boussinesq plume, it is common to use three fluxes: the buoyant flux, the mass flux and momentum flux (Carlotti & Hunt Reference Carlotti and Hunt2017). These are respectively defined as

(3.1) \begin{align} B(z) &= \int _0^\infty \left \langle g \left (1-\frac {\rho }{\rho _a}\right ) u_z \right \rangle 2\pi r\textrm{d}{r} , \end{align}

(3.2) \begin{align} G(z) &= \int _0^\infty \left \langle \rho u_z \right \rangle 2\pi r\textrm{d}{r} , \end{align}

(3.3) \begin{align} M(z) &= \int _0^\infty \left \langle \rho u_z^2 \right \rangle 2\pi r\textrm{d}{r} . \end{align}

The operator $\langle \boldsymbol{\cdot }\rangle$ represents a standard Reynolds average, so for any flow variable, $q = \langle q\rangle + q'$ , where $q'$ is the fluctuating variable. We will interchange the $\langle \boldsymbol{\cdot }\rangle$ with $\overline {\boldsymbol{\cdot }}$ when it is more convenient. Note that the infinite upper bound on the integral cannot be implemented exactly in the computations because the domain is finite; however, the outer limit of the domain is sufficient because $\rho$ and $u_z$ approach $\rho _a$ and 0, respectively, exponentially fast outside the plume core.

Most integral models of plumes require some assumptions on the relative contribution of the turbulence statistics compared with the mean statistics. van Reeuwijk et al. (Reference van Reeuwijk, Salizzoni, Hunt and Craske2016) provide a nice overview of this in the context of plumes. To perform this comparison, we decompose $B$ , $G$ and $M$ defined in (3.1), (3.2) and (3.3) into mean and turbulent components as

(3.4) \begin{align} {B}(z) &= \underbrace {\int _0^\infty g \left (1-\frac {\bar {\rho }}{\rho _a} \right ) \bar {u}_z 2\pi r\textrm{d}{r}}_{\overline {B}} \underbrace {-\int _0^\infty \frac {g}{\rho _a} \overline {\rho ' u_z'} 2\pi r\textrm{d}{r}}_{B'}, \end{align}

(3.5) \begin{align} {G}(z) &= \underbrace {\int _0^\infty \bar {\rho } \bar {u}_z 2\pi r\textrm{d}{r}}_{\overline {G}} +\underbrace {\int _0^\infty \overline {\rho ' u_z'} 2\pi r\textrm{d}{r}}_{G'}, \end{align}

(3.6) \begin{align} {M}(z) &= \underbrace {\int _0^\infty \bar {\rho } \bar {u}_z^2 2\pi r\textrm{d}{r}}_{\overline {M}} + \underbrace {\int _0^\infty \bar {\rho } \overline {u_z'^2} 2\pi r\textrm{d}{r}}_{M_1'} + \underbrace {\int _0^\infty \bar {u}_z \overline {\rho ' u_z'} 4\pi r\textrm{d}{r}}_{M_2'} + \underbrace {\int _0^\infty \overline {\rho ' u_z'^2} 2\pi r\textrm{d}{r}}_{M_3'}. \end{align}

Each component is a function of $z$ . For simplicity, we include the negative preceding the second integral in (3.4) into the definition of $B'$ .

An alternative approach to the decomposition shown in (3.4), (3.5) and (3.6) is to use a Favre decomposition as opposed to a Reynolds decomposition of the flow variables. From a modelling perspective, density weighting using Favre averages is appealing because it reduces the number of terms that need to be closed (Baltzer & Livescu Reference Baltzer and Livescu2020; Salizzoni et al. Reference Salizzoni, Vaux, Creyssels, Craske and van Reeuwijk2024). However, we elect to use Reynolds averages for two reasons. First, while Favre averages are mathematically convenient, they are more difficult to interpret since, in general, $\tilde {u}_i \neq \bar {u}_i$ , where $\tilde {u}_i = \overline {\rho u}_i/\bar {\rho }$ is a Favre-averaged velocity. A significant portion of this work is dedicated to analysis, making interpretability a priority. This is also consistent with what is most likely available from experiments since it is difficult to make simultaneous measurements of density and velocity. Second, while Reynolds averaging introduces additional terms, they provide more physical insight, as will be seen later in the text. Nonetheless, we plan to utilise Favre averages in future work when we model a broader physical parameter space.

We plot the fluxes in figure 3 for (a) buoyancy, $B$ ; (b) mass, $G$ ; and (c) momentum, $M$ , according to (3.1), (3.2) and (3.3), respectively, and the decompositions into average and turbulent fluxes as defined in (3.4), (3.5) and (3.6). We only show results from S1 since S2 shows very similar results. From this figure, we can make a number of important observations.

Figure 3.

Plots of the (a) buoyant fluxes, (b) mass fluxes and (c) momentum fluxes as a function of vertical distance $z$ for S1. Each of the total streamwise fluxes are decomposed into purely mean and turbulent contributions, respectively indicated with $\bar {\boldsymbol{\cdot }}$ and $\boldsymbol{\cdot }'$ . Definitions of each flux may be found in (3.4), (3.5) and (3.6). Note that, by definition, $G'$ and $B'$ are related by $-(\rho _a/g) \approx -15.1$ , making $G'$ a small negative number whose shape can be deduced from $B'$ in (a).

The total buoyant flux $B$ remains constant as a function of downstream distance since there is no mechanism to generate buoyancy except when helium is introduced at $z=0$ . Small spatial variations in $B$ are likely a result of the numerical procedure (e.g. spin-up time, averaging window, tolerancing, etc.) but are not expected to affect the physics. Interestingly, when $B$ is decomposed into $\overline {B}$ and $B'$ , the balance between the two varies with $z$ . Specifically, the puffing instability located at the base of the plume increases $B'$ rapidly due to the large coherent fluctuations. Once the flow has substantially mixed near $z/D\approx 5$ , $\textrm{d}{\overline {B}}/\textrm{d}{z} \approx 0$ and $\textrm{d}{B'}/\textrm{d}{z} \approx 0$ with $\overline {B}/B \approx 0.82-0.83$ . This finding needs to be taken into consideration when developing integrals models: using $B$ computed using inlet conditions is likely larger than the commonly used $\overline {B}$ .

The total mass flux $G$ is found to be almost identical to $\overline {G}$ . To demonstrate why $|G'| \ll |G|$ , a control volume analysis on the continuity equation yields $\textrm{d}{G}/\textrm{d}{z} = \rho _a f(z)$ , where $f(z)$ is a generic function. Since $\textrm{d}{B}/\textrm{d}{z} = 0$ , then $\textrm{d}{Q}/\textrm{d}{z} = f(z)$ , where $Q = \int _0^\infty \bar {u}_z 2 \pi r\, \textrm{d}{r}$ . Due to the presence of $\bar {u}_z$ in the integrand, we would expect $\overline {G}$ to be much closer to $\rho _a\, \textrm{d}{Q}/\textrm{d}{z}$ than $G'$ . This implies $\textrm{d}{G}/\textrm{d}{z} -\, \textrm{d}{\overline {G}}/\textrm{d}{z} \approx 0$ and $G'$ is roughly a constant. Since $G$ increases rapidly due to entrainment, $G'$ quickly becomes irrelevant. Note that $G'$ is not zero, but instead a small negative value since it differs by a constant negative factor with $B'$ , $G' = -(\rho _a/g) B' \approx -0.003$ . Despite the flow being non-Boussinesq, this finding from the simulations is consistent with Boussinesq plumes where $G'$ is exactly zero.

Similar to the buoyant fluxes, the mean momentum flux $\overline {M}$ makes up a significant amount of the total momentum flux $M$ , but $M_1'$ , the turbulent flux with velocity fluctuations and mean density in the integrand, acts as a secondary contributor. Generally, we find that $\overline {M}/M \approx 0.80$ for all $z$ , which is a slight decrease than the commonly observed value of 0.9 (van Reeuwijk et al. Reference van Reeuwijk, Salizzoni, Hunt and Craske2016). It is generally assumed that the fluctuations account for a small proportion of the total, allowing them to be neglected in most integral models.

3.3. Plume Richardson number

The source conditions lead to three distinct regimes of plumes: forced, pure and lazy plumes. These regimes can be distinguished by the plume Richardson number, defined as (Hunt & Kaye Reference Hunt and Kaye2001; Carlotti & Hunt Reference Carlotti and Hunt2005; Hunt & van den Bremer Reference Hunt and van den Bremer2010; Candelier & Vauquelin Reference Candelier and Vauquelin2012)

(3.7) \begin{equation} \varGamma = \rho _a^{1/2} \frac {B G^2}{M^{5/2}} . \end{equation}

The classic power-law solution of the plume equations corresponds to the pure-plume regime when $\varGamma$ equals $\varGamma _0$ , a constant value provided in (3.8). When $\varGamma \lt \varGamma _0$ , the plume has an excess of momentum and is said to be forced, and when $\varGamma \gt \varGamma _0$ , the plume has a deficit of momentum and is said to be lazy. For our source conditions of $B_s = 0.23$ m $^4$ s⁻ $^3$ , $G_s = 2.5\times 10^{-3}$ kg s⁻¹ and $M_s=7.6\times 10^{-4}$ kgm s⁻², this leads to $\varGamma _s \approx 95.2$ , placing us firmly in the category of a lazy plume.

Figure 4.

Plot of plume Richardson number as defined in (3.7) as a function of height. The left inset shows the same plot on a log–log scale and the right inset zooms into the far field to identify the asymptotic $\varGamma _0$ .

We plot $\varGamma$ as a function of downstream distance for S1 and S2 in figure 4 to better understand the spatial variability of $\varGamma$ . The two insets are the same lines but focus on different aspects: the left inset uses a log–log scale to better understand near-field variability while the right inset zooms-in on the asymptotic far-field value that Gamma approaches. For both S1 and S2, large values of $\varGamma$ are very short lived, with $\varGamma \gt 10$ only being present below $z/D = 0.1$ . The development of the boundary layer in S1 leads to slightly reduced values of $\varGamma$ compared with S2 since the boundary layer resists entrainment. Beyond the very near field, there is little difference between S1 and S2. Near  $z/D\approx 10$ , S1 achieves an asymptotic value of $\varGamma _0 \approx 0.30$ . This distance is consistent with the infinitely lazy plumes of van Reeuwijk et al. (Reference van Reeuwijk, Salizzoni, Hunt and Craske2016). In classical dynamical plume models (see § 4, Carlotti & Hunt (Reference Carlotti and Hunt2005), or Hunt & van den Bremer (Reference Hunt and van den Bremer2010) for examples), the asymptotic value of the plume Richardson number can be linked with the constant entrainment coefficient for a pure plume

(3.8) \begin{equation} \varGamma _0 = \frac {8 \sqrt {\pi }}{5} \alpha _0. \end{equation}

For the present simulation, the asymptotic value of $\varGamma _0$ leads to $\alpha _0 = 0.106$ , which is within the accepted range of $\alpha _0$ .

In accordance with theory (Hunt & Kaye Reference Hunt and Kaye2005; Hunt & van den Bremer Reference Hunt and van den Bremer2010; Candelier & Vauquelin Reference Candelier and Vauquelin2012) and experiments (Ezzamel, Salizzoni & Hunt Reference Ezzamel, Salizzoni and Hunt2015), we observe a monotonic decrease in $\varGamma$ as a function of height. This contrasts the simulations from van Reeuwijk et al. (Reference van Reeuwijk, Salizzoni, Hunt and Craske2016) and Marjanovic et al. (Reference Marjanovic, Taub and Balachandar2017) where some non-monotonicity is observed. A few reasons may have led to the difference including different source conditions, the Boussinesq approximation, spanwise dimensions of the computational domain and sufficient grid resolution near the source. Additional work is needed to identify how each of these factors affects $\varGamma$ .

3.4. Evolution of profiles with height

The centreline statistics in § 3.1 indicate there is a substantial region of the flow that does not exhibit power-law evolution expected for a self-similar pure plume. To understand the transition to self-similarity, we plot radial profiles of the average density deficiency,  $\overline {\xi }$ , and average vertical velocity, $\overline {u}_z$ , for various downstream locations in figure 5. Both quantities are normalised by their respective centreline values, and the radial coordinate is normalised by the plume radius as defined via the thresholding technique provided in § 3.5. In the self-similar region, all profiles should collapse to a self-similar profile.

Figure 5.

Radial profiles of the (a) density deficiency, $\xi$ , and (b) vertical velocity, $u_z$ , for S1 taken at five different downstream locations as indicated in the legend. Both quantities are normalised by their respective centreline values, $\xi _c$ and $u_c$ , as well as plume radius as defined via the thresholding technique defined in (3.14) using half the centreline value. We also plot a generic Gaussian profile using a black dash-dotted line and the top-hat profile used at the inflow using a black dashed line.

Very close to the source at $z/D=1/4$ , profiles of $\overline {\xi }$ and $\overline {u}_z$ do not monotonically decrease as a function of $r$ , and instead, we see peaks off the centreline, more noticeably for the vertical velocity near $r/b_{u,1/2} \approx 0.5$ but also near $r/b_{u,1/2} \approx 1.75$ when $\bar {u}_z/u_c$ dips slightly below 0. The reader is directed to Meehan & Hamlington (Reference Meehan and Hamlington2023) for more details on the near-field characteristics. The radial profiles further downstream begin trending towards the Gaussian profile quickly. By $z/D=2$ , both $\overline {\xi }/\xi _c$ and $\overline {u}_z/u_c$ appear very close to Gaussian although the flow does not exhibit power-law scaling at $z/D=2$ , as shown in § 3.1.

3.4.1. Profile integrals

To analyse the profiles further, we compute profile integrals defined as follows. A plume radius is defined as the radius measured from the plume centreline where the average velocity is equal to $u_c/e$ and denoted $b_*$ , where $e$ is the exponential constant. We will discuss the ambiguity surrounding this definition in § 3.5, but this is sufficient for the present purposes. We now define two generic functions $f$ and $h$ as follows:

(3.9) \begin{equation} f(r/b_*,z) = \frac {\bar {u}_z(r/b_*,z)}{u_c(z)}, \qquad h(r/b_*,z) = \frac {\rho _a - \bar {\rho }(r/b_*,z)}{\rho _a - \rho _c(z)}. \end{equation}

We use the centreline averages in the denominator of $f$ and $h$ instead of integral statistics as a convenient way to easily collapse the data at the centreline (i.e. figure 5) and readily evaluate non-Gaussian behaviour. For integer $k$ and using $s=r/b_*$ , the profile integrals are now defined as

(3.10) \begin{equation} I_k(z) = \int _0^\infty f^k(s,z) s\, \textrm{d}{s}, \qquad J_k(z) = \int _0^\infty f^k(s,z) h(s,z) s\, \textrm{d}{s} . \end{equation}

The threshold value of $e^{-1}$ is convenient because the flow is expected to be Gaussian in the far field making the integrals in (3.10) evaluate to simple fractions: $I_1\approx J_0 \approx 1/2$ , $I_2\approx J_1 \approx 1/4$ and $J_2\approx 1/6$ . Note that using $b$ does not imply self-similarity since $f$ and $h$ to depend explicitly on $z$ ; see Townsend (Reference Townsend1976).

Figure 6.

Plots of (a) integral coefficient values as defined in (3.10) and (b) ratios of the profile coefficients that appear in the integral equations computed for S1.

We plot integral profile coefficients $I_1$ , $I_2$ , $J_0$ , $J_1$ and $J_2$ as defined in (3.10) in figure 6(a) as a function of downstream distance. At the very base, all profiles are equal to 0.5 because the inflow is a top-hat profile. Immediately above the inflow, all profiles except $J_0$ rapidly increase as a result of the large off-centre peak in the vertical velocity. When $z/D \approx 0.5 - 5.0$ , i.e. the mid-field, there is some deviation from Gaussian for all coefficients. This result demonstrates the significance of minor deviations from Gaussian at $z/D=2$ in figure 5. Noticing $\bar {u}_z/u_c$ is below the Gaussian line near the tail of the radial profile and since the tail of the radial profile is weighted more heavily than near the centreline, this causes $I_1 \gt 0.5$ for the mid-field. As the flow reaches a few diameters downstream, the profiles nearly reach their asymptotic values, which are slightly below the expected Gaussian values. For $z/D \gt 5$ , $I_1$ and $J_0$ vary between 0.44 and 0.47.

Physically, we may expect non-Gaussian behaviour in $u_z$ near the tails since there is a transition from the turbulent plume to laminar ambient as described in the turbulent/non-turbulent interface (TNTI) (Khan & Rao Reference Khan and Rao2023). For $\xi$ , there must be an abrupt change from helium–air mixture to pure air, inconsistent with the Gaussian shape. Both considerations would lead to values of $I_1$ and $J_0$ being less than 0.5.

3.4.2. Arbitrariness of the choice of plume width and profile integral ratios

In the above analysis, we defined the plume width based on a threshold value connected with velocity profiles, but we easily may have selected a different threshold value or method to compute the plume width. While this choice is arbitrary, we show how different methods may be related to each other.

Consider a different method to compute the plume width and denote it $\hat {b}_*$ . The profile coefficients are now written as

(3.11) \begin{equation} \hat {f}(r/\hat {b}_*,z) = \frac {\bar {u}_z(r/\hat {b}_*,z)}{u_c(z)}, \qquad \hat {h}(r/\hat {b}_*,z) = \frac {\rho _a - \bar {\rho }(r/\hat {b}_*,z)}{\rho _a - \rho _c(z)}. \end{equation}

Letting $\hat {s}=r/\hat {b}_*$ , the profile integrals are now defined as

(3.12) \begin{equation} \hat {I}_k(z) = \int _0^\infty \hat {f}^k(\hat {s},z) \hat {s}\, \textrm{d}{s}, \qquad \hat {J}_k(z) = \int _0^\infty \hat {f}^k(\hat {s},z) \hat {h}(\hat {s},z) \hat {s}\, \textrm{d}{\hat {s}} . \end{equation}

A simple substitution shows that

(3.13) \begin{equation} \hat {I}_k(z) = \left ( \frac {b_*}{\hat {b}_*} \right )^2 I_k(z), \qquad \hat {J}_k(z) = \left ( \frac {b_*}{\hat {b}_*} \right )^2 J_k(z). \end{equation}

This shows that the choice of $b_*$ may always be related back to a different choice for the width in a direct manner. The ratio of profile integrals, however, are not affected by this arbitrariness and indeed represent real intrinsic properties of the plume dynamics, as will be seen in § 4.

In figure 6(b), we show ratios of these profile coefficients: $J_0/I_1$ , $J_1/I_2$ and $I_2 J_0 / I_1 J_1$ . These are convenient because deviations from 1 indicate a lack of self-similarity, and the last ratio is relevant in the momentum flux transport equation, as will be discussed later. As anticipated from figure 5(a), the profiles are not self-similar in the mid-field, but achieve a values very close to 1 for $z/D\gt 5$ . This finding emphasises that the profiles do exhibit similarity far from the plume source, but the profiles do not need to be strictly Gaussian. Future work to better estimate the profile shape may be useful but is not the primary concern of the present work.

3.5. Plume radius

In the previous section, a definition for the plume radius was required to evaluate the profile coefficients, but that definition was arbitrary. We chose a definition convenient to identify non-Gaussianity and consistent between $\xi$ and $u_c$ , but any other definition may have been used. For example, the instantaneous flow fields shown in figure 1 reveal a sharp transition between the turbulent and non-turbulent regions, and this could be used to identify the plume radius, as discussed in Tennekes & Lumley (Reference Tennekes and Lumley1972), Townsend (Reference Townsend1976) and Bisset, Hunt & Rogers (Reference Bisset, Hunt and Rogers2002).

We now examine how the plume radius depends on different common definitions using data from the present plume. We concentrate on using a given flow variable $q(x,y,z,t)$ , which may be any quantity measured in the flow (e.g. velocity or density defect). We also discuss some advantages and disadvantages associated with each. Beyond integral models, the relationship between definitions is practically useful for experiments or testing when the availability of diagnostics is limited.

3.5.1. Several definitions of plume width

In § 3.5, we used a threshold definition of plume width, given as

(3.14) \begin{equation} b_{q,s}(z) \mbox{ is the largest width where } \overline {q}\big(b_{q,s}, z\big) = s\overline {q}(0,z) , \end{equation}

where $0 \lt s \lt 1$ . Throughout the text, we will often use $b_{u,1/e}= b_*$ as the standard plume radius definition. The advantage of this definition is that it is simple to compute, but it requires an a priori specification of $s$ and may be ambiguous when the profile is not monotonic, potentially yielding multiple values of $b$ . For the same quantity $q$ we may also define an integral plume width as

(3.15) \begin{equation} b_{q,{int}} = \dfrac {\int _0^\infty \overline {q}(r,z) r^{k+1}\, \textrm{d}{r}}{\int _0^\infty \overline {q}(r,z) r^k\, \textrm{d}{r}} , \end{equation}

where $k \geqslant 0$ . Again, the advantage of this definition is its simplicity, but it requires a fine spatial resolution to be computed from simulation data and an a priori choice of $k$ . In general, $k=1$ is selected, as we do in the present work. A composite definition is possible by combining integral and centreline values as

(3.16) \begin{equation} b_{q,{comp}} = \left (\dfrac {1}{\overline {q}_c}\int _0^\infty \overline {q}(r,z) r\, \textrm{d}{r} \right )^{1/2}. \end{equation}

This is, again, simple but it is difficult to interpret. Returning to the sharp TNTI, the width may be derived by considering intermittency along the shear layer. The instantaneous enstrophy, $\varOmega = \omega ^2/2$ , may be used as a good indicator of the presence of turbulence. We use a threshold value of $\varOmega _{\textit{th}}=1.0$ s $^{-2}$ for the present simulations. This value captures the turbulent region as shown in the probability density functions of vorticity components from Appendix A of Meehan et al. (Reference Meehan, Wimer and Hamlington2022). The resulting instantaneous intermittency function is defined as

(3.17) \begin{equation} \chi (x,y,z,t) = \begin{cases} 1 & \text{if } \varOmega \gt \varOmega _{\textit{th}},\\ 0 & \text{otherwise}. \end{cases} \end{equation}

This function is relatively insensitive to the exact value of $\varOmega _{\textit{th}}$ for most plumes, provided the plume is fully turbulent and is a small fraction of the maximum $\varOmega$ (Townsend Reference Townsend1976; Bisset et al. Reference Bisset, Hunt and Rogers2002). We did not find any sensitivity of the present results on $\varOmega _{\textit{th}}$ . From this, a plume radius based on intermittency may be computed as

(3.18) \begin{equation} b_{\chi }(z) \mbox{ is the largest width where } \overline {\chi }\left (b_{\chi }, z\right ) = 0.5 \,. \end{equation}

This definition is convenient when using flow visualisations and insensitive to the choice of free parameters, but it is complicated to connect to integral plume models. Finally, using purely dimensional arguments, one may define a flux plume radius as

(3.19) \begin{equation} b_{\textit{flux}} = \frac {\overline {G}}{\left ( \rho _a \overline {M}\right )^{1/2}} . \end{equation}

This definition is easy to calculate with enough data and results in a number of simplifications in the plume equations, but it is difficult to interpret. We summarise each of these definitions in table 1.

Table 1.

List of common definitions used for computing the plume radius. From left to right, we list the colloquial name, mathematical symbol, mathematical definition and free parameters (if available). The last two columns refer to what $b$ evaluates to if the profile was a top hat or Gaussian as defined in (3.28) or (3.29). We use $q$ to represent a generic flow variable.

Many of the above definitions may be related to each other. For example, if $q$ is specified as the vertical velocity, we may use (3.9) to see that

(3.20) \begin{equation} b_{u,\textit{comp}} = I_{1}^{1/2} b_{u, 1/e}. \end{equation}

Additional examples may be derived, but an important relation is shown in § 3.5.2 between the threshold-based and flux-based definitions.

3.5.2. Connection between threshold and flux width

The threshold and flux widths may be directly connected using our integral profiles from (3.9) and profile coefficients from (3.10). Using the same plume definition from § 3.4, $b_* = b_{u, 1/e}$ according to (3.14), we may rewrite the three fluxes composed of mean variables as defined in (3.4), (3.5) and (3.6), as

(3.21) \begin{align} \overline {B} &= 2 \pi g {\xi }_c {u}_c \int _0^R f(r/b_*) h(r/b_*) r\, \textrm{d}{r} \nonumber \\ &= 2 \pi g {\xi }_c {u}_c b_*^2 J_1 , \end{align}

(3.22) \begin{align} \overline {G} &= 2 \pi \int _0^R [\rho _a - (\rho _a - {\rho }_c) h(r/b_*)] {u}_c f(r/b_*) r\, \textrm{d}{r} \nonumber \\ &= 2 \pi \rho _a {u}_c b_*^2 \left (I_1 - {\xi }_c J_1\right ) , \end{align}

(3.23) \begin{align} \overline {M} &= 2 \pi b^2 \int _0^R \left [\rho _a - (\rho _a - \rho _c) h(r/b_*)\right ] {u}_c^2 f^2(r/b_*) r\, \textrm{d}{r} \nonumber \\ &= 2 \pi \rho _a {u}_c^2 b_*^2 \left (I_2 - {\xi }_c J_2\right ) . \end{align}

With some algebra, we may rewrite ${u}_c$ , ${\xi }_c$ and $b_*$ in terms of $\overline {B}$ , $\overline {G}$ and $\overline {M}$ as

(3.24) \begin{align} {\xi }_c &= \frac {I_1/J_1}{1 + \frac {g}{\rho _a}\left (\overline {G}/\overline {B}\right )} , \end{align}

(3.25) \begin{align} b_* &= \frac {1}{(2 \pi \rho _a)^{1/2}} \frac {\left (I_2 - J_2 \xi _c \right )^{1/2}}{I_1 - J_1 \xi _c} \frac {\overline {G}}{\overline {M}^{1/2}} , \end{align}

(3.26) \begin{align} {u}_c &= \frac {I_1 - J_1 \xi _c}{I_2 - J_2 \xi _c} \frac {\overline {M}}{\overline {G}}. \\[9pt] \nonumber \end{align}

We include the expression for $u_c$ as it will be used later in § 4.2.2. We left $\xi _c$ in $\bar {u}_c$ and in $b_*$ as this leads to interesting physical insights into the equations. Equation (3.25) together with flux radius definition provided in (3.19) leads to

(3.27) \begin{equation} b_* = b_{u, 1/e} = \frac {1}{(2 \pi )^{1/2}} \frac {\left (I_2 - J_2 \xi _c \right )^{1/2}}{I_1 - J_1 \xi _c} b_{q, \textrm{flux}} . \end{equation}

3.5.3. Connection with top-hat and Gaussian widths

It is common to assume that the flow is self-similar using a top-hat profile or a Gaussian profile defined, respectively, as

(3.28) \begin{align} \text{Top-hat}\!&: \quad \bar {q}(r \leqslant b_*) = q_c \quad \textrm{and} \quad \bar {q}(r\gt b_*,z) = 0, \end{align}

(3.29) \begin{align} \textrm{Gaussian}\!&: \quad \overline {q}(r,z) = q_c \exp \big({-}(r/b_{q,1/e})^2\big). \end{align}

For the top-hat profile, $b_*$ represents the radial location where the abrupt transition occurs. Further details on the top-hat and Gaussian profiles may be found in Appendix B.

When the flow exhibits self-similarity of these forms, it is possible to directly evaluate the plume radius. These values are provided in table 1. For all definitions, the plume radius scales with the local length scale of the profile. In the flux definition, there is an additional density ratio factor that would not exist if the flow were Boussinesq. Despite the common application of $b_\chi$ when the flow is self-similar in the far field, it is not possible to directly calculate $b_\chi$ knowing the flow closely resembles a Gaussian profile. Lastly, it is possible to select a threshold value to force $b_{q,s}$ to be consistent with $b_{q,\textit{int}}$ or $b_{q,\textit{comp}}$ when the flow is Gaussian; those values are, respectively, $s_{\textit{int}} = \exp (-\pi /4) \approx 0.456$ and $s_{\textit{comp}} = \exp (-1/2) \approx 0.607$ .

3.5.4. Evaluation of plume width using DNS data

We now take these various definitions of the plume radius and apply them to S1 to assess their behaviour, particularly in the near field when the flow is not self-similar. We plot each as a function of downstream distance in figure 7, and when possible, we compare $b$ using $\xi _c$ and $u_c$ . Note that the plume radius $b_{u,1/e}$ indicated by the black line is the commonly used definition throughout the text and may be referenced as $b_*$ .

We make a few observations. First, each definition captures the ‘necking’ or contracting behaviour that is observed at the base of the plume but to varying degrees, with $b_{\textit{flux}}$ having the smallest neck. For definitions requiring a flow variable $q$ , there is a some dependence whether $q$ is $\xi$ or $u_c$ in the near field, but they essentially converge once the flow has propagated far enough downstream. Since the flow is nearly Gaussian beyond $z/D \approx 4$ , each definition appears be directly correlated with each other. This correlation, however, is not apparent in the near field where the flow is not self-similar, so caution must be used when interpreting plume radius data in the absence of self-similarity.

Figure 7.

Plot of plume radius as a function of downstream distance using various definitions provided in table 1. The inset provides a zoomed-in perspective of the near field. The definition used throughout the text is $b_{u,1/e}$ , or more simply referred to as $b_*$ .

4.1. Governing plume equations

The governing equations that fully describe the fluid flow of the helium plumes considered here are conservation equations for mass, momentum and helium mass fraction. From these, we will provide a brief derivation of the plume equations. Note that viscous and diffusive terms are present in the numerical simulations, but they do not have a major influence on the plume dynamics provided the flow is turbulent (Tennekes & Lumley Reference Tennekes and Lumley1972; Hunt & van den Bremer Reference Hunt and van den Bremer2010) so they are neglected in the present derivation. The results from the present simulations support this assumption.

The inviscid, immiscible governing equations are

(4.1) \begin{align} \frac {\partial {\rho }}{\partial {t}} + \frac {\partial {(\rho u_i)}}{\partial {x_i}} &= 0 , \end{align}

(4.2) \begin{align} \frac {\partial {(\rho u_i)}}{\partial {t}} + \frac {\partial {(\rho u_i u_j)}}{\partial {x_j}} &= -\frac {\partial {p}}{\partial {x_i}} + (\rho -\rho _a) g_i , \end{align}

(4.3) \begin{align} \frac {\partial {(\rho Y_{\textit{He}})}}{\partial {t}} + \frac {\partial {(\rho u_i Y_{\textit{He}})}}{\partial {x_i}} &= 0 . \end{align}

Notice that the pressure is separated into two parts: a mechanical part, $\partial p/\partial x_i$ , and a hydrostatic part, $(\rho -\rho _a) g_i$ , as is common in buoyant plumes (Townsend Reference Townsend1976; Bharadwaj & Das Reference Bharadwaj and Das2017; Meehan & Hamlington Reference Meehan and Hamlington2023). We will only consider scenarios where $g_i = [0,0,-g]$ , where $g$ is a positive constant.

To derive the plume equations, we consider a cylindrical coordinates system with a non-swirling flow (i.e. $u_\theta =0$ ) and perform a Reynolds average on (4.1), (4.2) and (4.3). We will ignore the radial momentum equation. We then perform a radial integration from 0 to $R$ where $R$ is a fixed radial distance far from the plume core. The conservation equation for mass yields the mass flux equation

(4.4) \begin{equation} \frac {\textrm{d}G}{\textrm{d}z} = 2 \pi \rho _a Ru_R, \end{equation}

where $u_R$ is the radial velocity at $r=R$ . The buoyant flux equation is a result of combining the mass and helium mass fraction conservation equations

(4.5) \begin{equation} \frac {\textrm{d}{B}}{\textrm{d}{z}} = 0. \end{equation}

A similar expression may be achieved if the plume was generated via temperature differences (assuming an ideal gas law) rather than differing fluid composition (Rooney & Linden Reference Rooney and Linden1996). Finally, evaluation of the conservation of momentum equation in the vertical direction yields

(4.6) \begin{equation} \frac {\textrm{d}{M}}{\textrm{d}{z}} = -\underbrace {\int _0^R \frac {\partial {\bar {p}}}{\partial {z}} 2 \pi r\, \textrm{d}{r}}_{\textrm{d}{P}/\textrm{d}{z}} +\underbrace {\int _0^R g(\rho _a - \bar {\rho }) 2 \pi r\textrm{d}{r}}_{D} . \end{equation}

We defined two additional fluxes, $\textrm{d}{P}/\textrm{d}{z}$ and $D$ , for subsequent discussion. Definitions for $G$ , $B$ and $M$ may be found in § 3.2. The equations (4.4), (4.5) and (4.6) are the basis for the governing plume equations.

Figure 8.

Term-by-term calculation of integral equations for (a) mass flux, (4.4), and (b) momentum flux, (4.6).

We now evaluate the production terms for (4.4) and (4.6) in figure 8. The buoyant flux was shown to be approximately constant in § 3.2 so we do not replot $B$ here. The terms for the balance equation for the mass flux are shown figure 8(a). The left- and right-hand sides balance perfectly: $\textrm{d}{G}/\textrm{d}{z} = -2\pi \rho _a Ru_R$ . This is not surprising since there are (essentially) no assumptions required to achieve this equation. Interestingly, we see $\textrm{d}{G}/\textrm{d}{z}$ increase approximately linearly, except for the very near field where boundary layer effects reduce the amount of mass entrained.

We next look at the terms for the momentum flux presented in figure 8(b). The momentum flux equation $\textrm{d}{M}/\textrm{d}{z} \approx -\textrm{d}{P}/\textrm{d}{z} + D$ does not include viscous contributions, and from our calculations, this is a well-supported approximation since the left- and right-hand sides are nearly equal. The commonly accepted boundary layer approximation assumes that $\textrm{d}{P}/\textrm{d}{z}\approx 0$ ; see Townsend (Reference Townsend1976, § 6.2) and Tennekes & Lumley (Reference Tennekes and Lumley1972, § 4.6). Our simulations show $\max (|\textrm{d}{P}/\textrm{d}{z}|/|\textrm{d}{M}/\textrm{d}{z}|) \approx 0.10-0.15$ , which is consistent with well-known experiments (Shabbir & George Reference Shabbir and George1994) and simulations (van Reeuwijk et al. Reference van Reeuwijk, Salizzoni, Hunt and Craske2016), and a result of the turbulent fluctuations that persist throughout the plume.

4.2. Closure without self-similarity

4.2.1. Closure from dimensional analysis

We need to close equations for mass and momentum fluxes presented in § 4.1 in terms of the three fluxes, $B$ , $M$ and $G$ , and the ambient density, $\rho _a$ ,

(4.7) \begin{align} \frac {\textrm{d}{G}}{\textrm{d}{z}} = \phi _1 (B, M, G, \rho _a, \ldots ) , \qquad \frac {\textrm{d}{M}}{\textrm{d}{z}} = \phi _2 (B, M, G, \rho _a, \ldots ), \end{align}

where $\phi _1$ and $\phi _2$ are functions to be determined. There are two principal strategies to do so: either by identifying relevant parameters and using dimensional analysis, or by using extra physical considerations. We start with dimensional analysis and leaving more physical arguments to later sections. Dimensional considerations lead to

(4.8) \begin{align} \frac {1}{2\left (\pi \rho _a M\right )^{1/2}} \phi _1 = \alpha \left ( \varGamma , \ldots \right ) ,\qquad \frac {1}{\rho _a} \frac {M}{BG} \phi _2 = \alpha ' \left ( \varGamma , \ldots \right ) , \end{align}

where $\alpha$ and $\alpha '$ are non-dimensional functions of the relevant non-dimensional parameters, including the Richardson number, $\varGamma = \rho _a^{1/2}B G^2 M^{-5/2}$ , defined previously in § 3.3. The numerical factor $2 \pi ^{1/2}$ has been introduced for convenience. This leads to the following closed-form equations:

(4.9) \begin{align} \frac {\textrm{d}{G}}{\textrm{d}{z}} &= 2\alpha\!\left (\pi \rho _a M\right )^{1/2} , \end{align}

(4.10) \begin{align} \frac {\textrm{d}{M}}{\textrm{d}{z}} &= \alpha '\rho _a \frac {BG}{M} . \end{align}

Here, $\alpha$ is the entrainment coefficient, and a closure model is needed at leading order, even when the flow is consistent with the standard assumptions (e.g. Boussinesq, self-similar, slender, etc.). On the other hand, $\alpha '$ is the momentum correction coefficient, and represents the small perturbation correction when the flow is inconsistent with the standard assumptions. In §§ 4.2.2 and 4.2.3, we will perform some analysis to gain physical insights into potential closure formulations, and in § 4.3, we will then use simulation data to estimate the actual evolution with respect to other parameters in § 4.3.

4.2.2. Closure of the momentum conservation equation

We now close the momentum flux equation provided in (4.6). We first make the boundary layer approximation such that $\textrm{d}{P}/\textrm{d}{z} = 0$ . This is a common assumption, but as shown in §§ 3.2 and 4.1, the error may not always be adequate. We next write the density deficit integral in (4.6) using the familiar profile coefficients and characteristic variables as

(4.11) \begin{equation} \int _0^R g(\rho _a - \bar {\rho }) 2\pi r\textrm{d}{r} = 2 \pi J_0 g(\rho _a - \rho _c) b^2 . \end{equation}

It is possible to use $\bar {u}_c$ , $\bar {\xi }_c$ and $b$ as derived in (3.24), (3.26) and (3.25) to rewrite (4.11) as

(4.12) \begin{eqnarray} 2 \pi J_0 g(\rho _a - \rho _c) b^2 &=& \rho _a \frac {J_0 I_2}{I_1 J_1} \boldsymbol{\cdot }\frac {1 - (J_2/I_2) \xi _c}{1 - (J_1/I_1) \xi _c} \frac {\overline {B} \,\overline {G}}{\overline {M}} . \end{eqnarray}

To obtain a closure for $\alpha '$ , we need to approximate the mean fluxes as approximately equal to the total fluxes, i.e. $B \approx \overline {B}$ , $G \approx \overline {G}$ and $M \approx \overline {M}$ . In doing so, we may compare (4.12) with (4.6) to obtain an approximate value for $\alpha '$

(4.13) \begin{equation} \alpha ' \approx \frac {J_0 I_2}{I_1 J_1} \boldsymbol{\cdot }\frac {1 - (J_2/I_2) \xi _c}{1 - (J_1/I_1) \xi _c} . \end{equation}

When the plume is Boussinesq (i.e. $\xi _c \ll 1$ ), this expression simplifies to $\alpha ' = J_0 I_2/I_1 J_1$ , and if the flow adheres to a top-hat or Gaussian profile, $\alpha ' \approx 1$ . This explains why it is generally assumed that $\alpha ' \approx 1$ expect Lamalle et al. (Reference Lamalle, Carlotti, Salizzoni and Perkins2013) who introduced it. For the present work, we defined $\alpha '$ using the total fluxes as shown in (4.6); from the discrepancies shown in § 3.2 and the boundary layer approximation in § 4.1, we expect $\alpha '$ to deviate from 1. We will use the simulation results to calculate $\alpha '$ in § 4.3.

4.2.3. Closure of the mass conservation equation

The closure of mass conservation (4.4) raises a different challenge because the value of $R u_R$ may not be expressed directly from plume quantities without further assumption. Assuming $u_z$ decays to 0 exponentially fast outside of the plume, we may then approximate $Ru_R$ as a constant via continuity when $R$ is much greater than the plume radius; we use a value of $R$ equivalent to the outer edge of the computational domain. This is an approximation from solving the irrotational flow problem outside of the plume, but it is good enough far from the plume. From this, we write $Ru_R$ using plume variables by equating the right-hand sides of (4.4) and (4.9) in order to understand the relationship of  $\alpha$ , yielding

(4.14) \begin{equation} 2\alpha \left (\pi \rho _a M\right )^{1/2} = 2 \pi \rho _a Ru_R . \end{equation}

We make the same approximation from § 4.2.2 that $M = \overline {M}$ , then using (3.23), we may rearrange (4.14) to arrive at an expression for $\alpha$

(4.15) \begin{eqnarray} \alpha = \frac {1}{\sqrt {2 \left (I_2 - {\xi }_c J_2\right )}} \frac {\textit{Ru}_R}{b_* u_c}. \end{eqnarray}

It is not possible to simplify this formulation further in a general sense.

We can show there is consistency between this formulation and classical plume theory after invoking similar assumptions (Morton et al. Reference Morton, Taylor and Turner1956; Rooney & Linden Reference Rooney and Linden1996; Hunt & van den Bremer Reference Hunt and van den Bremer2010; Carlotti & Hunt Reference Carlotti and Hunt2017). First, we invoke a top-hat profile consistent with the previous work, yielding $I_2 = J_2 = 1/2$ . Second, $Ru_R$ is independent of $r$ outside of the plume; therefore, with a sufficiently large $b_*$ , $R u_R \approx b_* u_e$ , where $u_e$ is horizontal velocity at the outer edge of the plume or the so-called entrainment velocity. Applying these results in

(4.16) \begin{equation} u_e = \alpha u_c \left (\frac {\rho _c}{\rho _a}\right )^{1/2}, \end{equation}

which is exactly the definition of entrainment coefficient of classical top-hat plume theory (Rooney & Linden Reference Rooney and Linden1996). This suggests that the present formulation retains the features of the plume equations for Boussinesq plumes under a top-hat assumption.

4.2.4. Final form of plume dynamical equations

Applying the closures presented above provides the governing plumes equations purely in flux form as

(4.17) \begin{align} \frac {\textrm{d}{B}}{\textrm{d}{z}} &= 0 , \end{align}

(4.18) \begin{align} \frac {\textrm{d}G}{\textrm{d}z} &= 2\alpha\!\left (\pi \rho _a \right )^{1/2}\!M^{1/2} , \end{align}

(4.19) \begin{align} \frac {\textrm{d}M}{\textrm{d}z} &= \alpha ' \rho _a \frac {BG}{M} . \end{align}

This set of equations constitutes a closed set of plume equations. We now need to estimate  $\alpha$ , which has a leading-order effect on the entrainment, and $\alpha '$ , where deviations from unity have a secondary effect on correcting the momentum flux. Both terms may have a functional dependence on $z$ or intrinsic non-dimensional numbers. We evaluate these in § 4.3. While a number of assumptions were made in the closure, we did not use a self-similarity hypothesis.

4.3. Estimation of plume coefficients

We now explore estimations of $\alpha$ and $\alpha '$ by directly computing $\alpha$ and $\alpha '$ from the numerical simulations. This may be done by computing $\alpha$ from (4.9) as

(4.20) \begin{equation} \alpha = \frac {1}{2}\left (\pi \rho _a M\right )^{-1/2}\frac {\textrm{d}G}{\textrm{d}z} . \end{equation}

In a similar fashion, we compute $\alpha '$ by rearranging (4.10)

(4.21) \begin{equation} \alpha ' = \frac {1}{\rho _a} \frac {M}{BG} \frac {\textrm{d}M}{\textrm{d}z} . \end{equation}

We plot the computed values of $\alpha$ and $\alpha '$ in figure 9 for both S1 and S2 and make a few observations.

Figure 9.

Evaluation of (a) the entrainment coefficient via (4.20) and (b) the momentum correction coefficient via (4.21). The inset in (a) is a zoomed-in region of the same plot to observe the near-field values. The asymptotic value for $\alpha$ is $\alpha _0 \approx 0.114$ , and the momentum correction coefficient is typically assumed to be $\alpha '=1$ .

The entrainment coefficient, $\alpha$ , is quite large near the base where the plume is very buoyant, with a maximum value of $\alpha \approx 0.6{-}0.7$ for S1 and $\alpha \approx 1.6{-}1.7$ for S2. This difference is a consequence of the no-slip boundary condition in S1, forcing horizontal velocity to 0 in the boundary layer close to the plane $z=0$ , while S2 has a boundary condition designed such that the horizontal velocity matches the value expected just above the boundary layer. The boundary layer also shifts the peak of $\alpha$ for S1 above the source. Interestingly, Webb, Wise & Hunt (Reference Webb, Wise and Hunt2023) estimated $\alpha \approx 1.2$ at the base of their experiments of heated plates. Since their source Richardson number is larger (e.g. $\varGamma =\infty$ ), we would expect increased entrainment, but our results suggest entrainment will increase further when the boundary layer is not present. Nonetheless, the large $\alpha$ near the base supports empirical evidence presented by simulations and experiments suggesting there is enhanced entrainment when $z/D$ is small (Pham et al. Reference Pham, Plourde and Kim2005; Lamalle et al. Reference Lamalle, Carlotti, Salizzoni and Perkins2013; Carlotti & Hunt Reference Carlotti and Hunt2017), and contradicts some theoretical propositions that there is little entrainment at the base of the plume (Hunt & Kaye Reference Hunt and Kaye2005; Diez & Dahm Reference Diez and Dahm2007). Theoretical work has been undertaken to improve the physical representation of the increased entrainment (Carazzo, Kaminski & Tait Reference Carazzo, Kaminski and Tait2006; van Reeuwijk & Craske Reference van Reeuwijk and Craske2015; Craske, Salizzoni & van Reeuwijk Reference Craske, Salizzoni and van Reeuwijk2017; Saeed et al. Reference Saeed, Weidner, Johnson and Mandel2022; Webb et al. Reference Webb, Wise and Hunt2023). Beyond the very near field, S1 and S2 are nearly identical. Farther downstream, $\alpha$ asymptotically approaches the self-similar Boussinesq pure-plume regime. Experimental and simulation data suggest an asymptotic value of $\alpha _0 \approx 0.10{-}0.13$ with a generally accepted value of $\alpha _0=0.12$ (Ricou & Spalding Reference Ricou and Spalding1961; Georges, Alpert & Tamanini Reference Georges, Alpert and Tamanini1977; Turner Reference Turner1979; List Reference List1982; Rooney & Linden Reference Rooney and Linden1996; Wang & Law Reference Wang and Law2002; Pham et al. Reference Pham, Plourde and Kim2005; Hunt & van den Bremer Reference Hunt and van den Bremer2010; Lamalle et al. Reference Lamalle, Carlotti, Salizzoni and Perkins2013; Ezzamel et al. Reference Ezzamel, Salizzoni and Hunt2015). From our data, we calculated a value of $\alpha _0\approx 0.115$ , conforming to the breadth of prior research. This asymptotic value, however, differs slightly from asymptotic behaviour of Richardson number described in § 3.3 due to the deviation from pure Gaussian profiles, as discussed in § 3.4.

The momentum coefficient, $\alpha '$ , is relatively constant for S1 throughout the domain and in the range of $1.0{-}1.2$ despite the changes in flow structure as the plume propagates downstream. This deviation from 1 is similar in magnitude with differences associated with approximations required to close the governing equations, as seen in §§ 3.2 and 4.1. The magnitude of $\alpha '$ for S2 is quite similar except for some statistical variability.

Expressions for $\alpha$ as function of non-dimensional plume parameters were proposed by many authors (see Carlotti & Hunt (Reference Carlotti and Hunt2017) and Ciriello & Hunt (Reference Ciriello and Hunt2020) and references therein). The two most popular expressions are

(4.22) \begin{equation} \alpha = \alpha _0 ( \varGamma /\varGamma _0)^\lambda , \qquad \alpha = \alpha _0 + \alpha _0 \lambda ({\varGamma }/{\varGamma _0} - 1). \end{equation}

The first equation was proposed by Carlotti & Hunt (Reference Carlotti and Hunt2017) while the second equation was introduced much earlier (see Ciriello & Hunt Reference Ciriello and Hunt2020) but is merely a linearisation of the first one when $0 \lt \varGamma / \varGamma _0 - 1 \ll 1$ (Ciriello & Hunt Reference Ciriello and Hunt2020). The constant $\lambda$ is unknown and needs to be determined from the data; physical considerations require $0 \leqslant \lambda \lt 1$ (Carlotti & Hunt Reference Carlotti and Hunt2017).

In figure 10, we present the variation of $\alpha$ with respect to $\varGamma$ by plotting $\tilde {\alpha } = \alpha /\alpha _0$ as a function of $\tilde {\varGamma } = \varGamma /\varGamma _0$ as computed from S1 and S2. We include coloured symbols where $z/D = 0.1$ and $z/D = 0.5$ to provide reference to the physical space. We also include lines corresponding to the nonlinear equation shown in (4.22) with $\lambda =1$ , $4/5$ and $3/5$ , and we include the linearisation using $\lambda =4/5$ . We see that $\lambda =4/5$ does a good job in predicting the dependence of $\alpha$ on $\varGamma$ for both the linear and nonlinear models for $\varGamma /\varGamma _0 \lesssim 3$ , but beyond that point, the entrainment is overpredicted by both formulations.

Figure 10.

Comparison of proposed models of $\alpha$ as a function of $\varGamma /\varGamma _0$ over the entire domain. We include coloured markers to provide reference to physical space.

We would like to propose an empirical correction for $\alpha$ that accounts for the inaccuracy of $\alpha$ at large $\varGamma$ . Based on the work of Webb et al. (Reference Webb, Wise and Hunt2023), we expect that $\lim _{\varGamma \to \infty }\alpha = \alpha _\infty \approx 1.2$ . The limiting behaviours of $\alpha$ would be

(4.23) \begin{equation} \alpha = \begin{cases} \alpha _0 ( \varGamma /\varGamma _0)^\lambda & \textrm{for} \quad \varGamma /\varGamma _0 \to 1^+ ,\\ \alpha _\infty & \textrm{for} \quad \varGamma /\varGamma _0 \to \infty .\end{cases} \end{equation}

We may incorporate this behaviour by multiplying our nonlinear expression in (4.22) by  $\phi$ , i.e. $\alpha = \alpha _0 (\varGamma /\varGamma _0)^\lambda \boldsymbol{\cdot }\phi (\varGamma /\varGamma _0)$ , where $\phi$ is a function of $\varGamma /\varGamma _0$ and is defined as

(4.24) \begin{equation} \phi \left (\frac {\varGamma }{\varGamma _0}\right ) = \left [ 1 + \frac {\alpha _0}{\alpha _\infty } \left (\frac {\varGamma }{\varGamma _0} - 1 \right ) \right ]^{-1} . \end{equation}

The proposed $\phi$ in (4.24) is a simple function that achieves the appropriate limiting behaviours provided in (4.23) without any additional tuning parameters. From this, our proposed entrainment function is

(4.25) \begin{equation} \frac {\alpha }{\alpha _0} = \left (\frac {\varGamma }{\varGamma _0}\right )^\lambda \left [ 1 + \frac {\alpha _0}{\alpha _\infty } \left (\frac {\varGamma }{\varGamma _0} - 1 \right ) \right ]^{-1} , \end{equation}

with $\lambda = 4/5$ . We plot this function using a green line in figure 10 and compare it with the simulation data. The function is in better agreement with the data at large $\varGamma$ compared with the expressions provided in (4.22), and it is consistent with known endpoint behaviour in the limit of large and small $\varGamma$ . There are, of course, many possible modifications to $\phi$ in (4.24) that could improve the agreement with the data, but we leave this exploration when additional simulation data with different buoyancy characteristics are available.

One additional consideration related to the entrainment model is the density ratio correction of $(\rho _c/\rho _a)^{1/2}$ that is typically used for non-Boussinesq plumes. For the present work, the density ratio correction is obfuscated by the lack of self-similarity, but when self-similarity is present, we are in alignment with Rooney & Linden (Reference Rooney and Linden1996), as shown in § 4.2.3. Of course, we do not achieve power-law scaling until far from this source when $\varGamma$ reaches its asymptotic value, but at that point, the flow is Boussinesq, as indicated by $\xi _c$ in figure 2(a). Moreover, it is difficult from the present singular set of plume source conditions to definitively identify this correction since $\alpha$ depends on $\varGamma$ , and $\varGamma$ depends on the density ratio. Nonetheless, the experiments of plumes by Lanzini et al. (Reference Lanzini, Marro, Creyssels and Salizzoni2024) and the simulations of jets by Salizzoni et al. (Reference Salizzoni, Vaux, Creyssels, Craske and van Reeuwijk2024) suggest this correction may need to be reconsidered.