Asymmetric behavior of surface waves induced by an underlying interfacial wave
Shixiao W. Jiang Gregor Kovačič Douglas Zhou
nlin.PS physics.flu-dyn 76B55, 76B07, 35L05, 65M22



          




         myheadings Asymmetric behavior of SWs induced by an underlying IWS.W. Jiang, G. Kovačič, and D. Zhou Asymmetric behavior of surface waves 

induced by an underlying interfacial waveReceived date, and accepted date (The correct dates will be entered by the editor).
    Shixiao W. JiangDepartment of Mathematics, the Pennsylvania State University, University Park, PA 16802-6400, USA, (suj235@psu.edu). 
Gregor KovačičMathematical Sciences Department, Rensselaer Polytechnic Institute, 110 8th Street, Troy, New York 12180, USA, (kovacg@rpi.edu).
Douglas ZhouSchool of Mathematical Sciences, MOE-LSC, and Institute of Natural Sciences, Shanghai Jiao Tong University, Shanghai 200240, China, (zdz@sjtu.edu.cn).
    Received: date / Accepted: date
===========================================================================================================================================================================================================================================================================================================================================================================================================================================================




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.
















76B55; 76B07; 35L05; 65M22






§ 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 <cit.>. 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
<cit.>. Generated by tidal flow, internal waves
usually can propagate thousands of kilometers from their source before
dissipation, sloping, and breaking extinguish them <cit.>.
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 <cit.>.












































































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 <cit.>.
The system can be formulated as a Hamiltonian system for a set of canonical variables <cit.>.
To evaluate the time derivatives of these canonical variables,
the high-order spectral (HOS) method <cit.> 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 <cit.>. 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 <cit.>.
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
<cit.>.
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 <cit.>.

[b] Ray-based theories: Many ray-based studies take a
statistical viewpoint of SWs modulated by a near-surface current induced by
IWs <cit.>.
These studies invoke phase-averaged models based on a wave-balance
equation and ray-based theory <cit.>, 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 <cit.>.
Class 3 triad resonance  is regarded as being responsible for the surface
signature of the underlying IWs <cit.>.
Based on a class 3 triad resonance condition,
many reduced models have been derived for the interfacial and surface waves
<cit.>
.
A detailed discussion of surface signature phenomena of IWs was presented in <cit.> 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 <cit.>.








































































































































































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. <ref>, we derive a reduced model for the two-layer fluid system.
In Sec. <ref>, we analyze the
basic properties of the model, including interfacial solitary-wave solutions
and  dispersion relations.
In Sec. <ref>, we present the numerical scheme and examine its numerical convergence.
In Sec. <ref>, 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. <ref>, and some mathematical details are presented in the
Appendix.







§ 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 ρ _1
and ρ _2, respectively. Here, ρ _2>ρ _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 ξ _2 between the two fluid layers, and their coupling with the
the overlaying free surface, ξ _1 [see Fig.
<ref>]. The velocity potential ϕ _i (i=1 for the
upper layer and i=2 for the lower layer) satisfies Laplace's equation,

    ϕ _ixx+ϕ _izz=0.

The kinematic equations for the continuity of the normal velocity at the
surface h_1+ξ _1, the interface ξ _2, and the flat topography -h_2 are given
in the form

    ξ _1t+ϕ _1xξ _1x=ϕ _1z,  at z=h_1+ξ _1,


    ξ _2t+ϕ _1xξ _2x=ϕ _1z,  at z=ξ _2,


    ξ _2t+ϕ _2xξ _2x=ϕ _2z,  at z=ξ _2,


    ϕ _2z=0,  at z=-h_2,

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,

    ϕ _1t+1/2( ϕ _1x^2+ϕ _1z^2) +gξ
    _1=0,  at z=h_1+ξ _1,


    ρ _1( ϕ _1t+1/2( ϕ _1x^2+ϕ
    _1z^2) +gξ _2)     =ρ _2( ϕ _2t+1/2
    ( ϕ _2x^2+ϕ _2z^2) +gξ _2) ,
    
    at   z    =ξ _2.

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,

    a/h_1=α≪ 1,   h_2/h_1=O(1).

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,

    h_1^2/L^2=β≪ 1,      h_2/h_1=O(1).

The two small parameters, α and β, control the nonlinear and
dispersive effects, respectively. Based on the scaling (<ref>)-(<ref>), we may nondimensionalize all
the physical variables by taking the original variables to be

    x=Lx^∗,  z=h_1z^∗,  t=(L/U_0)t^∗,


    (ϕ _1,ϕ _2)=( aLU_0/h_1)
    (ϕ _1^∗,ϕ _2^∗),  (ξ
    _1,ξ _2)=a(ξ _1^∗,ξ _2^∗),

where U_0=√(gh_1) is the characteristic speed of the gravity
waves. Here, all the variables with asterisks are assumed to be O(1) in α and β. In the dimensionless, starred variables, Laplace's
equation (<ref>) is formulated as

    βϕ _ix^∗x^∗^∗+ϕ _iz^∗z^∗^∗=0,    for  i=1,2.


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 ϕ _1^∗ in powers of β,

    ϕ _1^∗=ϕ _1^∗ (0)+βϕ _1^∗ (1)+O(β
    ^2),

which we use in the asymptotic analysis of the nondimensionalized problem of
equations (<ref>)-(<ref>) for small values of the
parameter β. From the leading O(1) term in Eqs. (<ref>) and (<ref>), ϕ _1^∗ (0) is found
to be independent of the height z,

    ϕ _1^∗ (0)=ϕ _1^∗ (0)(x^∗,t^∗).

The O(β ) terms in Eqs. (<ref>) and (<ref>) yield the equations

    ϕ _1x^∗x^∗^∗ (0)+ϕ _1z^∗z^∗^∗
    (1)=0,  αξ _2^∗<z^∗<h_1^∗+αξ _1^∗, with h_1^∗=1,


    ϕ _1z^∗^∗ (1)=ξ _1t^∗^∗+αϕ
    _1x^∗^∗ (0)ξ _1x^∗^∗,  at z^∗=h_1^∗+αξ _1^∗.

The expression for ϕ _1^∗ (1) is obtained as


    ϕ _1^∗ (1)=D_1ξ _1^∗Z_1-1/2ϕ _1x^∗x^∗^∗ (0)Z_1^2,

where Z_1=z^∗-h_1^∗-αξ _1^∗ and D_1ξ
_1^∗≡ξ _1t^∗^∗+αϕ _1x^∗^∗
(0)ξ _1x^∗^∗. Combining expressions (<ref>
) and (<ref>), to the first power of β, we can obtain
the solution ϕ _1^∗ as

    ϕ _1^∗=ϕ _1^∗ (0)+β( D_1ξ _1^∗Z_1-1/2ϕ _1x^∗x^∗^∗ (0)Z_1^2)
    +O(β ^2).

By integrating Eq. (<ref>) once from αξ
_2^∗ to h_1^∗+αξ _1^∗ with respect to z^∗, imposing the boundary conditions (<ref>) and (
<ref>), and substituting the expression (<ref>)
into equation (<ref>), we obtain the kinematic equation for the upper fluid
layer,

    η _1t^∗^∗+( η _1^∗u_1^∗) _x^∗-1/6β( h_1^∗) ^3
    u_1x^∗x^∗x^∗^∗-1/2β(
    h_1^∗) ^2ξ _1t^∗x^∗x^∗+O(αβ
    ,β ^2)=0,

where

    η _1^∗=h_1^∗+αξ _1^∗-αξ
    _2^∗,  u_1^∗=ϕ _1x^∗^∗
    (0).

Upon substitution of the velocity potential ϕ _1^∗ (<ref>) into the dynamical boundary condition (<ref>), we
obtain the dynamical equation governing the motion of the upper fluid layer,

    u_1t^∗^∗+αu_1^∗
    u_1x^∗^∗+gξ _1x^∗^∗+O(αβ
    ,β ^2)=0,

where the equation has been differentiated with respect to x^∗ once,
the terms in the first power of β are retained [the O(β) terms happen to
vanish in Eq. (<ref>)], and terms of O(αβ ,β ^2) are dropped. From the velocity potential ϕ
_1^∗, (<ref>), we obtain the horizontal velocity ϕ
_1x^∗^∗ as

    ϕ _1x^∗^∗=u_1^∗+βξ _1t^∗x^∗^∗Z_1-1/2βu_1x^∗x^∗^∗Z_1^2+O(αβ ,β ^2).

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

    u_1^∗=u_1^∗-1/2β
    h_1^∗ξ _1t^∗x^∗^∗-1/6β(
    h_1^∗) ^2u_1x^∗x^∗^∗+O(αβ ,β ^2);

the corresponding inverse is

    u_1^∗=u_1^∗+1/2β
    h_1^∗ξ _1t^∗x^∗^∗+1/6β(
    h_1^∗) ^2u_1x^∗x^∗^∗+O(αβ ,β ^2),

where the layer-mean horizontal velocity is defined as

    u
    _1^∗(x^∗,t^∗)=1η _1^∗
    ∫_αξ _2^∗^h_1^∗+αξ _1^∗ϕ
    _1x^∗^∗(x^∗,z^∗,t^∗)dz^∗.

After
substituting Eq. (<ref>) for the horizontal velocity u_1^∗, equations (<ref>) and (<ref>) 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 (ξ _1,ξ _2,
u_1,u_2), in the dimensional form, is

    η _1t+(η _1u_1)_x=0,  η _1=h_1+ξ
    _1-ξ _2,


    η _2t+(η _2u_2)_x=0,  η _2=h_2+ξ
    _2,


    u_1t+u_1u_1x+gξ _1x-1/3
    h_1^2u_1xxt+1/2h_1ξ _2ttx=0,


    u_2t+u_2u_2x+gξ _2x+ρ
    _rgη _1x
        -1/2ρ _rh_1^2u_1xxt+ρ _rh_1ξ
    _2ttx-1/3h_2^2u_2xxt=0,

where ρ _r is the density ratio ρ _1/ρ _2, and

    u_1(x,t)=1η _1∫_ξ _2^h_1+ξ
    _1ϕ _1x(x,z,t)dz,    u_2(x,t)=1η _2∫_-h_2^ξ _2ϕ
    _2x(x,z,t)dz,

are the layer-mean horizontal velocities.
We refer to the model (<ref>)-(<ref>) 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 <cit.>.
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,ρ_1,
ρ_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 √(gh_1)
, and the characteristic time is √(h_1/g) [see similar dimensionless forms used in
numerical simulations in <cit.>].
Note that, without loss of generality, our conclusions concerning broadening of IWs in Sec. <ref> and asymmetric behavior of SWs in Sec. <ref> also hold true for other parameter regimes of (h_1,h_2,g,ρ_1,
ρ_2).



§ 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.



 §.§ 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 (<ref>)-(<ref>). Then, we  use
these solitary-wave solutions as initial conditions for wave profiles and layer-mean velocities
for the subsequent numerical simulations in Sec. <ref>.
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, (ξ _i,u_i) [i=1,2], in the system (<ref>)-(<ref>):

    ξ _i( x,t) =Ξ _i( X) ,  u
    _i( x,t) =U_i( X) ,  
    X=x-C_st.

Substituting this ansatz into Eqs. (<ref>)-(<ref>) and integrating once with respect to X yields

    U_i=C_s( H_i-h_i) /H_i, with 
    H_1=h_1+Ξ _1-Ξ _2 and H_2=h_2+Ξ _2,

where we have assumed that H_i→ h_i as X→±∞, and h_1 (h_2) is the undisturbed thickness of the upper
(lower) fluid layer, respectively. Substituting the horizontal velocity (<ref>) for U_i into Eqs. (<ref>)-(<ref>) and integrating once with respect to X leads to the
equations

    -1/3C_s^2h_1^3H_1XX/H_1^2-1/2
    C_s^2h_1H_2XX   
    
    =C_s^2h_1^2/2( 1/H_1^2
    -1/h_1^2)     +g( H_1+H_2-h_1-h_2)
     -2/3C_s^2h_1^3H_1X^2/H_1^3,


    -1/2ρ _rC_s^2h_1^3H_1XX/H_1^2-ρ
    _rC_s^2h_1H_2XX-1/3C_s^2h_2^3H_2XX/
    H_2^2   
    
    =C_s^2h_2^2/2( 1/H_2^2-1/
    h_2^2) +g( H_2-h_2)
    +ρ _rg( H_1-h_1)     -ρ _rC_s^2h_1^3
    H_1X^2/H_1^3-2/3C_s^2h_2^3H_2X^2/
    H_2^3,

where we have assumed that H_iX, H_iXX→ 0 as X→±∞ for i=1,2. Since explicit solutions to Eqs. (<ref>)-(<ref>) are difficult to establish, we numerically
compute their solitary-wave solutions by applying the method in <cit.>.


In Fig. <ref>(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 <cit.>.
From Fig. <ref>(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 <cit.>. 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 
<cit.>, not only by the MCC-type models.










To quantify this broadening, we introduce a measure of the
effective width, Λ _s, for the interfacial solitary-wave solution <cit.>, defined as

    Λ _s=|1/Ξ _2(0)∫_0^∞Ξ
    _2(X)dX| .

Meanwhile, the effective width of MCC and KdV solutions are provided in the
reference <cit.>. Figure <ref>(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.










 §.§ Dispersion relations 




We now investigate the dispersion relation of the TWN model (<ref>)-(<ref>) using linear analysis. By substituting the
monochromatic solutions (ξ _i,u_i)∼exp [i (kx-μ
_kt)] into the system (<ref>)-(<ref>), the pure linear
dispersion relation in the absence of shear between the frequency μ _k
and the wavenumber k can be obtained as

    ( 1+ρ _rk^2h_1h_2+1/3k^2h_1^2+1/3
    k^2h_2^2+1/9k^4h_1^2h_2^2+1/12ρ
    _rk^4h_1^3h_2) μ _k^4
        -( gh_1+gh_2+1/3gk^2h_1h_2^2+1/3
    gk^2h_1^2h_2) k^2μ _k^2+( 1-ρ _r)
    g^2k^4h_1h_2=0.

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





Equation (<ref>) has 4 real roots  in the oceanic regime (the density
ratio ρ _r is close to 1). At the leading order in 1-ρ_r,  the dispersion
relations of the two-mode waves, denoted by Ω _k and ω _k [see Fig. <ref>],
can be approximated as

    Ω _k^2=( 1-ρ _r) g^2h_1h_2k^2/
    gh_1+gh_2+1/3gk^2h_1h_2^2+1/3gk^2h_1^2h_2
    ,

and

    ω _k^2=( gh_1+gh_2+1/3gk^2h_1h_2^2+
    1/3gk^2h_1^2h_2) k^2/1+ρ _rk^2h_1h_2+
    1/3k^2h_1^2+1/3k^2h_2^2+1/9
    k^4h_1^2h_2^2+1/12ρ _rk^4h_1^3h_2.

In the following, the two kinds of waves that correspond to the dispersion
relations Ω _k, (<ref>), and ω _k, (<ref>), are referred to as the slow-mode
waves and the fast-mode waves, respectively.














The modulated dispersion relation ω_k in the presence of
shear can be obtained by substituting (ξ _i,u
_i)∼ (Ξ _i,U_i)+exp [i(kx-ω_kt)]
into the system (<ref>)-(<ref>), where the shear is induced by an interfacial solitary wave
(Ξ _i,U_i) in above Sec. <ref>. The resulting equation is

    |[                 -ω_k+kU_1                 ω_k-k
    U_1          kh_1+kΞ _1-kΞ _2                         0;                         0                 -ω_k+kU_2                         0                kh_2+kΞ _2;                        gk             -1/2kh_1ω_k^2 -ω_k-
    1/3k^2h_1^2ω_k+kU_1                         0;                    ρ _rgk  (1-ρ _r)gk-ρ _rkh_1ω_k^2      -
    1/2ρ _rk^2h_1^2ω_k  -ω_k-1/3k^2h_2^2ω_k+kU_2 ]
    | =0,

where  |·| denotes the determinant of the enclosed matrix.
In the following Secs. <ref> and <ref>, 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_st. Note that the solitary-wave solutions (Ξ _i,U_i) are
steady in time T in this moving frame.
Then, the modulated dispersion
relation ν, corresponding to  the moving frame T and X,  is given by

    ν_k=ω_k-C_sk,

where
ω_k is the modulated dispersion relation in Eq. (<ref>) corresponding to the resting frame t and x.
Note that the dispersion relation ν_k is independent of time T since (Ξ _i,U_i) are steady in time T.
Moreover, the
wavelengths of (Ξ _i,U_i) are relatively
long with respect to the characteristic wavelengths of fast-mode waves, so (Ξ _i,U
_i) in the dispersion relation ν_k can be
locally treated as constant in X space.









































§ NUMERICAL SCHEME







For numerical computations, we cast Eqs. (<ref>)-(<ref>) in the conservation form in the right-moving frame with the
nonlinear phase velocity (soliton speed) C_s as follows:

    η _1T-C_sη _1X+(η _1u_1)_X=0,


    η _2T-C_sη _2X+(η _2u_2)_X=0,



    M_1T-C_sM_1X+( 1/2u_1^2+gξ _1)
    _X=0,


    M_2T-C_sM_2X+( 1/2u_2^2+gξ _2+ρ
    _rgη _1) _X=0,

where

    M_1=u_1-1/3h_1^2u_1XX-1/2
    h_1( η _2u_2) _XX,


    M_2=u_2-1/2ρ _rh_1^2u_1XX-
    1/3h_2^2u_2XX-ρ _rh_1( η _2u
    _2) _XX,

and T=t, X=x-C_st.
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 <cit.>.
The Kelvin-Helmholtz (KH) instability is suppressed by applying a low-pass filter 
<cit.>. (The wavenumbers for the KH instability are much larger than the wavenumbers of the trapped right-moving SWs as introduced in the following Sec. <ref>. Thus, these KH unstable wavenumbers are physically irrelevant in our computations.) In our simulations, we fix the parameter regime (ρ _1,ρ
_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 <ref>
shows the spatiotemporal evolution of the IWs' profile, ξ _2, for 0≤ T≤ 10000. We can see from Fig. <ref> 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 <cit.>.































Next, we examine the numerical convergence in time and space of our scheme for initially unperturbed
interfacial solitary-wave solutions [Eqs. (<ref>)-(<ref>)]. We compute the l_1
norm error for the SWs' profile defined as

    E_l_1=∑_j|ξ _1(X_j)-ξ
    _1^ ref(X_j)|Δ X,

where the reference solution ξ _1^ 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. <ref> that the scheme has fourth-order accuracy in
time and second-order accuracy in space.















§ ASYMMETRIC BEHAVIOR OF SWS IN THE PRESENCE OF AN UNDERLYING IW








































In this section, we present our numerical results for the system (<ref>)-(<ref>) 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 ξ _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. <ref>(b)]
for (ξ _1,ξ _2,u_1,u_2) is taken to be

    (Ξ _1+δ _1,Ξ _2,U_1,U_2),

where
(
Ξ _1,Ξ _2,U_1,U_2) is the
solitary-wave solutions described in  <ref>
and
the localized perturbation δ_1 is a narrow SW packet with a narrow band of
wavenumbers,

    δ _1=A_ε[tanh (X-X_0+x_0)-tanh
    (X-X_0-x_0)]cos (k_0(X-X_0)),

with A_ε=5× 10^-4, X_0=84, x_0=2, and k_0=5. For the interfacial solitary wave, Ξ _2, the
amplitude is -0.976, the
wavelength is ∼ 150, and the speed is C_s=0.0549.


Below, the group velocity v_g, the phase velocity v_p, and the
frequency ν 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 ω (=ν+C_sk) correspond to the resting frame (t,x).
The variable with the tilde, ·, 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. <ref>(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. <ref>(a)].
These trapped waves are all
right-moving waves, that is, their phase velocities v_p>0 are positive [dark stripes parallel to the black line in Fig. <ref>(c)
and to the magenta line in Fig. <ref>(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. <ref>, 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≤ T≤ 4300, the SW packets propagate towards the trailing edge with a relatively large wavenumber
k∼ 23 [Fig. <ref>(e)].
For 4300≤ T≤ 4700, the
SW packets propagate towards the leading edge with a relatively small wavenumber k∼ 1.8 [Fig. <ref>(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. <ref>(a)(c)(e), we can
see that for 0≤ T≤ 4300, the SW packets
propagate towards the trailing edge with a relatively small group
velocity  v_g∼ -0.045.
For 4300≤ T≤ 4700 [Figs. <ref>(a)(d)(f)], the SW packets propagate towards the leading edge
with a relatively large group velocity v_g∼ 0.37.
For 4800≤ T≤ 9100 [Figs. <ref>(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≤ T≤ 4300, the SW packets' amplitude increases to ∼
3× 10^-3 at the leading edge and then SWs propagate towards the trailing edge
[Fig. <ref>(e)]. For 4300≤ T≤ 4700, the
SW packets' amplitude decreases to ∼ 1× 10^-4 at the trailing edge
and then propagate  towards the leading edge
[Fig. <ref>(f)].























































To understand the dynamical
behavior of these right-moving SWs, we first quantify the dispersion relation of these waves.
Figure <ref>(a) shows the logarithmic modulus, log _10|ξ_1(k,ν )|^2, with its magnitude color-coded,
where ξ_1(k,ν ) is the spatiotemporal Fourier transform
of ξ _1(X,T). For comparison, also shown are the pure linear dispersion relation ν _k=ω _k-C_sk  [Eq. (<ref>), red dashed-dotted
curve in Fig. <ref>(a)]
and the modulated dispersion relation ν_k=
ω_k-C_sk [Eq. (<ref>), black solid curve in Fig. <ref>(a)].
For the modulated dispersion relation ν_k, we take
the amplitude of the interfacial solitary wave  to be (Ξ _1,Ξ _2,U_1,U
_2)=(0.001,-0.976,0.027,-0.026).
We can clearly see from
Fig. <ref>(a) that, for 15≤ k≤
25, the modulated dispersion relation ν_k can
capture the peak locations of the spectrum well, whereas the pure linear dispersion relation ν _k deviates greatly.
Therefore, these right-moving SWs
can be well characterized by the modulated dispersion relation ν_k and thereafter referred to as right-moving fast-mode SWs
[see the definition of fast-mode waves in Sec. <ref>].

Incidentally, there are two yellow spots on the pure linear dispersion
relation ν _k near k=5, as can be observed faintly in Fig.
<ref>(a). However,
the spectral power at these two yellow spots on the linear dispersion
relation ν _k is six orders of
magnitude lower than that at
the modulated dispersion relation ν_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 (<ref>)-(<ref>) in appendix <ref>.
Effective linearized equations (<ref>)-(<ref>) are obtained from the linearization of our TWN model (<ref>)-(<ref>) in the presence of the interfacial solitary wave. Their
mathematical details are presented in appendix <ref>.
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 <cit.>,

    dX/dT   =∂ν_k/∂ k,
    
    dk/dT   =-∂ν_k/∂ X,

where ν_k is the modulated dispersion relation  in
Eq. (<ref>).
Note that the modulated dispersion relation ν_k does not
explicitly depend on time T since the interfacial solitary wave  is
stationary in the moving frame. Clearly equations (<ref>) and (<ref>) constitute a Hamiltonian system with ν
_k as the Hamiltonian, X displacement and k momentum.
Equation (<ref>) 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 (<ref>), fast-mode waves
possess the peak locations X and wavenumbers k between the two green-level curves in Fig. <ref>(b).
For the minimal
wavenumber, the theoretical prediction k=1.7 can be attained at X=0
on the red-level curve  in Fig. <ref>(b).
This theoretical minimal
wavenumber is in good agreement with the measured one
k∼ 1.8 at (T,X)=(4500,0) in Fig. <ref>(f).
For the maximal wavenumber, the theoretical prediction ranges from 15.5 to 26.7 between
the two green-level curves in Fig. <ref>(b).
For the theoretical wavenumbers k=15.5 and k=26.7, the corresponding frequencies
are ν=1.895 and ν = 1.46, respectively.
This range of theoretical wavenumbers and frequencies [depicted by the green
rectangle in Fig. <ref>(a)] is in good agreement
with the the range of the measured ones in the spectrum in Fig. <ref>(a).
Therefore, the temporal evolution of the wavenumber can be characterized by the
ray-based theory (<ref>)-(<ref>) 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. <ref>(b)].













(ii) Next, we investigate the motion of the peak location X as a function of time T.
Figure <ref>(c) displays the temporal evolution of
numerically measured peak locations of the fast-mode waves for the TWN model (<ref>)-(<ref>) [green curve in Fig. <ref>(c)] as well as the prediction using the effective linearized equations (<ref>)-(<ref>) [blue curve in Fig. <ref>(c)].
For comparison, also displayed are the peak locations of
the wave packets predicted by the ray-based theory (<ref>)-(<ref>) [red curve in Fig. <ref>(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=∂ν
_k/∂ k.
As predicted by the ray-based theory, for 0≤ T≤ 4300, the group velocity
v_g∼ -0.045 is negative with a relatively small magnitude, whereas for 4300≤ T≤ 4700, the group velocity v_g∼ 0.37 is positive with a relatively large magnitude.
These two theoretical group velocities are in excellent agreement with the
measured ones, v_g∼ -0.045 and
v_g∼ 0.37, respectively.
As reflected in the zig-zag pattern in Fig. <ref>(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 <ref>(d) displays the maximal amplitude of
the fast-mode waves in our TWN model (<ref>)-(<ref>) [green curve] and that predicted using the effective
linearized equations (<ref>)-(<ref>) [blue curve]. One
can see very good agreement between them.
In addition, one
can observe from Fig. <ref>(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≤ T≤ 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 
<cit.>, the asymmetric behavior (iii) was earlier discovered in the reference <cit.>, 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 (<ref>)
is a small-amplitude, narrow-width SW packet with a narrow band of wavenumbers.










§ CONCLUSIONS AND DISCUSSION




























Using the long-wavelength  and
small-amplitude approximations, we have proposed a two-layer, weakly nonlinear (TWN) model (<ref>)-(<ref>) 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 <cit.>.
In Sec. <ref>, 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
(<ref>)-(<ref>) can always be well captured by
the predictions of the effective linearized equations
(<ref>)-(<ref>) and the ray-based theory
(<ref>)-(<ref>). 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. <ref>(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. <ref>(c)].

(iii) SWs become high in amplitude at the
leading edge and low at the trailing edge [Fig. <ref>(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 <cit.>. 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 <ref>].
The class 3 triad resonance
condition [Eq. (<ref>) in Appendix <ref>], 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 <cit.>.
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 <cit.>.
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.
















































§ 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.








§ EFFECTIVE LINEARIZED EQUATIONS




In this section, we present the effective linearized equations of
the TWN model (<ref>)-(<ref>).  The variables (ξ _i,u_i) in the TWN model are composed of two components,
one being the interfacial solitary wave (Ξ _i,U_i) and the other the perturbation of fast-mode
waves (ξ_i,u_i). By substituting (ξ _i,u_i)=(Ξ _i,U_i)+(ξ_i,
u_i) into the system (<ref>)-(<ref>) and collecting the linear terms with respect to (
ξ_i,u_i), we can obtain the effective
linearized equations for the fast-mode waves as follows,

    η_1T-C_sη_1X+[H_1u
    _1+U_1(ξ_1-ξ_2)]_X=0,


    η_2T-C_sη_2X+(H_2u
    _2+U_2ξ_2)_X=0,


    M_1T-C_sM_1X+( U_1
    u_1+gξ_1) _X=0,


    M_2T-C_sM_2X+( U_2
    u_2+gξ_2+ρ _rgη_1)
    _X=0,

where

    H_1=h_1+Ξ _1-Ξ _2,  H_2=h_2+Ξ _2,


    η_1=h_1+ξ_1-ξ_2,  
    η_2=h_2+ξ_2,


    M_1=u_1-1/3h_1^2
    u_1XX-1/2h_1( H_2u_2+
    U_2ξ_2) _XX,


    M_2=u_2-1/2ρ _rh_1^2
    u_1XX-1/3h_2^2u
    _2XX-ρ _rh_1( H_2u_2+U_2
    ξ_2) _XX.



































































§ 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 ω _k(k_1), ω _k(k_2), and Ω _k(k_1)  to satisfy the
conditions,

    k_1   -k_2 =k_3, 
    ω _k(k_1)    - ω _k(k_2) =Ω _k(k_3),

these three waves constitute a class 3 resonant triad. Moreover, if the wavenumbers
are specified as  k_1=k+Δ k/2, k_2=k-Δ k/2, k_3=Δ k, where Δ k≪ k and Δ k→ 0, then the
resonance condition reduces to

    c_g(k)=C_p(0),

where the group velocity c_g and the phase velocity C_p are given by
the equations

    c_g(k)≡∂ω _k( k) /∂ k,   C_p(0)≡. Ω _k(Δ k)/Δ k|
    _Δ k→ 0.

Here, ω _k corresponds to the dispersion relation of fast-mode waves
in Eq. (<ref>), and
Ω _k corresponds to the dispersion relation of slow-mode waves
in Eq. (<ref>).
Many early results
<cit.> have confirmed that there
exists a unique resonant wavenumber, denoted by k_res, satisfying Eq. (<ref>) in
the two-layer fluid system.
For the TWN model, one can observe from Fig. <ref> that
there exists a unique resonant wavenumber k_res, satisfying the resonance
condition (<ref>), c_g(k_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.








plain




