Internal waves generated by a plume impinging on a stratified fluid

Abstract

Inspired by laboratory experiments showing internal waves generated by a plume impinging upon a stratified fluid layer (Ansong & Sutherland. 2010 J. Fluid Mech. 648, 405–434), we perform large eddy simulations in three dimensions to examine the structure and source of internal waves emanating from the top of a plume that rises vertically into stratification whose strength ranges over two orders of magnitude between different simulations. Provided the plume is sufficiently energetic to penetrate into the stratified layer, internal waves are generated with frequencies in a relatively narrow band moderately smaller than the buoyancy frequency. Through adaptations of ray theory including viscosity and use of dynamic mode decomposition, we show that the waves originate from within the turbulent flow rather than at the turbulent/non-turbulent interface between the fountain top and the surrounding stratified fluid.

1. Introduction

Internal (gravity) waves propagate horizontally and vertically through stratified fluid driven by buoyancy forces. Because they vertically transport horizontal momentum, they can significantly influence atmospheric winds and, consequently, weather and climate. The waves can be generated by a variety of processes, including flow over topography, frontogenesis and convective storms (Fritts & Nastrom 1992). Of the last of these, internal waves can be excited when the top of a convective system impinges upon the base of the stratosphere. In the absence of mean winds and heating within the body of the plume, it has been proposed that internal waves can be excited by the vertically fluctuating motion of the cloud tops (Fritts & Alexander 2003). This perspective was drawn into question by laboratory experiments examining internal wave generation by a vertical plume impinging upon a stratified fluid layer (Ansong & Sutherland (2010), henceforth AS10). After penetrating into the stratified layer, the plume became negatively buoyant, transforming into a fountain. That study showed that the frequency spectrum of internal waves emanating from the plume/fountain was narrow and not related to the broad-banded spectrum associated with fluctuations of the turbulent/non-turbulent interface of the fountain top (corresponding with the cloud top in Fritts & Alexander (2003)).

Here we perform large eddy simulations (LES) of a plume impinging upon a stably stratified fluid layer in order to gain insight into the mechanism for generation of narrow-banded internal waves from broad-banded convective turbulent motion. Turbulent mixing and entrainment properties of this flow have previously been considered in Powell, Haynes & Taylor (2024). The numerical model is detailed in § 2. In § 3 we describe the flow and compare the generated wavefield observed in a set of simulations in which the squared buoyancy frequency of the stratified layer varies over two orders of magnitude. In § 4 we focus on a high-resolution simulation and present analyses which demonstrate that waves originate from within the fountain top rather than at the turbulent/non-turbulent interface.

2. Numerical model

We consider LES of a pure plume with source radius $r_0$ generated near the bottom ( $z=0$ ) and at the horizontal centre ( $x,y=L_H/2$ ) of a domain of height $L_z$ and horizontal extent $L_H$ , as shown in figure 1. Scales are chosen to be similar to those of the laboratory experiments of AS10, such that $r_0=0.005\,\rm m $ , $L_z=0.6\, \rm m$ and $L_H=1.0\, \rm m$ .

Figure 1. Schematic showing simulation set-up (a) and initial buoyancy profile $b_0(z)$ (b). The stratification begins at $H=0.2\, \rm m$ (black dashed line). The plume lies on the centreline $x=y=L_H/2$ of the domain of width $L_H = 1\, \rm m$ and height $L_z = 0.6\, \rm m$ . The forcing (i.e. where $f_w, f_b$ in (2.2), (2.3) are non-zero) occurs below the blue dashed line. Sponge layers are shaded in grey. Internal waves are indicated by blue wavy lines with wavevector $\boldsymbol {k}$ as shown. The flow within the plume is indicated by solid black arrows and structures referred to in the text are labelled in red. The maximum penetration height $z_{max }$ and quasi-steady-state height $z_{{ss}}$ , measured from the bottom of the stratified layer, are indicated by dotted lines.

Table 1. Simulation and time window parameters. Here $N_x$ , $N_y$ and $N_z$ are the number of grid cells in each direction.

At the source, the generated plume has integral buoyancy flux $F_0 = 2\int _0^\infty \left . \langle w\rangle \langle b \rangle \right |_{z=0}\,r\mathrm {d}r$ , where $\langle \cdot \rangle$ denotes a time and azimuthal average, $w$ is the vertical velocity and $b$ is the buoyancy. In all simulations we take $F_0=5\times 10^{-6}\,\mathrm {m}^4\,\mathrm {s}^{-3}$ . The plume carries a passive tracer which is used to mark the extent of the plume throughout each simulation. The density is uniform over the bottom $H=0.2\,\rm m$ of the domain and uniformly stratified above with squared buoyancy frequency, $N_0^2$ , that ranges from $1$ to $100\,\mathrm {s}^{-2}$ in moderate-resolution simulations and is $0.25\,\mathrm {s}^{-2}$ in a high-resolution simulation. Relevant simulation parameters are given in table 1. During penetration into the stably stratified layer, the plume transforms into a fountain that excites internal waves above, with the collapsing fluid from the fountain spreading laterally as an intrusion at its neutral buoyancy level. To inhibit reflection of internal waves from the top boundary, we include a sponge layer of depth $L_S = 0.1\, \rm m$ at the top of the domain (well above the fountain top in each of the simulations) where the velocity is damped towards zero. We also include a sponge layer of width $L_S$ on the four sides of the horizontally periodic domain to prevent low-frequency internal waves with a shallow angle from wrapping around the computational domain.

The simulations use the numerical method reported in Powell et al. (2024) which is briefly summarised here. The Boussinesq Navier–Stokes equations are integrated numerically using Fourier modes in the two periodic horizontal directions and second-order finite differences in the vertical. The boundary conditions on the top and bottom boundaries are no-stress and no-penetration. The full, filtered governing equations non-dimensionalised by $F_0$ and $N_0$ and including subgrid-scale (SGS) contributions are

(2.1) \begin{align} \nabla \cdot \widehat {\boldsymbol {u}} &= 0, \end{align}

(2.2) \begin{align} \frac {\mathrm {D} \widehat {\boldsymbol {u}}}{\mathrm {D}t} + \nabla \widehat {p} &= \frac {1}{{Re}} \nabla ^2 \widehat {\boldsymbol {u}} + \widehat {b} \hat {\boldsymbol {k}} - \nabla \cdot \unicode{x1D7BD} + f_w\hat {\boldsymbol {k}}, \end{align}

(2.3) \begin{align} \frac {\mathrm {D}\widehat {b}}{\mathrm {D}t} &= \frac {1}{{Re}Pr } \nabla ^2 \widehat {b} - \nabla \cdot \boldsymbol {\lambda }_b + f_b, \end{align}

where $\widehat {(\cdot )}$ indicates filtering at the resolved grid scale and $\hat {\boldsymbol {k}}$ is the unit vector in the vertical direction. The terms $f_w$ and $f_b$ represent the forcing applied to generate the buoyant plume, detailed in Powell et al. (2024). The SGS stress tensor $\unicode{x1D7BD}$ and SGS scalar flux $\boldsymbol {\lambda }_b$ are determined by the anisotropic minimum dissipation model, as detailed in Vreugdenhil & Taylor (2018). The two dimensionless parameters are the Reynolds number ${Re} = F_0^{1/2}/\left(\nu N_0^{1/2}\right)$ and the Prandtl number $Pr = \nu /\kappa =0.7$ , where $\nu$ is the molecular viscosity and $\kappa$ is the molecular diffusivity. Note that the wavelength of the dominant internal waves that we analyse here is always much larger than the grid spacing. In all regions of the domain, the LES captures the energy-containing scales, and the smallest resolved motions contain comparatively little energy. We therefore conclude that using LES to analyse the internal wavefield is justified. For details on verification of the numerical scheme, including a resolution sensitivity test that showed the turbulent flow is well represented at the resolution(s) used in this paper, we refer the reader to Powell et al. (2024).

All simulations are first run until the plume front reaches the base of the stratified layer, which we define to occur at time $t=0$ . Before beginning our analyses, the wavefield is allowed to develop until time $t_1 T_b$ , where $T_b = 2\pi /N_0$ is the buoyancy period. We then collect data at a fine temporal resolution $\Delta _t T_b$ until time $t_2 T_b$ . Values of the non-dimensional quantities $t_1$ , $t_2$ and $\Delta _t$ for various simulations are listed in table 1. Unless otherwise stated, analyses are performed on a spatial window excluding the sponge layers and focused on the stratified layer: $L_S \le x,y \le L_H-L_S$ and $H \le z \le L_z-L_S$ . Throughout this paper, perturbation components are calculated by subtracting a running mean of the azimuthally averaged field over one buoyancy period $T_b$ . The plume volume flux is small and the large-scale flow reaches a quasi-steady-state in which the background stratification and flow vary more slowly than the averaging period, motivating this choice. Practical constraints mean that azimuthal averages are stored for the buoyancy and velocity components only. Analyses are therefore restricted to using horizontal slices. Horizontal averages are weighted by the radial distance from the plume centreline. Finally, note that whilst (2.1)–(2.3) are stated in non-dimensional form, we state all values in dimensional units henceforth to avoid confusion when varying $N_0$ .

3. Comparison of plume and wavefield evolution

In each simulation, the plume rises through the uniform layer and penetrates into the stratified layer after approximately $3$ s. The plume overturns at its maximum penetration height above the bottom of the stratified layer, $z_{max } \propto F_0^{1/4} N_0^{-3/4}$ , once it becomes relatively less buoyant than its surroundings, before transforming into a fountain and forming a radially spreading intrusion at the neutral buoyancy height. The scaling for $z_{max }$ follows from dimensional analysis, has been verified in experiments (Briggs 1965), and fits our simulation data with a proportionality constant $4.4 \pm 0.2$ . A quasi-steady-state is reached in which fluid is continuously supplied by the plume to the stratified layer and spreads as an intrusion after mixing with the buoyant environment near the fountain top (Powell et al. 2024), with mean steady-state height $z_{ss}$ below $z_{max }$ . We ensure that this quasi-steady-state is reached before time $t_1 T_b$ when fine-resolution data collection begins. As noted in experimental studies, internal waves are not observed during rise to the maximum penetration height (AS10), instead appearing once plume fluid first overturns.

Figure 2. Instantaneous $x$ - $z$ slices of $w'$ showing the internal wavefield in the stratified layer in simulations HR, N0, N1 and N2 at $t=5\,T_b$ . Horizontal dotted and dash-dotted lines indicates the height at which spectra are calculated in figure 4. The passive tracer contour in black indicates the extent of the plume.

Figure 2 shows instantaneous $x$ - $z$ slices through the plume centreline of the perturbation vertical velocity $w'$ for all simulations and figure 3 shows horizontal slices at $z=0.5\,\rm m$ . Snapshots are taken at $t= 5\,T_b$ , during quasi-steady-state. A contour of the passive tracer field is shown to indicate the extent of the plume. The internal wavefield is evident, with coherent wave beams propagating upwards and outwards in the ambient fluid above the fountain top. The simulated internal waves are consistent with results presented in AS10. Figure 4(a) shows the time-averaged vertical energy flux, $F_{{wave}} = \int \langle w p \rangle \, \mathrm {d}A$ , computed from a Fourier–Bessel decomposition of $\langle w\rangle$ as detailed in AS10. Our results with $N_0^2 = 1$ (closest to the value used by AS10) show that $F_{{wave}} \approx O(10^{-7})$ throughout the stratified layer, consistent with AS10 results. The dotted line shows the theoretical scaling for $F_{{wave}}$ derived by Couston et al. (2018) for waves generated by Reynolds stresses due to eddies in an (isotropically) turbulent region below a very strongly stratified layer. Whilst there is poor agreement in the weakly stratified cases, the $N_0^2 = 100$ case matches the predicted scaling reasonably well in the far field.

Figure 3. Instantaneous horizontal slices of $w'$ showing concentric rings of the internal wavefield at $z=0.5\, \rm m$ in simulations HR, N0, N1 and N2 at $t=5\,T_b$ .

Figure 4. Analyses of simulations HR, N0, N1 and N2 showing (a) time-averaged vertical energy flux, $F_{{wave}}$ , compared with the theoretical prediction of Couston et al. (2018) for a strongly stratified regime with a stiff interface: $F_{{wave}} \sim (z-H)^{-13/8}$ (black dotted line), (b) total energy $\sum _{k_h, \omega } E \delta k_h \delta \omega$ at $(z-H)/z_{max } = 0.25, 1.4$ (crosses, circles), (c) normalised horizontal wavenumber spectrum $f_{k_h}$ and (d) normalised frequency spectrum $f_{\omega }$ at $(z-H)/z_{max }=1.4$ . Panel (e) compares the characteristic wave frequency $\omega _c$ and the plume forcing frequency $\omega _{{plume}}$ .

Frequency and horizontal wavenumber spectra are calculated from the (kinetic) energy density $E(k_h, \omega ; z) = \ ({1}/{2})\sum _i \lvert A^{(u^{\prime}_i)} \rvert ^2$ , where $A^{(u^{\prime}_i)}(\omega , k_h; z)$ are the 2-D fast Fourier transform (FFT) amplitudes of each perturbation velocity component $u^{\prime}_i$ at height $z$ . We apply an energy-corrected Hann window before computing the time FFT. The frequency spectrum is then $f_{\omega }(\omega ; z) = \sum _{k_h} E \delta k_h$ and the horizontal wavenumber spectrum is $f_{k_h}(k_h; z) = \sum _{\omega } E \delta \omega$ where $\delta k_h$ and $\delta \omega$ are the spacings in spectral space. Figure 4(b) shows the total energy $\sum _{k_h, \omega } E \delta k_h \delta \omega$ above (within) the plume, at $(z-H)/z_{max }=1.4$ ( $0.25)$ , indicated by crosses (circles). The two heights are indicated by the dotted and dot-dashed lines in figure 2. The energy at heights inside the turbulent plume is approximately two orders of magnitude larger than in the internal wavefield above the plume. The separation in energy scales is similar for $N_0^2 = 0.25, 1, 10$ and decreases when $N_0^2 = 100$ .

Figures 4(c) and 4(d) show $f_{k_h}$ and $f_\omega$ , respectively, above the plume at $(z-H)/z_{max }=1.4$ . The characteristic horizontal wavenumber $k_{h,c}$ , calculated as a power-weighted average from $f_{k_h}$ , scales with $N_0$ as $k_{h,c} \propto N_0^{0.5 \pm 0.05}$ . Thus the horizontal wavelength of the wave beams decreases with $N_0$ , as is evident in figure 2. The frequency spectra in figure 4(d) suggest that the characteristic wave frequency $\omega _c$ , calculated as a power-weighted average, remains approximately constant as a fraction of $N_0$ with $0.6 \lt \omega _c/N_0 \lt 1$ . This is consistent with experiments of the same set-up by AS10 which found $0.4 \lt \omega _c/N_0 \lt 0.9$ . The characteristic wave frequency is not related to the plume forcing frequency. This is shown by comparing $\omega _c$ with the characteristic frequency $\omega _{{plume}}$ of vertical fluctuations of the plume height around the quasi-steady-state height $z_{ss}$ (not shown). This comparison is shown in figure 4(e) and is consistent with AS10, showing no clear relation between the approximately fixed $\omega _c/N_0$ and the varying $\omega _{{plume}}$ . The question therefore remains: What determines the wave frequency spectrum?

Figure 5. Frequency spectrum $f_{\omega }(\omega ; z)$ at a range of heights in simulation HR shown (a) at all heights on a log–log scale, (b) within the plume and (c) above the plume. Dots indicate the raw spectrum which is smoothed to give the solid lines. Note the different scale between (b) and (c). The heights at which spectra are calculated are indicated by dashed coloured lines in (d), an instantaneous $x$ - $z$ snapshot of $w'$ at $t=5T_b$ with the same colour bar as in figure 2. In (a), (b) and (c) the vertical dashed line indicates $\omega /\overline {N} = 1$ .

4. Analyses

We now focus on analyses of the high-resolution simulation with $N_0^2 = 0.25\,\mathrm {s}^{-2}$ , for which the fountain top reaches a steady-state height of $z_{ss}\simeq 0.38\,\rm m$ .

4.1. Spectral analysis

Figure 5 shows the frequency spectrum $f_{\omega }$ at a range of heights within and above the plume. The raw spectrum is shown as coloured dots and smoothed in frequency space to give the coloured lines. The frequency axis is normalised by the time and horizontally averaged stratification strength $\overline {N}(z)$ , where $\overline {N}(z)=N_0$ sufficiently far above the fountain (see § 4.2). There is a sharp increase in the stratification strength at the fountain top where intense buoyancy gradients are established between the plume and the more buoyant environment (Powell et al. 2024), whilst the stratification is weaker deep inside the plume. Internal waves can propagate locally where $0 \lt \omega /\overline {N}(z) \lt 1$ . Within the plume ( $z \le 0.38\, \rm m$ , figure 5 b), the frequency spectrum $f_\omega (\omega ; z)$ is broad and decays with increasing frequency, as expected for turbulent flow. Above the plume (figure 5 c), the spectrum is narrow and peaks at approximately $0.7 \le \omega /\overline {N} \le 1$ . This is consistent with the spectra in figure 4(b) and results in AS10. In figure 5(b) there is evidence of the spectrum forming a peak at frequencies close to, but larger than, $\overline {N}(z)$ at heights $z=0.32, 0.36\, \rm m$ in the fountain top.

Figure 6(a) shows the wavenumber spectrum $f_{k_h}(k_h, z)$ at a range of heights. The dashed line shows an isotropic turbulence spectrum $f_{k_h} \sim k_h^{-5/3}$ and the dot-dashed line shows a 2-D (axisymmetric) turbulence spectrum $f_{k_h} \sim k_h^{-3}$ . The axisymmetric scaling most closely matches the observed spectra within the plume, consistent with direct numerical simulation studies of plumes (e.g. van Reeuwijk et al. (2016)). Figures 6(b) and 6(c) show the energy (density) spectrum $E$ at $(z-H)/z_{max }=0.25$ and $1.4$ , respectively, clearly showing the broad structure within the plume and restriction to small wavenumbers and a narrow frequency band above the plume.

Figure 6. (a) Horizontal wavenumber spectrum $f_{k_h}(k_h,z)$ at a range of heights in simulation HR, compared with an isotropic and axisymmetric turbulence scaling $k_h^{-5/3}, k_h^{-3}$ , shown as dashed and dot-dashed black lines, respectively. Line colours as in figure 5. (b,c) Energy spectrum $E$ at $(z-H)/z_{max }=0.25, 1.4$ . The black dashed line indicates $\omega /\overline {N}(z) = 1$ .

4.2. Viscous internal wave model

Taylor & Sarkar (2007) (henceforth TS07) present a simple linear model to explain the selection of a dominant range of frequencies in the spectrum of internal waves generated by a turbulent Ekman layer. The method starts with known wave amplitudes $A^{(w')}(\omega , k_h; z_0)$ computed from $w'$ at some initial height $z_0$ . Henceforth, we drop the superscript and assume amplitudes are computed from $w'$ . For § 4.2 only, we use an amplitude-corrected Hann window instead of energy-corrected. The amplitudes $\tilde {A}(\omega , k_h; z)$ at an arbitrary height $z$ are calculated based on the expected vertical propagation speed and viscous decay rate, assuming that the waves satisfy the linear dispersion relation and that the background fields are slowly varying in space and time. Here, following TS07, the spectrum $P(\omega ; z)$ is calculated from amplitudes $A$ as $P(\omega ; z) = \sqrt {\sum _{k_h} A(\omega , k_h; z)^2}$ . Compared with the original study, we move from a Cartesian to an axisymmetric wave geometry, implicitly assuming that the curvature in the waves is small enough to be approximated as plane waves. We account for the energy-conserving amplitude decrease of a spreading axisymmetric wave beam from a virtual point source at height $z_s\le z_0$ with a factor $\sqrt {r(z_0)/r(z)}=\sqrt {(z_0-z_s)/(z-z_s)}$ where $r(z)$ is the along-beam distance from the virtual source to the waves at height $z$ (Flynn, Onu & Sutherland 2003). In effect, $z_s$ is an unknown that may be estimated by choosing the value that minimises the error between predictions and observations. The predicted amplitude $\tilde {A}(\omega , k_h; z)$ for a given frequency $\omega$ and horizontal wavenumber $k_h$ at height $z$ is then

(4.1) \begin{equation} \!\frac {\tilde {A}(\omega , k_h; z)}{A(\omega , k_h; z_0)} = \frac {\overline {N}(z_0)}{\overline {N}(z)}\sqrt {\frac {z_0-z_s}{z-z_s}} \exp \!\left [\! -\frac {k_h^3}{\omega ^4}\int _{z_0}^z \overline {N}^4(z') \overline {\nu }_{{tot}}^{{plume}}(z') \left (\overline {N}^2(z') - \omega ^2\right )^{-1/2} \mathrm {d}z'\!\right ]_{_{}}\!\!, \end{equation}

where the total viscosity is the sum of the molecular and SGS viscosity, $\nu _{ {tot}} = \nu + \nu _{ {SGS}}$ , and $\overline {\nu }_{{tot}}^{{plume}}$ is the time and horizontal average of $\nu _{ {tot}}$ within the plume.

Figure 7. (a) Vertical profiles of time and horizontal average stratification strength $\overline {N}$ , time and horizontal average within the plume $\overline {N}^{{plume}}$ , time and azimuthal average $\langle N \rangle$ , background stratification $N_0$ , and two profiles of the time-average $\overline {N}^t(x, z)$ on the plume centreline $x=0.5$ and at $x=0.4$ . (b) Vertical profiles of time and horizontal average total viscosity $\overline {\nu }_{{tot}}$ , time and horizontal average within the plume $\overline {\nu }_{{tot}}^{{plume}}$ , and molecular viscosity $\nu$ .

Figure 7 compares the full horizontal average with the plume average of the stratification strength $N$ and total viscosity $\nu _{ {tot}}$ . Figure 7(a) also shows two profiles of the spatially varying but time-averaged stratification strength $\overline {N}^t(x, z)$ , one on the plume centreline and another at the edge of the plume cap. Points close to the plume centreline and near the bottom of the stratified layer where $(\overline {N}^{t})^2 \lt 0$ are set to zero and the profile is smoothed with a running mean. The amplitude prediction (4.1) uses the full horizontal average $\overline {N}(z)$ since this more faithfully represents the stratification in the region away from the centreline where the waves propagate outward – compare the dotted and dashed green lines in figure 7. At heights within the plume, $\nu _{\text {tot}}$ is between two and three orders of magnitude larger than the molecular viscosity $\nu$ , owing to strong turbulence, and limits to $\nu$ in the ambient. For the total viscosity, the plume average appears to better reflect the turbulent structure noted in Powell et al. (2024), with stronger/weaker turbulence in the plume cap/intrusion, respectively. We therefore use the plume average $\overline {\nu }_{{tot}}^{{plume}}$ in (4.1). The results are qualitatively the same with the optimal choice of $z_0$ and $z_s$ (see below) when using $\overline {\nu }_{{tot}}$ instead, but the plume average improves the prediction when $z_0$ lies within the plume.

Figure 8. Application of the viscous internal wave model from TS07 to simulation HR. Line colours as in figure 5. (a) Observed spectrum $P(\omega ; z_0)$ at initial height $z_0 = 0.36\, \rm m$ . (b) Comparison of predicted spectrum $\tilde {P}(\omega ; z)$ (dashed line) and observed spectrum $P(\omega ; z)$ (solid line) with virtual source height $z_s=0.32\, \rm m$ at $z = 0.4\, \rm m$ and $0.48\, \rm m$ . (c) The NMSE between the predicted and observed spectrum, averaged over $0.44\, \rm m$ $ \le z \le 0.48\, \rm m$ , as a function of initial height $z_0$ for a range of virtual source heights $z_s\le z_0$ (indicated in colour).

As an example of applying the TS07 model, figure 8(a) shows the observed spectrum $P(\omega ; z_0)$ at the initial height $z_0=0.36\, \rm m$ . Using this initial spectrum, the predicted spectrum $\tilde {P}(\omega ; z)$ at heights $z=0.4$ and $0.48\, \rm m$ given a virtual source at $z_s = 0.32\, \rm m$ are plotted as dashed lines in figure 8(b) and compared with the corresponding observed spectrum at the two heights, plotted as solid lines. We compare predictions with observations quantitatively by computing the mean squared error between them and normalising by the mean of the squared observed spectrum. This normalised mean squared error (NMSE) averaged over predictions at $0.44\, \rm m \le z\le 0.48\, \rm m$ every $\Delta z=0.01\, \rm m$ is plotted as a function of initial height, $z_0$ , in figure 8 for a range of virtual source heights $z_s\le z_0$ . For all values of $z_s$ , the error between predicted and observed waves well above the fountain top is minimised if the initial observation height is $z_0 \geq 0.36\, \rm m$ . This height is close to, but below, the fountain top. The lowest error is achieved with virtual source heights well within the fountain, around the height of the intrusion. Even with these optimal parameters, as in figure 8(b), the viscous decay model does not perform as well here as in TS07. In particular, whilst the model does capture the selection of high-frequency waves and the overall decay in power is well represented, the shape of the spectrum is poorly predicted, with the decay of intermediate frequencies underestimated. It is noted by TS07 that the predicted shape is sensitive to the shape of the initial spectrum. However, a key observation from (4.1) is that the maximum amplitude for a given height $z$ and horizontal wavenumber $k_h$ , assuming a fixed stratification and viscosity and a flat initial spectrum, occurs at $\omega /N_0 = \sqrt {4/5} \approx 0.9$ which is close to the peak in the observed spectra seen in figures 5(c) and 8(b). Overall, this analysis gives an indication that waves are generated within the body of the fountain, though their spectrum there is not clearly established.

4.3. Dynamic mode decomposition and ray tracing

Motivated by the implication that internal waves are generated within the fountain, we use the dynamic mode decomposition (DMD) method (Schmid 2010) to extract spatial structures associated with internal wave frequencies $0 \lt \omega /N_0 \lt 1$ . We then use ray tracing to identify coherent wave beams within these structures. Modal decomposition has previously been used in similar problems; for example, Nidhan, Schmidt & Sarkar (2022) use a flow decomposition technique to link wake-generated internal waves with coherent turbulent structures. Later, it was shown that properties of these waves are in close agreement with a linear theory similar to that used in § 4.2 (Gola, Nidhan & Sarkar 2024). To apply DMD we construct a ‘data matrix’ $\boldsymbol{X}$ from snapshots of $x$ - $z$ slices through the plume centreline, with four observables: the perturbation horizontal velocity $u'$ , buoyancy $b'$ , vertical velocity $w'$ and out-of-plane vorticity $\zeta _y = \partial _z u - \partial _x w$ . The capability to use multiple observables as input data lends itself to extracting wave modes which are spatiotemporally coherent across all observables and reduces sensitivity to noise. The snapshots are restricted to $0.3\, {\rm m}\le x \le 0.7\, {\rm m}$ and $0.24\, {\rm m}\le z \le 0.5\, {\rm m}$ to avoid any transient signal from the front of the spreading intrusion. Each column of the data matrix corresponds to a discrete time $t_k = k\Delta T$ in the range $t_1 \le t/T_b \le t_2$ with $t_1, t_2$ as given in table 1. A second matrix $\boldsymbol{X}'$ is constructed by advancing each column by one time step (and assuming periodicity, so that the last column becomes the first). The ‘exact DMD’ algorithm computes the eigendecomposition of the linear operator $\boldsymbol{A}$ which advances $\boldsymbol{X}$ to $\boldsymbol{X}'\approx \boldsymbol{A}\boldsymbol{X}$ . The decomposition yields eigenvalues whose imaginary part represents the temporal frequency $\omega _j$ of mode $j$ and whose real part represents the growth rate (which is close to zero here). The spatial structure of each observable $i$ associated with mode $j$ is given by the eigenvectors $\boldsymbol {\Phi }^{(i)}_j$ and amplitudes $\mathcal{A}_j$ , which come in complex conjugate pairs for real input data. The spatial structure of $w'$ (for example) in mode $j$ with conjugate mode $j^*$ is

(4.2) \begin{equation} w'_{ {DMD}}(\omega _j) = {Re} \left [ \mathcal{A}_j\Phi _j^{w'} + \mathcal{A}_{j^*}\Phi _{j^*}^{w'}\right ]. \end{equation}

Figure 9. Examples of DMD modes. Here, $w'_{ {DMD}}$ as defined in (4.2) is plotted. (a) Evanescent mode, (b) turbulent mode and (c–e) internal wave modes. In (c–e), green dotted lines indicate the wave beam angle $\theta = \arccos (\omega _j/N_0)$ derived from the mode frequency $\omega _j$ .

The spatial structure of the DMD modes determined from simulation HR is summarised in figure 9. The modes can be broadly categorised into three types: internal wave modes with $0 \lt \omega _j/N_0 \lt 1$ , evanescent wave modes with $\omega _j/N_0 \gtrsim 1$ and turbulent modes with $\omega _j/N_0 \gg 1$ . For the internal wave modes shown, we use the polarisation relation $\omega = N_0 \cos \theta$ to plot dotted lines with the wave beam angle expected from the mode frequency. The close agreement with the phase lines demonstrates that these DMD modes successfully capture internal waves matching the mode frequency.

We now use ray tracing to examine whether the time-periodic structures determined by the DMD analysis can be interpreted as internal wave beams originating from within the fountain. Linear ray theory implicitly assumes slowly varying background fields, meaning turbulent fluctuations are neglected here. We also neglect the mean flow when propagating rays, which is justified since the group velocity of waves with wavenumber $k_{h,c}$ and frequency $\omega _c$ is an order of magnitude larger than the mean velocities in the plume. To avoid interference from left- and right-moving waves, we first apply a Hilbert transform (Mercier, Garnier & Dauxois 2008) and filter each DMD mode into left-moving internal waves with $k_h \lt 0, k_z \lt 0$ and right-moving internal waves with $k_h \gt 0, k_z \lt 0$ . The filtering introduces artefacts close to the horizontal centreline of the mode, which are reduced by applying a low-pass filter $\left | k_h \right | \lt 250$ , but still present. From ray theory, the path of waves in the $x$ - $z$ plane is given by $\mathrm {d}x/\mathrm {d}z = \tan \theta$ in which $\theta$ is the angle formed between lines of constant phase and the vertical (Sutherland 2010). This angle is determined implicitly by the polarisation relation $\omega _j = \overline {N}^t\!(x, z) \cos \theta$ , giving $\theta$ as a function of $x$ and $z$ for fixed $\omega _j$ . Representative profiles of $\overline {N}^t\!(x,z)$ are shown in figure 7. For a particular DMD mode, we integrate starting from a height $z_0 = 0.28\, \rm m$ and at a range of horizontal starting positions $0.4\, {\rm m}\le x_0 \le 0.6\, {\rm m}$ . Along each ray, we calculate the phase $\varphi$ associated with the buoyancy and vertical velocity fields of the DMD modes. Under the plane wave assumption $b'_{ {DMD}} = \hat {b}e^{i\varphi }, \, w'_{ {DMD}} = \hat {w}e^{i\varphi }$ , with the amplitudes being related by the polarisation relation $\hat {b}\omega \cos \varphi = \hat {w}N_0^2\sin \varphi$ (Sutherland 2010). Since the phase is constant along an internal wave beam, we identify coherent sections of each ray where the phase is within $\pi /4$ of the mean phase along the ray, $\overline {\varphi }$ : $\lvert \varphi - \overline {\varphi } \rvert \le \pi /4$ . This method implicitly neglects interference from out-of-plane wave beams. Motivated by the coherent wave beam structure observed in horizontal cross-sections in figure 3, we assume that this interference is weak.

Figure 10 shows an example of the ray tracing results for modes with $\omega _j/N_0 \approx 0.5, 0.8$ . Rays are shown as thick solid black lines where coherent and dashed otherwise. The mean-subtracted phase along each ray is shown by coloured lines in figures 10(c) and 10(f). It is expected that there is some noise in the phase, especially where $\overline {N}^t(x,z)$ is noisy near the fountain boundary. The artefacts introduced by filtering the waves, as well as the imperfect nature of the modal decomposition, also contributes noise. Nonetheless there are several rays in each mode along which the phase is approximately constant. This demonstrates that internal waves apparent in the region above the fountain can be traced to a source within the fountain.

Figure 10. Ray tracing in DMD modes 6 and 7 with $\omega /N_0 = 0.6, 0.7$ from an initial height $z_0=0.28\, \rm m$ and horizontal starting positions $0.4\, {\rm m}\le x_0 \le 0.6\, {\rm m}$ shown by coloured dots. Panels (a,d) show filtered $w'_{ {DMD}}$ and (b,e) show filtered $b'_{ {DMD}}$ as described in the text. Colour bar is the same as in figure 9. Panels (c,f) show phase perturbation $\varphi - \overline {\varphi }$ as a function of along-beam distance $r-r(z_p)$ from the plume edge at height $z_p$ . Lines are coloured according to starting position and highlighted where 50 % of the ray within the plume is coherent. Rays are solid black (thin dashed) in (a,b,d,e) where coherent (incoherent), being coherent when the phase perturbation in (c,f) lies within the solid black lines of those plots.

5. Discussion and conclusions

Inspired by the laboratory experiments of AS10, we have performed LES of a buoyant plume penetrating into a stably stratified layer that then transforms into a fountain and excites internal waves that propagate horizontally and vertically away from the fountain top. In all simulations, the plume source conditions were identical with the plume encountering the base of the stratified layer $0.2\, \rm m$ above the bottom of the domain. Across a range of different simulations, the strength of the stratification, $N_0^2$ , varied over two orders of magnitude from $0.25$ to $100\,\mathrm {s}^{-2}$ . Although the depth of penetration of the plume into the stratified layer decreases significantly with increasing $N_0^2$ , the frequency spectrum above the fountain was found to be narrow-banded with energy concentrated in the range $0.6 \lt \omega /N_0 \lt 1$ and a peak frequency at a near constant fraction of $N_0$ at approximately $\omega /N_0 = 0.9$ . We applied a viscous internal wave model, introduced by Taylor & Sarkar (2007), to help understand the origin of the internal wave spectrum, concluding that the model can sufficiently capture the power decay and frequency selection of the far-field internal wave spectrum when initialised from a spectrum taken near the top of the fountain $z_0=0.36\, \rm m$ and assuming a virtual source height $z_s = 0.32\, \rm m$ within the fountain, near the intrusion height. Using DMD, we were able to extract flow structures associated with internal wave frequencies and used ray tracing to demonstrate that internal wave beams can be traced from inside the fountain. This method is subject to significant noise in the phase reconstruction owing to filtering of the wave modes; superior methods of identifying a wave signature in the turbulent flow within the plume may exist. Whilst our analyses do not elucidate the wave generation mechanism, the results imply waves are generated within the fountain and not at the turbulent/non-turbulent interface between the fountain top and ambient fluid.

Other numerical studies of convection interacting with a stratified layer have likewise shown the excitation of internal waves with their source originating within the convective region (Lecoanet & Quataert 2013; Lecoanet et al. 2015; Pinçon et al. 2016; Couston et al. 2018). However, in those studies convective cells scoured the underside of a strongly stratified layer having buoyancy frequency far exceeding the characteristic convective frequency. The convection excited horizontally long, hydrostatic waves with a relatively wide but low-frequency spectrum. The separation of spatial and temporal scales between the convective and stratified regions allowed them reasonably to adapt Lighthill’s theory for the generation of sound waves by turbulence (Lighthill 1952) to the problem of internal waves generated by convection. Thus they showed that internal waves were excited by Reynolds stresses within the convecting region with a spectrum that decayed rapidly with increasing frequency. By contrast, our study focuses on a different regime with internal waves generated by spatially localised penetrative convection, in which the frequencies associated with turbulence in the fountain top overlap with the observed frequencies of generated non-hydrostatic internal waves. The lack of spatial and temporal scale separation in our problem means that Lighthill’s theory cannot be applied, nor could it explain the observed narrow frequency band at which waves are excited. Indeed, comparison of energy spectra in figure 6 with those shown in Couston et al. (2018) clearly shows different structures in spectral space. However, our simulation with $N_0^2 = 100$ appears to mark a transition between the penetrative convection regime and the scouring regime: the energy scale separation between the turbulent plume and waves is much smaller (figure 4 b) and the vertical energy flux in figure 4(a) more closely matches the theoretical scaling derived by Couston et al. (2018) compared with weaker $N_0$ . This trend continues in simulations with $N_0^2=1000$ , not shown here. Whatever the strength of the stratification, internal waves are generated within the turbulent region, though we argue the excitation mechanism differs for penetrative convection.