Asymmetric behavior of surface waves induced by an underlying interfacial wave

TL;DR

This paper introduces a weakly nonlinear model that explains the asymmetric surface wave behaviors caused by underlying interfacial solitary waves in a two-layer fluid, relevant for understanding oceanic internal waves.

Contribution

The study presents a novel, tractable model capturing asymmetric surface wave behaviors induced by internal waves, enhancing theoretical and numerical understanding.

Findings

Surface waves shorten at the leading edge and lengthen at the trailing edge.

Surface wave group velocities differ significantly between edges.

Surface wave amplitudes are higher at the leading edge.

Abstract

We develop a weakly nonlinear model to study the spatiotemporal manifestation and the dynamical behavior of surface waves in the presence of an underlying interfacial solitary wave in a two-layer fluid system. We show that interfacial solitary-wave solutions of this model can capture the ubiquitous broadening of large-amplitude internal waves in the ocean. In addition, the model is capable of capturing three asymmetric behaviors of surface waves: (i) Surface waves become short in wavelength at the leading edge and long at the trailing edge of an underlying interfacial solitary wave. (ii) Surface waves propagate towards the trailing edge with a relatively small group velocity, and towards the leading edge with a relatively large group velocity. (iii) Surface waves become high in amplitude at the leading edge and low at the trailing edge. These asymmetric behaviors can be well quantified in the theoretical framework of ray-based theories. Our model is relatively easily tractable both theoretically and numerically, thus facilitating the understanding of the surface signature of the observed internal waves.

Equations

$${\varphi }_{ixx}+{\varphi }_{izz}=0.$$

$${\xi }_{1t}+{\varphi }_{1x}{\xi }_{1x}={\varphi }_{1z},\text{ at }z=h_1+{\xi }_1,$$

$${\xi }_{2t}+{\varphi }_{1x}{\xi }_{2x}={\varphi }_{1z},\text{ at }z={\xi }_2,$$

$${\xi }_{2t}+{\varphi }_{2x}{\xi }_{2x}={\varphi }_{2z},\text{ at }z={\xi }_2,$$

$${\varphi }_{2z}=0,\text{ at }z=-h_2,$$

$${\varphi }_{1t}+\frac{1}{2}\left({\varphi }_{1x}^2+{\varphi }_{1z}^2\right)+g{\xi }_1=0,\text{ at }z=h_1+{\xi }_1,$$

$${\rho }_1\left({\varphi }_{1t}+\frac{1}{2}\left({\varphi }_{1x}^2+{\varphi }_{1z}^2\right)+g{\xi }_2\right)$$

$$\text{at }z$$

$$a/h_1=\alpha \ll 1,h_2/h_1=O(1).$$

$$h_1^2/L^2=\beta \ll 1,h_2/h_1=O(1).$$

$$x=L{x}^{\ast },z=h_1{z}^{\ast },t=(L/U_0)t^{\ast },$$

$$\left({\varphi }_1,{\varphi }_2\right)=\left(aL{U}_0/h_1\right)\left({\varphi }_1^{\ast },{\varphi }_2^{\ast }\right),({\xi }_1,{\xi }_2)=a({\xi }_1^{\ast },{\xi }_2^{\ast }),$$

$$\beta {\varphi }_{i{x}^{\ast }{x}^{\ast }}^{\ast }+{\varphi }_{i{z}^{\ast }{z}^{\ast }}^{\ast }=0,\mathrm{for}i=1,2.$$

$${\varphi }_1^{\ast }={\varphi }_1^{\ast (0)}+\beta {\varphi }_1^{\ast (1)}+O({\beta }^2),$$

$${\varphi }_1^{\ast (0)}={\varphi }_1^{\ast (0)}(x^{\ast },t^{\ast }).$$

$${\varphi }_{1x^{\ast }{x}^{\ast }}^{\ast (0)}+{\varphi }_{1z^{\ast }{z}^{\ast }}^{\ast (1)}=0,\alpha {\xi }_2^{\ast }<z^{\ast }<h_1^{\ast }+\alpha {\xi }_1^{\ast },\text{ with }{h}_1^{\ast }=1,$$

$${\varphi }_{1z^{\ast }}^{\ast (1)}={\xi }_{1t^{\ast }}^{\ast }+\alpha {\varphi }_{1x^{\ast }}^{\ast (0)}{\xi }_{1x^{\ast }}^{\ast },\text{ at }{z}^{\ast }=h_1^{\ast }+\alpha {\xi }_1^{\ast }.$$

$${\varphi }_1^{\ast (1)}=D_1{\xi }_1^{\ast }{Z}_1-\frac{1}{2}{\varphi }_{1x^{\ast }{x}^{\ast }}^{\ast (0)}{Z}_1^2,$$

$${\varphi }_1^{\ast }={\varphi }_1^{\ast (0)}+\beta \left(D_1{\xi }_1^{\ast }{Z}_1-\frac{1}{2}{\varphi }_{1x^{\ast }{x}^{\ast }}^{\ast (0)}{Z}_1^2\right)+O({\beta }^2).$$

$${\eta }_{1t^{\ast }}^{\ast }+{\left({\eta }_1^{\ast }{\tilde{u}}_1^{\ast }\right)}_{x^{\ast }}-\frac{1}{6}\beta {\left(h_1^{\ast }\right)}^3{\tilde{u}}_{1x^{\ast }{x}^{\ast }{x}^{\ast }}^{\ast }-\frac{1}{2}\beta {\left(h_1^{\ast }\right)}^2{\xi }_{1t^{\ast }{x}^{\ast }{x}^{\ast }}+O(\alpha \beta ,{\beta }^2)=0,$$

$${\eta }_1^{\ast }=h_1^{\ast }+\alpha {\xi }_1^{\ast }-\alpha {\xi }_2^{\ast },{\tilde{u}}_1^{\ast }={\varphi }_{1x^{\ast }}^{\ast (0)}.$$

$${\tilde{u}}_{1t^{\ast }}^{\ast }+\alpha {\tilde{u}}_1^{\ast }{\tilde{u}}_{1x^{\ast }}^{\ast }+g{\xi }_{1x^{\ast }}^{\ast }+O(\alpha \beta ,{\beta }^2)=0,$$

$${\varphi }_{1x^{\ast }}^{\ast }={\tilde{u}}_1^{\ast }+\beta {\xi }_{1t^{\ast }{x}^{\ast }}^{\ast }{Z}_1-\frac{1}{2}\beta {\tilde{u}}_{1x^{\ast }{x}^{\ast }}^{\ast }{Z}_1^2+O(\alpha \beta ,{\beta }^2).$$

$${\overline{u}}_1^{\ast }={\tilde{u}}_1^{\ast }-\frac{1}{2}\beta h_1^{\ast }{\xi }_{1t^{\ast }{x}^{\ast }}^{\ast }-\frac{1}{6}\beta {\left(h_1^{\ast }\right)}^2{\tilde{u}}_{1x^{\ast }{x}^{\ast }}^{\ast }+O(\alpha \beta ,{\beta }^2);$$

$${\tilde{u}}_1^{\ast }={\overline{u}}_1^{\ast }+\frac{1}{2}\beta h_1^{\ast }{\xi }_{1t^{\ast }{x}^{\ast }}^{\ast }+\frac{1}{6}\beta {\left(h_1^{\ast }\right)}^2{\overline{u}}_{1x^{\ast }{x}^{\ast }}^{\ast }+O(\alpha \beta ,{\beta }^2),$$

$${\overline{u}}_1^{\ast }(x^{\ast },t^{\ast })=\frac{1}{{\eta }_1^{\ast }}{\int }_{\alpha {\xi }_2^{\ast }}^{h_1^{\ast }+\alpha {\xi }_1^{\ast }}{\varphi }_{1x^{\ast }}^{\ast }(x^{\ast },z^{\ast },t^{\ast })d{z}^{\ast }.$$

$${\eta }_{1t}+({\eta }_1{\overline{u}}_1{)}_x=0,{\eta }_1=h_1+{\xi }_1-{\xi }_2,$$

$${\eta }_{2t}+({\eta }_2{\overline{u}}_2{)}_x=0,{\eta }_2=h_2+{\xi }_2,$$

$${\overline{u}}_{1t}+{\overline{u}}_1{\overline{u}}_{1x}+g{\xi }_{1x}-\frac{1}{3}{h}_1^2{\overline{u}}_{1xxt}+\frac{1}{2}{h}_1{\xi }_{2ttx}=0,$$

$${\overline{u}}_{2t}+{\overline{u}}_2{\overline{u}}_{2x}+g{\xi }_{2x}+{\rho }_{r}g{\eta }_{1x}$$

$$-\frac{1}{2}{\rho }_r{h}_1^2{\overline{u}}_{1xxt}+{\rho }_r{h}_1{\xi }_{2ttx}-\frac{1}{3}{h}_2^2{\overline{u}}_{2xxt}=0,$$

Full text

Keywords. interfacial waves, surface waves, ray-based theory

Mathematics Subject Classification. 76B55; 76B07; 35L05; 65M22

1 Introduction

Internal waves with large amplitudes and long wavelengths are widely observed in coastal ocean regions, and are believed to be important for transferring momentum, heat, and energy in the ocean [30, 28]. Because of their strong turbulent mixing and breaking, they can influence many ocean processes, such as nutrient supply, sediment and pollutant transport, acoustic transmission, and interaction with man-made structures [20, 3]. Generated by tidal flow, internal waves usually can propagate thousands of kilometers from their source before dissipation, sloping, and breaking extinguish them [28, 5]. In recent decades, observation data for internal waves and their corresponding surface signature have been recorded using in situ measurements and Synthetic Aperture Radar (SAR) in many coastal seas worldwide [38, 20, 3].

Many previous studies have investigated the interaction between interfacial waves (IWs) and surface waves (SWs) in two-layer fluid systems. These works can be broadly subdivided into two classes. The first class focuses on the issue of addressing the direct numerical simulation of the Laplace equations for velocity potentials with the appropriate boundary conditions in a two-layer fluid system [2, 43]. The system can be formulated as a Hamiltonian system for a set of canonical variables [4, 2]. To evaluate the time derivatives of these canonical variables, the high-order spectral (HOS) method [2] is employed to solve the boundary value problems of Laplace equations. By applying the HOS method, it was found that energy can be transferred from SWs to IWs in the two-layer fluid system [43]. However, the direct numerical simulations are typically expensive, and it is also challenging to extract the mechanism underlying their results.

The second class focuses on the reduced models in the two-layer fluid system. To overcome the difficulty stemming from expensive computations, a common approach is to study a reduced model of the two-layer Euler system via multi-scale analysis [32, 27, 23, 19, 39, 6, 8, 29, 17]. These reduced models describing two-layer fluids mostly focus on three aspects:

[a] Traveling wave solutions: Traveling-wave solutions exhibiting oscillations are found in the reduced models. For example, generalized solitary waves with non-decaying oscillations along their tails in addition to the solitary pulse were found in a long-wave model [18, 39, 22]. Besides the generalized solitary waves, multi-humped solitary waves with a finite number of oscillations riding on the solitary pulse were found in a fully-nonlinear long-wave model [8].

[b] Ray-based theories: Many ray-based studies take a statistical viewpoint of SWs modulated by a near-surface current induced by IWs [24, 37, 9, 6, 5]. These studies invoke phase-averaged models based on a wave-balance equation and ray-based theory [9, 6], which can be applied to the remote-sensing observations of IWs via their surface manifestations.

[c] Resonant excitations: When two different modes coexist in a fluid system, a resonant interaction becomes possible between the modes to aid in transferring energy in the ocean [41]. Class 3 triad resonance is regarded as being responsible for the surface signature of the underlying IWs [41, 38, 27, 23, 16, 1]. Based on a class 3 triad resonance condition, many reduced models have been derived for the interfacial and surface waves[37, 32, 27, 23, 42, 36, 29, 17]. A detailed discussion of surface signature phenomena of IWs was presented in [16, 14, 15, 17] using a coupled Korteweg–de Vries (KdV) and a linear Schrödinger model. Narrow rough regions containing surface ripples were found and interpreted as a result of the energy accumulation in the localized bound states of the Schrödinger equation [17].

Different from previous works, we here develop a new reduced model to investigate the spatiotemporal manifestation of the small-amplitude SWs in the presence of an interfacial solitary wave in the two-layer setting. Based on model simulations, we demonstrate that our model is successful in characterizing many types of dynamical behavior of SWs, which can be well understood using the ray-based theories. In Sec. 2, we derive a reduced model for the two-layer fluid system. In Sec. 3, we analyze the basic properties of the model, including interfacial solitary-wave solutions and dispersion relations. In Sec. 4, we present the numerical scheme and examine its numerical convergence. In Sec. 5, we show the numerical results for the asymmetric behavior of SWs and quantify this asymmetric behavior using the ray-based theories. Conclusions and discussion are given in Sec. 6, and some mathematical details are presented in the Appendix.

2 The two-layer weakly nonlinear (TWN) model

We first introduce Euler equations for two immiscible layers of potential fluids with unequal densities. The two layers of fluids are assumed to be inviscid, irrotational, and incompressible. The unequal densities for the upper layer and for the lower layer are denoted by $\rho_{1}$ and $\rho_{2}$, respectively. Here, $\rho_{2}>\rho_{1}$ is assumed for the stable case. The horizontal and vertical coordinates are $x$ and $z$, respectively. We focus on the evolution of large-amplitude interfacial waves $\xi_{2}$ between the two fluid layers, and their coupling with the the overlaying free surface, $\xi_{1}$ [see Fig. 1]. The velocity potential $\phi_{i}$ ($i=1$ for the upper layer and $i=2$ for the lower layer) satisfies Laplace’s equation,

[TABLE]

The kinematic equations for the continuity of the normal velocity at the surface $h_{1}+\xi_{1}$, the interface $\xi_{2}$, and the flat topography $-h_{2}$ are given in the form

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $h_{1}$ and $h_{2}$ are the undisturbed thicknesses of the upper and lower layers, respectively. The dynamical equations for the continuity of pressure at the surface and the interface are the Bernoulli equations,

[TABLE]

[TABLE]

where $g$ is the gravitational acceleration.

For the small-amplitude approximation, we assume that the characteristic amplitude, $a$, of the IWs and SWs is much smaller than the thickness of the two fluid layers,

[TABLE]

For the long-wave approximation, we assume that the thickness of each fluid layer is much smaller than the characteristic wavelength, $L$, of the IWs and SWs,

[TABLE]

The two small parameters, $\alpha$ and $\beta$, control the nonlinear and dispersive effects, respectively. Based on the scaling (2.8)-(2.9), we may nondimensionalize all the physical variables by taking the original variables to be

[TABLE]

[TABLE]

where $U_{0}=\sqrt{gh_{1}}$ is the characteristic speed of the gravity waves. Here, all the variables with asterisks are assumed to be $O(1)$ in $\alpha$ and $\beta$. In the dimensionless, starred variables, Laplace’s equation (2.1) is formulated as

[TABLE]

We first focus on the upper fluid layer. The analogous derivation for the lower fluid layer follows a similar procedure. We seek an asymptotic expansion of $\phi_{1}^{\ast}$ in powers of $\beta$,

[TABLE]

which we use in the asymptotic analysis of the nondimensionalized problem of equations (2.1)-(2.7) for small values of the parameter $\beta$. From the leading $O(1)$ term in Eqs. (2.1) and (2.2), $\phi_{1}^{\ast(0)}$ is found to be independent of the height $z$,

[TABLE]

The $O(\beta)$ terms in Eqs. (2.1) and (2.2) yield the equations

[TABLE]

[TABLE]

The expression for $\phi_{1}^{\ast(1)}$ is obtained as

[TABLE]

where $Z_{1}=z^{\ast}-h_{1}^{\ast}-\alpha\xi_{1}^{\ast}$ and $D_{1}\xi_{1}^{\ast}\equiv\xi_{1t^{\ast}}^{\ast}+\alpha\phi_{1x^{\ast}}^{\ast(0)}\xi_{1x^{\ast}}^{\ast}$. Combining expressions (2.13) and (2.15), to the first power of $\beta$, we can obtain the solution $\phi_{1}^{\ast}$ as

[TABLE]

By integrating Eq. (2.1) once from $\alpha\xi_{2}^{\ast}$ to $h_{1}^{\ast}+\alpha\xi_{1}^{\ast}$ with respect to $z^{\ast}$, imposing the boundary conditions (2.2) and (2.3), and substituting the expression (2.16) into equation (2.1), we obtain the kinematic equation for the upper fluid layer,

[TABLE]

where

[TABLE]

Upon substitution of the velocity potential $\phi_{1}^{\ast}$ (2.16) into the dynamical boundary condition (2.6), we obtain the dynamical equation governing the motion of the upper fluid layer,

[TABLE]

where the equation has been differentiated with respect to $x^{\ast}$ once, the terms in the first power of $\beta$ are retained [the $O(\beta)$ terms happen to vanish in Eq. (2.18)], and terms of $O(\alpha\beta,\beta^{2})$ are dropped. From the velocity potential $\phi_{1}^{\ast}$, (2.16), we obtain the horizontal velocity $\phi_{1x^{\ast}}^{\ast}$ as

[TABLE]

By averaging (2.19) over the depth, we obtain the layer-mean horizontal velocity for the upper fluid layer,

[TABLE]

the corresponding inverse is

[TABLE]

where the layer-mean horizontal velocity is defined as

[TABLE]

After substituting Eq. (2.20) for the horizontal velocity $\widetilde{u}_{1}^{\ast}$, equations (2.17) and (2.18) provide Boussinesq-type equations governing the motion of the fluid in the upper layer.

Repeating a similar procedure, we can obtain the governing equations for the lower fluid layer. The final set of equations for the variables $(\xi_{1},\xi_{2},\overline{u}_{1},\overline{u}_{2})$, in the dimensional form, is

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\rho_{r}$ is the density ratio $\rho_{1}/\rho_{2}$, and

[TABLE]

are the layer-mean horizontal velocities. We refer to the model (2.21)-(2.24) as the two-layer weakly-nonlinear model or TWN model. The TWN model can also be obtained via a direct reduction from a fully nonlinear model (to which we refer as MCC model) given in [12, 8, 7]. For the numerical examples throughout the paper, all the variables and parameters are dimensionless and the parameters are fixed to be $(h_{1},h_{2},g,\rho_{1},\rho_{2})=(1,3,1,1,1.003)$. In particular, in our simulations, the characteristic length in both height and wavelength is $h_{1}$, the characteristic speed is $\sqrt{gh_{1}}$, and the characteristic time is $\sqrt{h_{1}/g}$ [see similar dimensionless forms used in numerical simulations in [13, 11]]. Note that, without loss of generality, our conclusions concerning broadening of IWs in Sec. 3.1 and asymmetric behavior of SWs in Sec. 5 also hold true for other parameter regimes of $(h_{1},h_{2},g,\rho_{1},\rho_{2})$.

3 Basic properties of the TWN model

In this section, we study some basic properties of the TWN model, including interfacial solitary wave solutions and dispersion relations.

3.1 Interfacial solitary-wave solution

To study the behavior of the overlaying SWs when an interfacial solitary wave moves beneath the surface, we first seek interfacial solitary-wave solutions of the TWN model (2.21)-(2.24). Then, we use these solitary-wave solutions as initial conditions for wave profiles and layer-mean velocities for the subsequent numerical simulations in Sec. 5. To look for the right-moving traveling waves that propagate with the nonlinear phase velocity $C_{s}$, we assume the following ansatz for the surface elevation, internal elevation, upper-layer velocity, and lower-layer velocity, $(\xi_{i},\overline{u}_{i})$ [$i=1,2$], in the system (2.21)-(2.24):

[TABLE]

Substituting this ansatz into Eqs. (2.21)-(2.22) and integrating once with respect to $X$ yields

[TABLE]

where we have assumed that $H_{i}\rightarrow h_{i}$ as $X\rightarrow\pm\infty$, and $h_{1}$ ($h_{2}$) is the undisturbed thickness of the upper (lower) fluid layer, respectively. Substituting the horizontal velocity (3.26) for $\overline{U}_{i}$ into Eqs. (2.23)-(2.24) and integrating once with respect to $X$ leads to the equations

[TABLE]

[TABLE]

where we have assumed that $H_{iX}$, $H_{iXX}\rightarrow 0$ as $X\rightarrow\pm\infty$ for $i=1,2$. Since explicit solutions to Eqs. (3.27)-(3.28) are difficult to establish, we numerically compute their solitary-wave solutions by applying the method in [26].

In Fig. 2(a), we show the numerical solutions of the TWN model for IWs with different amplitudes. For comparison, we also show the corresponding MCC and KdV solutions with the same amplitudes [13]. From Fig. 2(a), we can see that the TWN model and the MCC model can capture the broadening of internal waves that is often observed in the ocean. For instance, a single large internal wave in $340$ meters of water was observed in the northeastern South China Sea by [21]. The typical wavelength of the observed internal wave is longer than the KdV solution of the same amplitude that is used to fit this internal wave. It is worthwhile to mention that the broadening of interfacial solitary wave solutions can also be captured by other models [15, 25], not only by the MCC-type models.

To quantify this broadening, we introduce a measure of the effective width, $\Lambda_{s}$, for the interfacial solitary-wave solution [34], defined as

[TABLE]

Meanwhile, the effective width of MCC and KdV solutions are provided in the reference [13]. Figure 2(b) displays the comparison of effective width among the TWN solutions, MCC solutions, and KdV solutions. When the amplitude of waves is small, there is good agreement of the effective widths among all solutions. However, when the amplitude of the waves becomes large, discrepancy grows rapidly among these three solutions. When the amplitude increases to the limiting value, approximately $(h_{1}-h_{2})/2$, the TWN and MCC solutions become much broader than the KdV solutions. The maximum amplitude of our TWN model is approximately $(h_{1}-h_{2})/2$. Beyond the maximum amplitude, no solitary waves can exist for IWs.

3.2 Dispersion relations

We now investigate the dispersion relation of the TWN model (2.21)-(2.24) using linear analysis. By substituting the monochromatic solutions $(\xi_{i},\overline{u}_{i})\sim\exp[{\text{i}}(kx-\mu_{k}t)]$ into the system (2.21)-(2.24), the pure linear dispersion relation in the absence of shear between the frequency $\mu_{k}$ and the wavenumber $k$ can be obtained as

[TABLE]

Here, the shear is the interface and velocity jump induced by an interfacial solitary wave. The same dispersion relation can be found in [8].

Equation (3.30) has $4$ real roots in the oceanic regime (the density ratio $\rho_{r}$ is close to $1$). At the leading order in $1-\rho_{r}$, the dispersion relations of the two-mode waves, denoted by $\Omega_{k}$ and $\omega_{k}$ [see Fig. 3], can be approximated as

[TABLE]

and

[TABLE]

In the following, the two kinds of waves that correspond to the dispersion relations $\Omega_{k}$, (3.31), and $\omega_{k}$, (3.32), are referred to as the slow-mode waves and the fast-mode waves, respectively.

The modulated dispersion relation $\overline{\omega}_{k}$ in the presence of shear can be obtained by substituting $(\xi_{i},\overline{u}_{i})\sim(\Xi_{i},\overline{U}_{i})+\exp[{\text{i}}(kx-\overline{\omega}_{k}t)]$ into the system (2.21)-(2.24), where the shear is induced by an interfacial solitary wave $(\Xi_{i},\overline{U}_{i})$ in above Sec. 3.1. The resulting equation is

[TABLE]

where $\left|\cdot\right|$ denotes the determinant of the enclosed matrix. In the following Secs. 4 and 5, we will numerically study the TWN model in the right-moving frame with the nonlinear phase velocity $C_{s}$, that is, $T=t$ and $X=x-C_{s}t$. Note that the solitary-wave solutions $(\Xi_{i},\overline{U}_{i})$ are steady in time $T$ in this moving frame. Then, the modulated dispersion relation $\overline{\nu}$, corresponding to the moving frame $T$ and $X$, is given by

[TABLE]

where $\overline{\omega}_{k}$ is the modulated dispersion relation in Eq. (3.33) corresponding to the resting frame $t$ and $x$. Note that the dispersion relation $\overline{\nu}_{k}$ is independent of time $T$ since $(\Xi_{i},\overline{U}_{i})$ are steady in time $T$. Moreover, the wavelengths of $(\Xi_{i},\overline{U}_{i})$ are relatively long with respect to the characteristic wavelengths of fast-mode waves, so $(\Xi_{i},\overline{U}_{i})$ in the dispersion relation $\overline{\nu}_{k}$ can be locally treated as constant in $X$ space.

4 Numerical scheme

For numerical computations, we cast Eqs. (2.21)-(2.24) in the conservation form in the right-moving frame with the nonlinear phase velocity (soliton speed) $C_{s}$ as follows:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

and $T=t$, $X=x-C_{s}t$. The computational domain is set to be $[-M,M]$, with periodic boundary conditions. Even for an initially narrowly localized perturbation wave, radiation can be quickly emitted towards the two boundaries $x=-M$ and $x=M$. To eliminate possible reflected waves from these boundaries, two buffer zones in the regions $[-M,-M/2]$ and $[M/2,M]$ are established, and damping and diffusion terms are added to absorb the outgoing radiation. For numerical integration, we use the fourth-order Runge-Kutta method in time and the second-order collocation method in space [10]. The Kelvin-Helmholtz (KH) instability is suppressed by applying a low-pass filter [31]. (The wavenumbers for the KH instability are much larger than the wavenumbers of the trapped right-moving SWs as introduced in the following Sec. 5. Thus, these KH unstable wavenumbers are physically irrelevant in our computations.) In our simulations, we fix the parameter regime $(\rho_{1},\rho_{2},h_{1},h_{2},g)=(1,1.003,1,3,1)$ and the computational domain $M=300$. All the variables and parameters in our simulations are dimensionless.

We first focus on the evolution of initially unperturbed interfacial solitary-wave solutions. Figure 4 shows the spatiotemporal evolution of the IWs’ profile, $\xi_{2}$, for $0\leq T\leq 10000$. We can see from Fig. 4 that the IWs maintain their shape while traveling. This result is consistent with many experimental observations that large-amplitude internal waves typically can propagate over long distances with their shape virtually unchanged [28].

Next, we examine the numerical convergence in time and space of our scheme for initially unperturbed interfacial solitary-wave solutions [Eqs. (3.27)-(3.28)]. We compute the $l_{1}$ norm error for the SWs’ profile defined as

[TABLE]

where the reference solution $\xi_{1}^{\rm ref}$ is approximated by the numerical result obtained from a very small time step for the time accuracy test or from a very small spatial discretization for the spatial accuracy test. We can see from Fig. 5 that the scheme has fourth-order accuracy in time and second-order accuracy in space.

5 Asymmetric behavior of SWs in the presence of an underlying IW

In this section, we present our numerical results for the system (2.21)-(2.24) describing the behavior of a SW packet in the presence of an underlying interfacial solitary wave, and then compare them to the results of our theoretical analysis using ray-based theories. First, we initialize the SWs’ height $\xi_{1}$ to be a profile composed of a sufficiently-long-wavelength interfacial-solitary-wave solution and a localized perturbation, that is, the initial condition [Fig. 6(b)] for $(\xi_{1},\xi_{2},\overline{u}_{1},\overline{u}_{2})$ is taken to be

[TABLE]

where $\left(\Xi_{1},\Xi_{2},\overline{U}_{1},\overline{U}_{2}\right)$ is the solitary-wave solutions described in §§ 3.1 and the localized perturbation $\delta_{1}$ is a narrow SW packet with a narrow band of wavenumbers,

[TABLE]

with $A_{\varepsilon}=5\times 10^{-4}$, $X_{0}=84$, $x_{0}=2$, and $k_{0}=5$. For the interfacial solitary wave, $\Xi_{2}$, the amplitude is $-0.976$, the wavelength is $\sim 150$, and the speed is $C_{s}=0.0549$.

Below, the group velocity $v_{g}$, the phase velocity $v_{p}$, and the frequency $\nu$ correspond to the moving frame $(T,X)$. On the other hand, the group velocity $c_{g}=(v_{g}+C_{s})$, the phase velocity $c_{p}=(v_{p}+C_{s})$, and the frequency $\omega(=\nu+C_{s}k)$ correspond to the resting frame $(t,x)$. The variable with the tilde, $\widetilde{\cdot}$, stands for the numerical measurement of the corresponding quantity.

Initially, left-moving radiation is quickly emitted from the near field and eventually absorbed by our absorbing boundary condition in the buffer zones [dark stripe parallel to the yellow line in Fig. 6(a)]. After this initial transient, we can see that one SW packet is trapped in the near field [dark stripes parallel to the red line and blue line in Fig. 6(a)]. These trapped waves are all right-moving waves, that is, their phase velocities $\widetilde{v}_{p}>0$ are positive [dark stripes parallel to the black line in Fig. 6(c) and to the magenta line in Fig. 6(d)]. Thus, only the right-moving SWs that propagate in the same direction as the underlying interfacial solitary wave remain trapped in the near field.

We now study the spatiotemporal manifestation of these right-moving SWs in the near field. From Fig. 6, we can observe three features of these right-moving SWs:

(i) SWs become short in wavelength at the leading edge and long at the trailing edge. For $0\leq T\leq 4300$, the SW packets propagate towards the trailing edge with a relatively large wavenumber $\widetilde{k}\sim 23$ [Fig. 6(e)]. For $4300\leq T\leq 4700$, the SW packets propagate towards the leading edge with a relatively small wavenumber $\widetilde{k}\sim 1.8$ [Fig. 6(f)].

(ii) SW packets propagate towards the trailing edge with a relatively small group velocity, and towards the leading edge with a relatively large group velocity. From Figs. 6(a)(c)(e), we can see that for $0\leq T\leq 4300$, the SW packets propagate towards the trailing edge with a relatively small group velocity $\widetilde{v}_{g}\sim-0.045$. For $4300\leq T\leq 4700$ [Figs. 6(a)(d)(f)], the SW packets propagate towards the leading edge with a relatively large group velocity $\widetilde{v}_{g}\sim 0.37$. For $4800\leq T\leq 9100$ [Figs. 6(a)(g)], the SW packets again propagate towards the trailing edge with a relatively small group velocity.

(iii) SWs become high in amplitude at the leading edge and low at the trailing edge. For $0\leq T\leq 4300$, the SW packets’ amplitude increases to $\sim 3\times 10^{-3}$ at the leading edge and then SWs propagate towards the trailing edge [Fig. 6(e)]. For $4300\leq T\leq 4700$, the SW packets’ amplitude decreases to $\sim 1\times 10^{-4}$ at the trailing edge and then propagate towards the leading edge [Fig. 6(f)].

To understand the dynamical behavior of these right-moving SWs, we first quantify the dispersion relation of these waves. Figure 7(a) shows the logarithmic modulus, $\log_{10}|\widehat{\xi}_{1}(k,\nu)|^{2}$, with its magnitude color-coded, where $\widehat{\xi}_{1}(k,\nu)$ is the spatiotemporal Fourier transform of $\xi_{1}(X,T)$. For comparison, also shown are the pure linear dispersion relation $\nu_{k}=\omega_{k}-C_{s}k$ [Eq. (3.32), red dashed-dotted curve in Fig. 7(a)] and the modulated dispersion relation $\overline{\nu}_{k}=\overline{\omega}_{k}-C_{s}k$ [Eq. (3.33), black solid curve in Fig. 7(a)]. For the modulated dispersion relation $\overline{\nu}_{k}$, we take the amplitude of the interfacial solitary wave to be $(\Xi_{1},\Xi_{2},\overline{U}_{1},\overline{U}_{2})=(0.001,-0.976,0.027,-0.026)$. We can clearly see from Fig. 7(a) that, for $15\leq k\leq 25$, the modulated dispersion relation $\overline{\nu}_{k}$ can capture the peak locations of the spectrum well, whereas the pure linear dispersion relation $\nu_{k}$ deviates greatly. Therefore, these right-moving SWs can be well characterized by the modulated dispersion relation $\overline{\nu}_{k}$ and thereafter referred to as right-moving fast-mode SWs [see the definition of fast-mode waves in Sec. 3.2].

Incidentally, there are two yellow spots on the pure linear dispersion relation $\nu_{k}$ near $k=5$, as can be observed faintly in Fig. 7(a). However, the spectral power at these two yellow spots on the linear dispersion relation $\nu_{k}$ is six orders of magnitude lower than that at the modulated dispersion relation $\overline{\nu}_{k}$. These two blurry yellow spots correspond to the wave spectra of radiation waves in the far field, which is not the interest of this work.

To understand the asymmetric behavior (i)-(iii), we compare our numerical results of the TWN model with those of the effective linearized equations (1.46)-(1.49) in appendix A. Effective linearized equations (1.46)-(1.49) are obtained from the linearization of our TWN model (4.35)-(4.38) in the presence of the interfacial solitary wave. Their mathematical details are presented in appendix A. We also compare our TWN solutions with the theoretical predictions of ray-based theories. Due to the slow varying in space and time of the phase of fast-mode waves, the governing equations of space-time rays for the location $X$ and the wavenumber $k$ are given by [44, 6],

[TABLE]

where $\overline{\nu}_{k}$ is the modulated dispersion relation in Eq. (3.33). Note that the modulated dispersion relation $\overline{\nu}_{k}$ does not explicitly depend on time $T$ since the interfacial solitary wave is stationary in the moving frame. Clearly equations (5.44) and (5.45) constitute a Hamiltonian system with $\overline{\nu}_{k}$ as the Hamiltonian, $X$ displacement and $k$ momentum. Equation (5.44) states that the wave packet propagates at the group velocity.

We now quantify the asymmetric behavior of these right-moving fast-mode SWs by comparing the results of the TWN model, the results of the effective linearized equations, and theoretical predictions from ray-based theories:

(i) First, we provide a theoretical prediction for the temporal evolution of the wavenumber (equivalently the wavelength) of fast-mode waves. By the ray-based theory, when propagating with the initial perturbation (5.43), fast-mode waves possess the peak locations $X$ and wavenumbers $k$ between the two green-level curves in Fig. 7(b). For the minimal wavenumber, the theoretical prediction $k=1.7$ can be attained at $X=0$ on the red-level curve in Fig. 7(b). This theoretical minimal wavenumber is in good agreement with the measured one $\widetilde{k}\sim 1.8$ at $(T,X)=(4500,0)$ in Fig. 6(f). For the maximal wavenumber, the theoretical prediction ranges from $15.5$ to $26.7$ between the two green-level curves in Fig. 7(b). For the theoretical wavenumbers $k=15.5$ and $k=26.7$, the corresponding frequencies are $\nu=1.895$ and $\nu=1.46$, respectively. This range of theoretical wavenumbers and frequencies [depicted by the green rectangle in Fig. 7(a)] is in good agreement with the the range of the measured ones in the spectrum in Fig. 7(a). Therefore, the temporal evolution of the wavenumber can be characterized by the ray-based theory (5.44)-(5.45) for fast-mode waves. In particular, fast-mode SWs become short in wavelength at the leading edge [$X>0$] and long at the trailing edge [$X<0$] [see Fig. 7(b)].

(ii) Next, we investigate the motion of the peak location $X$ as a function of time $T$. Figure 7(c) displays the temporal evolution of numerically measured peak locations of the fast-mode waves for the TWN model (4.35)-(4.38) [green curve in Fig. 7(c)] as well as the prediction using the effective linearized equations (1.46)-(1.49) [blue curve in Fig. 7(c)]. For comparison, also displayed are the peak locations of the wave packets predicted by the ray-based theory (5.44)-(5.45) [red curve in Fig. 7(c)]. One can observe that there is excellent agreement between the numerical results and theoretical predictions for the motion of the peak locations. This confirms that the wave packet moves at the group velocity $v_{g}=\partial\overline{\nu}_{k}/\partial k$. As predicted by the ray-based theory, for $0\leq T\leq 4300$, the group velocity ${v}_{g}\sim-0.045$ is negative with a relatively small magnitude, whereas for $4300\leq T\leq 4700$, the group velocity ${v}_{g}\sim 0.37$ is positive with a relatively large magnitude. These two theoretical group velocities are in excellent agreement with the measured ones, $\widetilde{v}_{g}\sim-0.045$ and $\widetilde{v}_{g}\sim 0.37$, respectively. As reflected in the zig-zag pattern in Fig. 7(c), we can observe that SW packets propagate towards the trailing edge with a relatively small group velocity, and towards the leading edge with a relatively large group velocity.

(iii) Finally, we discuss the temporal evolution of the maximal amplitude, $A$, of fast-mode waves in the near field. Figure 7(d) displays the maximal amplitude of the fast-mode waves in our TWN model (4.35)-(4.38) [green curve] and that predicted using the effective linearized equations (1.46)-(1.49) [blue curve]. One can see very good agreement between them. In addition, one can observe from Fig. 7(d) that the amplitude $A$ is relatively large for the negative group velocity $v_{g}$ (interval marked by the red color), whereas the amplitude $A$ is relatively small for the positive group velocity $v_{g}$ (interval marked by the black color). Furthermore, for $4500\leq T\leq 9300$, the amplitude $A$ grows for positive $X$ (interval marked by the magenta color), whereas it decays for negative $X$ (interval marked by the cyan color). Therefore, SWs become high in amplitude at the leading edge ($X>0$) whereas low at the trailing edge ($X<0$).

To the best of our knowledge, the asymmetric behavior (i) was earlier discovered in references [6, 33], the asymmetric behavior (iii) was earlier discovered in the reference [37], but the asymmetric behavior (ii) was not reported in previous works. Here, we quantify these asymmetric types of behavior predicted by the ray-based theory for our TWN model when the initial perturbation (5.43) is a small-amplitude, narrow-width SW packet with a narrow band of wavenumbers.

6 Conclusions and discussion

Using the long-wavelength and small-amplitude approximations, we have proposed a two-layer, weakly nonlinear (TWN) model (2.21)-(2.24) describing the long-wave interactions between IWs and SWs. The TWN model captures the broadening of large-amplitude IWs that is a ubiquitous phenomenon in the ocean [40, 21]. In Sec. 5, we have investigated the spatiotemporal manifestation and the dynamical behavior of right-moving fast-mode SWs in the near field in the presence of an underlying IW. From our numerical results, the wavenumber, group velocity, and amplitude of fast-mode SW packets of our TWN model (4.35)-(4.38) can always be well captured by the predictions of the effective linearized equations (1.46)-(1.49) and the ray-based theory (5.44)-(5.45). The fast-mode waves behave as linear waves modulated by the underlying interfacial solitary wave. Importantly, the behavior of the right-moving fast-mode waves is asymmetric at the leading edge vs. the trailing edge when an underlying IW is present:

(i) *SWs become short in wavelength at the leading edge and long at the trailing edge [Fig. 7(b)]. *

(ii) SW packets propagate towards the trailing edge with a relatively small group velocity, and towards the leading edge with a relatively large group velocity [Fig. 7(c)].

(iii) *SWs become high in amplitude at the leading edge and low at the trailing edge [Fig. 7(d)]. *

In this work, we only focus on the spatiotemporal manifestation and dynamical behavior of SWs under a small-amplitude initial perturbation. As a natural extension of the above results, it is interesting to study the SWs when the amplitude of the perturbation is large, that is, the nonlinearity becomes prominent. In particular, it is important to understand how the nonlinearity and resonance affect the spatiotemporal manifestation and dynamical behavior of the right-moving fast-mode SWs in the presence of an underling IW.

Class 3 triad resonance is believed to be responsible for the surface signature of the underlying internal waves [38, 35, 29, 17]. The TWN model possesses two-mode waves, one slow and the other fast, and thus resonant interactions among different modes can occur during the energy exchange process [more details can be found in Appendix B]. The class 3 triad resonance condition [Eq. (2.50) in Appendix B], $c_{g}(k)=C_{p}(0)$, shows that, for resonantly-interacting waves, the group velocity of fast-mode waves $c_{g}(k)$ and the phase velocity of slow-mode waves $C_{p}(0)$ are equal [41, 27]. From many field observations, a narrow band of SWs with the resonant wavenumber was found to be located at the leading edge of an underlying IW and travel nearly at the same speed as the underlying IW [38, 35]. The surface phenomena may be related to both the triad resonance excitation and the three asymmetric types of behavior (i)-(iii). The allowance of triad resonance in the TWN model encourages us to investigate the spatiotemporal manifestation and dynamical behavior of SWs under large-amplitude initial perturbations in future work.

7 Acknowledgements

This work is supported by NYU Abu Dhabi Institute G1301, NSFC Grant No. 11671259, 11722107, and 91630208, and SJTU-UM Collaborative Research Program (D.Z.). We dedicate this paper to our late mentor David Cai.

Appendix A Effective linearized equations

In this section, we present the effective linearized equations of the TWN model (4.35)-(4.38). The variables $(\xi_{i},\overline{u}_{i})$ in the TWN model are composed of two components, one being the interfacial solitary wave $(\Xi_{i},\overline{U}_{i})$ and the other the perturbation of fast-mode waves $(\widehat{\xi}_{i},\widehat{\overline{u}}_{i})$. By substituting $(\xi_{i},\overline{u}_{i})=(\Xi_{i},\overline{U}_{i})+(\widehat{\xi}_{i},\widehat{\overline{u}}_{i})$ into the system (4.35)-(4.38) and collecting the linear terms with respect to $(\widehat{\xi}_{i},\widehat{\overline{u}}_{i})$, we can obtain the effective linearized equations for the fast-mode waves as follows,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Appendix B Class 3 triad resonance condition

In this section, we briefly verify the existence of solutions to the three-wave-resonance condition in the TWN model. If the dispersion relation allows the wavenumbers $k_{i}\ (i=1,2,3)$ and the corresponding frequencies $\omega_{k}(k_{1})$, $\omega_{k}(k_{2})$, and $\Omega_{k}(k_{1})$ to satisfy the conditions,

[TABLE]

these three waves constitute a class 3 resonant triad. Moreover, if the wavenumbers are specified as $k_{1}=k+\Delta k/2$, $k_{2}=k-\Delta k/2$, $k_{3}=\Delta k$, where $\Delta k\ll k$ and $\Delta k\rightarrow 0$, then the resonance condition reduces to

[TABLE]

where the group velocity $c_{g}$ and the phase velocity $C_{p}$ are given by the equations

[TABLE]

Here, $\omega_{k}$ corresponds to the dispersion relation of fast-mode waves in Eq. (3.32), and $\Omega_{k}$ corresponds to the dispersion relation of slow-mode waves in Eq. (3.31). Many early results [41, 27, 16, 17, 1] have confirmed that there exists a unique resonant wavenumber, denoted by $k_{\text{res}}$, satisfying Eq. (2.50) in the two-layer fluid system. For the TWN model, one can observe from Fig. 8 that there exists a unique resonant wavenumber $k_{\text{res}}$, satisfying the resonance condition (2.50), $c_{g}(k_{\text{res}})=C_{p}(k=0)$. Therefore, class 3 resonant triads exist among two fast-mode waves and one slow-mode wave for the TWN model.