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

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.


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 (1.1) where the wave phase, , is defined as (1. 2) The subscript = (0, 1, 2) is used throughout this paper to define the primary wave and the two secondary waves, respectively. Both (1.1) and (1.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 is defined as = ( , ) with magnitude | | = , where the components are given in Cartesian co-ordinates ( , ), marked in Figure 1, and denotes the frequency of the wave. In order to distinguish between the three wave beams in the triad and their corresponding parameters, we define where B indicates a wave beam with density and stream function Ψ fields, frequency , characteristic wavenumber vector and beam width Λ . For the primary beam 0 , 0 , and Λ 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 where is the angle between the lines of constant phase and the vertical and and are the characteristic wavenumber contributions from each beam. Here is the buoyancy frequency of the stratification given by where is the gravitational constant of acceleration. Under the assumptions of a Boussinesq, incompressible fluid, we decompose the total density as = 0 +¯( ) + ( , , ), where 0 is the reference density,¯is the background density stratification as a function of depth and is the perturbation density. We consider the density changes from perturbations and stratification to be small compared to the reference density, so that¯, 0 . Given the triadic resonant condition in (1.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 § 2 we detail the experimental set-up and processing procedure. In § 3 we then present the experimental results, looking first at the initial observations in § 3.1 and then at the long-term evolution of the experiments § 3.2. Based on these observations, we present the development of the two-dimensional weakly non-linear model in § 4. Here we outline the perturbation expansion used in § 4.1 and the subsequent development of the numerical solution in § 4.2 (comprising the two-dimensional advection scheme and the weakly non-linear interactions). In § 5 we present the results of the model. We start with § 5.1, where we examine the weakly non-linear interactions on their own before moving onto § 5.2, where the results of the weakly non-linear two-dimensional model are given. Conclusions are then drawn in § 6.

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(Dobra et al. , 2021Beckebanze 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 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.
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 saltstratification 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 = 0.45 ± 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 = 1.54 ± 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). 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.

Wave visualisation
A random dot pattern attached to an LED light bank was located 0.20 ± 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 ± 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 , which for convenience we non-dimensionalise according to where 0 is the horizontal wavelength of the primary wave beam given as 0 = 2 /| 0 |, where 0 is the horizontal component of the primary wave vector 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.

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.
where the locations , , , are respectively 7 /| 0 |, 9 /| 0 |, 13 /| 0 |, 15 /| 0 | and 0 is set to −0.05 mm −1 , giving a horizontal wavelength of 0 = 2 /| 0 | = 125.66 mm. As we restrict > 0 (for all ), having 0 < 0 means the primary wave beam is propagating to the left. The spatial structure of the forcing, described by (2.2), 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, − , 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 where = cos 2 ( /8 2 ) is the smoothing function on the outer flanks of the beam profile.
The long view of TRI in finite-width internal wave beams 7 The temporal forcing ( ) is then described as where 0 , 0 and 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 0 = 2 / 0 . The only two parameters to be varied in this experimental study are end and 0 . The run time, end , is either 90 or 180 minutes, while the non-dimensional amplitude 0 / 0 ranges between 0.028 -0.036. The amplitude threshold for the instability is achieved at 0 / 0 ≈ 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 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.

Initial observation and analysis
We start by examining one experiment from the set of 36 with an imposed amplitude displacement of 0 / 0 = 0.032. Figure 2(a) shows over the visualisation region at / 0 = 83. Here, B 0 , generated by the wavemaker, propagates energy up and to the left, at its respective group velocity 0 . The group velocity is defined for all wave beams by where again the subscript = (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 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 = 0 /( 0 ), where is the kinematic viscosity of 1 mm 2 s −1 . Here, ≈ 170. As the selected input amplitude displacement of 0 / 0 = 0.032 is above the instability threshold, B 0 becomes unstable, leading to the formation of two secondary beams. One of these beams, B 1 , is clearly visible in Figure 2(a). This beam emanates from the central region of 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, B 2 , that completes the triad is not visible. In order to understand the underlying modal structure of these beams, the flow field , is decomposed using Dynamic Mode Decomposition (DMD). DMD works by preforming an eigen-decomposition of a linearised representation of the shows the spatial averaging domain used for the primary beam, discussed in § 4.2.
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 . The maximum number of dynamic modes is given by the number of input frames in the sequence (in this case = 20 s, as we use a frame rate of 1 fps, which is just greater than the slowest period of the triad 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 / 0 86 are given in Figure 2(b)-(d).
As expected from the input forcing, Figure 2(b) corresponds to the input B 0 with nondimensional frequency 0 / = 0.62. Figure 2(c) then corresponds to B 1 with 1 / = 0.23 and (d) to the obscured B 2 , which propagates with a group velocity up and to the left, with non-dimensional frequency 2 / = 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, 0 = 1 + 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. The fourth (non-decaying mode) corresponds to / = 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 In this case, 0 = 0 + 0 − 0 , corresponds to the phase angle of B 0 and 0 = 0 − 0 − 0 to its reflection, B 0 , from the free surface. These wave beams will sum to produce a third wave beam with wavenumber vectorˇ= (0, 2 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 threedimensional 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 zerofrequency 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 B 0 and B 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 / 0 86 -the real parts of which are shown in Figures 2(b)-(d) -and calculating a twodimensional 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 / 0 = 0.005, given by the ratio of pixel resolution to region size. As the analysis preformed on the horizontal density gradient provides similar results to that of the vertical , 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 0 0 = (−6.28, −8.29), 0 1 = (3.73, 14.60), and 0 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 1 given 0 , based on both the dispersion relationship (1.4) and the TRI condition (1.1).
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 0 = 1 + 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 B 1 and B 2 are uniform over the field of view. In fact, as the experiment progresses, significant modulations to the structures of B 1 and 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 finiteamplitude 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.  intensity and width. This migratory behaviour persists for the full duration of the experiment, which lasts for over 800 periods of the primary beam.

Long-time development
Further quantitative analysis of this peculiar behaviour requires us to calculate the amplitude of the individual resonant wave beams B 0 , B 1 and B 2 . Decomposing by frequency into complex constituent fields using DMD, we find the inverse gradient (potential) field / 0 by integrating both the real and imaginary components of . 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 B 0 from its reflection from the free surface, B 0 . In order to have a singular value for amplitude that is independent of space, we then spatially mean average / 0 over the whole field of view denoted . 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 § 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 . 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 § 4.2. We then infer the amplitude from the measured displacements of the scalar twodimensional stream function Ψ = Ψ( , ) for each field, where the velocity vector = ∇ × (Ψˆ) and ∇ = ( / , / ). Assuming an oscillatory form Ψ =Ψ for each triadic wave beam, we define our spatially averaged amplitude for each constituent field as Ψ , 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 Ψ = Ψ ( ) , 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 =˜for each field.
Substituting the above forms for stream function and density into the inviscid linear conservation of mass equation / = − (¯/ ), and cancelling the fast time-scales, we find the reduced stream function amplitude for each wave field using where and are the frequency and horizontal wavenumber of the given beam , and |˜| 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 (3.3) for two experiments. In (a) we show the same experiment as Figure 4, while (b) corresponds to another experiment with the same amplitude ( 0 / 0 = 0.032) but with a much longer run time ( end / 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, 0 / 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 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 (3.3), as opposed to the imposed displacement from the wavemaker.
Another observable feature in Figure 5 is the gradual increase of the mean amplitude of 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Ψ using (3.3) -we were able to discount all extraneous sources that could contribute to these structural modulations.
In Figure 5, the amplitudes of B 1 (red) and B 2 (green) are positively correlated; their amplitudes are almost scaled values of each other. Meanwhile, the amplitude of B 0 is negatively correlated with B 1 and B 2 . When B 0 is at a local maximum, the amplitudes of B 1 and B 2 are concurrently at a local minimum and then versa when the amplitude of B 0 is at a minimum. This coupling of the modulations in amplitude between B 0 and the secondary B 1 and 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 B 1 /B 0 and B 2 /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 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 rate of 1 fps, the highest resolvable frequency (shortest time period) is / = 4.08. Several windowing functions were tested, and were not found to significantly affect the spectrogram results. The angled brackets, , again indicate that the results are averaged across the whole visualisation region. This underlying spectrogram, calculated using (3.4), therefore, reveals the details about the distribution of the frequency spectra for 1 and 2 . In contrast, as we are only selecting the three dominant modes obtained from the DMD over short time intervals ( / 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.
Both spectrograms in Figures 6(a) and (b) show a clear peak at 0 / = 0.62 for all time, consistent with the imposed displacement from ASWaM. Both secondary beams emerging from the instability become visible at approximately / 0 = 50, with peaks in the spectra around 1 / ≈ 0.23 and 2 / ≈ 0.39, though subsequently these modulate on a slow timescale 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 0 = 1 + 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) andBrouzet 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 / ≈ 0, which was also observed from the DMD and has already been discussed. The other two frequencies / ≈ 0.84 and / ≈ 1 correspond to two different TWIs, given in (3.2), between B 0 and either B 1 or 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 B 1 and B 2 associated with the instability. This is shown by the convergent 'wisps' on the 1 / and 2 / bands, where additional secondary beam pairs appear and merge with the continuous mode. This is highlighted for / 0 = 452 by the inset in Figure 6(b). Here we have plotted ln( ( )/ ) in cyan and ln( ( 0 − )/ ) in magenta against / 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 1 + 2 = 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 / 0 483, from the experiment presented in (3.6) where the widths of the Hamming windows are given by / 0 = 39 and = / 0 = 9.7, and the subscript on the angle brackets shows that ‡ ( , , ) 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, ‡ ( , , ), corresponding to 2 , around mid-height in the domain, there is a shift in both the amplitude and value of 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 / 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 1 and 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.
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 B 1 or 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 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 B 0 varies across the height of the domain. As 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 Λ 0 of the primary beam 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 M 2D . Details of its development are outlined in the following section. The goal here is to  Figure 4, showing the distribution of 1 and 2 over the non-dimensional height in the domain ( / ) calculated by (3.5). Here the image sequence is first temporally filtered in Fourier space to remove the signal from B 0 , and thus we do not see a peak at 0 . Both the surface plot colour and the height of the peaks, show the power spectral density ‡ ( , ). The background plot then shows 1 and 2 at the same instant in time in the Fourier plane of horizontal wavenumber component and of frequency. This contour plot is defined by (3.6). The timings of each image are given by the white dashed lines in Figure 6.
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 § 4.1 we examine the perturbation expansion used and in § 4.2 we look at the numerical solution to the obtained system of equations.

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 , the full Boussinesq non-linear equations of motion (comprised of the momentum, continuity and mass conser-vation equations) are given by where is the dynamic pressure and is the dynamic viscosity. Writing velocity in terms of the scalar stream function = ∇ × (Ψˆ), this system can be reduced to (4.5) Here, we focus on (4.4), 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 Ψ =Ψ( , , ) ( · − ) + c.c., (4.6) =˜( , , ) ( · − ) + c.c., (4.7) whereΨ and˜are the reduced forms of the stream function and density respectively and c.c. represents the complex conjugate. BothΨ and˜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 havẽ = − 2 0Ψ , (4.8) 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 = (|Ψ 00 | 2 0 ) −1 , (4.9) where, |Ψ 00 | is the characteristic magnitude of the stream function associated with the primary beam B 0 and = / is the kinematic viscosity. Here, can be viewed as a non-dimensional measure of the beam amplitude, which characterises the relative importance of the non-linear · ∇ terms in the momentum and conservation of mass equations. The spatial parameter defines the separation between the dominant wavelength and overall width of the primary beam. Finally, is the inverse of the Reynolds number defined earlier. In oceanographic settings, will be orders of magnitude smaller than and 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 ∼ 10 −1 1, broad beam ∼ 10 −1 1 and low viscosity ∼ 10 −2 1, and introduce re-scaled time and position through the five variables = , = , = , = , = . (4.12) As we shall see, accounts for the 'slow non-linear time' variations to the amplitude, while governs the 'slow advection time' scale. As, for the experiments, both 1 and 1, these scaled times account for change over long-time periods. The final time-scale given by represents the viscous decay of the wave beams. The spatial parameters and account for the gradual spatial variability of the wave beams in the domain. Utilising the dimensionless amplitude , we re-write the stream function in (4.6) as Ψ =Ψ( , , ) ( · − ) + c.c. =Ψ( , , , , ) ( · − ) + c.c., (4.13) whereΨ =Ψ. We therefore require |Ψ| ∼ 1 asΨ =Ψ = (|Ψ 00 | 2 0 / )Ψ 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 ( · − ) = . 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 ( , , , , , ) + c.c., (4.14) where the subscript indicates a locally plane-wave approximation to B 0 , B 1 or B 2 .
Each wave phase therefore represents the characteristic frequency and wavenumber contribution to B . A similar set of expressions can be written for as a sum of˘.
The superposition in (4.14) is then substituted into the non-linear equation (4.4) where, due to the separate space and time-scales, the partial derivatives in (4.4) with respect to ( , , ) become (4.15b) 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 the right-hand side (RHS) of (4.4) vanishes. Therefore, as the contributions from density on the RHS first appear at O ( 2 ), the linear relationship in (4.8) is valid up until order O ( 2 ). The resultant expression obtained provides the Boussinesq viscous equations of motion solely as a function of Ψ that can be used to examine both linear and non-linear dynamics between a triadic set of waves simply by collecting around orders of , and .
At O ( 0 0 ) and O ( 0 1 ) the expression will be zero as these orders correspond to a state of rest. At order O ( 1 0 ), we recover the linear wave solution, with the non-linearity in (4.4) vanishing and linear superposition applying such that the three waves propagate independently. Extracting terms with a common factor of , leaves which is the linear dispersion relationship for internal plane-waves. Indeed, for a single nonlinear plane-wave in the inviscid limit, the RHS of (4.4) 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 O ( 1 1 ). Again, as is still at first-order, the non-linear terms in (4.4) cancel. Looking at the terms, we obtain which after some re-arranging can be expressed as (4.18) where ∇ = ( / , / ) and ∇ = ( / , / ). 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 is the velocity at which energy is transported (e.g. Sutherland 2010). As energy scales with ∼ Ψ 2 , the fact that (4.18) shows thatΨ is also advected by , is not altogether surprising. Before examining the non-linear interaction terms at O ( 2 ), the viscous term at O ( ) needs to be considered in conjugation with the linear advection recovered at O ( 1 1 ). While these terms are of lower magnitude than the terms at O ( 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 O ( ), looking at terms with a common factor , we obtainΨ (4.19) which shows that the viscous decay of each beam scales as 2 . We choose to combine the evolution on time-scales and to obtain, at first order in , the advection equatioñ (4.20) This advection equation in (4.20) is solved in the 2D advection component of the M 2D model and its numerical implementation is addressed in § 4.2. We next consider terms at order O ( 2 0 ). Here, non-linearity enters the problem and terms are no longer only associated with , but are also comprised of cross-terms from · ∇ operator in (4.4), which have the form ( + ) , when expressed in Fourier modes (where , , are permutations of 0, 1, 2). We consider a triad of wave beams satisfying resonance conditions (1.2). Using (4.8) to eliminate˘from (4.4), at O ( 2 ) we recover the same non-linear interactions given in Bourget et al. (2013), specifically where an asterisk indicates the complex conjugate and the interaction term is given as It is possible to extend this expansion to examine the higher order O ( 2 1 ), however as this model sufficiently captures enough of the observed experimental behaviour, this further expansion was not required. Above O ( 2 1 ) (higher orders in ) the linear approximation in (4.8) is no longer valid for eliminating˘from (4.4). The evolution of these coupled ODEs in (4.21), are considered on their own and in § 5.1 and as part of the M 2D model in § 5.2.

Numerical implementation
This section outlines the development of the M 2D numerical model, built to solve the equations obtained at order O ( 1 1 ) and at order O ( 2 ). We start with the numerical scheme used to solve the advection equation (4.20), obtained at order O ( 1 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 in (4.20) is specific to each wave beam, , the number of numerical domains corresponds to the number of wave beams being considered, as, up to O ( 2 ), we only need to considerΨ and need not explicitly consider˜. 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 byB , as we are advecting the reduced stream functionΨ in (4.20). Being a second order scheme, the volume flux, , is calculated using the two upstream values of the reduced stream function. AsΨ 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 domainB 0 . To avoid confusion with indexing, compass coordinates are used for the nodes and the subscript = 0 is dropped fromΨ 0 in the equations below for clarity. Specifically, for domainB 0 , = 0Ψ and = 0Ψ , wherẽ the node indexing is shown in Figure 9 and Υ is the flux limiter which prevents the generation Figure 9: Sketch outlining the numerical advection domainB 0 . The bottom boundary forcing of the wavemaker is given by the sinusoid envelope in (2.2). The two outer grid layers required for the second-order scheme are shown in blue. In this domain 0 is directed to the north-west, meaning the advection of the complexΨ 0 is to the north-west.
For advection,Ψ 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.
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Ψ is then updated by the differences in the fluxes over the control volume using , (4.24) where the superscripts and † represent the value ofΨ before and after advection respectively and Δ is the time-step over which code is advanced. For the scheme to remain numerically stable, we limit the time step, Δ , to satisfy the CFL condition. Each domainB 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 Δ / 0 = Δ / 0 = 0.04, is much finer than the smallest wavelength considered. At every time step, random complex background noise of magnitude 10 −7 (2 ) is added to each cell in every domain. Both and are independent random numbers spanning the range [0, 1], giving a mean magnitude of the order 5 × 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 isB 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Ψ, the boundary condition is given by the envelope of the experimental forcing, which is obtained by removing the fast time forcing, 0 , from (2.2) and (2.4). This boundary condition forB 0 is shown in Figure 9. We implement a prescribed displacement condition by computing the corresponding complex valued stream function amplitudeΨ in , as opposed to 0 , whereΨ in = |Ψ in | (2 ) . Here, is a constant value within the range [0, 1]. This ensures that the phase angle of the slowly evolving complex amplitudeΨ 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 (4.21) are calculated.
The non-linear interactions are included by converting (4.21) into the numerical format whereΨ corresponds to every cell in domainB and † and + 1 respectively represent the values ofΨ in each domain after advection (4.24) and at time + Δ 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 § 5.1 we first consider the non-linear interactions on their own, before examining the results from the M 2D model in § 5.2.

Zero-dimensional model results
We note that the coupled non-linear equations in(4.21) can be recovered in zero-dimensional space by considering the form of the stream functionΨ ( ) in (4.4), 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 (4.21) 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 and length 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 (4.21) becomẽ where is given in (4.22), 0 is a unit vector in the direction of 0 and the forcing term , represents the energy flux for the primary wave beam through the control area with an incoming amplitude ofΨ in . We convert this input amplitude to the non-dimensional measure in = 2 0 |Ψ in |/ . We recognise that in is comparable to the experimental input amplitude = 2 0 |Ψ 0 | / , where |Ψ 0 | is magnitude of the reduced stream function averaged over the black domain in Figure 2(b). The terms on the end of (5.1b) and (5.1c) 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 (5.1), which we use to describe the energy exchange in TRI in the context of a finite-width beam, as the zero-dimensional model M 0D .
The M 0D model is numerically integrated to examine its prediction for the development of the triad. To match the experimental set-up, the width and length of the interaction region are set to Λ 0 (defined in (2.3)) and 2Λ  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 M 0D . This is achieved by defining a resonant triad as T Φ = {B 0 , B 1 Φ , B 2 Φ }, where the subscript Φ 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.1) and dispersion (1.4), meaning all of the configurations lie exactly on the green curve in Figure 3. For the M 0D model we consider configurations T and T , marked by the blue circle and star on Figure 3, respectively, with parameters shown in Table 1.
The results of the numerical integration for the M 0D model are given in Figure 10 for T and T , respectively, across a range of five non-dimensional input forcing amplitudes for the primary beam 0.092 in 0.184. The resultant amplitudes of the triadic beams are given in terms of = 2 0 |Ψ |/ . 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 in order to de-correlate the signal with the position of a beam in space. As the M 0D model is not a function of space, no spatial averaging is required for . While and are therefore not quantitatively comparable, they are qualitatively. For details on the different measures of amplitude, see Table 2.
We first consider the general behaviour of the M 0D model, consistent between both Figure  10 Table 1. growth of |Ψ in | and then remains at a constant amplitude. There then exists a critical forcing amplitude, above which, TRI is observed. For T (Figure 10(a)) this occurs when in 0.111, while for T (Figure 10(b)), instability occurs when in 0.097. While both of these limits are lower than the experimental limit of 0.166, the M 0D 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 |Ψ 0 | always decays to the same asymptotic value, while the amplitudes of |Ψ 1 | and |Ψ 2 | asymptote at higher values as in is increased. Upon investigation, the asymptotic amplitude of |Ψ 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 T and T 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 T , a result that agrees with T 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 3 , this lower value of the secondary beams might be due to greater viscous dissipation in T , due to > .
Comparing Figure 10 with the experimental results in Figure 5, we see remarkably different behaviour in the evolution of the instability. For the M 0D 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 M 0D 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 . 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 M 0D 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 M 2D model to see how it describes the instability when considered in a two-dimensional domain.

Two-dimensional model results
We now consider the results from the M 2D model, where in each simulation we provide three domains -one for each triadic beam -with input parameters chosen to satisfy the triadic condition (1.2) and dispersion (1.4). Six triad configurations T Φ = {B 0 , B 1 Φ , B 2 Φ } are considered, with wavenumber vectors distributed across the solid green loci branch in Figure 3 and parameters given in Table 1. The parameters for B 0 (common to all six triads), match the experiments and are given in (5.2). Figure 11 shows the results of 36 simulations, where sub-plots (a) to (f) correspond to triad configurations T to T respectively, each shown with six input amplitudes, in . Again, in is comparable to the experimental values . 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 T and T (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 T ) 1 / 0 1 and both correspond to the smallest values of 1 considered. Indeed, for input amplitudes that caused the other triadic configurations to become unstable, we found that no pairs with 1 / 0 < 1 generated TRI for the M 2D 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 1 / 0 > 1. When 1 / 0 < 1, we have the condition 1 < 0 < 2 . Thus energy transfers to both the larger length scale ( 1 ) and the smaller ( 2 ). Bourget et al. (2014) show that the linear growth rate from the M 0D model is larger for this triad configuration when Λ 0 < 7 0 . However, in our experiments where Λ 0 ≈ 3 0 , it is still a triad configuration with 1 / 0 > 1 (and therefore 0 < 1 , 2 , located on the outer solid green branch of Figure 3), that is selected. The effect of having a smaller value of 1 is discussed further below. 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 T to T . 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 M 2D and the experimental values. For the triad configurations T to T , instability occurred when in 0.161, while experimentally 0.166 triggered instability. We also notice that for all of these simulations that become unstable, as we increase in , not only do the secondary wave beam amplitudes increase, their growth also occurs at earlier times. This observation is consistent with the M 0D model. We can measure this initial linear growth using the growth rate , of the form , 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 T , forced at a non-dimensional input amplitude of 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 .
Here we see the initial development of B 1 occurring at the top of B 0 , where it grows in strength. (We note that calculations preformed in a larger domain demonstrated that the generation region of 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 0 and 2 , B 2 is not obvious in this visualisation region as it propagates within the confines of B 0 . Over time, 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 B 0 . By construction, 1 and 2 also remain fixed. The final behavioural evolution seen in Figure 11, most obvious for the largest forcing amplitude using T , is when the triadic system reaches a stable equilibrium, closely resembling the steady state results from M 0D 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 (1.4), a greater value of Figure 13: Schematic showing the effect of different 1 and 2 combinations on 1 (red) and 2 (green), the distances over which the secondary wave beams can extract energy with the primary beam before exiting the boundary. spatial region over which B 1 can extract energy from B 0 and consequently an increased distance for B 1 to grow. This is shown schematically in Figure 13. We define the lengths of these interaction regions for B 1 and B 2 as 1 and 2 , respectively. Figure 14(a) shows the relationship between 1 and 1 , together with the relationship between 2 and 2 . The grey shading in the background marks the range of frequencies obtained experimentally. As 1 increases, not only does 1 increase, | 1 | also decreases, as shown in Figure 14(b). A decrease in | 1 | causes B 1 to remain within the spatial confines of 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 | 1 | and | 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 1 at each 1 that can satisfy both dispersion (1.4) and the triadic resonant condition (1.2), (shown by the loci on Figure 3), both B 1 and 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 1 and 2 were observed, changing in both the physical location and duration of the experiment. As | 1 |, defined in (3.1), 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 Λ 0 ), enclosed in a Gaussian envelope. Various windowing functions were tested and were not found to significantly alter the range of wavenumbers obtained.
The result of these two changing factors, and | | (where = 1 or 2), shown in Figure 14(a) and (b), respectively, are combined to form a residence time = /| | that characterises how long each secondary wave beam spends within B 0 . The non-dimensional residence time, given as a function of frequency, is shown in Figure 14(c). As 1 increases, 1 also increases. While an increase in 1 results in a decrease to 2 and therefore a shallower angle for B 2 , as B 2 and B 0 are propagating the same direction, the residence time 2 is always greater than 1 and it is never the limiting factor in the interaction. This is shown by the consistently larger values of 2 compared to 1 . Overlaid on Figure 14(c) are the inverse of the linear growth rates (given by the form ), which are marked for all 36 simulations shown in Figure 11 by the same T Φ marker style as Figure 3, along with many others for different simulations, marked with a hexagon. of both B 1 (red) and B 2 (green) as a function of non-dimensional frequency. The red and green shaded regions corresponds to the range of | | 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 1 / 0 (red) and 2 / 0 (green) as a function of non-dimensional frequency, along with the non-dimensional inverse linear growth rates ( 0 ) −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 | | shown in (b).
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 M 2D 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 B 1 propagates out of 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 T and T , where no growth of the secondary wave beams is observed. In this case, 1 < 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 B 1 . For these cases, therefore, the development time of the secondary wave beams is comparable with the time taken for B 1 to propagate across 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, 1 ≈ 1. Using this 1 measure allows us to account for the forcing amplitude of the primary beam (which affects ). For example, for the simulations marked with hexagons at low values of 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 1 > 1, the development times (marked by the blue triangle and cross) are sufficiently short compared to the range of residence times 1 . This means 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 ) before exiting the underlying 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 M 2D 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 M 2D 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 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 M 2D model simulations presented here, the interaction region of the triad did not move in physical space, as, by construction, B 1 and B 2 have fixed frequencies and wavenumbers, causing growth in a specific location.

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 (M 2D ) 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 B 0 . By comparing the linear growth rate of the M 2D model simulations with the residence time of the secondary beams , we have identified the conditions under which these modulations appear to occur. When the these two time-scales are comparable (i.e. ∼ 1), the secondary beams are able to grow, yet the system is unable to reach an equilibrium state as B 1 and B 2 are unable to extract sufficient energy from B 0 . When either the amplitude of B 0 (increasing ) or the residence time 2 is sufficiently increased, the system is able to extract sufficient energy to reach equilibrium.
Through the M 2D 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 M 2D 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 B 1 and 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 M 2D model did not exhibit movement of the triad in physical space, since, for the results presented in this paper, B 1 and 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 M 2D 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.