# The long view of triadic resonance instability in finite-width internal
gravity wave beams

K. M. Grayson1<EMAIL_ADDRESS>Stuart B. Dalziel1  Andrew G. W.
Lawrie2  1Department of Applied Mathematics and Theoretical Physics,
University of Cambridge, Wilberforce Road, Cambridge, CB3 0WA, UK 2Department
of Engineering, University of Bristol, Academic Street, Camford CF3 5QL, UK

###### Abstract

This paper presents our investigation into the modification of a finite-width
internal gravity wave beam arising from triadic resonance instability. We
present both experimental and weakly non-linear modelling to examine this
instability mechanism, in which a primary wave beam generates two secondary
wave beams of lower frequencies and shorter length scales. Through a versatile
experimental set-up, we examine how this instability develops over hundreds of
buoyancy periods. Unlike predictions from previous zero-dimensional weakly
non-linear theory, we find that the approach to a saturated equilibrium state
for the triadic interactions is not monotonic; rather, the amplitudes and
structures of the constituent beams continue to modulate without ever reaching
a steady equilibrium. To understand this behaviour we develop a weakly non-
linear approach to account for the spatio-temporal evolution of the amplitudes
and structures of the beams over slow time-scales and long distances, and
explore the consequences using a numerical scheme. Through this approach, we
establish that the evolution of the instability is remarkably sensitive to the
spatio-temporal triadic configuration for the system and how part of the
observed modulations can be attributed to a competition between the linear
growth rate of the secondary wave beams and the finite residence time of the
triadic perturbations within the underlying primary beam.

###### keywords:

Authors should not enter keywords on the manuscript, as these must be chosen
by the author during the online submission process and will then be added
during the typesetting process (see Keyword PDF for the full list). Other
classifications will be added at the same time.

##  1 Introduction

The meridional overturning circulation is critical in the regulation of the
earth’s climate, and understanding the processes essential for maintaining
this circulation is of key importance in global climate models. Munk (1966)
was amongst the first to suggest that internal gravity waves play a
significant part in the deep water vertical mixing of the density
stratification within the open ocean, and hence the maintenance of these
currents. It is now well established that the breaking of internal waves
contributes to the turbulent mixing in the ocean (Staquet & Sommeria, 2002;
Wunsch & Ferrari, 2004), yet only recently have the pathways by which internal
waves transfer energy to smaller scales and the eventual breaking events been
examined in more detail. As noted by Dauxois et al. (2018), our understanding
of these dissipative mechanisms, as opposed to internal wave generation,
leaves several open questions.

Various key mechanisms have been cited for how large-scale internal waves
cascade energy to smaller scales. These include internal wave reflection off
sloping boundaries (Nash et al., 2004), critical angle reflection (Dauxois &
Young, 1999) and scattering due to small scale topography (Peacock et al.,
2009). A review by Sarkar & Scotti (2017) suggests that no single mechanism is
responsible for the internal wave contribution to the energy cascade, rather
it is a combination of multiple linear and eventually non-linear processes.
MacKinnon & Winters (2005) and Alford et al. (2007) suggested that,
equatorward of a critical latitude (Richet et al., 2018), parametric sub-
harmonic instability (PSI), plays a dominant role in the energy transformation
of the internal tide into higher-mode near-inertial waves. Indeed, Sutherland
(2013) argues that away from sea-floor boundaries, and neglecting the
distorting influence of ocean currents, PSI is one of the primary mechanisms
for the energy cascade in the abyssal ocean. PSI can be viewed a special case
of triadic resonance instability (TRI), which is a weakly non-linear, slowly-
growing resonant mechanism whereby a primary wave becomes unstable due to
infinitesimal perturbations within the flow. As the instability grows, a
resonant triad interaction forms whereby the primary wave transfers energy to
two secondary waves of lower frequency and shorter length scale (Staquet &
Sommeria, 2002).

In the inviscid limit and under the assumption of an infinite plane-wave, the
frequencies of the secondary waves in the triad are equal to half of the
primary wave, motivating the traditional terminology of PSI (Fan & Akylas,
2019). While one often makes the appropriate assumption of oceanic scales
being inviscid, in the laboratory setting (where scales are smaller), viscous
effects cannot be neglected and resonant wave frequencies deviate away from
this sub-harmonic relationship. Moreover, for certain beam widths, the finite-
amplitude manifestation of this instability is unable to access these sub-
harmonic frequencies (Bourget et al., 2014). In the context of a viscous
finite-width beam it is therefore more appropriate to refer to TRI as opposed
to PSI.

The first reported experimental evidence of TRI for internal and interfacial
waves was approximately 50 years ago by Davis & Acrivos (1967), McEwan (1971)
and McEwan & Plumb (1977), who showed that for finite-width beams there exists
an amplitude threshold that must be surpassed for instability to occur. This
threshold is not found in the limiting case of an infinite plane-wave, where
infinitesimal perturbations may induce the development of the instability
(Koudella & Staquet, 2006). In fact, in the special case of a linearly
stratified Boussinesq fluid, a plane-waveform holds the peculiar property of
being an exact solution to the full non-linear equations at any amplitude
(e.g. Thorpe, 1968; Thorpe & Haines, 1986; Sutherland, 2006), albeit not a
linearly stable one. However, while single monochromatic plane-waves are
convenient mathematically, in nature waves will never take this form.
Realistically, oceanic waves are generated from baroclinic tides across ocean
ridges and will manifest as beams confined locally in space and therefore
broadly distributed over the wavenumber spectrum (Lamb, 2004; Gostiaux et al.,
2007). The focus of analyses using plane-wave solutions has been highlighted
in the review by Dauxois et al. (2018), who argue (correctly in our view) that
the effects of finite-width and envelope shape play an important, but
generally overlooked role, when considering the non-linearities of internal
waves.

In attempting to address these concerns, researchers have turned towards
exploring the dynamics of TRI in spatially localised internal wave beams.
Building on the work of Bourget et al. (2013), Bourget et al. (2014) calculate
a growth rate for the instability based on a energy balance that accounts for
the role of a finite-width beam. Using direct numerical simulations they also
show that the amplitude threshold for instability decreases as the beam width
is increased. This decrease is due to any perturbations having a larger
spatial field (and hence a longer time) in which to interact with the
underlying primary beam, thereby increasing the spatial region (and time
interval) over which energy can be transferred. These findings align with the
theoretical work of Karimi & Akylas (2014), who show how the form of the
carrier envelope for a finite-width wave beam has a significant influence on
its ability to become unstable based on the wavenumber spectrum produced from
the windowing. These works highlight the duality of interpretation for finite-
width beams in terms of both the physical parameters and the spectrum in
Fourier space.

Triadic resonance can arise due to the sustained spatio-temporal interactions
that occur when

${\phi_{0}}={\phi_{1}}+{\phi_{2}},$ (1)

where the wave phase, $\phi_{p}$, is defined as

$\phi_{p}=\boldsymbol{k}_{p}\cdot\boldsymbol{x}-\omega_{p}t.$ (2)

The subscript $p=(0,1,2)$ is used throughout this paper to define the primary
wave and the two secondary waves, respectively. Both (1) and (2) are true for
three-dimensions, but, without loss of generality we can rotate to a two-
dimensional (2D) co-ordinate system. The 2D wave vector of wave $p$ is defined
as $\boldsymbol{k}_{p}=(l_{p},m_{p})$ with magnitude
$|\boldsymbol{k}_{p}|=\kappa_{p}$, where the components are given in Cartesian
co-ordinates $(x,z)$, marked in Figure 1, and $\omega_{p}$ denotes the
frequency of the wave. In order to distinguish between the three wave beams in
the triad and their corresponding parameters, we define

$\mathbbm{B}_{p}=\\{\rho_{p},\Psi_{p};\omega_{p},\boldsymbol{k}_{p},\Lambda_{p}...\\},$
(3)

where $\mathbbm{B}_{p}$ indicates a wave beam with density $\rho_{p}$ and
stream function $\Psi_{p}$ fields, frequency $\omega_{p}$, characteristic
wavenumber vector $\boldsymbol{k}_{p}$ and beam width $\Lambda_{p}$. For the
primary beam $\omega_{0}$, $\boldsymbol{k}_{0}$, and $\Lambda_{0}$ are imposed
control parameters, whereas for the secondary beams they arise from the
triadic conditions and (weakly) non-linear dynamics. All triadic wave beams
must also satisfy the dispersion relationship for internal waves given as

$\frac{\omega_{p}}{N}=\pm\cos{\theta_{p}}=\pm\frac{\left|l_{p}\right|}{\sqrt{{l_{p}}^{2}+{m_{p}}^{2}}},$
(4)

where $\theta_{p}$ is the angle between the lines of constant phase and the
vertical and $l_{p}$ and $m_{p}$ are the characteristic wavenumber
contributions from each beam. Here $N$ is the buoyancy frequency of the
stratification given by

$N=\sqrt{-\frac{g}{\varrho_{0}}\frac{\partial\bar{\rho}}{\partial z}},$ (5)

where $g$ is the gravitational constant of acceleration. Under the assumptions
of a Boussinesq, incompressible fluid, we decompose the total density
$\varrho$ as $\varrho=\varrho_{0}+\bar{\rho}(z)+\rho(x,z,t)$, where
$\varrho_{0}$ is the reference density, $\bar{\rho}$ is the background density
stratification as a function of depth and $\rho$ is the perturbation density.
We consider the density changes from perturbations and stratification to be
small compared to the reference density, so that $\bar{\rho}$,
$\rho\ll\varrho_{0}$.

Given the triadic resonant condition in (1), it is easy to assume that the
instability selects one particular triad, comprised of three distinct
frequencies and wavenumbers for all time. More recently, our understanding of
triad selection is evolving for finite-width beams. Indeed, while examining
the transient start up of the instability, Koudella & Staquet (2006) note that
not just one triad is responsible for the initial instability, rather, a
number of triads form around the maximum linear growth rate. In addition,
recent work by Fan & Akylas (2020) shows how classic TRI theory is unable to
explain the instability in the context of a thin beam due to the broadband
wavenumber spectrum corresponding to the primary beam.

The novelty of the present paper lies in the examination of the long-term
evolution of the instability. Due to the 11 m long tank used in the
experimental set-up, we are able to observe the experiment for hours without
interference from side wall reflections or significant changes to the
stratification. We show experimentally that, over long time-scales, the
constituent triadic waves synchronously modulate in amplitude and in the
physical location of the two secondary wave beams. Further investigation shows
that part of these modulations are coincident with the growth and decay of
separate triads, all linked through the primary wave beam. Through two-
dimensional weakly non-linear modelling, we are then able to show how the
evolution of the instability in a finite-width beam is remarkably sensitive to
these separate triads. This sensitivity is due to their affect on the
residence time of the secondary wave beams with the underlying primary beam.

The outline of the remainder of this paper is as follows. In $\S$ 2 we detail
the experimental set-up and processing procedure. In $\S$ 3 we then present
the experimental results, looking first at the initial observations in $\S$
3.1 and then at the long-term evolution of the experiments $\S$ 3.2. Based on
these observations, we present the development of the two-dimensional weakly
non-linear model in $\S$ 4. Here we outline the perturbation expansion used in
$\S$ 4.1 and the subsequent development of the numerical solution in $\S$ 4.2
(comprising the two-dimensional advection scheme and the weakly non-linear
interactions). In $\S$ 5 we present the results of the model. We start with
$\S$ 5.1, where we examine the weakly non-linear interactions on their own
before moving onto $\S$ 5.2, where the results of the weakly non-linear two-
dimensional model are given. Conclusions are then drawn in $\S$ 6.

##  2 Experimental Procedure

###  2.1 Experimental setup

Experiments were undertaken in an 11 m long, 0.48 m deep, 0.29 m wide Perspex
(acrylic) tank. Along a 1 m section of the tank floor, 2.5 m away from one
end, sits the Arbitrary Spectrum Wave Maker (ASWaM), also known as the magic
carpet. This flexible horizontal boundary can generate sinusoidal forcing
(Dobra et al., 2022, 2021; Beckebanze et al., 2021) (as well as aperiodic
configurations (Dobra et al., 2019)), with the ability to vary amplitude,
frequency and wavenumber in both the temporal and spatial domain. The
wavemaker is comprised of a series of 96 computer-controlled linear actuators
that sit below the tank. Each actuator is mounted to a vertical drive rod that
passes through the base of the tank and connects to a 0.28 m long horizontal
rod that spans the tank width. These rods are spaced at 10 mm intervals along
the wavemaker. A 3 mm thick neoprene foam sheet covers the full length and
width of the wavemaker thus interpolating between the horizontal rods to allow
smooth forcing. The lengthwise edges of the neoprene slide against the tank
walls and beneath there is a 80 mm cavity into which glycerol is added to help
prevent salt crystallising and causing leakage around the seals that enable
the drive rods to pass into the tank from the bank of actuators beneath.
Provided the chosen waveforms preserve a zero-mean displacement across the
length of the flexible surface, the pressure gradient available to drive flow
around the edges of the neoprene foam is negligible. Thus flow in either
direction between the cavity and the visualisation region may be considered
negligible. When submerged in a stratified fluid, the wavemaker can generate
quasi-two-dimensional, internal wave beams at amplitudes sufficient to permit
wave breaking at distances away from the source. For full details of ASWaM’s
construction see Dobra (2018) and Dobra et al. (2019).

The procedure for filling the tank is as follows. First, glycerol is gravity
fed into the wavemaker cavity. The tank is then filled over the course of 8
hours with a linear salt-stratification using two computer controlled gear
pumps, operated via the software DigiFlow (Dalziel et al., 2007). Each pump
draws from either a fully saturated salt water or fresh-water reservoir. This
filling method allows for more precise control of the density stratification
compared with the traditional double bucket technique (Oster & Yamamoto,
1963), enabling the density gradient and fluid depth to be pre-determined. The
pumps used are Coleparmer Ismatec BVP-Z Analog gear pump drives mounted with
two magnetically driven Coleparmer Micropump L20562 A-Mount Suction Shoe pump
heads. The depth of the stratification $H$ = 0.45 $\pm$ 0.01 m.

To measure the density profile created by the pumps, an aspirating
conductivity probe is mounted to a linear traverse above the tank. One minute
before the start and one minute after the end of an experiment, the probe is
traversed downwards through the stratification to measure the conductivity of
the saline solution passing through the probe tip. For the experimental
campaign reported here, a linear density stratification with a buoyancy
frequency $N$ = 1.54 $\pm$ 0.04 rad s-1 is used. The variation in the buoyancy
frequency is attributed to the evolution of the density stratification over
the course of the week that the experiments were undertaken for. A schematic
of the tank, as viewed from the front, can be seen in Figure 1(a).

Figure 1: (a) A schematic showing the front view of the tank as would be seen
by the camera. The wavemaker is located along a 1 m section of the tank floor,
sitting below a 0.45 m density stratification. A conductivity probe is mounted
above the tank which is used to measure the density profile. (b) A schematic
showing the side view of the tank in order to visualise the optical
arrangement for Synthetic Schlieren. The thermal tunnel is not shown for
clarity.

###  2.2 Wave visualisation

Synthetic Schlieren (Dalziel et al., 2000, 2007) was used to visualise our
experiments. This non-intrusive technique takes advantage of the differential
refraction of light in a refractive index gradient and the Gladstone-Dale
relationship between refractive index and fluid density, such that light rays
curve towards regions of higher density. Internal waves make local
perturbations to the density field and thus the direction of light rays
passing through them will also be perturbed. The resulting distortion of a
textured background image yields a measurable signal associated with the
density perturbations. To minimise the effects of convective thermal
fluctuations on the Synthetic Schlieren measurements in the air between the
tank and the camera, a ‘thermal tunnel’ ran from the camera lens to the
perimeter of the visualisation region on the tank.

A random dot pattern attached to an LED light bank was located 0.20 $\pm$ 0.04
m behind the tank, while a 12 MPixel ISVI IC-X12CXP camera with a Nikkor
35-135 mm zoom lens was located 3.50 $\pm$ 0.10 m from the front. This optical
arrangement is shown in the side-view schematic in Figure 1(b). The large
distance between the camera and the tank was chosen in order to reduce
parallax (Thomas et al., 2009).

We compute the line-of-sight mean of the gradient vector of the density
perturbation field $\rho$, which for convenience we non-dimensionalise
according to

$\boldsymbol{\beta}=(\beta_{x},\beta_{z})=\frac{\lambda_{x_{0}}}{\varrho_{0}}\bigg{(}\frac{\partial{\rho}}{\partial
x},\frac{\partial{\rho}}{\partial z}\bigg{)},$ (6)

where $\lambda_{x_{0}}$ is the horizontal wavelength of the primary wave beam
given as $\lambda_{x_{0}}$ = 2$\pi/|l_{0}|$, where $l_{0}$ is the horizontal
component of the primary wave vector $\boldsymbol{k}_{0}$. Our experiment is
configured to generate and diagnose quasi-two-dimensional internal wave
structures, up to the limit of wave breaking, the point at which the mapping
of ray paths to density perturbations is no longer an aim.

###  2.3 Internal wave forcing

The experimental campaign presented in this paper comprises of 36 experiments.
To reduce uncertainties associated with the test conditions both within the
tank and in the laboratory ambient, the campaign was run within a seven-day
period without refilling the tank, allowing a period of 3 hours between each
experiment for any residual motion to dissipate. We note that other
experimental campaigns were also undertaken over the course of a year that
exhibited the same behaviour detailed below; for simplicity, we focus here on
this one campaign.

Following arguments laid out by Dobra et al. (2019), for all the experiments
in this campaign the vertical displacement $z=h(x,t)$ imposed on the neoprene
foam to generate the primary beam, $\mathbbm{B}_{0}$, is

$z=h(x,t)=\begin{cases}\Re\big{(}f(t)\hskip
1.70717pte^{il_{0}x}\cos^{2}\big{(}\frac{x-B}{8\pi^{2}}\big{)}\big{)},&\qquad\qquad
A<x<B,\\\ \Re\big{(}f(t)\hskip 1.70717pte^{il_{0}x}\big{)},&\qquad\qquad
B<x<C,\\\ \Re\big{(}f(t)\hskip
1.70717pte^{il_{0}x}\cos^{2}\big{(}\frac{x-C}{8\pi^{2}}\big{)}\big{)},&\qquad\qquad
C<x<D,\\\ 0,&\qquad\qquad\textmd{elsewhere},\\\ \end{cases}$ (7)

where the locations $A,B,C,D$ are respectively $7\pi/|l_{0}|$, $9\pi/|l_{0}|$,
$13\pi/|l_{0}|$, $15\pi/|l_{0}|$ and $l_{0}$ is set to $-$0.05 mm-1, giving a
horizontal wavelength of $\lambda_{x_{0}}$ = 2$\pi/|l_{0}|$ $=$ 125.66 mm. As
we restrict $\omega_{p}>0$ (for all $p$), having $l_{0}<0$ means the primary
wave beam is propagating to the left. The spatial structure of the forcing,
described by (7), takes the form of a beam with the inner two wavelengths
reaching maximum amplitude and the outer wavelengths being smoothed by a
cosine-squared envelope. Due to this cosine squared smoothing on the edges of
the beam profile, we do not consider the full width, $D-A$, for energy
transfer. Rather, we estimate the contribution from one of the smoothed edges
using the integral measure employed by Dalziel et al. (1999), giving a
horizontal beam width of

$\Lambda_{x_{0}}=\Lambda_{0}/\cos\theta=2\lambda_{x_{0}}+2\int_{0}^{\lambda_{x_{0}}}\alpha\big{(}1-\alpha\big{)}dx=277.41\textrm{
mm},$ (8)

where $\alpha=\cos^{2}(x/8\pi^{2})$ is the smoothing function on the outer
flanks of the beam profile. The temporal forcing $f(t)$ is then described as

$f(t)=\begin{cases}0,&\qquad\qquad\hskip 11.38109ptt\leq 0\>\mathrm{s},\\\
\eta_{0}\hskip
1.70717pt\big{(}\frac{t}{30}\big{)}e^{-i\omega_{0}t},&\qquad\qquad\hskip
11.38109pt0\leq t\leq 30\>\mathrm{s},\\\ \eta_{0}\hskip
1.70717pte^{-i\omega_{0}t},&\qquad\qquad\hskip 11.38109pt30\>\mathrm{s}\leq
t\leq t_{\textrm{end}}-30\>\mathrm{s}\>,\\\ \eta_{0}\hskip
1.70717pt\big{(}\frac{t_{\textrm{end}}-t}{30}\big{)}e^{-i\omega_{0}t},&\qquad\qquad\hskip
11.38109ptt_{\textrm{end}}-30\>\mathrm{s}\leq t\leq
t_{\textrm{end}}\>,\end{cases}\quad$ (9)

where $\omega_{0}$, $\eta_{0}$ and $t_{\textrm{end}}$ are respectively, the
forcing frequency of 0.95 rad s-1, the nominal forcing amplitude in mm of the
primary beam and the end time of the experiment in seconds. Experiments are
captured at 1 frame per second (fps), which is more than sufficient to capture
the fast time evolution of the wave field given by the primary beam period
$T_{0}=2\pi/\omega_{0}$.

The only two parameters to be varied in this experimental study are
$t_{\textrm{end}}$ and $\eta_{0}$. The run time, $t_{\textrm{end}}$, is either
90 or 180 minutes, while the non-dimensional amplitude
$\eta_{0}/\lambda_{x_{0}}$ ranges between 0.028 - 0.036. The amplitude
threshold for the instability is achieved at $\eta_{0}/\lambda_{x_{0}}\approx$
0.031 (reducing by $0.002$ throughout the week due a slow evolution of the
stratification). Our focus is on the weakly non-linear regime, so we seek to
minimise unnecessary mixing induced by wave actuation and limit our amplitudes
to those just sufficient to exceed the instability threshold calculated by
(Davis & Acrivos, 1967).

Since the tank extends well beyond the field of view in both directions,
internal wave beams with typical dominant wavenumbers of
$\boldsymbol{k}_{0}=(-0.05,-0.06)$ mm-1 reflecting off the far boundary wall
return to the viewing region with only 2$\%$ of their original amplitude, due
to viscous dissipation over a beam length exceeding 4 m (the horizontal travel
distance of the beam before re-interaction). We thus consider wave-wave
interactions involving these reflected beams to be negligible.

##  3 Experimental results

###  3.1 Initial observation and analysis

We start by examining one experiment from the set of 36 with an imposed
amplitude displacement of $\eta_{0}/\lambda_{x_{0}}$ = 0.032. Figure 2(a)
shows $\beta_{z}$ over the visualisation region at $t/T_{0}$ = 83. Here,
$\mathbbm{B}_{0}$, generated by the wavemaker, propagates energy up and to the
left, at its respective group velocity $\boldsymbol{c}_{g_{0}}$. The group
velocity is defined for all wave beams by

$\boldsymbol{c}_{g_{p}}=\bigg{(}\frac{\partial}{\partial
l_{p}},\frac{\partial}{\partial
m_{p}}\bigg{)}\omega_{p}=\textrm{sgn}(l_{p})\frac{Nm_{p}}{\kappa_{p}^{3}}(m_{p},-l_{p}),$
(10)

where again the subscript $p=(0,1,2)$ corresponds to the primary beam and the
two secondary beams, respectively, and the broadband wavenumber spectrum of
each beam is approximated with a characteristic wavenumber. The primary beam
$\mathbbm{B}_{0}$ reflects off the free surface, causing the vertical
component of its group velocity to change sign and the wave-packets
subsequently move down and to the left. An appropriate Reynolds number for the
flow is given by $Re=\boldsymbol{c}_{g_{0}}/(\nu\kappa_{0})$, where $\nu$ is
the kinematic viscosity of 1 mm2 s-1. Here, $Re\approx 170$. As the selected
input amplitude displacement of $\eta_{0}/\lambda_{x_{0}}$ = 0.032 is above
the instability threshold, $\mathbbm{B}_{0}$ becomes unstable, leading to the
formation of two secondary beams. One of these beams, $\mathbbm{B}_{1}$, is
clearly visible in Figure 2(a). This beam emanates from the central region of
$\mathbbm{B}_{0}$ but moves in nearly the opposite direction, with a group
velocity down and to the right. From Figure 2(a), the third beam,
$\mathbbm{B}_{2}$, that completes the triad is not visible. In order to
understand the underlying modal structure of these beams, the flow field
$\boldsymbol{\beta}$, is decomposed using Dynamic Mode Decomposition (DMD).

DMD works by preforming an eigen-decomposition of a linearised representation
of the underlying evolution operator for a given flow field (Schmid, 2010).
The ‘dynamic modes’ are the recurrent spatial structures that accurately
describe the dominant behaviour captured in the data sequence. Where DMD
excels is in determining the frequencies and structure of the modes from short
time series where there is a discrete spectrum that can reasonably be
approximated by a combination of delta functions at slowly evolving
frequencies. The ability to extract the modes from short time series allows
exploration of the slowly evolving structure and frequency of the modes. This
linear approximation for the evolution operator is valid for the experiments
shown here due to the two discrete time-scales, whereby the slow time
evolution of the beam amplitude is much less than the fast time-scales
$\omega_{p}t$.

The maximum number of dynamic modes is given by the number of input frames in
the sequence $\delta t$ (in this case $\delta t$ = 20 s, as we use a frame
rate of 1 fps, which is just greater than the slowest period of the triad
$T_{1}$) and if the obtained mode is complex then it is coupled as a conjugate
pair. Here, however, we are only interested in those modes with an eigenvalue
modulus very close to one, as they represent the steady, non-decaying modes of
the system. When instability occurs experimentally, four non-decaying modes
are obtained, three of which are conjugate pairs of eigen-values. The real
part of these three modes produced over the temporal window $83\leq
t/T_{0}\leq$ 86 are given in Figure 2(b)-(d). As expected from the input
forcing, Figure 2(b) corresponds to the input $\mathbbm{B}_{0}$ with non-
dimensional frequency $\omega_{0}/N$ = 0.62. Figure 2(c) then corresponds to
$\mathbbm{B}_{1}$ with $\omega_{1}/N=$ 0.23 and (d) to the obscured
$\mathbbm{B}_{2}$, which propagates with a group velocity up and to the left,
with non-dimensional frequency $\omega_{2}/N=$ 0.39. To understand if these
additional frequencies are the result of TRI, we sum the frequencies of the
secondary beams and see that the temporal condition for triadic resonance,
$\omega_{0}=\omega_{1}+\omega_{2}$, is satisfied. We remark that this
frequency relationship is not enforced at any stage of experimental post-
processing, but arises naturally from prominent signals found in the temporal
spectrum.

$\boldsymbol{k}_{0}$$\boldsymbol{k}_{1}$$\boldsymbol{k}_{2}$ Figure 2: (a)
$\beta_{z}$ of the full flow field at $t/T_{0}=$ 83 into an experiment forced
at $\eta_{0}/\lambda_{x_{0}}$ = 0.032. (b)–(d) The real part of three of the
dominant frequencies produced from the DMD over $83\leq t/T_{0}\leq$ 86\. The
black arrows overlaid indicate the orientation of the respective wavenumber
vectors $\boldsymbol{k}_{p}$. In panel (b) we see solely the wave field
$\mathbbm{B}_{0}$ with $\omega_{0}/N$ = 0.62. In (c) we see $\mathbbm{B}_{1}$
with $\omega_{1}/N=$ 0.23 and in (d) $\mathbbm{B}_{2}$ with $\omega_{2}/N=$
0.39. The black box in (b) shows the spatial averaging domain
$\langle\rangle_{r}$ used for the primary beam, discussed in $\S$ 4.2.

The fourth (non-decaying mode) corresponds to $\omega/N$ = 0 and is not shown
here. This is generated from a two-wave interaction (TWI), in which two wave
beams interact to produce a third wave beam, with a phase angle relationship

$\check{\phi}=\pm\phi_{0}\mp\phi^{\prime}_{0}.$ (11)

In this case, $\phi_{0}=l_{0}x+m_{0}z-\omega_{0}t$, corresponds to the phase
angle of $\mathbbm{B}_{0}$ and $\phi^{\prime}_{0}=l_{0}x-m_{0}z-\omega_{0}t$
to its reflection, $\mathbbm{B}^{\prime}_{0}$, from the free surface. These
wave beams will sum to produce a third wave beam with wavenumber vector
$\check{\boldsymbol{k}}=(0,2m_{0})$ aligned with the vertical and with zero
frequency. This non-propagating disturbance can not be classed as a wave, but
instead should be treated as a forced oscillatory structure that is confined
to the interaction region of the primary beam with its reflection. If
considered analytically (Thorpe & Haines, 1986) or numerically (Grisouarda et
al., 2013) in a two-dimensional setting, only weak horizontal vorticity is
generated, which is partially suppressed by the background stratification
(Beckebanze et al., 2019). When considered in a three-dimensional setting,
however, Grisouarda et al. (2013) show, both experimentally and numerically,
how a stronger slowly evolving three-dimensional horizontal mean flow develops
from the interaction region of the primary beam with its reflection. This flow
has a vertical component to its vorticity field. Indeed, if viscous
attenuation and cross-beam variations are present, it is possible for a three-
dimensional mean flow to be generated from the wave beam interacting with
itself, as shown analytically by Kataoka & Akylas (2015) and experimentally by
Bordes et al. (2012). In all of the three-dimensional cases cited above,
however, the wave beam is propagating in a tank wider than the beam width.
This allows for a recirculating mean flow to develop in the transverse
direction, outside of the spatial extent of the beam. As noted by Sutherland
(2006) in experiments where wave beams are confined laterally by tank side
walls, as is the case in the experiments presented here, horizontal mean flow
of this type is unable to develop. The observed zero-frequency mode in our
experiments, closely resembles the two-dimensional simulations of Grisouarda
et al. (2013) and, while the disturbance does slowly exit the interaction
region of $\mathbbm{B}_{0}$ and $\mathbbm{B}^{\prime}_{0}$, no strong
recirculating mean flow is seen to develop and as such does not impact the
evolution of TRI described here.

We proceed to determine the wave vectors corresponding to the primary and
secondary wave beams by taking our frequency-decomposed gradient field over
the temporal window $83\leq t/T_{0}\leq$ 86 – the real parts of which are
shown in Figures 2(b)–(d) – and calculating a two-dimensional power spectra on
each constituent field separately. Each image is embedded in a zero filled
matrix in order to improve resolution and limit spatial aliasing. The
wavenumber is determined by fitting a quadratic curve to the peak of the
resultant power spectra and finding the wavenumber corresponding to the peak
of the curve. This procedure is preformed on every row and column of the
domain and subsequently mean averaged over both spatial dimensions. The
smallest resolvable length scale is 2 pixels, equivalent to the non-
dimensional length $x/\lambda_{x_{0}}=0.005$, given by the ratio of pixel
resolution to region size. As the analysis preformed on the horizontal density
gradient $\beta_{x}$ provides similar results to that of the vertical
$\beta_{z}$, we use only the results from the vertical gradient for
simplicity. The non-dimensional characteristic wavenumbers for the vertical
gradient fields shown in Figure 2 are $\lambda_{x_{0}}\boldsymbol{k}_{0}$ =
($-$6.28, $-$8.29), $\lambda_{x_{0}}\boldsymbol{k}_{1}$ = (3.73, 14.60), and
$\lambda_{x_{0}}\boldsymbol{k}_{2}$ = ($-$9.93, $-$24.88). These wave vectors
are shown by the blue arrows on Figure 3, where the underlying black, red and
green curves provides the locus of all possible solutions for
$\boldsymbol{k}_{1}$ given $\boldsymbol{k}_{0}$, based on both the dispersion
relationship (4) and the TRI condition (1).

Figure 3: The underlying solid and dashed black, red and green curves give all
of the possible locations for the tip of $\boldsymbol{k}_{1}$ that satisfy
both the dispersion relationship (4) and the TRI condition (1) for a given
$\mathbbm{B}_{0}$. The dark blue arrows show the experimentally produced,
characteristic, wavenumber vectors of the resonant triad shown in Figure 2,
obtained from taking the Fourier transform in $(x,z)$ of the gradient field.
The shaded grey region then indicates the range of wavenumber vectors obtained
over the course of the experiment. The six dark blue marks correspond to the
different triad wave vector configurations used in the weakly non-linear
modelling and are discussed in $\S$ 5.1. The panel in the bottom right corner
shows an enlarged view of the region enclosed by the black rectangle.

While the calculated characteristic wave vectors shown in Figure 3 lie almost
in a closed triangle, their alignment is not perfect, potentially indicating
that the spatial triadic resonance condition
$\boldsymbol{k}_{0}=\boldsymbol{k}_{1}+\boldsymbol{k}_{2}$ is not exactly
satisfied. The reason for this slight misalignment is due to three factors.
Firstly, there is the impact of inevitable experimental noise. Secondly, as we
are considering finite-width beams as opposed to plane-waves, each beam is
comprised of a broadband wavenumber spectrum. By defining a single
characteristic wavenumber for the beam – taken from the peak of the Fourier
spectrum – we are therefore approximating this wavenumber distribution.
Thirdly, we are assuming that the spatial structures of $\mathbbm{B}_{1}$ and
$\mathbbm{B}_{2}$ are uniform over the field of view. In fact, as the
experiment progresses, significant modulations to the structures of
$\mathbbm{B}_{1}$ and $\mathbbm{B}_{2}$ are observed, revealing that this
assumption of spatial uniformity is inappropriate.

While it is clear that TRI was indeed being witnessed experimentally in a
finite-width beam, this in itself is not novel. In an experiment actuated by
an oscillating cylinder, Clark & Sutherland (2010) attribute the breakdown of
a wave beam due to TRI, showing how the instability evolves from infinitesimal
perturbations in the flow. This work has recently been developed by Fan &
Akylas (2020), who discuss the validity of TRI theory in thin wave beams.
Moreover Joubaud et al. (2012) and Bourget et al. (2013) clearly show the
growth of the instability for a finite-width beam in experiments using their
sidewall wavemaker. In our work, the regime of interest is not the initial
growth of the instability, but rather the finite-amplitude unsteady
modulations that occur afterwards. As noted, the expected saturated
equilibrium state for the weakly non-linear instability is not observed,
rather we witness slow modulations of the amplitudes and structures of the
constituent beams in the triad, revealing much more dynamical behaviour than
anticipated. We investigate the long-term evolution of this unsteady behaviour
for the remainder of the paper.

###  3.2 Long-time development

Figure 4 shows 8 instantaneous images of the experiment shown in Figure 2.
Figure 4(a) is captured at $t/T_{0}$ = 53 into the experiment, just as
$\mathbbm{B}_{0}$ becomes visibly unstable. By the instant shown in 4(b) (the
same image shown in Figure 2(a)) $\mathbbm{B}_{1}$ has clearly developed, with
a group velocity propagating down and to the right. A particularly interesting
feature of the subsequent time frames is the modulation of $\mathbbm{B}_{1}$
over time. Not only is its region of generation not constant – it migrates
across the full height of $\mathbbm{B}_{0}$ – the beam itself also varies in
both intensity and width. This migratory behaviour persists for the full
duration of the experiment, which lasts for over 800 periods of the primary
beam.

Figure 4: Sequence of images showing the vertical gradient of the density
perturbation of an experiment with $\eta_{0}/\lambda_{x_{0}}$ = 0.032. (a)
$t/T_{0}$ = 53, (b) $t/T_{0}$ = 83, (c) $t/T_{0}$ = 113, (d) $t/T_{0}$ = 144,
(e) $t/T_{0}$ = 174, (f) $t/T_{0}$ = 213, (g) $t/T_{0}$ = 432, (h) $t/T_{0}$ =
582. The black lines in (g) indicate where the $\mathbbm{B}_{1}$ beam changes
frequency, evidenced by the subtle change in angle in between the two lines.

Further quantitative analysis of this peculiar behaviour requires us to
calculate the amplitude of the individual resonant wave beams
$\mathbbm{B}_{0}$, $\mathbbm{B}_{1}$ and $\mathbbm{B}_{2}$. Decomposing by
frequency into complex constituent fields using DMD, we find the inverse
gradient (potential) field $\rho_{p}/\varrho_{0}$ by integrating both the real
and imaginary components of $\boldsymbol{\beta}$. We then isolate the wave
beam of interest further using a Hilbert Transform, first used for internal
waves by Mercier et al. (2008). This filtering technique is applied to isolate
the quadrant of Fourier space containing the wave vectors of the beam of
interest from other signals of the same temporal frequency e.g. separating
$\mathbbm{B}_{0}$ from its reflection from the free surface,
$\mathbbm{B}^{\prime}_{0}$.

In order to have a singular value for amplitude that is independent of space,
we then spatially mean average $\rho_{p}/\varrho_{0}$ over the whole field of
view denoted $\langle\rangle_{w}$. This choice of spatial averaging ensures
that our measure of amplitude is decorrelated with the position of a beam in
space, a topic that will warrant further discussion in $\S$ 5.1. An
unavoidable consequence of this choice, however, is that this average measure
no longer represents the local amplitude within a beam. To account for this,
the other region used for spatial averaging is shown by the black box in
Figure 2(b), which we donate as $\langle\rangle_{r}$. This region is only ever
used for the primary beam and is used to compare the experimental input
amplitude with the two-dimensional and zero-dimensional modelling discussed in
$\S$ 4.2.

We then infer the amplitude from the measured displacements of the scalar two-
dimensional stream function $\Psi=\Psi(\boldsymbol{x},t)$ for each field,
where the velocity vector
$\boldsymbol{u}=\boldsymbol{\nabla}\times(\Psi\hat{\boldsymbol{y}})$ and
$\boldsymbol{\nabla}=(\partial/\partial x,\partial/\partial z)$. Assuming an
oscillatory form $\Psi_{p}=\tilde{\Psi}_{p}e^{i\phi_{p}}$ for each triadic
wave beam, we define our spatially averaged amplitude for each constituent
field as $\langle\tilde{\Psi}_{p}\rangle_{w}$, a quantity that is independent
of the oscillatory time-scale. As our focus here is on the slow time-scale
evolution of the field, we allow
$\langle\tilde{\Psi}_{p}\rangle_{w}=\langle\tilde{\Psi}_{p}(t)\rangle_{w}$, on
a slow time-scale well-separated from the oscillation period of all wave beams
in the system. In the same way, we define the density perturbation
$\rho_{p}=\tilde{\rho}_{p}e^{i\phi_{p}}$ for each field.

Substituting the above forms for stream function and density into the inviscid
linear conservation of mass equation $\partial\rho/\partial
t=-w(d\bar{\rho}/dz)$, and cancelling the fast time-scales, we find the
reduced stream function amplitude for each wave field using

$\langle|\tilde{\Psi}_{p}|\rangle_{w}=\bigg{|}\frac{\omega_{p}}{l_{p}}\frac{g}{N^{2}}\bigg{|}\frac{\langle|\tilde{\rho}_{p}|\rangle_{w}}{\varrho_{0}},$
(12)

where $\omega_{p}$ and $l_{p}$ are the frequency and horizontal wavenumber of
the given beam $p$, and $\langle|\tilde{\rho}_{p}|\rangle_{w}$ is the root
mean square (magnitude) of the complex output of the wave field after being
spatially filtered by the Hilbert Transform.

Figure 5 shows the amplitude calculated from (12) for two experiments. In (a)
we show the same experiment as Figure 4, while (b) corresponds to another
experiment with the same amplitude ($\eta_{0}/\lambda_{x_{0}}$ = 0.032) but
with a much longer run time ($t_{\textrm{end}}/T_{0}=1633$). We first note
that the growth of the secondary wave beams appears earlier in (a) than (b)
and that the maximum amplitude of the primary wave beam in (a) is larger,
despite both experiments having the same amplitude displacement,
$\eta_{0}/\lambda_{x_{0}}$, from the wavemaker. This is due to the
deterioration of the stratification throughout the week of experiments, which
results in decreased transmission from the wavemaker to $\mathbbm{B}_{0}$ and
hence a slight reduction in the instability threshold. This change emphasises
the need to use the measured wave beam amplitude, calculated using (12), as
opposed to the imposed displacement from the wavemaker.

Figure 5: The non-dimensional amplitude of the reduced stream function
$\epsilon_{w}=\kappa_{0}^{2}\langle|\tilde{\Psi}_{p}|\rangle_{w}/N$, for
$\mathbbm{B}_{0}$ (blue), $\mathbbm{B}_{1}$ (red), $\mathbbm{B}_{2}$ (green)
for two experiments. Here $\langle|\tilde{\Psi}|\rangle_{w}$ is defined in
(12), where the spatial averaging of each signal is performed over the whole
domain $\langle\rangle_{w}$. The dimensional amplitude of each wave field is
given on the right axis. (a) $t_{\textrm{end}}/T_{0}=816$,
$\eta_{0}/\lambda_{x_{0}}$ = 0.032. (b) $t_{\textrm{end}}/T_{0}=1633$,
$\eta_{0}/\lambda_{x_{0}}$ = 0.032.

Another observable feature in Figure 5 is the gradual increase of the mean
amplitude of $\mathbbm{B}_{0}$ over time. This behaviour is also seen in lower
amplitude forcing experiments that did not become unstable to TRI (not shown
here). This increase can not be directly due to the instability, as the TRI
mechanism transfers energy from the primary beam to the two secondary wave
beams, as opposed to injecting energy into the primary wave beam. Rather, the
amplitude increase is believed to be due to the peristaltic motion of the
wavemaker leading to a sharpening of the stratification directly above the
wavemaker, resulting in an increased transmission efficiency between the
energy transfer from the wavemaker to the internal waves.

The most prominent feature in Figure 5 is the amplitude modulations of all the
triadic wave beams, observed in every experiment that became unstable. While
these modulations were anticipated from qualitatively observing the
experiments, quantitatively they are found to be unexpectedly large and
without obvious periodicity. This behaviour was so striking that we initially
sought explanations unrelated to the physics of the system, such as
measurement errors in converting raw video footage to density gradient fields
or discrepancies that might be introduced by frequency-decomposition into
constituent fields. After careful examination of both the raw data and the
tool chain, including replicating the harmonic analysis of Mercier et al.
(2008) – a technique that relies solely on Fourier transforms to isolate waves
before calculating $\tilde{\Psi}$ using (12) – we were able to discount all
extraneous sources that could contribute to these structural modulations.

In Figure 5, the amplitudes of $\mathbbm{B}_{1}$ (red) and $\mathbbm{B}_{2}$
(green) are positively correlated; their amplitudes are almost scaled values
of each other. Meanwhile, the amplitude of $\mathbbm{B}_{0}$ is negatively
correlated with $\mathbbm{B}_{1}$ and $\mathbbm{B}_{2}$. When
$\mathbbm{B}_{0}$ is at a local maximum, the amplitudes of $\mathbbm{B}_{1}$
and $\mathbbm{B}_{2}$ are concurrently at a local minimum and then versa when
the amplitude of $\mathbbm{B}_{0}$ is at a minimum. This coupling of the
modulations in amplitude between $\mathbbm{B}_{0}$ and the secondary
$\mathbbm{B}_{1}$ and $\mathbbm{B}_{2}$, reveals a continuous energy exchange
flux between the wave beams in the triad that does not saturate to a steady
equilibrium. For these experiments, the pattern of slow modulation appears to
be independent of the primary wave beam amplitude, as, when normalised, the
amplitude ratios $\mathbbm{B}_{1}/\mathbbm{B}_{0}$ and
$\mathbbm{B}_{2}/\mathbbm{B}_{0}$ are similar across all experiments that
become unstable independent of the amplitude of the forcing. Despite the clear
pattern of modulations shown in Figure 5, there is sufficient randomness that
the signal does not have a clear dominant frequency in Fourier space. This
observation is common to all experiments where instability develops.

Both the physical positioning of the secondary wave beams (seen in Figure 4)
and their amplitudes (shown in Figure 5) undergo slow modulation. Less obvious
is that the beam frequencies also simultaneously modulate. This is evidenced
in Figure 4(g), where the angle of $\mathbbm{B}_{1}$ is noticeably closer to
the horizontal in the lower part of the domain in comparison to the upper part
of the domain (in between the two black lines). To further understand the slow
evolution of beam frequencies, Figure 6 shows the temporal-frequency spectra
computed using a Fourier-transform for both experiments presented in Figure 5,
along with the corresponding DMD estimates of the triadic frequencies overlaid
in white. The amplitude of the spectra is determined by

$S_{\beta_{z}}(\omega,t)=\bigg{\langle}\bigg{|}\frac{1}{T_{T}}\int_{-\infty}^{+\infty}\beta_{z}(x,z,t^{\prime})e^{-i2\pi(\omega
t^{\prime})}W(t^{\prime}-t;T_{T})dt^{\prime}\bigg{|}^{2}\bigg{\rangle}_{w},$
(13)

where $W(t^{\prime};T_{T})$ is a Hamming window of non-dimensional width
$T_{T}/T_{0}$ = 39. For the frame rate of 1 fps, the highest resolvable
frequency (shortest time period) is $\omega/N=4.08$. Several windowing
functions were tested, and were not found to significantly affect the
spectrogram results. The angled brackets, $\langle\rangle_{w}$, again indicate
that the results are averaged across the whole visualisation region. This
underlying spectrogram, calculated using (13), therefore, reveals the details
about the distribution of the frequency spectra for $\omega_{1}$ and
$\omega_{2}$. In contrast, as we are only selecting the three dominant modes
obtained from the DMD over short time intervals ($\delta t/T_{0}=3$), this
methods approximates the underlying energy spectrum by a series of delta
functions, allowing us to clearly see the slow-time evolution of these
dominant modes.

$\frac{\omega_{2}}{N}$$\frac{\omega_{1}}{N}$$\frac{\omega_{0}}{N}$$\big{\\}}$$\big{\\}}$$\big{\\}}$$\big{\\{}$$\big{\\{}$$\big{\\{}$
Figure 6: Time-frequency spectra computed for the experiments in Figure 5. The
spectral density is computed by (13) and is normalised by the total energy
$S_{E}$ = $\sum_{z=0}^{H}S_{\beta_{z}}(\omega)^{2}$ for each instant in time.
The dominant frequencies for each experiment obtained from the DMD frequency-
decomposition are overlaid in white. The white dashed lines in (a) indicate
the times of the instantaneous images shown in Figure 8. The subplot overlaid
on (b) shows a transect in time at $t/T_{0}$ = 452, marked by the black and
magenta arrows. Here we have plotted $\ln(S_{\beta_{z}}(\omega)/S_{E})$ in
cyan and $\ln(S_{\beta_{z}}(\omega_{0}-\omega)/S_{E})$ in magenta against
$\omega/\omega_{0}$.

Both spectrograms in Figures 6(a) and (b) show a clear peak at $\omega_{0}/N$
= 0.62 for all time, consistent with the imposed displacement from ASWaM. Both
secondary beams emerging from the instability become visible at approximately
$t/T_{0}$ = 50, with peaks in the spectra around $\omega_{1}/N$ $\approx$ 0.23
and $\omega_{2}/N$ $\approx$ 0.39, though subsequently these modulate on a
slow time-scale throughout the duration of an experiment. The overlaid DMD
frequency estimates match almost perfectly the three frequency peaks on the
spectrogram, following the same pattern of slow modulations. Despite this
modulation, the temporal triadic relationship
$\omega_{0}=\omega_{1}+\omega_{2}$, is satisfied at all times for the
frequencies obtained from the DMD. As noted previously, the triadic
requirement is not built into the DMD analysis. Interestingly, a similar
variation in frequency has been witnessed by both Bourget et al. (2013) and
Brouzet et al. (2016) in their experimental studies, however the phenomenon
was not the focus of their work.

In addition to the triadic frequencies, there are three other distinct
frequency bands found in the time-frequency spectrograph. The band with the
lowest frequency corresponds to $\omega/N$ $\approx$ 0, which was also
observed from the DMD and has already been discussed. The other two
frequencies $\omega/N$ $\approx$ 0.84 and $\omega/N$ $\approx$ 1 correspond to
two different TWIs, given in (11), between $\mathbbm{B}_{0}$ and either
$\mathbbm{B}_{1}$ or $\mathbbm{B}_{2}$, respectively.

What is perhaps most striking from these time-frequency spectra is how, at
certain points in time, there are multiple sets of $\mathbbm{B}_{1}$ and
$\mathbbm{B}_{2}$ associated with the instability. This is shown by the
convergent ‘wisps’ on the $\omega_{1}/N$ and $\omega_{2}/N$ bands, where
additional secondary beam pairs appear and merge with the continuous mode.
This is highlighted for $t/T_{0}$ = 452 by the inset in Figure 6(b). Here we
have plotted $\ln(S_{\beta_{z}}(\omega)/S_{E})$ in cyan and
$\ln(S_{\beta_{z}}(\omega_{0}-\omega)/S_{E})$ in magenta against
$\omega/\omega_{0}$. The presence of a spectrum of triadic relations here is
evidenced by the strong correlation between the two traces, indicating that
the triadic requirement $\omega_{1}+\omega_{2}=\omega_{0}$ persists across all
the spectrum. To analyse these frequency modulations further, Figure 7 shows
the real part of the dynamic modes associated with the three dominant pairs of
frequencies from the DMD over frames 481 $\leq t/T_{0}\leq$ 483, from the
experiment presented in Figure 4. Unlike its earlier counterpart in Figures
2(c) and (d), where there was one distinct frequency and wavenumber pair for
both $\mathbbm{B}_{1}$ and $\mathbbm{B}_{2}$, Figures 7(c) and (d) show that,
at this instant in time, TRI is occurring at two different locations over the
height of the primary beam. For both of these modes, the signal is
discontinuous across a transition region where the two out-of-phase wave beams
de-constructively meet, highlighted by the black dashed rectangle in Figure
7(d). Indeed, the presence of these separate beams is confirmed by splitting
the domain in half horizontally and preforming the DMD analysis separately on
the two halves. For the upper half of the domain $\omega_{1}/N$ = 0.211 and
$\omega_{2}/N$ = 0.406, while in the bottom half of the domain $\omega_{1}/N$
= 0.214 and $\omega_{2}/N$ = 0.403.

Figure 7: (a) $\beta_{z}$ of the density perturbation field for the experiment
presented in Figure 4 at $t/T_{0}$ = 481. (b)–(d) The real part of the three
pairs of dominant modes extracted using DMD over a time-window from 481 $\leq
t/T_{0}\leq$ 483\. Specifically, (b) $\mathbbm{B}_{0}$ with $\omega_{0}/N$=
0.617 (c) $\mathbbm{B}_{1}$ with $\omega_{1}/N$ = 0.211 (d) $\mathbbm{B}_{2}$
with $\omega_{2}/N$ = 0.405. The black dashed box in (d) indicates the region
of discontinuity in the $\mathbbm{B}_{2}$ beam.

The spatial dependence of the instability is highlighted further in Figure 8,
where the three-dimensional surface plots shows the horizontal components of
the wave vectors ${l_{1}}$ and ${l_{2}}$ as a function of height in the domain
for 8 different instances in time, given by the white dashed lines on Figure
6(a). The surface is defined by

$S_{\beta_{z}^{\ddagger}}(l,z,t)=\bigg{|}\int_{-\infty}^{+\infty}\beta_{z}^{\ddagger}(x^{\prime},z,t)e^{-i2\pi
lx}W(x^{\prime}-x;X_{X})dx^{\prime}\bigg{|}^{2},$ (14)

where $W(x^{\prime};X_{X})$ is a Hamming window of width
$X_{X}=x/\lambda_{x_{0}}=$ 9.7, spanning the full width of the domain. Here,
$\beta_{z}^{\ddagger}(x,z,t)$ corresponds to the instantaneous vertical
density perturbation gradient that has already been temporally filtered in
Fourier space to remove the signal from $\omega_{0}$. The surface plots show
(14) evaluated at each height in the domain to obtain the horizontal component
of wavenumber $l_{1}$ and $l_{2}$. The contour plots behind show the
corresponding frequency-wavenumber spectrogram. This is obtained by the two-
dimensional Fourier transform (in $x$ and $t$)

$S_{\beta_{z}^{\ddagger}}(\omega,l,t)=\bigg{\langle}\bigg{|}\frac{1}{T_{T}}\int\int_{-\infty}^{+\infty}\beta_{z}^{\ddagger}(x^{\prime},z,t^{\prime})e^{-i2\pi(lx^{\prime}+\omega
t^{\prime})}W(x^{\prime}-x;X_{X})W(t^{\prime}-t;T_{T})dx^{\prime}dt^{\prime}\bigg{|}^{2}\bigg{\rangle}_{z},$
(15)

where the widths of the Hamming windows are given by $T_{T}/T_{0}$ = 39 and
$X_{X}=x/\lambda_{x_{0}}=$ 9.7, and the subscript $z$ on the angle brackets
shows that $S_{\beta_{z}^{\ddagger}}(\omega,l,t)$ is averaged over the height
of the domain. The region of spatio-temporal discontinuity shown in physical
space by the black dashed rectangle in Figure 7(d) is clearly visible in
wavenumber space in Figure 8(e). Examining the peak of the spectral
isosurface, $S_{\beta_{z}^{\ddagger}}(l,z,t)$, corresponding to $l_{2}$,
around mid-height in the domain, there is a shift in both the amplitude and
value of $l_{2}$ where the peak occurs. The presence of this discontinuous
region indicates that two wave beams, of slightly different frequency and
wavenumber, are destructively interfering with each other. Later, at $t/T_{0}$
= 547 in (f), the triadic interaction in the lower part of the domain has
decayed (as there is only a very low amplitude signal for both $l_{1}$ and
$l_{2}$), while the interaction occurring in the upper region of the domain is
still present. This continuously varying range of wavenumbers explains why the
grey region of experimentally obtained characteristic wavenumbers on Figure 3
does not exactly fit the spatial triadic conditions of the underlying green
branch of the loci.

Figure 8: Eight surface plots for eight moments in time of the experiment
presented in Figure 4, showing the distribution of $l_{1}$ and $l_{2}$ over
the non-dimensional height in the domain ($z/H$) calculated by (14). Here the
$\beta_{z}$ image sequence is first temporally filtered in Fourier space to
remove the signal from $\mathbbm{B}_{0}$, and thus we do not see a peak at
$l_{0}$. Both the surface plot colour and the height of the peaks, show the
power spectral density $S_{\beta_{z}^{\ddagger}}(l,z)$. The background plot
then shows $l_{1}$ and $l_{2}$ at the same instant in time in the Fourier
plane of horizontal wavenumber component and of frequency. This contour plot
is defined by (15). The timings of each image are given by the white dashed
lines in Figure 6.

We speculate that the reason for these modulations – observed in both real and
Fourier space – is due to the finite-width of the primary wave beam. As a
packet of energy in $\mathbbm{B}_{1}$ or $\mathbbm{B}_{2}$ exits the
underlying primary beam, the energy exchange between the triad is broken. The
time taken for both these secondary beams to exit the spatial confines of
$\mathbbm{B}_{0}$ is dependant on the group velocities of the beams, which are
functions of their wavenumbers, and the relative orientation of the beams
determined by their frequencies. If the secondary beams are unable to extract
sufficient energy before propagating out of the primary beam, the triad system
will not be able to form a stable equilibrium and another triadic perturbation
will grow in another location. Moreover, all the triadic beams are comprised
of a broadband wavenumber spectrum due to their finite-width, as indicated in
Figure 8. This introduces a range of group velocities in the secondary beams
which will exit the underlying beam at different times, enhancing the unsteady
transfer of energy.

Additionally, the structure of the underlying $\mathbbm{B}_{0}$ varies across
the height of the domain. As $\mathbbm{B}_{0}$ propagates upwards through the
tank, it decays due to viscosity, resulting in a broadening in spatial extent
and reduction in amplitude (Fan & Akylas, 2020). These combine to give
considerable variation in both real and Fourier space over the height of the
domain, where different locations will favour slightly different triadic
perturbations. Indeed, as different perturbations grow, the secondary wave
beams with very similar frequencies could interact with each other non-
linearly via the primary wave beam, generating a slow ‘beating’ effect. This
interaction could cause the secondary beams to decay in some locations, while
in others it causes a growth in amplitude, amplifying the effects of the
modulations. Making the approximation that there is a single discrete set of
parameters corresponding to the secondary wave beam for the whole domain is
therefore an oversimplification that ignores the spatial variation of the
instability.

From the experimental results presented above, we believe that the unsteady
behaviour of the instability is a function of the finite-width $\Lambda_{0}$
of the primary beam $\mathbbm{B}_{0}$. We therefore seek to understand this
interaction in a two-dimensional context. We pursue this through the
development of a two-dimensional weakly non-linear model, which we will refer
to as $\mathcal{M}_{2\textrm{D}}$. Details of its development are outlined in
the following section. The goal here is to dissect the experiments and to
isolate the dynamics that are observed experimentally, in order to improve the
understanding of the system. A computational fluid dynamics (CFD) code would
be an inappropriate choice to achieve this, as little would be learnt about
the physical mechanisms governing the behaviour. In $\S$ 4.1 we examine the
perturbation expansion used and in $\S$ 4.2 we look at the numerical solution
to the obtained system of equations.

##  4 Weakly non-linear model construction

###  4.1 Perturbation expansion

Assuming a two-dimensional incompressible continuously stratified Boussinesq
fluid in background hydrostatic balance with constant buoyancy frequency $N$,
the full Boussinesq non-linear equations of motion (comprised of the momentum,
continuity and mass conservation equations) are given by

$\displaystyle\varrho\frac{D\boldsymbol{u}}{Dt}$
$\displaystyle=-\boldsymbol{\nabla}p-g\varrho\mathbf{\hat{z}}+\mu\boldsymbol{\nabla}^{2}\boldsymbol{u},$
(16) $\displaystyle\frac{D\rho}{Dt}$ $\displaystyle=-w\frac{d\bar{\rho}}{dz},$
(17) $\displaystyle\boldsymbol{\nabla}\cdot\boldsymbol{u}$ $\displaystyle=0,$
(18)

where $p$ is the dynamic pressure and $\mu$ is the dynamic viscosity. Writing
velocity in terms of the scalar stream function
$\boldsymbol{u}=\boldsymbol{\nabla}\times(\Psi\hat{\boldsymbol{y}})$, this
system can be reduced to

$\displaystyle\frac{D^{2}}{Dt^{2}}(\boldsymbol{\nabla}^{2}\Psi)-\frac{D}{Dt}\nu\boldsymbol{\nabla}^{2}(\boldsymbol{\nabla}^{2}\Psi)$
$\displaystyle+N^{2}\frac{\partial^{2}\Psi}{\partial
x^{2}}=\frac{g}{\varrho_{0}}\bigg{(}\frac{\partial}{\partial
x}\frac{D}{Dt}-\frac{D}{Dt}\frac{\partial}{\partial{x}}\bigg{)}\rho,$ (19)
$\displaystyle\frac{D\rho}{Dt}$ $\displaystyle=-\frac{\partial\Psi}{\partial
x}\frac{d\bar{\rho}}{dz}.$ (20)

Here, we focus on (19), which represents the non-linear momentum balance in
terms of stream function and density. We seek the simplest form of the stream
function and density that will describe the behaviour of TRI in a finite-width
beam. Specifically we define

$\Psi=\tilde{\Psi}(x,z,t)e^{i(\boldsymbol{k}\cdot\boldsymbol{x}-\omega
t)}+\textrm{c.c.},$ (21)
$\rho=\tilde{\rho}(x,z,t)e^{i(\boldsymbol{k}\cdot\boldsymbol{x}-\omega
t)}+\textrm{c.c.},$ (22)

where $\tilde{\Psi}$ and $\tilde{\rho}$ are the reduced forms of the stream
function and density respectively and c.c. represents the complex conjugate.
Both $\tilde{\Psi}$ and $\tilde{\rho}$ are given as functions of space and
time because, based on the experimental results, the behaviour of TRI in a
finite-width beam is spatio-temporally dependant. We note that in the limit of
linear waves we have

$\tilde{\rho}=-\frac{l}{\omega}\frac{N^{2}\varrho_{0}}{g}\tilde{\Psi},$ (23)

a relationship that also holds for non-linear plane-waves. This relationship
will prove useful later in our exploration of the weakly non-linear
interactions.

We define the three non-dimensional parameters

$\displaystyle\epsilon$
$\displaystyle=(|\tilde{\Psi}_{00}|\kappa_{0}^{2})N^{-1},$ (24)
$\displaystyle\gamma$ $\displaystyle=(\kappa_{0}\Lambda_{0})^{-1},$ (25)
$\displaystyle\chi$
$\displaystyle=(\nu\kappa_{0})\boldsymbol{c}_{g_{0}}^{-1},$ (26)

where, $|\tilde{\Psi}_{00}|$ is the characteristic magnitude of the stream
function associated with the primary beam $\mathbbm{B}_{0}$ and $\nu=\mu/\rho$
is the kinematic viscosity. Here, $\epsilon$ can be viewed as a non-
dimensional measure of the beam amplitude, which characterises the relative
importance of the non-linear $\boldsymbol{u}\cdot\boldsymbol{\nabla}$ terms in
the momentum and conservation of mass equations. The spatial parameter
$\gamma$ defines the separation between the dominant wavelength and overall
width of the primary beam. Finally, $\chi$ is the inverse of the Reynolds
number defined earlier. In oceanographic settings, $\chi$ will be orders of
magnitude smaller than $\epsilon$ and $\gamma$ as, at the scales of interest,
the role of viscosity in the ocean can be considered negligible. Here,
however, as experiments are inherently more viscous than similar motion
patterns at oceanic scales, we retain the leading order effects of viscosity.
We shall restrict our attention to low amplitude $\epsilon\sim 10^{-1}\ll 1$,
broad beam $\gamma\sim 10^{-1}\ll 1$ and low viscosity $\chi\sim 10^{-2}\ll
1$, and introduce re-scaled time and position through the five variables

$\begin{split}\tau_{\gamma}=\gamma t,\qquad\xi&=\gamma x,\qquad\zeta=\gamma
z,\\\ \qquad\tau_{\epsilon}=\epsilon t,&\qquad\tau_{\chi}=\chi t.\\\
\end{split}$ (27)

As we shall see, $\tau_{\epsilon}$ accounts for the ‘slow non-linear time’
variations to the amplitude, while $\tau_{\gamma}$ governs the ‘slow advection
time’ scale. As, for the experiments, both $\epsilon\ll 1$ and $\gamma\ll 1$,
these scaled times account for change over long-time periods. The final time-
scale given by $\tau_{\chi}$ represents the viscous decay of the wave beams.
The spatial parameters $\xi$ and $\zeta$ account for the gradual spatial
variability of the wave beams in the domain. Utilising the dimensionless
amplitude $\epsilon$, we re-write the stream function in (21) as

$\Psi=\tilde{\Psi}(x,z,t)e^{i(\boldsymbol{k}\cdot\boldsymbol{x}-\omega
t)}+\textrm{c.c.}=\epsilon\breve{\Psi}(\tau_{\gamma},\xi,\zeta,\tau_{\epsilon},\tau_{\chi})e^{i(\boldsymbol{k}\cdot\boldsymbol{x}-\omega
t)}+\textrm{c.c.},$ (28)

where $\tilde{\Psi}=\epsilon\breve{\Psi}$. We therefore require
$|\breve{\Psi}|\sim$ 1 as
$\tilde{\Psi}=\epsilon\breve{\Psi}=(|\tilde{\Psi}_{00}|\kappa_{0}^{2}/N)\breve{\Psi}$
and we are interested in small perturbations to the flow. The ‘fast time’ is
associated with the phase variations of the waves given by the wave frequency
and is captured by the complex exponential wave form
$e^{i(\boldsymbol{k}\cdot\boldsymbol{x}-\omega t)}=e^{i\phi}$. As we are
interested in the wave triad, we express the stream function as the summation
in the same way as McEwan & Plumb (1977), by

$\Psi=\sum_{p=0}^{2}{\Psi}_{p}=\sum_{p=0}^{2}\tilde{\Psi}_{p}(x,z,t)e^{i\phi_{p}}+\textrm{c.c.}=\sum_{p=0}^{2}\epsilon\breve{\Psi}_{p}(\tau_{\gamma},\xi,\zeta,\tau_{\epsilon},\tau_{\epsilon\gamma},\tau_{\chi})e^{i\phi_{p}}+\textrm{c.c.},$
(29)

where the subscript $p$ indicates a locally plane-wave approximation to
$\mathbbm{B}_{0}$, $\mathbbm{B}_{1}$ or $\mathbbm{B}_{2}$. Each wave phase
$e^{i\phi_{p}}$ therefore represents the characteristic frequency and
wavenumber contribution to $\mathbbm{B}_{p}$. A similar set of expressions can
be written for $\rho$ as a sum of $\epsilon\breve{\rho}_{p}e^{i\phi_{p}}$.

The superposition in (29) is then substituted into the non-linear equation
(19) where, due to the separate space and time-scales, the partial derivatives
in (19) with respect to $(x,z,t)$ become

$\displaystyle\begin{split}\frac{\partial\Psi_{p}}{\partial
x}&=\epsilon\bigg{(}il_{p}+\gamma\frac{\partial}{\partial\xi}\bigg{)}\breve{\Psi}_{p}e^{i\phi_{p}},\end{split}$
(30a) $\displaystyle\begin{split}\frac{\partial\Psi_{p}}{\partial
z}&=\epsilon\bigg{(}im_{p}+\gamma\frac{\partial}{\partial\zeta}\bigg{)}\breve{\Psi}_{p}e^{i\phi_{p}}.\end{split}$
(30b) $\displaystyle\begin{split}\frac{\partial\Psi_{p}}{\partial
t}&=\epsilon\bigg{(}-i\omega_{p}+\gamma\frac{\partial}{\partial\tau_{\gamma}}+\epsilon\frac{\partial}{\partial\tau_{\epsilon}}+\chi\frac{\partial}{\partial\tau_{\chi}}\bigg{)}\breve{\Psi}_{p}e^{i\phi_{p}},\end{split}$
(30c)

This substitution and following manipulations were performed with the aid of
Mathematica (Wolfram Research, 2021) to ensure reliability of the lengthy
algebraic manipulations required.

We note that at first order in $\epsilon$ the right-hand side (RHS) of (19)
vanishes. Therefore, as the contributions from density on the RHS first appear
at $\mathcal{O}(\epsilon^{2})$, the linear relationship in (23) is valid up
until order $\mathcal{O}(\epsilon^{2})$. The resultant expression obtained
provides the Boussinesq viscous equations of motion solely as a function of
$\Psi$ that can be used to examine both linear and non-linear dynamics between
a triadic set of waves simply by collecting around orders of $\epsilon$,
$\gamma$ and $\chi$.

At $\mathcal{O}(\epsilon^{0}\gamma^{0})$ and
$\mathcal{O}(\epsilon^{0}\gamma^{1})$ the expression will be zero as these
orders correspond to a state of rest. At order
$\mathcal{O}(\epsilon^{1}\gamma^{0})$, we recover the linear wave solution,
with the non-linearity in (19) vanishing and linear superposition applying
such that the three waves propagate independently. Extracting terms with a
common factor of $e^{i\phi_{p}}$, leaves

$\frac{\omega_{p}}{N}=\frac{l_{p}}{\sqrt{(l_{p}^{2}+m_{p}^{2})}},$ (31)

which is the linear dispersion relationship for internal plane-waves. Indeed,
for a single non-linear plane-wave in the inviscid limit, the RHS of (19)
vanishes while the left hand side returns the dispersion relationship for non-
trivial solutions, regardless of the wave amplitude.

We then collect around the next order $\mathcal{O}(\epsilon^{1}\gamma^{1})$.
Again, as $\epsilon$ is still at first-order, the non-linear terms in (19)
cancel. Looking at the $e^{i\phi_{p}}$ terms, we obtain

$\gamma\frac{\partial\breve{\Psi}_{p}}{\partial\tau_{\gamma}}=\gamma\bigg{(}\frac{l_{p}(\omega_{p}^{2}-N^{2})}{\omega_{p}\kappa_{p}^{2}}\frac{\partial\breve{\Psi}_{p}}{\partial\xi}+\frac{\omega_{p}m_{p}}{\kappa_{p}^{2}}\frac{\partial\breve{\Psi}_{p}}{\partial\zeta}\bigg{)},$
(32)

which after some re-arranging can be expressed as

$\gamma\frac{\partial\breve{\Psi}_{p}}{\partial\tau_{\gamma}}=-\gamma\big{(}\boldsymbol{c}_{g_{p}}\cdot\boldsymbol{\nabla}_{\boldsymbol{\xi}}\big{)}\breve{\Psi}_{p}\rightarrow\frac{\partial\tilde{\Psi}_{p}}{\partial
t}=-\big{(}\boldsymbol{c}_{g_{p}}\cdot\boldsymbol{\nabla}\big{)}\tilde{\Psi}_{p},$
(33)

where
$\boldsymbol{\nabla}_{\boldsymbol{\xi}}=(\partial/\partial\xi,\partial/\partial\zeta)$
and $\boldsymbol{\nabla}=(\partial/\partial x,\partial/\partial z)$. This
linear advection equation shows that the stream function of each wave beam in
the triad is advected at its respective group velocity. It is already well
known that, for small amplitude internal waves, the group velocity
$\boldsymbol{c}_{g}$ is the velocity at which energy is transported (e.g.
Sutherland, 2010). As energy scales with $\sim\Psi^{2}$, the fact that (33)
shows that $\breve{\Psi}$ is also advected by $\boldsymbol{c}_{g}$, is not
altogether surprising.

Before examining the non-linear interaction terms at
$\mathcal{O}(\epsilon^{2})$, the viscous term at $\mathcal{O}(\epsilon\chi)$
needs to be considered in conjugation with the linear advection recovered at
$\mathcal{O}(\epsilon^{1}\gamma^{1})$. While these terms are of lower
magnitude than the terms at $\mathcal{O}(\epsilon^{2})$, based on the
experimental and oceanographic parameters, they govern the viscous decay of
individual wave beams irrespective of the non-linear interactions, making it
appropriate to consider their role in conjunction with the advection. At
$\mathcal{O}(\epsilon\chi)$, looking at terms with a common factor
$e^{i\phi_{p}}$, we obtain

$\chi\frac{\partial\breve{\Psi}_{p}}{\partial\tau_{\chi}}=-\frac{\boldsymbol{c}_{g_{0}}}{\kappa_{0}}\frac{\chi}{2}\kappa_{p}^{2}\breve{\Psi}_{p}\rightarrow\frac{\partial\tilde{\Psi}_{p}}{\partial
t}=-\frac{\nu}{2}\kappa_{p}^{2}\tilde{\Psi}_{p},$ (34)

which shows that the viscous decay of each beam scales as $\kappa_{p}^{2}$. We
choose to combine the evolution on time-scales $\tau_{\gamma}$ and
$\tau_{\chi}$ to obtain, at first order in $\epsilon$, the advection equation

$\frac{\partial\tilde{\Psi}_{p}}{\partial
t}=-\big{(}\boldsymbol{c}_{g_{p}}\cdot\boldsymbol{\nabla}\big{)}\tilde{\Psi}_{p}-\frac{\nu}{2}\kappa_{p}^{2}\tilde{\Psi}_{p}.$
(35)

This advection equation in (35) is solved in the 2D advection component of the
$\mathcal{M}_{2\textrm{D}}$ model and its numerical implementation is
addressed in $\S$ 4.2.

We next consider terms at order $\mathcal{O}(\epsilon^{2}\gamma^{0})$. Here,
non-linearity enters the problem and terms are no longer only associated with
$e^{i\phi_{p}}$, but are also comprised of cross-terms from
$\boldsymbol{u}\cdot\boldsymbol{\nabla}$ operator in (19), which have the form
$e^{i(\phi_{q}+\phi_{r})}$, when expressed in Fourier modes (where $p,q,r$ are
permutations of 0, 1, 2). We consider a triad of wave beams satisfying
resonance conditions (2). Using (23) to eliminate $\breve{\rho}$ from (19), at
$\mathcal{O}(\epsilon^{2})$ we recover the same non-linear interactions given
in Bourget et al. (2013), specifically

$\displaystyle\epsilon^{2}\frac{\partial\breve{\Psi}_{0}}{\partial\tau_{\epsilon}}=\epsilon^{2}I_{0}\breve{\Psi}_{1}\breve{\Psi}_{2}\rightarrow\frac{\partial\tilde{\Psi}_{0}}{\partial
t}=I_{0}\tilde{\Psi}_{1}\tilde{\Psi}_{2},$ (36a)
$\displaystyle\epsilon^{2}\frac{\partial\breve{\Psi}_{1}}{\partial\tau_{\epsilon}}=\epsilon^{2}I_{1}\breve{\Psi}_{0}\breve{\Psi}_{2}^{*}\rightarrow\frac{\partial\tilde{\Psi}_{1}}{\partial
t}=I_{1}\tilde{\Psi}_{0}\tilde{\Psi}_{2}^{*},$ (36b)
$\displaystyle\epsilon^{2}\frac{\partial\breve{\Psi}_{2}}{\partial\tau_{\epsilon}}=\epsilon^{2}I_{2}\breve{\Psi}_{0}\breve{\Psi}_{1}^{*}\rightarrow\frac{\partial\tilde{\Psi}_{2}}{\partial
t}=I_{2}\tilde{\Psi}_{0}\tilde{\Psi}_{1}^{*},$ (36c)

where an asterisk indicates the complex conjugate and the interaction term is
given as

$I_{p}=\frac{l_{q}m_{r}-m_{q}l_{r}}{2\omega_{p}\kappa_{p}^{2}}\bigg{[}\omega_{p}(\kappa_{q}^{2}-\kappa_{r}^{2})+l_{p}N^{2}\bigg{(}\frac{l_{q}}{\omega_{q}}-\frac{l_{r}}{\omega_{r}}\bigg{)}\bigg{]}.$
(37)

It is possible to extend this expansion to examine the higher order
$\mathcal{O}(\epsilon^{2}\gamma^{1})$, however as this model sufficiently
captures enough of the observed experimental behaviour, this further expansion
was not required. Above $\mathcal{O}(\epsilon^{2}\gamma^{1})$ (higher orders
in $\epsilon$) the linear approximation in (23) is no longer valid for
eliminating $\breve{\rho}$ from (19). The evolution of these coupled ODEs in
(36), are considered on their own and in $\S$ 5.1 and as part of the
$\mathcal{M}_{2\textrm{D}}$ model in $\S$ 5.2.

###  4.2 Numerical implementation

This section outlines the development of the $\mathcal{M}_{2\textrm{D}}$
numerical model, built to solve the equations obtained at order
$\mathcal{O}(\epsilon^{1}\gamma^{1})$ and at order
$\mathcal{O}(\epsilon^{2})$. We start with the numerical scheme used to solve
the advection equation (35), obtained at order
$\mathcal{O}(\epsilon^{1}\gamma^{1})$. Specifically, we use a monotonic
second-order upwind finite volume scheme to advect the complex stream
function. The finite volume method discretizes the governing equations into
arbitrary control volumes around each node and the advective fluxes are then
evaluated across the upwind faces of each control volume. For internal waves,
their advection velocity is determined by their relative group velocity. As
$\boldsymbol{c}_{g_{p}}$ in (35) is specific to each wave beam, $p$, the
number of numerical domains corresponds to the number of wave beams being
considered, as, up to $\mathcal{O}(\epsilon^{2})$, we only need to consider
$\tilde{\Psi}$ and need not explicitly consider $\tilde{\rho}$. There is no
limit to the number of superposed domains that can be evolved simultaneously,
provided there are suitable interaction terms that couple them. For the
results in this paper, we only consider three domains, allowing us to model a
three-beam system. Limiting the number of domains to three means we fix the
Fourier components of the triad. This turns out to be sufficient to capture
the amplitude modulations observed in our experiments, though it excludes
frequency and wavenumber fluctuations seen experimentally in the secondary
wave beams.

We define each domain by $\tilde{\mathbbm{B}}_{p}$, as we are advecting the
reduced stream function $\tilde{\Psi}_{p}$ in (35). Being a second order
scheme, the volume flux, $F$, is calculated using the two upstream values of
the reduced stream function. As $\tilde{\Psi}_{p}$ is complex, the imaginary
and real parts are advected separately and combined back into a complex value
after being advected. Figure 9 shows a schematic of the domain
$\tilde{\mathbbm{B}}_{0}$. To avoid confusion with indexing, compass co-
ordinates are used for the nodes and the subscript $p=0$ is dropped from
$\tilde{\Psi}_{0}$ in the equations below for clarity. Specifically, for
domain $\tilde{\mathbbm{B}}_{0}$, $F_{w}=c_{g_{0_{x}}}\tilde{\Psi}_{w}$ and
$F_{n}=c_{g_{0_{z}}}\tilde{\Psi}_{n}$, where

$\displaystyle\tilde{\Psi}_{w}=\tilde{\Psi}_{M}+\frac{\Upsilon_{W}}{2}\big{(}\tilde{\Psi}_{M}-\tilde{\Psi}_{E}\big{)}\quad\quad$
$\displaystyle\textrm{for}\quad c_{g_{0_{x}}}<0,\leftarrow$ (38a)
$\displaystyle\tilde{\Psi}_{n}=\tilde{\Psi}_{M}+\frac{\Upsilon_{N}}{2}\big{(}\tilde{\Psi}_{M}-\tilde{\Psi}_{S}\big{)}\quad\quad$
$\displaystyle\textrm{for}\quad c_{g_{0_{z}}}>0,\uparrow$ (38b)

the node indexing is shown in Figure 9 and $\Upsilon$ is the flux limiter
which prevents the generation of oscillations inherent in second-order schemes
(Roe, 1986). The flux limiter used is the ‘min-mod’ function which takes the
ratio of the downstream to upstream gradient. The flux limiter is calculated
for both the real and imaginary parts and the more restrictive value (i.e. the
one closest to enforcing a first order scheme) is used for both components.
Accounting for viscous attenuation, the value of the stream function at the
central node $\tilde{\Psi}_{M}$ is then updated by the differences in the
fluxes over the control volume using

$\tilde{\Psi}_{M}^{\dagger}=\tilde{\Psi}_{M}^{i}+\Delta
t\bigg{(}-\boldsymbol{\nabla}\cdot\boldsymbol{c}_{g_{0}}\tilde{\Psi}_{M}^{i}-\frac{\nu}{2}\kappa_{0}^{2}\tilde{\Psi}_{M}^{i}\bigg{)}=\tilde{\Psi}_{M}^{i}+\Delta
t\Bigg{(}\frac{F_{w}^{i}-F_{e}^{i}}{\Delta
x}+\frac{F_{s}^{i}-F_{n}^{i}}{\Delta
z}-\frac{\nu}{2}\kappa_{0}^{2}\tilde{\Psi}_{M}^{i}\Bigg{)},$ (39)

where the superscripts $i$ and $\dagger$ represent the value of $\tilde{\Psi}$
before and after advection respectively and $\Delta t$ is the time-step over
which code is advanced. For the scheme to remain numerically stable, we limit
the time step, $\Delta t$, to satisfy the CFL condition.

Figure 9: Sketch outlining the numerical advection domain
$\tilde{\mathbbm{B}}_{0}$. The bottom boundary forcing of the wavemaker is
given by the sinusoid envelope in (7). The two outer grid layers required for
the second-order scheme are shown in blue. In this domain
$\boldsymbol{c}_{g_{0}}$ is directed to the north-west, meaning the advection
of the complex $\tilde{\Psi}_{0}$ is to the north-west. For advection,
$\tilde{\Psi}_{0}$ is spilt into its real and imaginary parts and converted
back into a complex value after advection is complete. The right panel then
show the close up of the real flux (front domain) and imaginary flux (behind
domain) over a control volume.

Each domain $\tilde{\mathbbm{B}}_{p}$ has 242 x 82 cells in the horizontal and
vertical respectively, matching the aspect ratio of the experimental
visualisation window. The non-dimensional grid spacing $\Delta
x/\lambda_{x_{0}}=\Delta z/\lambda_{x_{0}}=$ 0.04, is much finer than the
smallest wavelength considered. At every time step, random complex background
noise of magnitude $10^{-7}\delta e^{i(2\pi\vartheta)}$ is added to each cell
in every domain. Both $\delta$ and $\vartheta$ are independent random numbers
spanning the range $[0,\hskip 2.84526pt1]$, giving a mean magnitude of the
order $5\times 10^{-8}$ with a uniformly distributed phase angle. This
background noise is randomly selected at every time step. Without the addition
of a perturbation, the instability would not be triggered. Each domain in the
system hosts a boundary condition tailored for the wave beam it contains. For
the boundaries corresponding to outgoing waves, standard non-reflecting
boundary conditions are used. For the incoming boundaries, small amplitude
complex noise, of the same structure as the background noise is advected into
the domain.

The only domain containing a different inflow boundary is
$\tilde{\mathbbm{B}}_{0}$, where it is necessary to impose bottom boundary
forcing to simulate the wavemaker. As the model considers the advection of the
complex valued reduced stream function $\tilde{\Psi}$, the boundary condition
is given by the envelope of the experimental forcing, which is obtained by
removing the fast time forcing, $e^{i\phi_{0}}$, from (7) and (9). This
boundary condition for $\tilde{\mathbbm{B}}_{0}$ is shown in Figure 9. We
implement a prescribed displacement condition by computing the corresponding
complex valued stream function amplitude $\tilde{\Psi}_{\textrm{in}}$, as
opposed to $\eta_{0}$, where $\tilde{\Psi}_{\textrm{in}}$ =
$|\tilde{\Psi}_{\textrm{in}}|e^{i(2\pi\varsigma)}$. Here, $\varsigma$ is a
constant value within the range $[0,\hskip 2.84526pt1]$. This ensures that the
phase angle of the slowly evolving complex amplitude $\tilde{\Psi}_{0}$ is
fixed at the boundary, in the same way as the experimental wavemaker. At each
time step, after each domain is advected at its respective group velocity, the
non-linear interactions (36) are calculated.

The non-linear interactions are included by converting (36) into the numerical
format

$\displaystyle\tilde{\Psi}_{0}^{i+1}$
$\displaystyle=\tilde{\Psi}_{0}^{\dagger}+\Delta
tI_{0}\tilde{\Psi}_{1}^{\dagger}\tilde{\Psi}_{2}^{\dagger},$ (40a)
$\displaystyle\tilde{\Psi}_{1}^{i+1}$
$\displaystyle=\tilde{\Psi}_{1}^{\dagger}+\Delta
tI_{1}\tilde{\Psi}_{0}^{\dagger}\tilde{\Psi}_{2}^{\dagger*},$ (40b)
$\displaystyle\tilde{\Psi}_{2}^{i+1}$
$\displaystyle=\tilde{\Psi}_{2}^{\dagger}+\Delta
tI_{2}\tilde{\Psi}_{0}^{\dagger}\tilde{\Psi}_{1}^{\dagger*},$ (40c)

where $\tilde{\Psi}_{p}$ corresponds to every cell in domain
$\tilde{\mathbbm{B}}_{p}$ and $\dagger$ and $i+1$ respectively represent the
values of $\tilde{\Psi}_{p}$ in each domain after advection (39) and at time
$t+\Delta t$ after the non-linear interactions are calculated. As the
interaction co-efficient is applied to the whole domain, there is no need to
distinguish here the different cells using compass indexing. In the following
section $\S$ 5.1 we first consider the non-linear interactions on their own,
before examining the results from the $\mathcal{M}_{2\textrm{D}}$ model in
$\S$ 5.2.

##  5 Weakly non-linear behaviour

###  5.1 Zero-dimensional model results

We note that the coupled non-linear equations in(36) can be recovered in zero-
dimensional space by considering the form of the stream function
$\tilde{\Psi}_{p}(t)e^{i\phi_{p}}$ in (19), where the slowly evolving reduced
stream function is considered solely as a function of time. Indeed, this was
first proposed by McEwan & Plumb (1977) and further developed by Koudella &
Staquet (2006) and Bourget et al. (2013), who obtain these coupled ODEs that
govern the development of the triad.

While this theory considers the ‘slow-time’ development of the amplitude of
the beams, as noted by Sutherland (2013), it is still based upon the
assumption that the waves are monochromatic in space and time. In the case of
a finite-width beam becoming unstable to TRI, the secondary wave beams have a
finite time with which to interact with the underlying primary beam. This
limitation was addressed by Bourget et al. (2014), who adapted the equations
in (36) to examine the energy flux across a finite region of the primary wave
beam. Using an energy balance, they define a two-dimensional control area of
width $W$ and length $L$ over which the resonant beams can interact with the
primary. Accounting for the energy flux through this control area via non-
linear interactions, viscous attenuation and incoming and outgoing energy flux
of the primary beam, the ODEs in (36) become

$\displaystyle\frac{d\tilde{\Psi}_{0}}{dt}$
$\displaystyle=I_{0}\tilde{\Psi}_{1}\tilde{\Psi}_{2}-\nu\bigg{(}\frac{\kappa_{0}^{2}}{2}\bigg{)}\tilde{\Psi}_{0}+T,$
(41a) $\displaystyle\frac{d\tilde{\Psi}_{1}}{dt}$
$\displaystyle=I_{1}\tilde{\Psi}_{0}\tilde{\Psi}_{2}^{*}-\bigg{(}\frac{\nu\kappa_{1}^{2}}{2}+\frac{|\boldsymbol{c}_{g_{1}}\cdot\boldsymbol{e}_{k_{0}}|}{2W}\bigg{)}\tilde{\Psi}_{1},$
(41b) $\displaystyle\frac{d\tilde{\Psi}_{2}}{dt}$
$\displaystyle=I_{2}\tilde{\Psi}_{0}\tilde{\Psi}_{1}^{*}-\bigg{(}\frac{\nu\kappa_{2}^{2}}{2}+\frac{|\boldsymbol{c}_{g_{2}}\cdot\boldsymbol{e}_{k_{0}}|}{2W}\bigg{)}\tilde{\Psi}_{2},$
(41c)

where $I_{p}$ is given in (37), $\boldsymbol{e}_{k_{0}}$ is a unit vector in
the direction of $\boldsymbol{k}_{0}$ and the forcing term
$T=|\boldsymbol{c}_{g_{0}}|(\tilde{\Psi}_{\textrm{in}}^{*}\tilde{\Psi}_{\textrm{in}}-\tilde{\Psi}_{0}^{*}\tilde{\Psi}_{0})/(2L\tilde{\Psi}_{0}^{*})$
in (41a), represents the energy flux for the primary wave beam through the
control area with an incoming amplitude of $\tilde{\Psi}_{\textrm{in}}$. We
convert this input amplitude to the non-dimensional measure
$\epsilon_{\textrm{in}}=\kappa_{0}^{2}|\tilde{\Psi}_{\textrm{in}}|/N$. We
recognise that $\epsilon_{\textrm{in}}$ is comparable to the experimental
input amplitude
$\epsilon_{r}=\kappa_{0}^{2}\langle|\tilde{\Psi}_{0}|\rangle_{r}/N$, where
$\langle|\tilde{\Psi}_{0}|\rangle_{r}$ is magnitude of the reduced stream
function averaged over the black domain in Figure 2(b). The terms on the end
of (41b) and (41c) represent the viscous decay within, and flux of energy out
of, the control area. For the remainder of this paper, we will refer to the
above set of spatially zero-dimensional ODEs in (41), which we use to describe
the energy exchange in TRI in the context of a finite-width beam, as the zero-
dimensional model $\mathcal{M}_{0\textrm{D}}$.

The $\mathcal{M}_{0\textrm{D}}$ model is numerically integrated to examine its
prediction for the development of the triad. To match the experimental set-up,
the width $W$ and length $L$ of the interaction region are set to
$\Lambda_{0}$ (defined in (8)) and $2\Lambda_{0}$ respectively (see Bourget et
al. (2014) for details). The parameters for $\mathbbm{B}_{0}$ are also kept
consistent with the experimental ones. In the limit of
$\mathcal{M}_{0\textrm{D}}$, $\mathbbm{B}_{0}$ reduces to

$\mathbbm{B}_{0}=\\{\omega_{0}/N=0.62,\hskip
5.69054pt\lambda_{x_{0}}\boldsymbol{k}_{0}=(-6.28,-8.29)\\}.$ (42)

We are curious to see how varying the wavenumbers and frequencies of the
secondary wave beams impact the evolution of the instability described by the
$\mathcal{M}_{0\textrm{D}}$. This is achieved by defining a resonant triad as
$\mathbbm{T}_{\Phi}=\\{\mathbbm{B}_{0},\mathbbm{B}_{1_{\Phi}},\mathbbm{B}_{2_{\Phi}}\\}$,
where the subscript $\Phi$ corresponds to a specific triadic configuration,
obtained by changing the characteristic frequencies and wavenumbers of the
secondary wave beams. We require that all the triadic configurations satisfy
both the resonant condition (1) and dispersion (4), meaning all of the
configurations lie exactly on the green curve in Figure 3. For the
$\mathcal{M}_{0\textrm{D}}$ model we consider configurations $\mathbbm{T}_{a}$
and $\mathbbm{T}_{d}$, marked by the blue circle and star on Figure 3,
respectively, with parameters shown in Table 1.

Triad configuration $\mathbbm{T}_{\Phi}$ |  $\omega_{1}/N$ |  $\omega_{2}/N$ |  $\lambda_{x_{0}}l_{1}$ |  $\lambda_{x_{0}}l_{2}$ |  $\lambda_{x_{0}}m_{1}$ |  $\lambda_{x_{0}}m_{2}$
---|---|---|---|---|---|---
$\mathbbm{T}_{a}$ | 0.198 | 0.419 | 2.136 |  $-$ 8.671 | 10.304 |  $-$18.598
$\mathbbm{T}_{b}$ | 0.206 | 0.411 | 2.388 |  $-$8.922 | 11.561 |  $-$19.855
$\mathbbm{T}_{c}$ | 0.222 | 0.395 | 3.267 |  $-$9.802 | 14.577 |  $-$22.871
$\mathbbm{T}_{d}$ | 0.227 | 0.390 | 3.644 |  $-$10.179 | 15.708 |  $-$24.127
$\mathbbm{T}_{e}$ | 0.231 | 0.386 | 4.021 |  $-$10.556 | 16.965 |  $-$25.258
$\mathbbm{T}_{f}$ | 0.239 | 0.378 | 4.775 |  $-$11.310 | 19.227 |  $-$27.520
Table 1: The input parameters of $\mathbbm{B}_{1_{\Phi}}$ and
$\mathbbm{B}_{2_{\Phi}}$ for each
$\mathbbm{T}_{\Phi}=\\{\mathbbm{B}_{0},\mathbbm{B}_{1_{\Phi}},\mathbbm{B}_{2_{\Phi}}\\}$,
used in the $\mathcal{M}_{0\textrm{D}}$ and $\mathcal{M}_{2\textrm{D}}$ model.
The wave vector locations of each secondary wave beam pair can be seen by the
blue marks on Figure 3.

The results of the numerical integration for the $\mathcal{M}_{0\textrm{D}}$
model are given in Figure 10(a) and (b) for $\mathbbm{T}_{a}$ and
$\mathbbm{T}_{d}$, respectively, across a range of five non-dimensional input
forcing amplitudes for the primary beam
$0.092\leq\epsilon_{\textrm{in}}\leq$0.184. The resultant amplitudes of the
triadic beams are given in terms of
$\epsilon_{n}=\kappa_{0}^{2}|\tilde{\Psi}_{p}|/N$. In order to draw the most
meaningful comparison between the spatially 0D model and 2D experiments, we
average the experimental results over the whole visualisation window
$\langle\rangle_{w}$ in order to de-correlate the signal with the position of
a beam in space. As the $\mathcal{M}_{0\textrm{D}}$ model is not a function of
space, no spatial averaging is required for $\epsilon_{n}$. While
$\epsilon_{n}$ and $\epsilon_{w}$ are therefore not quantitatively comparable,
they are qualitatively. For details on the different measures of amplitude,
see Table 5.1.

Figure 10: Results of the $\mathcal{M}_{0\textrm{D}}$ model by (Bourget et
al., 2014) given in (41) across a range of forcing amplitudes
$0.092\leq\epsilon_{\textrm{in}}\leq 0.184$. (a) corresponds to the triad
configuration $\mathbbm{T}_{a}$, while (b) corresponds to $\mathbbm{T}_{d}$.
The parameters for these configurations are given in Table 1.

We first consider the general behaviour of the $\mathcal{M}_{0\textrm{D}}$
model, consistent between both Figure 10(a) and (b). For the lowest forcing
amplitude of $\epsilon_{\textrm{in}}$ = 0.092, no instability occurs, so
$|\tilde{\Psi}_{1}|$ (red line) and $|\tilde{\Psi}_{2}|$ (green line) remain
zero while $|\tilde{\Psi}_{0}|$ (blue line) increases with the initial growth
of $|\tilde{\Psi}_{\textrm{in}}|$ and then remains at a constant amplitude.
There then exists a critical forcing amplitude, above which, TRI is observed.
For $\mathbbm{T}_{a}$ (Figure 10(a)) this occurs when
$\epsilon_{\textrm{in}}\geq$ 0.111, while for $\mathbbm{T}_{d}$ (Figure
10(b)), instability occurs when $\epsilon_{\textrm{in}}\geq$ 0.097. While both
of these limits are lower than the experimental limit of $\epsilon_{r}\geq$
0.166, the $\mathcal{M}_{0\textrm{D}}$ still captures the fact that there
exists an amplitude threshold that must be surpassed before a finite-width
beam exhibits TRI. We next note, that for all solutions where TRI is observed,
the amplitude of $|\tilde{\Psi}_{0}|$ always decays to the same asymptotic
value, while the amplitudes of $|\tilde{\Psi}_{1}|$ and $|\tilde{\Psi}_{2}|$
asymptote at higher values as $\epsilon_{\textrm{in}}$ is increased. Upon
investigation, the asymptotic amplitude of $|\tilde{\Psi}_{0}|$ after the
initial decay is equal to the amplitude threshold for instability.
Interestingly, for the larger forcing amplitudes, the approach to the
equilbiurm state takes the form of an under-damped non-linear oscillator,
shown by a small oscillation after the initial onset of the instability.

The difference between how the triad configurations $\mathbbm{T}_{a}$ and
$\mathbbm{T}_{d}$ affect the evolution of the instability is subtle. On
inspection, Figure 10(b) shows a quicker growth of the secondary beams for the
same input amplitude as $\mathbbm{T}_{a}$, a result that agrees with
$\mathbbm{T}_{d}$ having the lower amplitude threshold for instability.
Strangely, despite this quicker growth, the resultant amplitudes of the
secondary wave beams are smaller than in Figure 10(a). This is unintuitive, as
one would expect that a larger decay of the primary beam would result in
larger amplitudes of the secondary beams. As the viscous decay along a beam
scales with $\kappa^{3}$, this lower value of the secondary beams might be due
to greater viscous dissipation in $\mathbbm{T}_{d}$, due to
$\kappa_{d}>\kappa_{a}$.

non-dimensional |  |  |  |
---|---|---|---|---
amplitude |  Experimental (E) Modelling ($\mathcal{M}_{2\textrm{D}}$ or $\mathcal{M}_{0\textrm{D}}$) |  dimensional |  |
|  units |  Description |  |
$\eta_{0}\lambda_{x_{0}}$ | E | $\eta_{0}$ | mm |  Imposed half peak to peak displacement
|  |  |  | from the wavemaker
$\epsilon_{r}$ | E | $\langle|\tilde{\Psi}_{0}|\rangle_{r}$ | mm s-2 |  Spatially averaged measure of $\mathbbm{B}_{0}$ over
|  |  |  | region shown in Figure 2(b)
$\epsilon_{w}$ | E $\&$ $\mathcal{M}_{2\textrm{D}}$ | $\langle|\tilde{\Psi}_{p}|\rangle_{w}$ | mm s-2 |  Spatially averaged measure of all wave beam fields over their respective domains
$\epsilon$ | $\mathcal{M}_{2\textrm{D}}$ | $|\tilde{\Psi}_{00}|$ | mm s-2 |  Characteristic measure of primary wave amplitude for perturbation expansion
$\epsilon_{\textrm{in}}$ |  $\mathcal{M}_{2\textrm{D}}$ $\&$ $\mathcal{M}_{0\textrm{D}}$ | $|\tilde{\Psi}_{\textrm{in}}|$ | mm s-2 |  Imposed amplitude input for $\tilde{\mathbbm{B}}_{0}$ ($\mathcal{M}_{2\textrm{D}}$) and $\mathbbm{B}_{0}$ ($\mathcal{M}_{0\textrm{D}}$)
$\epsilon_{n}$ | $\mathcal{M}_{0\textrm{D}}$ | $|\tilde{\Psi}_{p}|$ | mm s-2 |  Measure of amplitude for all wave fields
|  |  |  | from $\mathcal{M}_{0\textrm{D}}$

Table 2: Details of the different non-dimensional measures of amplitude for
clarity.

Comparing Figure 10 with the experimental results in Figure 5, we see
remarkably different behaviour in the evolution of the instability. For the
$\mathcal{M}_{0\textrm{D}}$ model results, larger forcing engenders a
significant decay in amplitude of the primary wave beam. This large decay of
the primary beam is not observed experimentally, where the amplitude of the
primary beam oscillates around a mean value comparable with that set by the
forcing from the wavemaker. The most obvious difference between the two then
arises in their description of the long-term development. The
$\mathcal{M}_{0\textrm{D}}$ model results predict that after the initial
instability, the energy exchange between the triad saturates to a steady
equilibrium, set by the non-linear interaction term $I$. This is clearly not
the case for the experimental results. The slow synchronous amplitude
modulations seen experimentally reveal a continuous fluctuation in the energy
exchange between the primary and the two secondary beams.

As the $\mathcal{M}_{0\textrm{D}}$ model does not consider the triadic
interaction as a function of space it is unable to describe the modulations
witnessed experimentally. We next consider the results of the
$\mathcal{M}_{2\textrm{D}}$ model to see how it describes the instability when
considered in a two-dimensional domain.

###  5.2 Two-dimensional model results

We now consider the results from the $\mathcal{M}_{2\textrm{D}}$ model, where
in each simulation we provide three domains – one for each triadic beam – with
input parameters chosen to satisfy the triadic condition (2) and dispersion
(4). Six triad configurations
$\mathbbm{T}_{\Phi}=\\{\mathbbm{B}_{0},\mathbbm{B}_{1_{\Phi}},\mathbbm{B}_{2_{\Phi}}\\}$
are considered, with wavenumber vectors distributed across the solid green
loci branch in Figure 3 and parameters given in Table 1. The parameters for
$\mathbbm{B}_{0}$ (common to all six triads), match the experiments and are
given in (42).

Figure 11 shows the results of 36 simulations, where sub-plots (a) to (f)
correspond to triad configurations $\mathbbm{T}_{a}$ to $\mathbbm{T}_{f}$
respectively, each shown with six input amplitudes, $\epsilon_{\textrm{in}}$.
Again, $\epsilon_{\textrm{in}}$ is comparable to the experimental values
$\epsilon_{r}$. Looking at all of the plots in Figure 11, it is relatively
easy to categorise three different behavioural evolutions of the triad
simulations. The first behavioural evolution, observed for all the simulations
using triad configurations $\mathbbm{T}_{a}$ and $\mathbbm{T}_{b}$ (shown by
the blue circle and cross on Figure 3), is when no growth of the secondary
wave beams occur and the system remains as a single stable primary beam. The
reason for this is twofold: both triad configurations satisfy (or nearly
satisfy in the case of $\mathbbm{T}_{b}$) $\kappa_{1}/\kappa_{0}\leq 1$ and
both correspond to the smallest values of $\omega_{1}$ considered. Indeed, for
input amplitudes that caused the other triadic configurations to become
unstable, we found that no pairs with $\kappa_{1}/\kappa_{0}<1$ generated TRI
for the $\mathcal{M}_{2\textrm{D}}$ model. This observation is further
reinforced by our experiments: the grey shaded region on Figure 3 shows that
all observed triadic combinations of wave vectors satisfy
$\kappa_{1}/\kappa_{0}>1$. When $\kappa_{1}/\kappa_{0}<1$, we have the
condition $\kappa_{1}<\kappa_{0}<\kappa_{2}$. Thus energy transfers to both
the larger length scale ($\kappa_{1}$) and the smaller ($\kappa_{2}$). Bourget
et al. (2014) show that the linear growth rate from the
$\mathcal{M}_{0\textrm{D}}$ model is larger for this triad configuration when
$\Lambda_{0}<7\lambda_{0}$. However, in our experiments where
$\Lambda_{0}\approx 3\lambda_{0}$, it is still a triad configuration with
$\kappa_{1}/\kappa_{0}>1$ (and therefore $\kappa_{0}<\kappa_{1},\kappa_{2}$,
located on the outer solid green branch of Figure 3), that is selected. The
effect of having a smaller value of $\omega_{1}$ is discussed further below.

Figure 11: Amplitude plots generated using the two-dimensional weakly non-
linear $\mathcal{M}_{2\textrm{D}}$ model. Each sub-plot shows a range forcing
amplitudes for $\mathbbm{B}_{0}$ between $0.092\leq\epsilon_{\textrm{in}}\leq
0.184$. The difference between each sub-plot is the parameters for the
secondary wave beams in each triad, given in Table 1. (a) $\mathbbm{T}_{a}$
(b) $\mathbbm{T}_{b}$ (c) $\mathbbm{T}_{c}$ (d) $\mathbbm{T}_{d}$ (e)
$\mathbbm{T}_{e}$ (f) $\mathbbm{T}_{f}$. The run time is
$t_{\textrm{end}}/T_{0}=816$. The resultant non-dimensional amplitudes are
calculated over the whole visualisation window
$\epsilon_{w}=\kappa_{0}^{2}\langle|\tilde{\Psi}|\rangle_{w}/N$.

The second behavioural evolution seen in Figure 11, is when TRI occurs and the
amplitudes of all beams undergo coupled modulations, similar to those seen
experimentally in Figure 5. This behaviour is observed for many of the
simulations using triad configurations $\mathbbm{T}_{c}$ to $\mathbbm{T}_{f}$.
For these configurations, we see that an amplitude threshold must be surpassed
before instability can occur. We note a very close agreement in the amplitude
threshold required for instability between $\mathcal{M}_{2\textrm{D}}$ and the
experimental values. For the triad configurations $\mathbbm{T}_{c}$ to
$\mathbbm{T}_{f}$, instability occurred when $\epsilon_{\textrm{in}}\geq$
0.161, while experimentally $\epsilon_{r}\geq$ 0.166 triggered instability. We
also notice that for all of these simulations that become unstable, as we
increase $\epsilon_{\textrm{in}}$, not only do the secondary wave beam
amplitudes increase, their growth also occurs at earlier times. This
observation is consistent with the $\mathcal{M}_{0\textrm{D}}$ model. We can
measure this initial linear growth using the growth rate $\sigma$, of the form
$e^{\sigma t}$, to characterise how quickly the instability develops. This
growth rate term will prove useful later.

Focusing on a specific simulation with coupled amplitude modulations, we
consider triadic configuration $\mathbbm{T}_{d}$, forced at a non-dimensional
input amplitude of $\epsilon_{\textrm{in}}$ = 0.161 (Figure 11(d)). Figure 12
shows six instantaneous images from this simulation obtained by the
superposition of the three domains multiplied by their respective fast time
and short length scales $e^{i\phi_{p}}$. Here we see the initial development
of $\mathbbm{B}_{1}$ occurring at the top of $\mathbbm{B}_{0}$, where it grows
in strength. (We note that calculations preformed in a larger domain
demonstrated that the generation region of $\mathbbm{B}_{1}$ occurs at the
same location, showing that it is not an effect from the boundary.) As with
the experimental images in Figure 4, due to the similar alignment and
direction of $\boldsymbol{k}_{0}$ and $\boldsymbol{k}_{2}$, $\mathbbm{B}_{2}$
is not obvious in this visualisation region as it propagates within the
confines of $\mathbbm{B}_{0}$. Over time, $\mathbbm{B}_{1}$ grows in both
amplitude and width before decaying in a quasi-periodic manner. Unlike the
experimental results, however, its generation region remains approximately
fixed and does not traverse the height of $\mathbbm{B}_{0}$. By construction,
$\omega_{1}$ and $\omega_{2}$ also remain fixed.

$\delta_{1}$ Figure 12: Sequence of images at $t/T_{0}=61$ apart, showing the
superposition of the three domains in the model, using triad configuration
$\mathbbm{T}_{d}$. Each domain is multiplied by its respective short length
and fast time-scales $e^{i\phi_{p}}$. The forcing amplitude for
$\mathbbm{B}_{0}$ is $\epsilon_{\textrm{in}}$ = 0.161. The spatially averaged
amplitude of each wave beam over time is shown in Figure 11(d). The timing of
each image is (a) $t/T_{0}$ = 121, (b) $t/T_{0}$ = 151, (c) $t/T_{0}$ = 272,
(d) $t/T_{0}$ = 333, (e) $t/T_{0}$ = 424, (f) $t/T_{0}$ = 484. The black line
in (f) marks the length $\delta_{1}$, which defines the spatial distance over
which $\mathbbm{B}_{1}$ can extract energy from $\mathbbm{B}_{0}$.

The final behavioural evolution seen in Figure 11, most obvious for the
largest forcing amplitude using $\mathbbm{T}_{f}$, is when the triadic system
reaches a stable equilibrium, closely resembling the steady state results from
$\mathcal{M}_{0\textrm{D}}$ shown in Figure 10. Here, the amplitudes of the
triadic beams do not exhibit any modulations and the triad quickly reaches a
stable equilibrium, after a large, smooth decay in amplitude of the primary
beam.

The only parameters being varied across these simulations in Figure 11 are the
frequencies and wavenumbers of the secondary beams in the triad and the
amplitude of the primary beam. This varying behavioural evolution of the
triadic beams must, therefore, be due to these changing parameters. Based on
the dispersion relationship (4), a greater value of $\omega_{1}$ results in a
$\mathbbm{B}_{1}$ beam with closer alignment to $\mathbbm{B}_{0}$. This
steeper angle leads to a greater spatial region over which $\mathbbm{B}_{1}$
can extract energy from $\mathbbm{B}_{0}$ and consequently an increased
distance for $\mathbbm{B}_{1}$ to grow. This is shown schematically in Figure
13. We define the lengths of these interaction regions for $\mathbbm{B}_{1}$
and $\mathbbm{B}_{2}$ as $\delta_{1}$ and $\delta_{2}$, respectively. Figure
14(a) shows the relationship between $\delta_{1}$ and $\omega_{1}$, together
with the relationship between $\delta_{2}$ and $\omega_{2}$. The grey shading
in the background marks the range of frequencies obtained experimentally.

As $\omega_{1}$ increases, not only does $\delta_{1}$ increase,
$|\boldsymbol{c}_{g_{1}}|$ also decreases, as shown in Figure 14(b). A
decrease in $|\boldsymbol{c}_{g_{1}}|$ causes $\mathbbm{B}_{1}$ to remain
within the spatial confines of $\mathbbm{B}_{0}$ for longer and hence
increases the time in which it can extract energy. The red and green shaded
regions on Figure 14(b), associated with $|\boldsymbol{c}_{g_{1}}|$ and
$|\boldsymbol{c}_{g_{2}}|$, respectively, mark the range of non-dimensional
group velocities contained within the secondary beams due to their broadband
spectrum. While there is only one wavenumber vector $\boldsymbol{k}_{1}$ at
each $\omega_{1}$ that can satisfy both dispersion (4) and the triadic
resonant condition (2), (shown by the loci on Figure 3), both
$\mathbbm{B}_{1}$ and $\mathbbm{B}_{2}$ are beams that are broadly distributed
over the wavenumber spectrum due to their finite-width and so there will be
exact triads selected from this distribution. A wavenumber distribution was
clearly shown experimentally in Figure 8, where a spectrum of $l_{1}$ and
$l_{2}$ were observed, changing in both the physical location and duration of
the experiment. As $|\boldsymbol{c}_{g_{1}}|$, defined in (10), is a function
of wavenumber, this spectrum strongly impacts the range of group velocities
present in the beam. The strength of the shading in Figure 14(b) corresponds
to the amplitude of the power spectrum for each wavenumber, obtained from
Fourier transforming each secondary wave beam profile. The analytically
calculated profile assumes a sinusoid with characteristic wavenumber given by
the solid branch of the loci in Figure 3 and width given by the geometry of
the triad (assuming a primary beam width $\Lambda_{0}$), enclosed in a
Gaussian envelope. Various windowing functions were tested and were not found
to significantly alter the range of wavenumbers obtained.

Figure 13: Schematic showing the effect of different $\omega_{1}$ and
$\omega_{2}$ combinations on $\delta_{1}$ (red) and $\delta_{2}$ (green), the
distances over which the secondary wave beams can extract energy with the
primary beam before exiting the boundary.

The result of these two changing factors, $\delta_{q}$ and
$|\boldsymbol{c}_{g_{q}}|$ (where $q=1$ or 2), shown in Figure 14(a) and (b),
respectively, are combined to form a residence time
$R_{q}=\delta_{q}/|\boldsymbol{c}_{g_{q}}|$ that characterises how long each
secondary wave beam spends within $\mathbbm{B}_{0}$. The non-dimensional
residence time, given as a function of frequency, is shown in Figure 14(c). As
$\omega_{1}$ increases, $R_{1}$ also increases. While an increase in
$\omega_{1}$ results in a decrease to $\omega_{2}$ and therefore a shallower
angle for $\mathbbm{B}_{2}$, as $\mathbbm{B}_{2}$ and $\mathbbm{B}_{0}$ are
propagating the same direction, the residence time $R_{2}$ is always greater
than $R_{1}$ and it is never the limiting factor in the interaction. This is
shown by the consistently larger values of $R_{2}$ compared to $R_{1}$.

Figure 14: (a) The non-dimensional interaction length of both
$\mathbbm{B}_{1}$ and $\mathbbm{B}_{2}$ with $\mathbbm{B}_{0}$, again as a
function of non-dimensional frequency. The grey region in the background
highlights the range of experimentally obtained frequencies. (b) The non-
dimensional group velocity of both $\mathbbm{B}_{1}$ (red) and
$\mathbbm{B}_{2}$ (green) as a function of non-dimensional frequency. The red
and green shaded regions corresponds to the range of $|\boldsymbol{c}_{g}|$
possible for a fixed frequency due to a broadband wavenumber spectrum. Details
of how these ranges are calculated are provided in the text. (c) The non-
dimensional residence time $R_{1}/T_{0}$ (red) and $R_{2}/T_{0}$ (green) as a
function of non-dimensional frequency, along with the non-dimensional inverse
linear growth rates $(T_{0}\sigma)^{-1}$, for the simulations shown in Figure
11 and many others. The style of the growth rate marker indicates the triad
configuration tested (given in Figure 3), while smaller hexagons are for other
simulations not shown in Figure 11). The colour of the marker indicates the
behaviour of the simulation, characterised as no growth of secondary waves
(magenta), amplitude modulations to the triadic beams (black) or steady
equilbiurm of all amplitudes in the triad (blue). The red and green shaded
regions associated with the residence time are given from the range of
$|\boldsymbol{c}_{g}|$ shown in (b).

Overlaid on Figure 14(c) are the inverse of the linear growth rates $\sigma$
(given by the form $e^{\sigma t}$), which are marked for all 36 simulations
shown in Figure 11 by the same $\mathbbm{T}_{\Phi}$ marker style as Figure 3,
along with many others for different simulations, marked with a hexagon. The
linear growth rates are obtained by a linear fit on a logarithmic-linear plot
to the initial growth of each simulation. The inverse growth rate can be
viewed as a ‘development time’, that characterises how long the secondary
beams take to grow. The colour of each mark indicates which behavioural
evolution the simulation corresponds to. The magenta is used when no
instability arose, the black markers represent those simulations with
amplitude modulations and the blue marks indicate the simulations that
achieved steady state. The behaviours of the simulations (shown with different
style markers and larger size) can be verified from Figure 11.

Figure 14 suggests why the $\mathcal{M}_{2\textrm{D}}$ model is so sensitive
to the secondary wave parameters and why three different behavioural
evolutions are observed across the triadic configurations tested. For the
cases where no instability occurs, marked in magenta on Figure 14(c), the
development time of the secondary beams is greater than the residence time,
meaning $\mathbbm{B}_{1}$ propagates out of $\mathbbm{B}_{0}$ before
sufficient energy transfer can occur. This is case for all the input
amplitudes shown in Figure 11(a) and (b) using the triad configurations
$\mathbbm{T}_{a}$ and $\mathbbm{T}_{b}$, where no growth of the secondary wave
beams is observed. In this case, $R_{1}\sigma<1$.

For the triad configurations and amplitudes that exhibited the quasi-periodic
modulations in Figure 11, their development times are marked in black. The
majority of these points lie within the range of residence times for
$\mathbbm{B}_{1}$. For these cases, therefore, the development time of the
secondary wave beams is comparable with the time taken for $\mathbbm{B}_{1}$
to propagate across $\mathbbm{B}_{0}$. The secondary beams are able to grow
but have insufficient time to saturate to a stable equilibrium, as wave
perturbations will have moved out of the interaction region before this can
occur. Moreover, due to the range of group velocities present in the beam,
energy will be leaving the primary beam at different times, enhancing the
modulations. For these cases, $R_{1}\sigma\approx 1$. Using this $R_{1}\sigma$
measure allows us to account for the forcing amplitude of the primary beam
(which affects $\sigma$). For example, for the simulations marked with
hexagons at low values of $\omega_{1}$, the forcing amplitude was
significantly increased in order to get the secondary beams to grow.

The final behavioural evolution, observed for the higher amplitude forcing in
Figure 11(e) and (f), is when the secondary beams grow and no modulations are
observed; rather we see a smooth, rapid decay of the primary beam and the
system reaching a steady state. For these cases, where $R_{1}\sigma>1$, the
development times (marked by the blue triangle and cross) are sufficiently
short compared to the range of residence times $R_{1}$. This means
$\mathbbm{B}_{1}$ is able to extract sufficient energy from the primary beam
to reach a steady equilibrium (set by the value of the non-linear interaction
term $I_{p}$) before exiting the underlying $\mathbbm{B}_{0}$. This reflects
how the system would act in the limit of a plane-wave, where the triadic
interactions occur infinitely over space and time and the residence time is
always greater than the development time.

Through only considering discrete triadic configurations in the
$\mathcal{M}_{2\textrm{D}}$ model we have been able to understand the
sensitivity of the instability to the secondary wave parameters and how, in
the context of a finite-width beam, the spatio-temporal configuration of the
triad plays a fundamental role in the evolution of the instability. The
$\mathcal{M}_{2\textrm{D}}$ model shows that when the development time has a
comparable time-scale to the duration of residence of a wave-packet, the
secondary beams are able to grow but are unable to extract sufficient energy
to reach a saturated equilibrium state, resulting is continuous amplitude
modulations. This phenomenon is also seen experimentally, yet here there is a
whole range of perturbations present in the underlying flow. This means that
at different locations in physical space, separate triad configurations will
be selected, based on the varying structure of $\mathbbm{B}_{0}$ across the
domain and a range of background perturbations. As one triadic interaction
decays, another forms, but this time at a different physical location with
modified secondary beam parameters. This explains why, experimentally,
modulations were observed not only in physical space but also in Fourier
space. In the $\mathcal{M}_{2\textrm{D}}$ model simulations presented here,
the interaction region of the triad did not move in physical space, as, by
construction, $\mathbbm{B}_{1}$ and $\mathbbm{B}_{2}$ have fixed frequencies
and wavenumbers, causing growth in a specific location.

##  6 Conclusions

Novel experimental results have shown that, when TRI arises in a finite-width
internal gravity wave beam, the physical regions containing the secondary wave
beams modulate over long time-scales without reaching a steady equilibrium.
Analysis using Dynamic Mode Decomposition and Fourier methods show that these
modulations are present also in the amplitude of the triadic system and in the
Fourier space parameters of the secondary beams.

Through the development and implementation of a two-dimensional weakly non-
linear ($\mathcal{M}_{2\textrm{D}}$) model, we have then been able to dissect
the experimental set-up by analysing how individual triad configurations
effect the evolution of the instability. This model highlights the importance
of considering the instability as a function of space, showing how different
frequencies of the secondary wave beams alter both their orientation and their
group velocity, hence changing their residence time within the underlying
$\mathbbm{B}_{0}$. By comparing the linear growth rate $\sigma$ of the
$\mathcal{M}_{2\textrm{D}}$ model simulations with the residence time of the
secondary beams $R_{q}$, we have identified the conditions under which these
modulations appear to occur. When the these two time-scales are comparable
(i.e. $\sigma R_{q}\sim 1$), the secondary beams are able to grow, yet the
system is unable to reach an equilibrium state as $\mathbbm{B}_{1}$ and
$\mathbbm{B}_{2}$ are unable to extract sufficient energy from
$\mathbbm{B}_{0}$. When either the amplitude of $\mathbbm{B}_{0}$ (increasing
$\sigma$) or the residence time $R_{2}$ is sufficiently increased, the system
is able to extract sufficient energy to reach equilibrium.

Through the $\mathcal{M}_{2\textrm{D}}$ model, we have been able to isolate
the weakly non-linear dynamics of a single triad in a two-dimensional
framework and explain the conditions under which certain triad configurations
result in amplitude modulations. Yet, in the $\mathcal{M}_{2\textrm{D}}$
model, when these modulations in amplitude occurred, the spatial location of
the instability was fixed in one region of the domain. Experimentally, the
reasons we also see modulations in the physical location of the triad and in
the Fourier space parameters of $\mathbbm{B}_{1}$ and $\mathbbm{B}_{2}$ is due
to a whole range of perturbations being present in the underlying flow. A
range of perturbations means that at different locations in physical space,
different triad configurations with similar linear growth rates will be
selected. As one triadic interaction decays another forms, but this time at a
different physical location with modified secondary beam parameters. This is
why the $\mathcal{M}_{2\textrm{D}}$ model did not exhibit movement of the
triad in physical space, since, for the results presented in this paper,
$\mathbbm{B}_{1}$ and $\mathbbm{B}_{2}$ only correspond to discrete
parameters. Further work not reported here (Grayson, 2021), suggests that when
a host of triad configurations are present in the $\mathcal{M}_{2\textrm{D}}$
model, modulations of the secondary beams are also witnessed in both physical
and Fourier space.

While there are many other mechanisms that need to be considered in
oceanographic data, understanding the evolution of freely evolving finite-
width internal wave beams in unbounded domains is fundamental. Here we have
elucidated new developments to our knowledge of how the triadic resonance
instability mechanism may manifest in scenarios more akin to those found in
the ocean as opposed to monochromatic plane waves.

##  7 Acknowledgements

K. M. Grayson acknowledges support from an Natural Environmental Research
Council (NERC) studentship (grant no. NE/L002507/1) and from an Engineering
and Physical Sciences Research Council (EPSRC) fellowship (grant no.
EP/W522600/1).

For the purpose of open access, the authors have applied a Creative Commons
Attribution (CC BY) licence to any Author Accepted Manuscript version arising
from this submission.

Declaration of Interests. The authors report no conflicts of interest.

## References

  * Alford et al. (2007) Alford, M. H., MacKinnon, J. A., Zhao, Z., Pinkel, R., Klymak, J. & Peacock, T. 2007 Internal waves across the pacific. Geophysical Research Letters 34, 2–7.
  * Beckebanze et al. (2021) Beckebanze, F., Grayson, K. M., Maas, L. R.M. & Dalziel, S. B. 2021 Experimental evidence of internal wave attractor signatures hidden in large-amplitude multi-frequency wave fields. Journal of Fluid Mechanics 915.
  * Beckebanze et al. (2019) Beckebanze, F., Raja, K. J. & Maas, L. R. M. 2019 Mean flow generation by three-dimensional nonlinear internal wave beams. Journal of Fluid Mechanics 864, 303–326.
  * Bordes et al. (2012) Bordes, G., Venaille, A., Joubaud, S., Odier, P. & Dauxois, T. 2012 Experimental observation of a strong mean flow induced by internal gravity waves. Physics of Fluids 24.
  * Bourget et al. (2013) Bourget, B., Dauxois, T., Joubaud, S. & Odier, P. 2013 Experimental study of parametric subharmonic instability for internal plane waves. Journal of Fluid Mechanics 723, 1–20.
  * Bourget et al. (2014) Bourget, B., Scolan, H., Dauxois, T., Le Bars, M., Odier, P. & Joubaud, S. 2014 Finite-size effects in parametric subharmonic instability. Journal of Fluid Mechanics 759, 739–750.
  * Brouzet et al. (2016) Brouzet, C., Ermanyuk, E. V., Joubaud, S., Sibgatullin, I. & Dauxois, T. 2016 Energy cascade in internal-wave attractors. Epl 113 (44001).
  * Clark & Sutherland (2010) Clark, H. A. & Sutherland, B. R. 2010 Generation, propagation, and breaking of an internal wave beam. Physics of Fluids 22 (7), 1–16.
  * Dalziel et al. (2007) Dalziel, S. B., Carr, M., Sveen, J. K & Davies, P. A. 2007 Simultaneous synthetic schlieren and piv measurements for internal solitary waves. Measurement Science and Technology 18 (3), 533.
  * Dalziel et al. (2000) Dalziel, S. B., Hughes, G. O. & Sutherland, B. R. 2000 Whole field density measurements by ‘synthetic schlieren’. Experiments in Fluids 28, 322–335.
  * Dalziel et al. (1999) Dalziel, S. B., Linden, P. & Youngs, D. L. 1999 Self-similarity and internal structure of turbulence induced by rayleigh-taylor instability. Journal of Fluid Mechanics 399, 1–48.
  * Dauxois et al. (2018) Dauxois, T., Joubaud, S., Odier, P. & Venaille, A. 2018 Instabilities of internal gravity wave beams. Annual Review of Fluid Mechanics 50, 131–156.
  * Dauxois & Young (1999) Dauxois, T. & Young, W. R. 1999 Near critical reflection of internal waves. Journal of Fluid Mechanics 390, 271–295.
  * Davis & Acrivos (1967) Davis, R. E. & Acrivos, A. 1967 The stability of oscillatory internal waves. Journal of Fluid Mechanics 30 (4), 723–736.
  * Dobra (2018) Dobra, T. E. 2018 Nonlinear interactions of internal gravity waves. PhD thesis, University of Bristol.
  * Dobra et al. (2019) Dobra, T. E., Lawrie, A. & Dalziel, S. B. 2019 The magic carpet: an arbitrary spectrum wave maker for internal waves. Experiments in Fluids .
  * Dobra et al. (2021) Dobra, T. E., Lawrie, A. G.W. & Dalziel, S. B. 2021 Harmonics from a magic carpet. Journal of Fluid Mechanics 911.
  * Dobra et al. (2022) Dobra, T. E., Lawrie, A. G.W. & Dalziel, S. B. 2022 A hierarchical decomposition of internal wave fields. Journal of Fluid Mechanics 934, 1–40.
  * Fan & Akylas (2019) Fan, B. & Akylas, T. R. 2019 Effect of background mean flow on psi of internal wave beams. Journal of Fluid Mechanics 869, R1.
  * Fan & Akylas (2020) Fan, B. & Akylas, T. R. 2020 Finite-amplitude instabilities of thin internal wave beams: Experiments and theory. Journal of Fluid Mechanics 904.
  * Gostiaux et al. (2007) Gostiaux, L., Didelle, H., Mercier, S. & Dauxois, T. 2007 A novel internal waves generator. Experiments in Fluids 42, 123–130.
  * Grayson (2021) Grayson, K. M. 2021 Triadic resonance instability in finite-width internal gravity wave beams. PhD thesis, University of Cambridge.
  * Grisouarda et al. (2013) Grisouarda, N., Leclair, M., Gostiaux, L. & Staquet, C. 2013 Large scale energy transfer from an internal gravity wave reflecting on a simple slope. Procedia IUTAM 8, 119–128.
  * Joubaud et al. (2012) Joubaud, S., Munroe, J., Odier, P. & Dauxois, T. 2012 Experimental parametric subharmonic instability in stratified fluids. Physics of Fluids 24.
  * Karimi & Akylas (2014) Karimi, H. H. & Akylas, T. R. 2014 Parametric subharmonic instability of internal waves: Locally confined beams versus monochromatic wavetrains. Journal of Fluid Mechanics 757, 381–402.
  * Kataoka & Akylas (2015) Kataoka, T. & Akylas, T. R. 2015 On three-dimensional internal gravity wave beams and induced large-scale mean flows. Journal of Fluid Mechanics 769, 621–634.
  * Koudella & Staquet (2006) Koudella, C. R. & Staquet, C. 2006 Instability mechanisms of a two-dimensional progressive internal gravity wave. Journal of Fluid Mechanics 548, 165–196.
  * Lamb (2004) Lamb, K. G. 2004 Nonlinear interaction among internal wave beams generated by tidal flow over supercritical topography. Geophysical Research Letters 31 (9), 1–4.
  * MacKinnon & Winters (2005) MacKinnon, J. A. & Winters, K. B. 2005 Subtropical catastrophe: Significant loss of low-mode tidal energy at 28.9 degrees. Geophysical Research Letters 32, 1–5.
  * McEwan (1971) McEwan, A. D. 1971 Degeneration of resonantly-excited standing internal gravity waves. Journal of Fluid Mechanics 50, 431–448.
  * McEwan & Plumb (1977) McEwan, A. D. & Plumb, R. A. 1977 Off-resonant amplification of finite internal wave packets. Dynamics of Atmospheres and Oceans 2, 83–105.
  * Mercier et al. (2008) Mercier, M. J., Garnier, N. B. & Dauxois, T. 2008 Reflection and diffraction of internal waves analyzed with the hilbert transform. Physics of Fluids 20.
  * Munk (1966) Munk, W. 1966 Abyssal recipes. Deep-Sea Research and Oceanographic Abstracts 13, 707–730.
  * Nash et al. (2004) Nash, J. D., Kunze, E., Toole, J. M. & Schmitt, R. W. 2004 Internal tide reflection and turbulent mixing on the continental slope. Journal of Physical Oceanography 34 (5), 1117–1134.
  * Oster & Yamamoto (1963) Oster, G. & Yamamoto, M. 1963 Density gradient techniques. Chemical Reviews 63, 257–268.
  * Peacock et al. (2009) Peacock, T., Mercier, M., Didelle, H., Viboud, S. & Dauxois, T. 2009 A laboratory study of low-mode internal tide scattering by finite-amplitude topography. Physics of Fluids 21, 1–4.
  * Richet et al. (2018) Richet, O., Chomaz, J. M. & Muller, C. 2018 Internal tide dissipation at topography: Triadic resonant instability equatorward and evanescent waves poleward of the critical latitude. Journal of Geophysical Research: Oceans 123 (9), 6136–6155.
  * Roe (1986) Roe, P. L. 1986 Characteristic-based schemes for the euler equations. Annual Review of Fluid Mechanics 1, 337–365.
  * Sarkar & Scotti (2017) Sarkar, S. & Scotti, A. 2017 From topographic internal gravity waves to turbulence. Annual Review of Fluid Mechanics 49, 195–220.
  * Schmid (2010) Schmid, P. J. 2010 Dynamic mode decomposition of numerical and experimental data. Journal of Fluid Mechanics 656, 5–28.
  * Staquet & Sommeria (2002) Staquet, C. & Sommeria, J. 2002 Internal gravity waves: From instabilities to turbulence. Annual Review of Fluid Mechanics 34, 559–593.
  * Sutherland (2006) Sutherland, B. R. 2006 Internal wave instability: Wave-wave versus wave-induced mean flow interactions. Physics of Fluids 18.
  * Sutherland (2010) Sutherland, B. R. 2010 Internal Gravity Waves. Cambridge University Press.
  * Sutherland (2013) Sutherland, B. R. 2013 The wave instability pathway to turbulence. Focus on Fluids 724, 1–4.
  * Thomas et al. (2009) Thomas, L. P., Marino, B. M. & Dalziel, S. B. 2009 Synthetic schlieren: Determination of the density gradient generated by internal waves propagating in a stratified fluid. Journal of Physics: Conference Series 166.
  * Thorpe (1968) Thorpe, S. A. 1968 On the shape of progressive internal waves. Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences 263, 563–614.
  * Thorpe & Haines (1986) Thorpe, S. A. & Haines, A. P. 1986 On the reflection of a train of finite-amplitude internal waves from a uniform slope. Journal of Fluid Mechanics 178, 279–302.
  * Wolfram Research (2021) Wolfram Research, Inc. 2021 Mathematica, Version 12.3. Wolfram Research, Inc.
  * Wunsch & Ferrari (2004) Wunsch, C. & Ferrari, R. 2004 Vetical mixing, energy, and the general circulation of the oceans. Annual Review of Fluid Mechanics 36, 281–314.