2.1. Integral equations

In this section, equations for the integral volume, momentum and buoyancy fluxes are derived. The derivation is consistent with that of Weil (Reference Weil1988) and Schatzmann (Reference Schatzmann1978), but is performed here without any simplifications, such as the assumption of axisymmetry (Schatzmann Reference Schatzmann1978) or neglecting turbulence and pressure terms. A schematic detailing the nomenclature used here is given in figure 1. Cartesian coordinates $(x, y, z)$ are used which represent the streamwise, transverse and vertical directions, respectively. The crossflow speed in the $x$ direction is denoted $U$ and is uniform in $z$. It is convenient to also introduce a plume-following coordinate system $(s,y,\eta )$, where $s$ is the distance from the source along the plume centreline, $y$ is the transverse direction (which remains unchanged) and $\eta$ is the coordinate perpendicular to $s$ in the $(x, z)$ plane. The plume centreline makes an angle $\varphi$ with the $x$ axis. The plume source is positioned at (0, 0, 0) and the release is momentumless with a constant integral buoyancy flux $F_0 = r_0^2 \phi _b$, where $r_0$ is the source radius and $\phi _b$ is the (diffusive) buoyancy flux. This infinitely lazy plume condition is chosen so as not to introduce a further parameter, namely the source momentum flux, and to manage the complexity of the problem.

Figure 1. A schematic of the plume with the Cartesian $(x,z)$ coordinate system and the curvilinear $(s, \eta )$ coordinate system. The $y$ direction is perpendicular to this plane. Values for $L_x$, $L_z$ and $L_n$ for the five simulations can be found in table 1.

The incompressible Navier–Stokes equations with the Boussinesq approximation are

(2.1a)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} =0, \end{gather}

(2.1b)\begin{gather}\frac{\partial\boldsymbol{u}}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u} ={-}\boldsymbol{\nabla} p + \nu \nabla^2 \boldsymbol{u} + b \boldsymbol{\hat{k}}, \end{gather}

(2.1c)\begin{gather}\frac{\partial b}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} b = \kappa \nabla^2 b, \end{gather}

where $\boldsymbol {u} = (u,v,w)$ is the velocity of the fluid, $b=g (\rho _0 - \rho ) / \rho _0$ is the buoyancy and $p=\tilde p/\rho _0 + gz$ is the kinematic pressure perturbation, where $g$ is the gravitational acceleration, $\rho _0$ is a constant reference density and $\tilde p$ is the standard pressure. The kinematic viscosity and scalar diffusivity are denoted by $\nu$ and $\kappa$, respectively. Assuming high-Reynolds-number and high-Péclet-number flow, applying Reynolds averaging to (2.1) and making use of the fact that the flow is statistically steady results in

(2.2a)\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \bar{\boldsymbol{u}} = 0, \end{gather}

(2.2b)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} (\bar{\boldsymbol{u}}\overline{u_i} +\overline{\boldsymbol{u}'u_i'}) ={-}\frac{\partial\bar{p}}{\partial x_i} + \bar{b} \delta_{i3}, \end{gather}

(2.2c)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} (\bar{\boldsymbol{u}}\bar{b} +\overline{\boldsymbol{u}'b'}) = 0, \end{gather}

where the overbars denote time averaging and fluctuating quantities are defined as $\phi '=\phi -\bar {\phi }$ for any variable $\phi$.

Using that $\overline {\boldsymbol {u}} \bar {\phi } + \overline {\boldsymbol {u}'\phi '}=\overline {\boldsymbol {u} \phi }$, the equations above are all of the form (e.g. Weil Reference Weil1988)

(2.3)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \left(\overline{\boldsymbol{u}\phi}\right) = G, \end{equation}

which can be integrated over the plume cross-section $A$, whose domain is denoted $\varOmega$, to obtain (van Reeuwijk, Vassilicos & Craske Reference van Reeuwijk, Vassilicos and Craske2021) (see Appendix A)

(2.4)\begin{equation} \frac{\text{d}}{\text{d} s} \displaystyle\iint_\varOmega \bar{u}_s \bar{\phi} + \overline{{u}'_s\phi'} \,\text{d} A = \displaystyle\iint_\varOmega G\,\text{d} A + {\rm \pi} E[\phi] + {\rm \pi} \mathcal{L}[\phi], \end{equation}

where $\bar {u}_s=\bar {u} \cos \varphi + \bar {w} \sin \varphi$ is the mean velocity in the $s$ direction and $\varOmega (s)$ represents the domain occupied by the plume. See § 3.2 for a description of how $\varOmega$ is determined practically. In deriving the equation above it was assumed that the curvature of the plume trajectory $\text {d} \varphi / \text {d} s \ll 1$ (see Appendix A for more details). The boundary terms $E[\phi ]$ and $\mathcal {L}[\phi ]$ are line integrals of the form

(2.5a,b)\begin{equation} E[\phi] = \frac{1}{{\rm \pi}}\oint_{\partial \varOmega} \overline{\boldsymbol{u}_\perp \phi} \boldsymbol{\cdot} \boldsymbol{n}\, \text{d} \ell, \quad \mathcal{L}[\phi] = \frac{1}{{\rm \pi}}\oint_{\partial \varOmega} \overline{u_s \phi} \frac{N_s}{| \boldsymbol{N}_\perp|}\, \text{d} \ell, \end{equation}

where $\partial \varOmega$ denotes the plume boundary, $\boldsymbol {N} = (N_s, N_y, N_\eta )^\textrm {T}$ denotes the three-dimensional inward-pointing normal along $\partial \varOmega$, $\boldsymbol {N}_\perp = (N_y, N_\eta )^\textrm {T}$ and $\boldsymbol {u}_\perp = (u_y, u_\eta )^\textrm {T}$ are the vector components in the $(y, \eta )$ plane and $\boldsymbol {n} = \boldsymbol {N}_\perp / |\boldsymbol {N}_\perp |$ is the inward-pointing two-dimensional normal in the $(y, \eta )$ plane. As the plume boundary will be detected via the average buoyancy field $\bar {b}$ (see § 4), $\boldsymbol {N}$ is defined as $\boldsymbol {N} = \boldsymbol {\nabla } \bar {b} / | \boldsymbol {\nabla } \bar {b}|$ (see also van Reeuwijk et al. Reference van Reeuwijk, Vassilicos and Craske2021). The factor $1/{\rm \pi}$ in the definition of $\mathcal {L}[\phi ]$ and $E[\phi ]$ is introduced for consistency with the convention in plume theory of dividing the integral quantities through by ${\rm \pi}$ (see (2.7)).

The terms $E[\phi ]$ and $\mathcal {L}[\phi ]$ represent distinct physical processes by which the quantity $\phi$ enters the plume. The term $E[\phi ]$ represents the advection of $\phi$ into the plume; the term $\mathcal {L}[\phi ]$ represents incorporation of the quantity $\phi$ via expansion of the plume into a non-quiescent atmosphere – if the plume surface area $A$ over which the governing equations are integrated changes as a function of $s$, and the value of $\phi$ at the boundary of the plume is non-zero, then the surface integral will be increased by a quantity equal to $[{\overline {u_s\phi }}]_b \,\text {d} A/\text {d} s$, where ${[\overline {u_s\phi }]}_b$ is the average value of $\overline {u_s\phi }$ at the boundary.

Integrating the Navier–Stokes equations over the plume cross-section and applying (2.4) with $\phi = 1, u-U, w$ and $b$, respectively, results in

(2.6a)\begin{gather} \frac{\text{d} Q}{\text{d} s} = E[1] + \mathcal{L}[1]= E + \mathcal{L}, \end{gather}

(2.6b)\begin{gather}\frac{\text{d} M_x}{\text{d} s} + \frac{\text{d} M'_x}{\text{d} s} ={-}P_x + E[u-U] + \mathcal{L}[u-U], \end{gather}

(2.6c)\begin{gather}\frac{\text{d} M_z}{\text{d} s} + \frac{\text{d} M'_z}{\text{d} s} ={-}P_z + B + E[w] + \mathcal{L}[w], \end{gather}

(2.6d)\begin{gather}\frac{\text{d} F}{\text{d} s} + \frac{\text{d} F'}{\text{d} s} = E[b] + \mathcal{L}[b], \end{gather}

where the integral volume flux $Q$, the horizontal momentum excess flux $M_x$, the vertical momentum flux $M_z$, the buoyancy flux $F$ and the integral horizontal and vertical pressures $P_x$ and $P_z$ are defined as

(2.7)\begin{equation} \left. \begin{array}{ll@{}} Q = \dfrac{1}{{\rm \pi}} \displaystyle\iint_\varOmega \bar{u}_s \,\text{d} A, & F = \dfrac{1}{{\rm \pi}} \displaystyle\iint_\varOmega \bar{u}_s \bar{b} \,\text{d} A,\\M_x = \dfrac{1}{{\rm \pi}} \displaystyle\iint_\varOmega \bar{u}_s (\bar{u}-U) \,\text{d} A, & P_x = \dfrac{1}{{\rm \pi}} \displaystyle\iint_\varOmega \dfrac{\partial \bar{p}}{\partial x} \,\text{d} A,\\M_z = \dfrac{1}{{\rm \pi}} \displaystyle\iint_\varOmega \bar{u}_s \bar{w} \text{d} A, & P_z = \dfrac{1}{{\rm \pi}} \displaystyle\iint_\varOmega \dfrac{\partial \bar{p}}{\partial z} \,\text{d} A,\\M'_x = \dfrac{1}{{\rm \pi}} \displaystyle\iint_\varOmega \overline{u'_s u'} \,\text{d} A, & F' =\dfrac{1}{{\rm \pi}} \displaystyle\iint_\varOmega \overline{u'_sb'} \,\text{d} A,\\M'_z = \dfrac{1}{{\rm \pi}} \displaystyle\iint_\varOmega \overline{u'_s w'} \text{d} A, & B = \dfrac{1}{{\rm \pi}} \displaystyle\iint_\varOmega \bar{b} \,\text{d} A. \end{array} \right\} \end{equation}

The scalar $\bar {u}-U$ is chosen for the $x$-momentum equation so as to be consistent with the work of Weil (Reference Weil1988) in the definition of $M_x$. The integral therefore represents the momentum flux in the $x$ direction in excess of the $x$-momentum flux of the ambient atmosphere. A consequence of this choice of variable is that the Leibniz and entrainment terms for the $x$-momentum equation are significantly reduced, due to the fact that to first order $u \approx U$ at the boundary of the plume.

2.2. Plume models

In conventional plume theory $\bar {u} - U$, $\bar {w}$ and $\bar {b}$ are assumed to vanish at the plume boundary. In addition $\mathcal {L}[1]$ is assumed to be small. It thus follows that all entrainment and Leibniz terms are negligible except for $E \equiv E[1]$. Furthermore, it is commonly assumed that the pressure at the plume boundary is equal to the ambient pressure, which implies that $P_x$ and $P_z$ are zero (although vertical pressure gradients are taken into account implicitly via an added-mass coefficient). The buoyancy integral $B = {\rm \pi} ^{-1} \int \bar {b} \text {d} A$ is generally modelled as $B=F/U_s$ (Weil Reference Weil1988), where $U_s$ is the characteristic velocity of the plume, defined as

(2.8)\begin{equation} U_s = (U + u_m) \cos \varphi + w_m \sin \varphi. \end{equation}

Here, $u_m$ is is the characteristic horizontal plume velocity relative to the crossflow (horizontal velocity excess), $w_m$ is the characteristic vertical velocity and $\varphi$ is the angle of the plume. These characteristic quantities are linked to the integral quantities by

(2.9a–d)\begin{equation} Q = r_m^2 U_s, \quad M_x = r_m^2 U_s u_m, \quad M_z = r_m^2 U_s w_m, \quad F = r_m^2 U_s b_m. \end{equation}

Here, $b_m$ is the characteristic buoyancy or reduced gravity and $r_m$ is the characteristic plume radius. Equation (2.9a–d) implicitly provides the definition of $u_m$, $w_m$ and $b_m$ as $M_x/Q$, $M_z/Q$ and $F/Q$, respectively. The characteristic radius $r_m$ is defined as $r_m=(Q/U_s)^{1/2}$ (i.e. using the definition of $Q$; (2.9a)), which encapsulates both the pure plume definition $r_m = Q/M_z^{1/2}$ for $U\ll w_m$, $u_m \ll w_m$, and the bent-over solution $r_m = (Q/U)^{1/2}$ for $U\gg u_m, U \gg w_m$.

The entrainment term $E$ is typically parameterised as (Weil Reference Weil1988; Devenish et al. Reference Devenish, Rooney, Webster and Thomson2010)

(2.10)\begin{equation} E = 2 \gamma r_m w_m, \end{equation}

where $\gamma$ is the entrainment coefficient and $r_m$ and $w_m$ are a characteristic radius and vertical velocity, respectively. Invoking the assumptions stated at the start of this section ($\mathcal {L} \ll 1$ and $E \equiv E[1]$), neglecting turbulent fluxes and substituting (2.10) into the integral equations (2.6a–d) result in

(2.11a–d)\begin{equation} \frac{\text{d} Q}{\text{d} s} = 2 \gamma \frac{M_z}{r_m U_s}, \quad \frac{\text{d} M_x}{\text{d} s} = 0, \quad (1+k)\frac{\text{d} M_z}{\text{d} s} = \frac{F}{U_s}, \quad \frac{\text{d} F}{\text{d} s} = 0. \end{equation}

Here $k$ is the added-mass constant which takes into account the effects of an adverse pressure gradient on the vertical motion (e.g. Weil Reference Weil1988), which compensates for the fact that no explicit integral pressure terms are present. It is introduced by analogy with the motion of a solid body and is equivalent to using a drag term (see Briggs Reference Briggs1982; Schatzmann Reference Schatzmann1978; Davidson Reference Davidson1989). Kinematic relations for the plume trajectory are required to evolve the plume in space:

(2.12a,b)\begin{equation} \frac{\text{d} x}{\text{d} s} = \cos\varphi, \quad \frac{\text{d} z}{\text{d} s} = \sin\varphi. \end{equation}

To close the set of equations, the plume angle $\varphi$ is conventionally assumed to be the angle of the local streamlines $\varphi _u$, defined via

(2.13)\begin{equation} \tan \varphi_u = \frac{w_m}{U+u_m}. \end{equation}

Under the assumption that $\varphi = \varphi _u$, $U_s$ represents the velocity magnitude $\sqrt {(U\!+\!u_m)^2 \!+\! w_m^2}$, which can then be equivalently expressed in terms of the integral quantities as $Q^{-1} \sqrt {(Q U+M_x)^2 + M_z^2}$. Here, we make an explicit difference between $\varphi$ and $\varphi _u$ since it will turn out that these quantities are not identical for the flow cases under consideration (see § 4). When $\varphi \ne \varphi _u$, there is a component of velocity in the $\eta$ direction:

(2.14)\begin{equation} U_\eta ={-}(U + u_m) \sin \varphi + w_m \cos \varphi. \end{equation}

2.3. Bent-over plumes

In this section we briefly provide an overview of analytical results for bent-over plumes in a neutral stratification and a uniform crossflow. Bent-over plumes are defined to have their centreline inclined at a very small angle to the horizontal (in which case the equations are formally equal to those for a line thermal; Briggs Reference Briggs1984; Lee & Chu Reference Lee and Chu2003). This implies that $w_m \ll U$, $u_m \ll U$ and $\text {d} s \approx \text {d} x$. Assuming that in this case the entrainment rate is a constant, $\gamma = \beta '$, the following relations can be derived (e.g. Weil Reference Weil1988). Integrating the volume flux equation (2.11a) yields

(2.15)\begin{equation} r_m = \beta' z_c + r_0, \end{equation}

where $r_0$ is the source radius. Integrating the vertical momentum flux equation (2.11c) twice and making use of (2.12a,b), (2.13) and (2.15), assuming $\varphi =\varphi _u$, results in

(2.16)\begin{equation} \frac{z_c}{\ell_b} = \left(\frac{3}{2\beta^2} \right)^{1/3} \left( \frac{x - x_v}{\ell_b} \right)^{2/3} \end{equation}

for $r_0 \ll \beta ' z_c$, where $\ell _b = F/U^3$ is the buoyancy length scale and $x_v$ is a virtual origin correction, which is a constant subtracted from the $x$ coordinate of the plume trajectory that compensates for any near-field effects. It is necessary because the region close to the source (where $r_0 \ll \beta ' z_c$ no longer holds) is not expected to obey the two-thirds law. The value of $\beta$ is not equal to $\beta '$, because of added-mass effects (Lee & Chu Reference Lee and Chu2003), and the two can be related as $\beta = \beta '\sqrt {1+k}$ (Weil Reference Weil1988, p. 164; see also § 8), where $k$ is the same added-mass coefficient as in (2.11a–d). Experiments indicate that $0.3 \lesssim k \lesssim 1.3$ (e.g. Briggs Reference Briggs1982; Weil Reference Weil1988; Contini et al. Reference Contini, Donateo, Cesari and Robins2011) and $0.45 \lesssim \beta ' \lesssim 0.75$ (e.g. Hewett, Fay & Hoult Reference Hewett, Fay and Hoult1971; Hoult & Weil Reference Hoult and Weil1972; Briggs Reference Briggs1984; Weil Reference Weil1988; Huq & Stewart Reference Huq and Stewart1996). Although $\beta '$ is a constant coefficient, it may vary between simulations or experiments dependent on factors such as the strength of the crossflow and the initial buoyancy flux $F_0$ (Briggs Reference Briggs1984).

Several entrainment models are in use for plumes in a crossflow (Costa et al. Reference Costa2016). The classical model for a bent-over plume (e.g. Briggs Reference Briggs1984; Weil Reference Weil1988) assumes $\gamma = \beta '$. An entrainment model which is also valid for weak-crossflow plumes is given by (Devenish et al. Reference Devenish, Rooney, Webster and Thomson2010, DRWT)

(2.17)\begin{equation} \gamma = \frac{\sqrt[n]{(\alpha w_m)^n + (\beta' U )^n}}{U_s}. \end{equation}

The DRWT model uses two parameters, $\alpha$ and $\beta '$, to account for entrainment due to vertical and horizontal motion, respectively, and assumes that the added-mass coefficient $k=0$. Here, $\beta '$ refers to the same bent-over entrainment coefficient as the one that is used in (2.15). A further parameter $n$ determines the function used to combine the entrainment resulting from the vertical and horizontal momentum into a single entrainment quantity. In their paper, Devenish et al. (Reference Devenish, Rooney, Webster and Thomson2010) test the values $n = \{1, 1.5, 2\}$ and conclude that $n=1.5$ provides the best match between the model and their LES data.

3.1. Simulation details

The simulations concern the evolution of a plume in an unstratified environment with a uniform crossflow speed $U$. A constant integral buoyancy flux $F_0$ is imposed on a circular area source of diameter $2 r_0$, which is positioned at $x=0$. The source has no initial velocity, which implies the plume is infinitely lazy at the source (Hunt & Kaye Reference Hunt and Kaye2005). The buoyancy flux is therefore imposed as a diffusive flux using a Neumann boundary condition, which will be transported through a very thin thermal boundary layer before it becomes convectively unstable and starts rising as a conventional plume.

The relevant non-dimensional parameters for buoyant plumes in a crossflow can be formed from the exit velocity $w_0=M_{0}^{1/2}/r_0$, the buoyancy velocity scale $U_b=(F_0 / r_0)^{1/3}$, the crossflow speed $U$, and for completeness the viscous and thermal velocity scales $\nu / r_0$ and $\kappa / r_0$, respectively. All quantities with subscript zero represent their respective quantities at the source (e.g. $M_{0}$ represents $M$ evaluated at the source). From these five velocity scales, four dimensionless quantities can be constructed:

(3.1a–d)\begin{equation} R_0=\frac{w_0}{U}, \quad Ri_0 = \frac{F_0}{r_0 w_{0}^3} = \frac{b_{0} r_{0}}{w_{0}^2}, \quad Re_0 = \frac{2F_0^{1/3}r_0^{2/3}}{\nu}, \quad Pr = \frac{\nu}{\kappa}, \end{equation}

where $R_0$ is the jet-to-crossflow speed ratio, $Ri_0$ is the source Richardson number, $Re_0$ is the source Reynolds number and $Pr$ is the Prandtl number. Note that the Froude number, which is also commonly used, is the square root of the reciprocal of $Ri_0$. These dimensionless groups are not ideal for the current simulation set-up since $w_0=0$ which implies that $w_0/U=0$ and $Ri_0 = \infty$. In this study, we thus use the crossflow Richardson number $Ri_U$, defined as

(3.2)\begin{equation} Ri_U = Ri_0 R_0^3 = \frac{F_0}{U^3 r_0} = \frac{U_b^3}{U^3}, \end{equation}

which is the cube of the ratio of the buoyancy velocity scale $U_b$ to the crossflow speed $U$. Equivalently, using the buoyancy length scale $\ell _b=F_0/U^3$, $Ri_U$ can be interpreted as a ratio of length scales, $Ri_U = \ell _b / r_0$.

Five simulations S1–S5 have been performed at different values of $Ri_U$, details of which are presented in table 1. Simulation S1 reproduces a case characterised by a weak crossflow, with $Ri_U=8$, which corresponds to a flow in which the crossflow speed $U$ is equal to half the buoyancy velocity scale $U_b$. Simulation S2, with $Ri_U=2.4$, also represents a weak crossflow. Simulations S4 and S5, with $Ri_U=0.5$ and $Ri_U=0.3$, respectively, correspond to strong crossflow cases for which $U> U_b$. Simulation S3 has $Ri_U=1$, which implies that $U = U_b$. All simulations use an identical source Reynolds number $Re_0 = 1000$ and Prandtl number $Pr=1$. This value is chosen because it is the highest value at which the simulations may currently be conducted, due to computational cost. This value of $Re_0$ does ensure that the flow is fully turbulent.

Table 1. Simulation details. The source Reynolds number $Re_0=1000$ and Prandtl number $Pr=1$ for all simulations.

Due to the fact that the strength of the crossflow affects the evolution of the plume, the domain sizes were chosen large enough such that each plume can evolve inside the domain without interacting with its boundaries, but tight enough not to waste computational resources on flow regions where only ambient flow is present. The domain dimensions are given in table 1. Note that the $x$ coordinate ranges from $x/r_0 = -5$ to $x/r_0 = L_x/r_0-5$. For all simulations, the percentage of the domain occupied by the plume in the $(y, z)$ plane is less than $9\,\%$. This is sufficient to ensure that the interactions with the boundaries are negligible. For a row of plumes in a quiescent environment, the mean plume boundaries plotted by Rooney (Reference Rooney2015) show that a plume which spans ${\sim }25\,\%$ of the domain is not significantly affected by the competing entrainment of the other plumes. The grid resolution was chosen such that $\Delta x/\eta _K$ peaks at just below 3 close to the source for all simulations, where $\eta _K = (\nu ^3/\varepsilon )^{1/4}$ is the Kolmogorov length scale and $\Delta x = L_x / N_x$ is the grid resolution. The same grid resolution is used in the transverse and vertical directions, so $\Delta x$ is also equal to $L_y/N_y$ and $L_z/N_z$. The ratio $\Delta x/\eta _K$ inside the plume decays rapidly as a function of $x$; it is below 2 within $4r_0$ downstream of the source for all simulations, indicating that the grid resolution is sufficiently high to be considered DNS. Indeed, the peak of the dissipation rate takes place at a scale of ${\sim }24\eta _K$ for isotropic turbulence (Pope Reference Pope2000).

The simulations were performed using the DNS code SPARKLE (Craske & van Reeuwijk Reference Craske and van Reeuwijk2015), which integrates (2.1) on a cuboidal domain and is fully parallelised using domain decomposition in two directions. The spatial differential operator employs a fourth-order symmetry-preserving central difference scheme. Incompressibility was enforced by taking the divergence of (2.1) and solving the resulting Poisson equation for $p$ by performing fast Fourier transforms in the lateral directions. A third-order Adams–Bashforth method is used for the time integration.

Periodic boundary conditions were applied in the lateral directions and free-slip boundary conditions were applied at $z = 0$ and $z = L_z$. Zero-flux Neumann boundary conditions were used for buoyancy at $z = 0$ and $z = L_z$, with the exception of a circular region of radius $r_0$, centred at $(0,0,0)$, through which a constant buoyancy flux was imposed. In order to enforce a uniform inflow boundary condition at $x/r_0 = -5$ on our periodic domain we used a nudging region of length $L_n$ at the end of the domain over which the fluid velocity and buoyancy were gradually adjusted to become those of the environment (Stevens, Graham & Meneveau Reference Stevens, Graham and Meneveau2014). Denoting the distance from the beginning of the nudging region as $x^*$, we adjusted the velocity according to $\boldsymbol {u}^* = (1 - x^*/L_n) \boldsymbol {u} + (x^*/L_n) \boldsymbol {u}_a$, where $\boldsymbol {u}$ is the original DNS velocity and $\boldsymbol {u}_a = (U,0,0)$ is the ambient environmental velocity. The same linear transition method is used for the temperature. This procedure was implemented immediately prior to applying the pressure correction to maintain incompressibility. For all simulations the nudging length was set to $L_n/r_0=4$. Careful analysis of the simulation results was carried out to ensure that the nudging does not affect the statistics upstream, which for an incompressible flow could occur via pressure gradients generated by the nudging. No pressure gradients were observed near the nudging region, implying that it does not contaminate the upstream statistics.

3.2. Plume identification

The determination of all integral quantities requires identification of the plume fluid. This was performed by applying a threshold on the average buoyancy $\bar {b}$ with a value of 1 % of the maximum mean buoyancy at fixed $x$. Examples of the identified plume boundary are shown with a red line in figure 3.

3.3. Budget calculations

For the presentation of the budgets (2.6a–c) that are discussed in §§ 5 and 6, the entrainment terms $E$ and $\mathcal {L}$ need to be determined. Calculating these directly from the definition in (2.5a,b) can be challenging, since it requires a boundary integral in the plane perpendicular to the plume, which in turn requires defining the computational cells which constitute the boundary, taking account of any disconnected regions, defining normal vectors and taking into account the fact that the underlying grid is staggered. It is computationally more convenient to use the divergence theorem in reverse to calculate $E[\phi ]$ via

(3.3)\begin{equation} E[\phi] ={-}\displaystyle\iint_\varOmega \frac{\partial \overline{v \phi}}{\partial y} + \frac{\partial \overline{u_\eta \phi}}{\partial \eta}\,\text{d} A \end{equation}

and to use the original definition of $\mathcal {L}[\phi ]$ (van Reeuwijk et al. Reference van Reeuwijk, Vassilicos and Craske2021):

(3.4)\begin{equation} \mathcal{L[\phi]} = \frac{\text{d}}{\text{d} s}\displaystyle\iint_\varOmega \overline{u_s\phi} \,\text{d} A -\displaystyle\iint_\varOmega\frac{\partial \overline{u_s\phi}}{\partial s}\,\text{d} A. \end{equation}

4.1. General plume behaviour

Instantaneous and time-averaged snapshots of the buoyancy field $b$ in the $(x, z)$ plane are shown in figure 2, normalised by $b_{max}(x)=\max _{y,z} \bar {b}$. The plume, which is infinitely lazy at the source, must first diffuse through the thermal boundary layer at the wall, during which time it is advected downstream by the crossflow. Once the plume has traversed the thermal boundary layer (which typically occurs between 0.5 and 2 source radii downstream depending on the crossflow speed), the plume accelerates upwards, causing some necking in the $y$ direction. Following this, the plume angle $\varphi$ gradually decreases as the plume rises through the ambient, and the plume radius increases due to turbulent entrainment. This increase is clearly dependent on the flow speed: the plumes are observed to be narrower for higher crossflow speeds. Another notable feature of these snapshots is that there is a tendency for a small amount of buoyancy to be swept into the wake of the plume. The buoyancy detrains from the plume close to the source as a result of the initial transition instabilities, and remains close to the $z=0$ boundary. This requires some careful attention when it comes to determining the plume boundary because in general it is inappropriate to regard this region, which may have a small amount of buoyancy but is otherwise quiescent, as being within the plume. For clarity we have removed the buoyancy concentration that would be in these quiescent regions from the figure. Cross-sections in the $(y, \eta )$ plane at $s/r_0 = 5$ are shown in figure 3. These clearly show a buoyant core structure with twin vortices, which becomes narrower and less diffuse as the crossflow speed increases.

Figure 2. Cross-section at the $y=0$ centre plane of the instantaneous (a–c) and time-averaged (d–f) buoyancy field. (a,d) Simulation S1. (b,e) Simulation S3. (c,f) Simulation S5.

Figure 3. The buoyancy field in the $(y, \eta )$ plane at $s/r_0 = 5$, for the (a–c) instantaneous buoyancy and (d–f) time-averaged buoyancy. (a,d) Simulation S1. (b,e) Simulation S3. (c,f) Simulation S5. The red line in (d–f) is the applied threshold on buoyancy at 1 % of the maximum of the time-averaged value.

4.2. Plume centreline

The plume centreline is a fundamental quantity in integral plume models. Theoretically, it is typically represented as the central streamline of the plume (e.g. Weil Reference Weil1988). This is relatively straightforward to do using numerical simulation when the full velocity field is known (Yuan & Street Reference Yuan and Street1998; Yuan et al. Reference Yuan, Street and Ferziger1999; Muppidi & Mahesh Reference Muppidi and Mahesh2005; De Wit et al. Reference De Wit, Van Rhee and Keetels2014; Cintolesi et al. Reference Cintolesi, Petronio and Armenio2019) via a streamline that starts at the centre of the source. Most laboratory studies, but also some numerical studies (Devenish et al. Reference Devenish, Rooney, Webster and Thomson2010), use either the maximum velocity or buoyancy or use the centre of mass of the buoyancy/passive scalar or velocity (Contini et al. Reference Contini, Donateo, Cesari and Robins2011).

Shown in figure 4 are a number of centreline definitions as a function of $x$, plotted together with isolines of the buoyancy averaged over the plume width (in $y$) for simulation S3. The centre of mass was calculated according to

(4.1)\begin{equation} z_c(x) = \dfrac{\displaystyle\iint_\varOmega z \bar{b}\, \text{d} y \,\text{d} z}{\displaystyle\iint_\varOmega \bar{b}\, \text{d} y\, \text{d} z}. \end{equation}

Here, integration is performed over the plume fluid area $\varOmega$. Extracting the central streamline $z_U$ required some careful consideration as the plume is infinitely lazy at its source, implying that the velocities are zero on the boundary. This makes it impossible to start a streamline in the centre of the source. Instead, we let the streamline start at $z_c(x=2r_0)$, which is where the plume lifts off from the surface. As can be seen, all the other plume trajectory indicators are practically identical at that stage. The maximum buoyancy, velocity and turbulent kinetic energy were determined by taking the mean value across the plume width (in the $y$ direction), and then finding the maximum value for each $x$. By averaging over the plume width, the centreline represents the entire plume area. Finally, a trajectory is shown that is based on the centre of mass of the velocity $V=\sqrt {\bar {u}^2+\bar {w}^2}$, using (4.1) but substituting $\bar {b}$ with $V$.

Figure 4. Various estimates for the plume centreline, based on simulation S3 (tke, turbulent kinetic energy). The isocontours denote the average buoyancy across the plume width.

Near $x=0$, some discontinuities can be seen in the buoyancy isocontours which are associated with the plume lift-off from the ground and that will not be considered in the analysis in the remainder of the paper. The central streamline $z_U$ and centre of mass $z_c$ evolve very similarly, although $z_U$ has a slightly higher trajectory initially and then follows a slightly smaller slope further downwind. The maxima of the buoyancy and velocity are very similar, and follow a trajectory which is lower than the centres of mass and the maximum turbulent kinetic energy. Figure 2 suggests that the maximum of buoyancy will roughly align with the vortex centres of the double-roll structures, where buoyancy accumulates. Unsurprisingly, the centres of mass change much more smoothly than the centrelines based on finding a maximum. The centre of mass trajectory based on $V$ is slightly higher than that of $\bar {b}$, but this is a small difference compared with the difference between $z_c$ and $z_{U}$.

4.3. Plume angle

Having calculated the integral plume quantities, we verify whether the plume slope $\varphi$ is identical to the mean velocity streamline angle $\varphi _u$ (2.13), which would imply that

(4.2)\begin{equation} \frac{\text{d} z_c}{\text{d} x} = \frac{w_m}{U + u_m} = \frac{M_z}{QU + M_x}. \end{equation}

This relation is plotted in figure 5 for simulation S3, together with the gradient of $z_c$ and $z_U$. Figure 5 shows that $w_m/(U + u_m) = M_z/(QU + M_x)$ is nearly half of $\text {d} z_c / \text {d} x$ and $\text {d} z_U / \text {d} x$, and therefore implies that the ratio of the integral momentum fluxes is a poor estimate of the plume slope. Note that $z_c$ was approximated by two cubic splines in order to avoid numerical errors resulting from repeated differentiation (see Jordan (Reference Jordan2021) for details). Simulation S3 is representative of the other simulations. Figure 5 shows that the streamline angle $w_m/(U + u_m)$ is not identical to the plume angle for $z_c$ or $z_{U}$. The assumption that the plume centre of mass follows a streamline dictated by $w_m$ and $U+u_m$ is central to the integral theory of plumes in a crossflow. Therefore, the implications of the observed discrepancy are substantial.

Figure 5. Slope of the plume centreline for simulation S3, together with several integral-quantity estimates.

Figure 6(a) shows the velocity streamlines through the plume for simulation S3. The streamlines were constructed from the velocity inside the plume, by integrating the velocity over the $y$ direction within the plume and then dividing through by the local plume width, thereby obtaining a velocity field representative across the width of the plume. Figure 6(a) is consistent with figure 5 in demonstrating that the angle of the streamlines is substantially different from the plume centreline (thick black line), for both $z_c$ and $z_{U}$. The reason that $z_U$ is not aligned with the mean streamlines across the plume is the double-roll structure which causes positive vertical velocities on the centreplane and negative velocities away from the centreplane. For the case under consideration (S3), the streamlines enter the top of the plume (indicated by a dashed line) at nearly 45$^{\circ }$, but have much smaller angles near the bottom boundary. However, neither the centre of mass of the plume $z_c$ (thick black line) nor the central streamline $z_{U}$ (thick dashed black line) is parallel to the streamlines. Thus, the plume angle $\varphi$ is different from the velocity angle $\varphi _u$, which implies that the conventional assumption of the plume aligning with the central velocity streamline ($\varphi =\varphi _u$) is violated.

Figure 6. Streamlines (blue) for simulation S3, together with $z_c$ (thick black line), $z_{U}$ (thick dotted line), the plume boundaries (dashed line) and the $y$-averaged buoyancy field (in grey). Streamline quantities are averaged in $y$ across the local plume width. (a) Streamlines of $\boldsymbol {u} = (\bar {u},\bar {w})$. (b) Field lines of $\bar {\boldsymbol {f}} = (\bar {u}\bar {b}, \bar {w}\bar {b})$. (c) Field lines of $\boldsymbol {f} = (\bar {u}\bar {b}+\overline {u'b'}, \bar {w}\bar {b}+\overline {w'b'})$.

Since $z_c$ was determined from the centre of mass of the buoyancy of the plume, it stands to reason to investigate whether it is more suitable to consider field lines based on the local mean buoyancy flux $\bar {\boldsymbol {f}} = (\bar {u}\bar {b}, \bar {w}\bar {b})$, which are shown in figure 6(b). Similar to the mean velocity, the mean buoyancy flux was determined by integrating the individual components over $y$ and dividing by the local plume width. Although the mean buoyancy field lines follow the bottom boundary almost perfectly, it is clear that the behaviour is nearly identical to that of figure 6(a).

It is only when turning to the field lines of the total buoyancy flux $\boldsymbol {f} = (\bar {u}\bar {b}+\overline {u'b'}, \bar {w}\bar {b}+\overline {w'b'})$ that the field lines are found to be fully representative of $z_c$. The central streamline $z_{U}$ also correlates well with these field lines, even when $z_U$ is based on the velocity field only. A secondary conclusion that can be drawn from comparing figures 6(b) and 6(c) is that the region below the plume centreline is dominated by mean buoyancy transport and the top by a combination of mean and turbulent buoyancy fluxes. This is consistent with previous observations (Huq & Stewart Reference Huq and Stewart1997). Furthermore this image suggests that the plume is entraining strongly at the top. It is likely that the high degree of turbulence towards the top of the plume and the impinging crossflow contribute to the strong entrainment there. The fact that the total buoyancy transport determines the angle of the plume is verified in figure 5, which is shown to follow $\text {d} z_c / \text {d} x$ very closely.

Similar conclusions can be drawn from the integral fluxes by defining decomposed fluxes for $F$ and $F'$:

(4.3a,b)\begin{gather} F_x = \frac{1}{{\rm \pi}} \displaystyle\iint_\varOmega \bar{u} \bar{b}\, \text{d} A, \quad F_z = \frac{1}{{\rm \pi}} \displaystyle\iint_\varOmega \bar{w} \bar{b}\, \text{d} A, \end{gather}

(4.4a,b)\begin{gather}F'_x = \frac{1}{{\rm \pi}} \displaystyle\iint_\varOmega \overline{u'b'}\, \text{d} A, \quad F'_z = \frac{1}{{\rm \pi}} \displaystyle\iint_\varOmega \overline{w'b'} \,\text{d} A. \end{gather}

Shown in figure 5 is $F_z/F_x$, which is the approximation of $\text {d} z_c / \text {d} x$ based on mean buoyancy fluxes in the horizontal and vertical directions. It is clear that this result is very similar to that given by the ${w_m}/({U + u_m})$ estimate for $\text {d} z_c / \text {d} x$, consistent with figure 6(a,b). Defining the slope based on the total buoyancy flux as $\text {d} z_c / \text {d} x = (F_z + F_z')/(F_x+F_x')$, consistent with figure 6(c), produces an accurate estimate of the plume slope $\text {d} z_c / \text {d} x$. For this reason, $z_c$ is used to represent the plume trajectory from this point onwards. From the observation that $M_z/(UQ+M_x) \approx F_z/F_x$ (figure 5), the plume angle can be approximated as follows:

(4.5)\begin{equation} \frac{\text{d} z_c}{\text{d} x} \approx \frac{F_z + F'_z}{F_x + F_x'} \approx \tan \varphi_u + \frac{F'_z}{F_x} = \frac{\sin \varphi_u + \theta_f}{\cos \varphi_u}, \end{equation}

where $\theta _f=F'_z/F$. In arriving at the result above, use was made of (2.13), $F_x'\ll F_x$ (justified from analysing the DNS data) and using that $F_x/F \approx \cos \varphi _u$ since $F=\sqrt {F_x^2 + F_z^2}$. The ratio of the turbulent buoyancy flux to the mean buoyancy flux $F_z'/F$ was denoted $\theta _f$ in van Reeuwijk & Craske (Reference van Reeuwijk and Craske2015), which motivates its name in the current paper. An a priori calculation of $\text {d} z_c / \text {d} x$ using (4.5) is shown in figure 5 which shows good agreement with the observed plume angle. The comparison of $\text {d} z_c / \text {d} x$ with parameterised $\theta _f$ for all cases is shown in § 5.

4.4. Integral fluxes

The plume centrelines $z_c$ are shown as a function of $x/r_0$ in figure 7(a) for all simulations. In order to validate the simulations, we compare these centrelines with the analytical prediction (2.16) of a bent-over plume. Shown in figure 7(b) is $(z_c/\ell _b)^{3/2}$ against $x$. As can be observed, all plume trajectories evolve in a linear fashion, confirming scaling consistent with the expected bent-over plume behaviour. The advantage of plotting $(z_c/\ell _b)^{3/2}$ against $x$ in a linear plot, rather than $(z_c/\ell _b)$ against $x$ in a log–log plot, is that the power-law exponent extracted from a log–log plot can be strongly influenced by the virtual origin correction. Comparison with experimental and numerical data of plumes in crossflow is carried out based on the fit to (2.16), rather than displaying them in figure 7 (as the main difference would be the slope, which depends on $R_0$; Briggs Reference Briggs1984).

Figure 7. Plume trajectories. (a) Centreline $z_c$ as a function of $x$. (b) $z_c^{3/2}$ as a function of $x$. (c) Plume radius $r_m$ as a function of $z_c$.

By fitting a straight line through the last part of the trajectory for each simulation shown in figure 7(b), $\beta '$ and the virtual origin $x_v$ can be identified. Figure 7(c) shows the linear dependence between $r_m$ and $z_c$, conforming with the bent-over plume predictions (2.15). The values of $\beta$, $\beta '$, $k$ and $x_v/r_0$ are presented in table 1 for each simulation. They are calculated from (2.15) and (2.16) and from the observations in figure 7.

We find values of $\beta '$ in the range $0.49 \le \beta '\le 0.88$ which is higher than in the experiments of De Wit et al. (Reference De Wit, Van Rhee and Keetels2014), and the LES of Cintolesi et al. (Reference Cintolesi, Petronio and Armenio2019), suggesting the plumes considered here ascend less rapidly. We suspect these higher values of $\beta '$ are due to the fact that the plume under consideration is infinitely lazy at the source. Indeed, Briggs (Reference Briggs1984) developed an empirical prediction $\beta ' = 0.4 + {1.2}/{R_0}$ for the dependence of $\beta '$ on $R_0 = {w_0}/{U}$, demonstrating that this parameter depends on conditions at the source, particularly at low values of $R_0$.

Once the plume trajectory and the plume fluid have been identified, the integral volume flux $Q$, horizontal momentum excess flux $M_x$ and vertical momentum flux $M_z$ can be calculated from their definitions (2.9a–d). These quantities are shown in figure 8 as a function of $s$. The volume flux $Q$ grows increasingly fast as it evolves due to the plume's increasing size (see § 5). Contrary to the conventional assumptions, the horizontal momentum excess flux $M_x$ is not conserved for the plumes. Instead, it increases with $s$. This is caused by the pressure gradient term $P_x$ in (2.6b), which is non-zero. Simulation S1 is clearly of a different character from the other simulations, as the second derivative of $M_x$ is of opposite sign. The vertical momentum flux $M_z$ also grows as a function of $s$. Unlike $M_x$, this quantity is expected to be non-zero due to the influence of buoyancy in (2.6c). However, the influence of pressure is non-negligible for this quantity also. Consistent with figure 8(a,b), the larger the crossflow speed $U$, the smaller the gradient, but this can be partially explained by the way in which the figures are normalised. The physical interpretation of these effects on the momentum fluxes is that the pressure field, which is generated by the interaction of the crossflow with the leading edge of the plume, transfers some of the upward momentum generated by the buoyancy into horizontal momentum. This leads to $M_x$ being higher than expected, while $M_z$ is lower than expected. The role of pressure is investigated in further detail in § 6.

Figure 8. Integral plume quantities as a function of $s$ for each simulation. (a) Volume flux $Q$. (b) Horizontal momentum excess flux $M_x$. (c) Vertical momentum flux $M_z$. (d) Buoyancy flux $F$.

The mean buoyancy flux $F$ increases rapidly as measured from the centre of the source as it is advected by the mean flow, and the initially diffusive buoyancy flux is transferred to an advective flux, which is displayed here. The buoyancy flux is practically constant across all simulations within two source diameters of the source, and remains close to its initial value $F_0$. This is consistent with the buoyancy flux equation (2.11d). The value of the mean buoyancy flux $F_0$ differs from unity because a part of the buoyancy transport is due to turbulence. Turbulence transports a higher percentage of the buoyancy for weaker crossflows, reaching a limiting value of approximately 20 % in the absence of wind (van Reeuwijk et al. Reference van Reeuwijk, Salizzoni, Hunt and Craske2016).