On the interaction of short linear internal waves with internal solitary waves

We study the interaction of small-scale internal wave packets with a large-scale internal solitary wave using high-resolution direct numerical simulations in two dimensions. A key finding is that for wave packets whose constituent waves are short in comparison to the solitary wave width, the interaction leads to an almost complete destruction of the short waves. For mode-1 short waves in the packet, as the wavelength increases, a cutoff is reached, and for larger wavelengths the waves in the packet are able to maintain their structure after the interaction. This cutoff corresponds to the wavelength at which the phase speed of the short waves upstream of the solitary wave exceeds the maximum current induced by the solitary wave. For mode-2 waves in the packet, however, no corresponding cutoff is found. Analysis based on linear theory suggests that the destruction of short waves occurs primarily due to the velocity shear induced by the solitary wave, which alters the vertical structure of the waves so that significant wave activity is found only above (below) the deformed pycnocline for overtaking (head-on) collisions. The deformation of vertical structure is more significant for waves with a smaller wavelength. Consequently, it is more difficult for these waves to adjust to the new solitary-waveinduced background environment. These results suggest that through the interaction with relatively smaller length scale waves, internal solitary waves can provide a means to decrease the power observed in the short-wave band in the coastal ocean.

between near-inertial waves and ocean fronts, and found that inertial waves could travel on two distinct characteristics at a front, one flat and one tilted upward, implying the existence of critical reflections from the ocean surface.
Interaction between internal waves of different length scales also occurs naturally (Sun and Pinkel, 2012). When the disparity in length scales between the participating waves is large, the relatively smaller length scale wave essentially plays the role of the disturbance to the "background flow" induced by the relatively larger length scale wave, as they interact with each 5 other. Previous literature has considered wave-wave interaction in a variety of contexts. For example, using ray theory for linearized waves and the principle of wave action conservation, Broutman and Young (1986) studied the interaction of short high-frequency progressive internal waves and long progressive near-inertial waves, and found that there was a net energy transfer from the inertial wave field to the short internal waves.  investigated the interaction between two mode-1 internal solitary waves (ISWs), and showed that interaction of solitary waves did not correspond to soliton dynamics, since 10 energy exchange was observed and small-amplitude trailing waves of possibly higher modes were generated. More recently, Stastna et al. (2015) examined the interaction between mode-1 and mode-2 internal solitary (solitary-like in cases when the mode-2 wave was breaking) waves, and demonstrated that the interaction yielded a nearly complete disintegration of the relatively smaller mode-2 wave. In particular, the majority of kinetic energy carried by the mode-2 wave was lost and the disturbance to the flow field after the collision no longer had a mode-2 structure. When the length scales of participating waves 15 are similar, Sutherland (2016) found that nonlinear self-interaction might occur, which resulted in energy being transferred to superharmonic disturbances. These disturbances were a superposition of modes such that the amplitude was largest where the change in background buoyancy frequency with depth was largest.
In this work, we study the interaction of small-scale mode-1 internal waves initialized from linear waves with an ISW initialized from the exact Dubreil-Jacotin-Long equation, using high-resolution direct numerical simulations in two dimensions. 20 Internal waves that are short in terms of wavelength compared to the fluid depth are generally less documented in the nonlinear wave literature. In fact, the derivation of the model equations of most weakly nonlinear theories, such as the Korteweg-de Vries equation and its variations, assumes large horizontal scales and thus filters out short waves (Lamb and Yan, 1996). Nevertheless, such waves occupy a non-negligible portion of the Garrett-Munk spectrum of internal waves in the oceans (Thorpe, 2005), and it is important to understand their behavior in order to fully describe internal wave dynamics. 25 The remainder of the paper is organized as follows. Theoretical descriptions of internal waves are introduced in section 2.
The problem is formulated in section 3. The simulation results are presented in section 4. A key finding is that for waves that are short in comparison to the ISW width, the interaction leads to an almost complete destruction of the short waves, but that for mode-1 short waves there is a cutoff determined by the wavelength of short waves, and waves longer than this cutoff maintain their structure after interaction. We show that this is a key difference from mode-1-mode-2 interaction, which is examined 30 in Stastna et al. (2015). The energy transfer during the interaction is discussed in section 5, and a summary concluding the findings of this study is given in section 6.

Internal wave theories
In the classical linear theory, the horizontal structure of internal waves is usually described by the travelling wave ansatz where k is the horizontal wave number and is related to the wavelength λ by the formula k = 2π/λ, and c p is the phase speed. The vertical structure is described by solutions of the eigenvalue problem often referred to as the Taylor-Goldstein (TG) equation (Kundu et al., 2012), which is given by where U = U (z) is the background horizontal velocity, N is the buoyancy frequency defined by withρ being the (dimensionless) undisturbed density profile in the background, and H is the height of the water column. If there are no critical layers (i.e. c p − U = 0 for all z), for physically relevant N (z), the TG equation has an infinite set of discrete Due to the nonlinear nature of fluid flows, purely linear waves are a mathematical idealization. For large-amplitude waves or on time scales long enough for nonlinear effects to manifest themselves, results predicted by the linear theory do not agree with measurements. Weakly nonlinear theory attempts to better describe internal wave dynamics by expanding flow variables asymptotically and retaining corrections that correspond to finite amplitude (nonlinearity) and wavelength (dispersion). The most famous weakly nonlinear model is probably the Korteweg-de Vries (KdV) equation, given by (Lamb, 2005), where B is the horizontal structure function describing the propagation and evolution of the wave in the x-direction, c lw is the speed of linear long wave (i.e. waves with k = 0), and the parameters α and β measures the nonlinearity and dispersion of the wave, respectively. Analysis of the KdV equation shows that it has solitary wave solutions of the form 10 where a measures the wave amplitude and V is the nonlinear wave propagation speed. Internal solitary waves (ISWs) are translating waves of permanent form, which consist of a single wave crest. They are one of the most commonly observed types of internal waves in the field (e.g. Klymak and Moum (2003); Scotti and Pineda (2004); also see Helfrich and Melville (2006) for a more complete review).
While the KdV theory correctly predicts some properties of internal waves, it can only be expected to perform well within 15 certain asymptotic limits (e.g. small amplitude and long wave). For large amplitude waves, solutions of the KdV equation and its variations have been shown to be different from wave forms predicted by the fully nonlinear theory (Lamb and Yan, 1996;Lamb, 1999). Fully nonlinear ISWs can be computed by solving the nonlinear eigenvalue problem known as the Dubreil-Jacotin-Long (DJL) equation, which, in a zero background current, takes the form (Stastna and Lamb, 2002;Lamb, 2005) 20 In this equation, c isw is the solitary wave propagation speed (equivalent to V in the KdV theory), and η = η(x, z) is the vertical displacement of the isopycnal relative to its far-upstream depth. The DJL equation is equivalent to the full set of stratified Euler equations in a frame moving with the wave, where no assumptions are made with respect to the nonlinearity of the fluid flow.
Hence, its solutions are exact solitary wave solutions. For non-constant N , the DJL equation has no analytical solutions. In this 25 work, the DJL equation is solved numerically using the method described in Dunphy et al. (2011). The algorithm for solving the DJL equation is based on the variational scheme developed in Turkington et al. (1991), which seeks a solution iteratively that minimizes the kinetic energy, subject to the constraint that the scaled available potential energy (APE) is held fixed.

Governing equations and numerical method
We consider a right-handed Cartesian coordinate system, with the origin fixed at the lower left corner of the domain. The position vector is expressed as x = (x, z), with the x-axis directed to the right along the flat bottom and the z-axis pointing up towards the surface. The governing equations for the present work are the incompressible Navier-Stokes equations under the 5 rigid lid and Boussinesq approximations, given by (Kundu et al., 2012) where D/Dt is the material derivative defined by and In these equations, u = (u, w) describes the velocity field with u being the horizontal velocity and w being the vertical velocity, ρ describes the density field, p is the pressure, g is the gravitational acceleration,k is the unit vector in the vertical direction 15 (positive upwards), ν is the kinematic viscosity, and κ is the molecular diffusivity. As is the common practice under the Boussinesq approximation, the equations are in dimensional form, except that the density ρ and pressure p are scaled by the reference density ρ 0 . For all simulations, we fix the viscosity at ν = 10 −6 m 2 s −1 and the diffusivity at κ = 2 × 10 −7 m 2 s −1 .
This gives a Schmidt number Sc = ν/κ = 5. Periodic boundary conditions are used in the horizontal direction, and free-slip boundary conditions are used at the top and bottom boundaries. We note that no-slip conditions have also been tested but the 20 difference is insignificant, since the majority of the velocity perturbations are found near the pycnocline. The effect of the Earth's rotation is neglected, and hence the simulations are performed in an inertial frame of reference.
A complete description of the numerical model used in this study can be found in Subich et al. (2013), where a detailed validation and accuracy analysis through several test cases is also given. The model employs a spectral collocation method, which yields highly accurate results at moderate grid resolutions. From a purely numerical point of view, the spectral method 25 requires a minimum of two points in the horizontal direction in order to completely resolve a wave (Trefethen, 2000). Nevertheless, in this study we employ high resolution in order to better resolve the thin pycnocline, with at least ten grid points in the pycnocline and across the short waves. For spatial discretization, equally spaced grid points are used in both horizontal and vertical directions. As appropriate for the boundary conditions, the Fourier transform is applied in the x-direction, whereas the Fourier sine or cosine transform is applied in the z-direction depending on the variable of interest. For time stepping, the We focus on flows in a quasi two-layer stratification with a dimensionless density difference ∆ρ = 0.01, for which the Boussinesq approximation can be safely adopted. The background density profile, non-dimensionalized by the reference density ρ 0 , is given bȳ where z 0 is the location of the pycnocline and d is the half-width of the pycnocline. The specific location of the pycnocline does not affect dynamics of the interaction between the ISW and the linear waves in general, except for the case where the 15 pycnocline is close to the surface such that the ISW could be breaking (Lamb, 2002(Lamb, , 2003. In this work, we set z 0 = 0.4 m in order to avoid this situation. The thickness of the pycnocline can affect the gradient Richardson number through the buoyancy frequency profile it determines, which may have an impact on the interaction. However, this topic is not the focus of the present work (see section 6). In this work, we simply set d = 0.01 m for all simulations.
The initial solitary wave is specified by interpolating a solution of the DJL equation (5) onto the ISW subdomain. We Table 1. Solitary wave parameters. Here, the amplitude is measured by the maximum isopycnal displacement ηmax, and the wavelength is measured by the horizontal velocity profile along the inviscid upper boundary according to the formula λisw = 2(xR − xL), where xR and xL satisfy the equation u(xR, Lz) = u(xL, Lz) = 0.5umax.  Re based on the amplitude and maximum wave-induced current as

Amplitude Propagation speed Maximum current Minimum current Wavelength Reynolds number Richardson number
While there are a variety of Reynolds number estimates available in the literature, this simple estimate is more relevant to the length and velocity scales set by the ISW. The gradient Richardson number Ri is defined by It measures the ratio between the strength of the stratification and the shear stress in a parallel shear flow. The horizontal velocity profile of the ISW, together with the Richardson number contours, are shown in figure 3. The figure shows that the Richardson number has a local minimum in the high shear region near the pycnocline along the wave crest, and is very large outside this region. We note that while Ri min given in table 1 is slightly smaller than the critical Richardson number 0.25, the Richardson number criterion Ri < 0.25 is a necessary but not sufficient condition for linear stability in a parallel shear flow. Table 2. Linear wave parameters. In the case labels, O indicates an "overtaking" collision and H indicates a "head-on" collision, and the proceeding digits correspond to the wavelength of the linear waves.
Case label Wavelength Wave number Phase speed Group speed Time scale only when Ri is considerably smaller than 0.25 over a region long enough for perturbations to amplify in space (Lamb and Farmer, 2011).
We perform a suite of simulations in which the solitary wave propagates to the right and interacts with a small-scale wave packet initialized from linear waves. The linear waves are specified by solving the TG equation (1) numerically using a pseudospectral technique (Trefethen, 2000) in the linear wave subdomain. In order to ensure a smooth transition across the boundaries 5 between the ISW subdomain and the linear wave subdomain, an envelope function is applied to the amplitude of linear waves.
The particular form of the envelope function used here is given by although by testing other forms we found that results are not sensitive to the exact shape of the envelope.
We consider linear waves of wavelengths ranging from 0.2 to 0.6 m, whose parameters are given in table 2. For waves with 10 a wavelength less than 0.2 m, nonlinear self-interaction (Sutherland, 2016) becomes evident and may affect the interaction.
For waves longer than 0.6 m, they may no longer be considered "short" in comparison to the ISW width and hence will not be discussed. We examine two types of interaction in particular. An "overtaking collision" means that the ISW and the linear waves propagate in the same direction, whereas a "head-on" collision means that the two propagate in the opposite direction.
The amplitude of linear waves is set to be 1 cm for all cases. According to the linear theory, the propagation of linear internal 15 waves is independent of their amplitude, at least in the limit of small-amplitude waves. In fact, simulations with an amplitude of 2 cm have produced quantitatively similar results, and thus will not be discussed in this paper. For each experiment, we measure the time, τ , over which the solitary wave (which moves at the speed c isw ) and the linear wave packet (which moves at the speed c g ) experience a full collision cycle by defining At t = τ , the location of the solitary wave relative to the linear wave packet is approximately the same as it was in the initial field. In the figures, reported time T is scaled by this quantity such that T = t/τ .
Finally, we would like to mention that these linear waves are in fact not purely linear during the simulations. However, by scaling the relevant terms (i.e. B t and BB x ) in the KdV equation (3) using the amplitudes and wavelengths of these waves, we found that the time scale at which the nonlinearity becomes important is on the order of 1000 seconds, at least for waves with 10 an amplitude of 1 cm and a wavelength larger than 0.2 m. In contrast, the time scale of the interaction, as indicated in table 2, is on the order of 100 seconds. Hence, for clarity of notation, the small-scale waves will still be referred to as "linear waves", as opposed to the "solitary wave" or the "ISW", in the remainder of this paper. behind the solitary wave has a spatial structure that is completely different from the initial linear waves. To demonstrate that such deformation of linear waves does not occur naturally but is a result of the collision, we performed an additional simulation with the same linear wave packet but without the solitary wave. The resulting density field at T = 1 is shown in panel (d).
The density profiles of the case O6 are shown in figure 5. The initial density profile, visible in panel (a), shows again the disparity between the solitary wave and the linear waves, though in this case the wavelength of the linear waves is three times 25 larger than that in the previously discussed case (or 22.5 % of the wavelength of the solitary wave). Panels (b) and (c) show, however, that unlike in the previously discussed case, the linear waves are able to retain their spatial structure throughout the collision. The amplitude is also maintained, suggesting that energy loss during the collision is small. Instead, comparison to panel (d), the corresponding density profile obtained from the simulation with linear waves only, suggests that the primary effect of the collision on the linear waves is a phase shift, as indicated by the vertical lines. We will revisit the energy loss in 30 these cases below in section 5. To understand what causes the deformation of short waves in these cases, we performed an analysis similar to Stastna et al. (2015), in particular their figure 9. We first extracted the background horizontal velocity and buoyancy frequency profiles at  (head-on) collision case have a wave-like structure above (below) the pycnocline. We also note that if there is no velocity shear in the background, the vertical structure of linear waves of a given wavelength (e.g. λ = 0.2 m in this case, as indicated by a dotted curve) depends only on the stratification, regardless of the direction they propagate. Moreover, the vertical structure with respect to the pycnocline center is essentially unchanged. This suggests that the velocity shear in the background alters the vertical structure of the short waves in a nonlinear manner and leads to the observation that a head-on collision manifests In figure 7 (b), we made the same plot for linear waves with a wavelength of 0.6 m. The figure shows that the key difference in the initial structure function is that it has a non-negligible value over a much larger vertical extend. As a result, at some given depth, in particular that near the deformed pycnocline center, changes of the vertical structure functions from their initial state are much less dramatic in both overtaking and head-on collision cases. Therefore, longer waves are able to adapt to the ISW-induced background environment more easily and hence are more likely to survive the collision with the solitary wave. 5 We also note that formally changing the amplitude of linear waves does not change their vertical structures and thus does not affect the dynamics of the collision process, though in practice larger amplitude waves are expected to have a different (i.e. Stokes wave) structure.

Comparison to mode-1-mode-2 collision
The above analysis suggests that as the linear waves enter into the solitary wave-induced background state, they are subject 10 to a modified stratification and a velocity shear due to solitary wave-induced current, and it is this velocity shear across the deformed pycnocline that leads to the deformation of short waves. This process is in many ways similar to that found in Stastna et al. (2015). However, a key difference is that the disintegration of mode-2 waves due to the collision is much less dependent on their wavelength. To compare and contrast with their results, we performed an additional simulation, with mode-2 waves of amplitude of 1 cm and wavelength of 0.6 m interacting with the same ISW with an overtaking collision. The phase and 15 group speeds of the mode-2 waves are c p = 1.43 cm s −1 and c g = 1.35 cm s −1 , respectively, much smaller than their mode-1 counterparts. Figure 8 shows that after the collision with the ISW, the mode-2 waves are almost completely destroyed, except for some mode-1 like disturbances found near x = 8 m in panel (b).
In figure 9, we plotted the vertical structure functions for mode-2 waves in the ISW-induced background environment. The show characteristics of mode-1 waves, with essentially no perturbation below (above) the pycnocline for the overtaking (headon) case. This is similar to figure 9 (b) in Stastna et al. (2015), but fundamentally different from our figure 7 (b), implying that mode-1-mode-2 collisions are different from mode-1-mode-1 collisions. In fact, mode-2 waves were unable to maintain their coherent structure after the collision with mode-1 waves in all simulations in Stastna et al. (2015). Recent experiments (M. Carr, personal communication) suggest that the situation is more complex when the mode-1 wave amplitude is comparable to 5 the mode-2 wave amplitude, though it is unclear if such a situation had relevance to situations in the ocean.

Change of phase speed
Recall from figure 5 that a secondary effect of the interaction is a phase shift of the linear waves. To explain this observation, consider the linear long wave speed c lw in a two-layer stratification, defined by 10 where h 1 is the upper layer depth, h 2 is the lower layer depth and H is the total depth. The long wave speed sets the limit of the phase speed of linear waves in a two-layer stratification such that c p approaches c lw as the wavelength approaches infinity.
Thus c lw provides a good estimate of the maximum phase speed in a quasi-two layer stratification. Using the long wave speed as a guide, we note that the phase speed reaches its maximum value when h 1 = h 2 (i.e. when the two layers are equal in depth), provided other parameters (e.g. wavelength) remain constant. In our simulations, since we consider an ISW of depression, the 15 pycnocline at the wave crest is close to the mid-depth. Hence, the linear waves will experience an increase in phase speed as they propagate through the ISW.
In addition to the stratification, the presence of background current will also modify the phase speed. In figure 10, we explore the change of phase speed due to the presence of ISW-induced shear current for mode-1 linear waves. For overtaking collisions shown in panel (a), in the long wave limit, the phase speed in the shear background current is very close to that in a zero background current. However, in the short wave limit, the figure shows that the phase speed in the shear background current approaches the maximum ISW-induced current, whereas the phase speed in a zero background current approaches zero instead. This again suggests that short mode-1 waves are more likely to be influenced by the nonlinear interaction with ISW.
In particular, the critical wavelength that determines whether the phase speed is significantly influenced by the shear current is approximately 0.5 m, where the phase speed in a zero background current intersects the maximum velocity of the shear spectrum of wavelength, since the phase speed in a zero background current is less dependent on the wavelength and is always much smaller than the maximum velocity of the shear current. This is also consistent with the fact that mode-2 waves are less persistent after nonlinear interactions with the ISW. The fact that the ISW-induced maximum current essentially sets the lower limit for the phase speed of short waves implies that a critical layer does not exist for the ISW used in our simulations (as well as those with smaller amplitude). While the above analysis is performed for overtaking collisions (i.e. for linear waves 5 propagating to the right), we also examined head-on collisions. As shown in panel (b), in the short wave limit the behavior of the phase speed as a function of wavelength is very similar to that in the cases of an overtaking collision, except that now the phase speed is approaching the minimum current induced by the ISW.
We would like to note that given the nonlinear nature of the ISW, the interaction is indeed a nonlinear process, and thus the linear theory can only provide some rough guide for the flow behavior. To measure the nonlinearity of the fluid flows, we shall 10 introduce the Froude number which, in the context of internal wave dynamics, is usually defined as where U is the background current and c is the phase speed of the linear waves. The flow is said to be critical if Fr = 1, in which case the nonlinear effects are dominant. In our simulations, for any x the vertically integrated U is essentially zero since the flow is non-divergent in the simulation domain. Thus, a better estimation of U would be the effective horizontal velocity in 15 a reference frame moving with the ISW, which is essentially −c isw . In this reference frame, the estimated c would be −c isw +c p where c p > 0 for an overtaking collision and c p < 0 for a head-on collision. Hence, we can define the Froude number in a reference frame moving with the ISW as (16) Figure 11 shows the Froude number of the linear waves in the ISW-induced background current as a function of wavelength.

20
The figure shows that Fr < 1 for an overtaking collision and Fr > 1 for a head-on collision. In both cases, Fr approaches 1 in the short wave limit since short waves propagate slower. This implies that in the interaction of ISWs with short waves, the nonlinear effects become more important as the wavelength of short waves becomes smaller.

25
A function that describes a physical process can be represented either in the physical space or in the Fourier space. The two different representations are connected through the Fourier transform. Suppose f is a function of position x in the physical space, then the corresponding Fourier transformed variable F is a function of the horizontal wave number k and is given by If x is bounded, then k takes discrete values Parseval's theorem states that the total power in a signal is the same whether it is computed in the physical space or in the Fourier space (Press et al., 2002). That is, 5 From this theorem, we can define the power spectral density (PSD) of the function f as The PSD is a function of the wave number k. It can be interpreted as the strength of the signal at each wave number. For this reason, it provides a powerful tool for analyzing physical processes.
In the remainder of this section, we compute the PSD of horizontal velocity in the layer above the pycnocline and use it 10 to estimate the amount of wave energy being transferred during the collisions. The location of horizontal layer chosen for the analysis is z = 0.43 m (i.e. 3 cm above the pycnocline), though we have also calculated the PSD at other depths and found that results are not sensitive to the particular choice of horizontal layer. The PSD profiles of the initial horizontal velocity fields for some of the overtaking collision cases (O2, O3, O4 and O6) are plotted in figure 12. The figure clearly shows the wave number peaks due to the small-scale waves, which occur at considerably larger wave numbers than those associate with the 15 ISW spectrum (the peak near k = 0). This suggests that these small-scale waves are indeed "short" in comparison to the ISW width. For each simulation, we scale the PSD computed at the scaled time T = 1 by the maximum PSD of the considered linear waves observed in the initial field. According to Parseval's theorem, this ratio remains the same when mapped back into the physical space. Although only the horizontal velocity is used here, computation of the PSD of the vertical velocity yields quantitatively similar results, as it usually decays in a way similar to that of the horizontal velocity due to the interaction. Thus, 20 the scaled PSD of horizontal velocity represents the relative strength of horizontal current at T = 1 and hence provides an estimate of the percentage of kinetic energy remaining after one full collision cycle.

Reduction of wave energy
In figure 13 we examine the energy reduction of linear waves due to collision by plotting the PSD of horizontal velocity in the wave number domain. The figure shows that in all cases, there is a net loss of wave energy due to the collision. It also shows that for a given solitary wave, the wavelength of the linear waves (which remains unchanged after the collision) is the single most important factor that determines the amount of PSD (and hence wave energy) remaining after the collision. While the longest 5 waves may retain as much as 85% of the kinetic energy they had initially, the shortest waves lose almost all of their initial energy such that the peaks of the PSD can hardly be distinguished from background noise. Among other factors, a head-on collision is slightly more efficient in destroying the linear waves than an overtaking collision, except for the small wavelength limit. This may be explained by the fact that during a head-on collision, the structure function (especially its peak) shifts further away from its initial state than during an overtaking collision, as shown in figure 7, such that the new, ISW-induced background 10 environment is more difficult for the linear waves to adjust to. In contrast, the initial amplitude of the linear waves has very little impact on the net energy transfer due to collision, since curves produced from simulations in which the linear waves have an amplitude of 2 cm (not shown) are almost exactly the same as their smaller amplitude counterparts shown in panel (a).
Though we did not consider large amplitude short waves, since these will have their own complex dynamics.
The maximum value of the scaled PSD as a function of wavelength is plotted in figure 14, along with a quantitative mea-  suggest that the maximum values of the scale PSD at T = 1 are at the level of 90 % but are never larger than 100 %, implying that very little wave energy is being transferred from the short waves during the collision and that no energy is transferred from the ISW to the short waves. For the longest waves, the slight decrease in PSD is at a similar level to viscous dissipation. We note that this observation is consistent with the result shown in figure 10, since above the critical wavelength λ = 0.5 m, very little energy exchange occurs due to the interaction.

5
For the cases O6 and H6, simulations were performed for an extended period of time in order to allow for repeated collisions between the solitary wave and the linear wave packet. For each of these cases, four complete collision cycles were observed, and the scaled PSD has been computed at T = 1, 2, 3 and 4, as shown in figure 15. The corresponding measurement of the scaled PSD at each peak is given in table 4. The figure and table suggest that for both cases, the scaled PSD is reduced after each subsequent collision, down to 60.20 % in the case O6 and 33.61 % in the case H6 at T = 4. Nevertheless, they are still larger

Influence of the interaction on the ISW
During the collision with the linear wave packet, the solitary wave is also affected by the linear waves that pass through it.
We note however, that the kinetic energy carried by the linear waves is much smaller than that carried by the solitary wave, and hence the impact of linear waves on the solitary wave is also small. Here, we define the kinetic energy (KE) per unit mass following standard convention (which drops the reference density and hence changes the dimensions of the quantity) by We found that when measured in terms of vertically integrated kinetic energy at the wave crest, the linear waves are about 1 % as energetic as the solitary wave.
To analyze changes in the solitary wave and determine if they are results of the collision, we performed an additional simulation with the same solitary wave but without the linear waves. We estimated the vertically integrated KE at the crest of T) in figure 16 over one complete collision cycle. Mathematically, this quantity is computed as where A is the normalization factor defined as the maximum vertically integrated KE of the initial solitary wave. The subscript has increased by at least 1 %. We are able to confirm that such an energy increase in the solitary wave is robust since we have also performed additional simulations with a longer linear wave packet (not shown), and found that the maximum vertically integrated KE increases approximately linearly with respect to the length of the wave packet. On the other hand, for waves with a wavelength λ = 0.6 m shown in panel (b), energy increase in the solitary waves after the collision is insignificant. In all 10 cases, the curves shown demonstrate periodicity associated with their particular wavelengths.
We have also attempted to detect the phase shift of the solitary waves from the locations of maximum vertically integrated KE. However, we found that such a phase shift, if it exists at all, is on the order of millimeters. In other words, the detected phase shift is on the grid scale and is subject to numerical error. For this reason, the results are not shown here.

15
In this work we performed two-dimensional direct numerical simulations to study the interaction between a large-scale, fully nonlinear ISW and small-scale linear internal waves. We demonstrated that there was a net energy transfer from the small-scale linear waves to the large-scale solitary waves. This contrasts the conclusion in Broutman and Young (1986), made for a different type of internal wave interaction, that energy is transferred from large-scale waves to small-scale waves. Our simulation results suggest that during the interaction, the solitary wave essentially acts as a filter through which only long waves may pass. For waves with a smaller wavelength, the interaction leads to a reduction of their initial energy and a destruction of their spatial structure. These processes occur in a background state set by the solitary wave-induced stratification and current. During the 5 interaction, adjustment of the short waves to this new background environment extracts their wave energy and modifies the wave structure. The fact that short waves may not survive the interaction with a solitary wave, or more generally, any localized nonlinear background environment which both deforms the pycnocline and induces shear, implies that the observed spectrum of wavelengths of internal waves in locations with large amplitude ISWs (such as Straits) is likely to be poor in short waves.
At the time of writing we are unaware of measurements to support or contradict this hypothesis. 10 We performed analysis based on linear wave theory and showed that during the nonlinear interaction with the ISW, the destruction of short linear waves occurs primarily due to the presence of ISW-induced velocity shear, which alters the vertical structure of the short waves in a nonlinear manner, leading to significant wave amplitudes on only one side of the deformed pycnocline center. On the other hand, a shift of the location of pycnocline plays a secondary role during the collision, as its main effect is to alter the propagation speed of the linear waves, and shift the location of the maximum of the vertical 15 structure downward. However, the vertical structure is unchanged with respect to the pycnocline center. We also demonstrated that a critical layer is not present during the collision, regardless of the wavelength of the linear waves, since the phase speed approaches the maximum ISW-induced current asymptotically as the wavelength approaches zero.
A clear avenue of future research is to explore the parameter space, in particular the Richardson number effect, of the solitary wave. In the present work we studied the ISW whose minimum Richardson number is 0.246. While none of the 20 simulations show evidence of the generation of shear instability, this does not necessarily mean that the wave-wave interaction considered in the present work is Richardson number independent. Moreover, Lamb and Farmer (2011) showed that it is not only the minimum Richardson number, but also the length of the unstable region with a low Richardson number relative to the wavelength of ISW, that is jointly responsible for the generation of shear instability in an ISW. It is thus reasonable to assume that the relative length of the region with a low-Richardson number in the ISW also has an influence on the wave-wave 25 interaction. Future research will explore these effects in detail.
We note that our findings are in many ways similar to those in Stastna et al. (2015). Their study also concluded that the direction of energy transfer during the interaction is from the small-scale, weakly nonlinear wave (i.e. the mode-2 wave) to the large-scale solitary wave (i.e. the mode-1 wave), and that such energy transfer is more efficient when a head-on, instead of overtaking, collision is involved. The main difference is that in mode-1-mode-1 interaction, there is a cutoff determined by 30 the wavelength of short waves, above which the small-scale waves maintain their structure after interaction, whereas mode-1mode-2 interaction is much less dependent on the wavelength. In mode-1-mode-1 interaction, this cutoff corresponds to the wavelength at which the phase speed of the short waves upstream of the solitary wave exceeds the maximum ISW-induced current. In mode-1-mode-2 interaction, however, this cutoff does not exist since the maximum ISW-induced current is always larger than the phase speed for any given wavelengths.

35
While all of the simulations discussed in this work are performed on the laboratory scale, the scaling-up of the current experiments to the field scale is left as a topic for future work. When the field scale is considered, waves with a much larger range of wavelengths can be expected to breakdown, including short waves affected by self-interaction (Sutherland, 2016).
Also, a higher Reynolds number implies that the overturning seen in figure 6 may eventually lead to significant overturns. The three-dimensionalization of the flow field should also be examined. As shown in Andreassen et al. (1994), two-dimensional 5 models are unable to describe properly the physics or the consequences of the wave breaking process, in particular that induced by the presence of a critical layer. We also note that in two dimensions, the only possible form of wave-wave interaction is either an overtaking collision or a head-on collision. However, observational evidences (e.g. Quaresma et al. (2007); in particular, see their figures 2 and 8) suggest that internal waves do not generally propagate parallel to each other but may interact at different angles. The effects of directionality of wave propagation is another topic that can be considered in forthcoming studies.