Nonlinear Traveling Internal Waves in Depth-Varying Currents

# Nonlinear Traveling Internal Waves in Depth-Varying Currents

K. L. Oliveras Mathematics Department, Seattle University, Seattle, WA 98122  and  C. W. Curtis Department of Mathematics and Statistics, San Diego State University, San Diego, CA 92182
###### Abstract.

In this work, we study the nonlinear traveling waves in density stratified fluids with depth varying shear currents. Beginning the formulation of the water-wave problem due to [1], we extend the work of [4] and [18] to examine the interface between two fluids of differing densities and varying linear shear. We derive as systems of equations depending only on variables at the interface, and numerically solve for periodic traveling wave solutions using numerical continuation. Here we consider only branches which bifurcate from solutions where there is no slip in the tangential velocity at the interface for the trivial flow. The spectral stability of these solutions is then determined using a numerical Fourier-Floquet technique. We find that the strength of the linear shear in each fluid impacts the stability of the corresponding traveling wave solutions. Specifically, opposing shears may amplify or suppress instabilities.

## 1. Introduction

As noted in [3] and [19], internal waves are known to be a common feature in coastal oceanic flows. Resulting from instabilities along pynoclines and thermoclines, internal waves with amplitudes from 5-100 m [3, 19, 26] have been observed. While several generation mechanisms have been proposed and studied, it appears that in most cases internal waves are generated in coastal regions due to the presence of tide induced shearing forces in stratified fluids moving over varying bathymetry [15, 19]. What is interesting though is that while depth varying shear currents should be anticipated in such environments, most modeling attempts have traditionally ignored the impact of vorticity on internal wave dynamics. This work will address this shortcoming in part by studying the impact of depth varying shear currents on the existence and flow properties of internal traveling wave solutions (TWSs) in density stratified environments.

While the impact of vorticity on internal waves is not well understood, a great deal of work on the influence of shear in free surface flows has appeared. The simplest case of free surfaces moving over depth varying currents is when the vorticity is assumed to be a constant throughout the fluid. The literature on this subject is too vast to address here, though we point out the seminal studies [10, 30, 33] and the extensive bibliography in [6]. We also note that a key part of the approach we take in this paper follows the work in [4], which in part makes use of the Unified Transform Method (UTM) pioneered in [16].

Perhaps surprisingly, there are still outstanding and significant issues even for constant vorticity flows. In particular, the impact of vorticity on the stability of traveling wave solutions is not fully understood. As shown in [24], the impact of vorticity modifies the amplitudes of Benjamin–Feir (BF) instabilities depending on the sign of the vorticity. Similar results have analytically been found in approximate models [5, 20, 32]. In addition, constant vorticity introduces “high-frequency” instabilities similar to those observed in irrotational waves [22, 17, 13]. These instabilities may be more dominant than the BF instability depending on the amplitude of the TWS and the strength of the linear shear [24]. Furthermore, for a TWS with a region of lower pressure trapped under the wave [10, 34], new instabilities, perhaps related to a Rayleigh-Taylor instability, are visible in the spectrum. As shown in [34], the existence of these low pressure regions is directly related to the existence of stagnation points within the fluid. While usually avoided in much of the existing literature due to the analytic difficulties they introduce, it is becoming clearer that stagnation points can have significant impacts on wave shape and stability. Thus, in part due to stagnation point induced low pressure regions, and in part due to high-frequency instabilities, no TWSs have been found which are stable to all perturbations despite the ability of a constant shear current to suppress the onset of BF instabilities.

Given these results for constant vorticity, it is interesting then to ask how depth variation in the vorticity affects the existence of these instabilities, and whether more complicated shear profiles will support the existence of stable internal or even free surface traveling waves. The case of non-constant vorticity flows has received far less attention. For the case of non-steady waves, this was first examined in [11] for a stratified fluid with two differing densities and and two differing constant vortices and . We refer to such fluids as bistratified; see Figure 1 for reference. More recently, such systems have been looked at experimentally and numerically in [31, 21]. Assuming a rigid lid and zero vorticity, i.e. , in the lower fluid region, nonlinear, large amplitude traveling interfacial waves for bistratified fluids were studied in [27, 29]. A complete characterization of the linear stability of the small amplitude limit of these waves was presented in [28]. Generalizing this work, the evolution of a free surface and a free internal layer for a bistratified fluid in a shallow-water approximation was studied in [9]. Significant nonlinear phenomena like internal dispersive shock waves as a result of bistratification was reported for the first time.

Focusing then on the case of single-valued internal layers, by assuming a rigid lid at the fluid-air interface and an impermeable sea bed, using the UTM techniques in [18, 4] and incorporating techniques developed for traveling wave solutions as found in [25], we derive a closed system describing the evolution of the free internal surface in a bistratified fluid. Using this formulation, we can analytically determine at which wave speeds bifurcations from flat states will occur. By developing a Stokes expansion at this bifurcation point, we can analytically determine the behavior of the bifurcation curve in terms of wave height and speed in a small amplitude limit, thereby allowing for a better understanding of when stagnation points in the fluid velocity appear and when turning points in bifurcation curves appear.

We then use numerical continuation to generate bifurcation curves which provide the speed-amplitude relationship for the traveling nonlinear interfacial waves. Of particular interest is the role different choices of vorticities play in determining the number of stagnation points within the fluid regions, and how these stagnation points influence the shape and properties of the interfacial TWSs. As we show, the position and number of stagnation points determines whether waves are elevated, depressed, symmetric in nature, or especially peaked. To complement these results, we also present plots of the streamlines to help illustrate the different impacts on dynamics that the stagnation points have throughout the bulk of the fluid as well as on the interface.

Finally, we study the stability of various TWSs. One of the most striking features associated with bistratification is the complicated relationship between the appearance of modulational instabilities (MIs) and the choices of the shear strength and depth ratio of the layers. As we show, by making the ratio of the depths of the two layers large enough, MIs appear to be completely suppressed. However, we also show results in which MIs are introduced and suppressed by shear, thereby illustrating the complicated balances mediated by nonlinearity between these parameter choices. These results generalize and expand on those for constant vorticity in single layer fluids [32, 24], in which vorticity is shown to strongly influence whether MIs exist.

However, at least for the solutions examined in this paper, all traveling solutions are unstable with respect to higher-frequency perturbations. Given the more complicated nature of the physical problem studied in this paper, it is a non-trivial question to determine if all traveling wave solutions are unstable, especially given the stabilizing mechanisms we illustrate via our numerical result. Likewise, the complex nature of instabilities for multi-valued interfaces remain unknown and will be explored in future work.

The structure of the paper is as follows. The derivation of our model is given in Section 2. The computation of the TWSs and analytic derivations of the nonlinear impact on wave speed and the presence of stagnation points are given in Section 3. Section 4 explains the means by which we numerically determine the stability of the TWSs computed in this paper, and Section 5 contains the results of these stability computations. Concluding remarks appear in Section 6, and an Appendix collections technical details used throughout the paper.

## 2. Boundary Conditions and AFM Formulation

We assume the flow is incompressible in each region in a frame of reference moving to the right with speed . The resulting equations must satisfy Euler’s equations in the bulk:

 ∇⋅uj=0, (1) ∂tuj+(uj⋅∇)uj=−1ρj∇pj−g^z, (2)

where and corresponds to the fluid velocities in the top and bottom layer of the fluid respectively, is the pressure at any point within each fluid, is the density of each fluid, and is the constant acceleration due to gravity.

At the interface between each fluid, we require continuity in pressure, as well as continuity in the normal velocities of each fluid at the interface. This provides the boundary conditions

 p1=p2, and u1⋅N=u2⋅Nat z=η(x,t), (3)

where represents a unit normal vector to the interface defined by . Furthermore, the interface is evolved by the fluid velocities normal to the interface. That is,

 ηt=uj⋅N,

for or as they are equivalent. Finally, the rigid lid and solid boundary conditions require

 v1=0 at z=h1,v2=0 at z=−h2. (4)

### 2.1. A Non-local Formulation

Equations (1-4) represent the equations we wish to solve provided that there is a constant linear shear within each fluid domain. That is, that

 ∂xwj−∂zuj=ωj

within each fluid where we have chosen a particular sign convention for the vorticity consistent with conventions used elsewhere in the literature (see [10, 21] for example). A disadvantage of the above formulation is that it requires one to solve for the velocities and pressure within the bulk of each fluid as well as the interface separating the two fluids. In our work, we choose to cast the equations of motion using a nonlocal reformulation that allows us to work completely within the context of variable defined at the interface. Here we following the work of [4, 18] in order to recast the problem into variables defined at the interface. However, to take advantage of this formulation, we first pose (1-4) in terms of new bulk variables.

Let for represent the stream functions in the top and bottom layer of the fluid respectively. Furthermore, we introduce the functions for to represent the contributions to the velocity from potential flow. Thus, we can write the fluid velocities as

 (uj−cwj)=(∂xϕj−ωjz∂zϕj)=(∂zψj−∂xψj).

Using (1)-(4), the stream- and potential-functions satisfy

 Δϕj=0, (x,z)∈Dj (5) Δψj=−ωj (x,z)∈Dj, (6) ∂tϕj+ωjψj+12|∇ϕj−ωjz^x|2+pjρj+gz=0, (x,z)∈Dj (7)

within the bulk of each fluid.

###### Remark 1.

In integrating (2) within the fluid domain to find (7), we arrive at an equation of the form

 ∂tϕj+ωjψj+12|∇ϕ−ωjz^x|2+pjρj+gz=Bj(t)

where we call the time varying function the Bernoulli function. However, by introducing the change of variables

 ϕj=~ϕj+∫t0Bj(s)ds,

we can eliminate the Bernoulli function, and thus we set it to zero throughout the remainder of this text.

#### 2.1.1. Boundary Conditions

Examining the conditions prescribed at the interface given by (3), we find

 ∂zϕj−(∂xϕj−ωjη)∂xη=∂tη.

Likewise, we can readily show that

 ddx(ψj(x,η(x,t),t))=−∂tη,

so that

 ψj(x,η(x,t),t)=−∂−1x∂tη

where we have picked the arbitrary constant for the domain. Invoking the pressure condition that at the free interface, we then find

 ∂tϕ1−ω1∂−1x∂tη+12((∂xϕ1−ω1η)2+(∂zϕ1)2)+gη=ρ(∂tϕ2−ω2∂−1x∂tη+12((∂xϕ2−ω2η)2+(∂zϕ2)2)+gη)

where . Finally, at the rigid lid and the solid bottom boundary, we have the zero flux conditions which require

 ∂zϕ1=0 at z=h1, and  ∂zϕ2=0 at −h2.

### 2.2. Nonlocal Formulation in Interface Quantities

Following the work of [1, 4, 18], we introduce the harmonic functions such that

 φj=e−ikxcosh(k(z+(−1)jhj)).

Noting that

 (∂zφj)(Δψj+ωj)+∂zψj(Δφj)=0,

where solves (6), we integrate the above quantity over the respective fluid domain. By using Green’s second identity and integration by parts, we can readily show that for .

 ∫2πL0e−ikx((ik(∂tη)−ωj)Cj−k(∂xqj−ωjη)Sj)dx=0,∀k∈Λ (8)

where , is given by

 qj(x,t)=ϕj(x,η(x,t),t),

and the quantities , are defined by

 Cj=cosh(k(η+(−1)jhj)),Sj=sinh(k(η+(−1)jhj)).

Using the definition for , we can rewrite both and both evaluated at in terms of our new variables and as

 ∂xϕj−ωjη=∂xqj−ωjη−∂xη∂tη1+(∂xη)2,∂zϕj=∂tη+∂xη(∂xqj−ωjη)1+(∂xη)2.

Similarly, we can also write as

 ∂tϕj=∂tqj−∂tη(∂tη+∂xη(∂xqj−ωjη))1+(∂xη)2

Using the boundary condition

 p1(x,η,t)=p2(x,η,t),

along with our new variables for , we have

 ∂tq1−ω1∂−1x∂tη+12(∂xq1−ω1η)2−12(∂tη+∂xη(∂xq1−ω1η))21+(∂xη)2+gη=ρ(∂tq2−ω2∂−1x∂tη+12(∂xq2−ω2η)2−12(∂tη+∂xη(∂xq2−ω2η))21+(∂xη)2+gη). (9)

Thus, (8) for along with (9) represent a closed system of three equations for the three unknown quantities .

### 2.3. Non-Dimensionalization

Introducing the non-dimensional parameters given by

 ~x=xL,~z=zh1,~t=c0Lt,qj=ϵc0L~qj,η=a~η,~k=Lk,

and

 ~ωj=ωjh1c0,~c=cc0,μ=h1L,ϵ=ah1,δ=h2h1,

where , we have the following system of three equations

 ∂t(q1−ω1∂−1xη)+ϵ2(∂xq1−ω1η)2+η−ϵμ22(∂tη+(ϵ∂xη)(∂xq1−ω1η))21+(ϵμ∂xη)2=ρ(∂t(q2−ω2∂−1xη)+ϵ2(∂xq2−ω2η)2+η−ϵμ22(∂tη+(ϵ∂xη)(∂xq2−ω2η))21+(ϵμ∂xη)2), (10)
 ∫2π0e−ikx((ikϵμ2∂tη−ω1)C1−ϵμk(∂xq1−ω1η)S1)dx=0,∀k∈Z/{0}, (11)

and

 ∫2π0e−ikx((ikϵμ2∂tη−ω2)C2−ϵμk(∂xq2−ω2η)S2)dx=0,∀k∈Z/{0}, (12)

## 3. Traveling Wave Equations

Working with the nondimensional equations in a traveling coordinate frame, letting as well as for , the equations for stationary solutions in a traveling frame become

 (12(ϵ(∂xq1−ω1η)−c)21+(ϵμ∂xη)2+ϵη−12c2)=ρ(12(ϵ(∂xq2−ω2η)−c)21+(ϵμ∂xη)2+ϵη−12c2)
 ∫2π0e−ikx(C1ω1+μkS1(ϵ(∂xq1−ω1η)−c))dx=0.
 ∫2π0e−ikx(C2ω2+μkS2(ϵ(∂xq2−ω2η)−c))dx=0.

where and with and . In the above fomulation, we are required to solve three equations for three unknown functions , and . However, as in [13, 4], we can reduce the dimensionality of the problem by formally solving the Bernoulli equation for either the quantity or . Without loss of generality, we solve for the tangential velocity approaching the interface from above in terms of to find

 ϵ(∂xq1−ω1η)−c=σ√ρ(ϵ∂xq2−ϵω2η−c)2+(c2−2ϵη)(1+(ϵμ∂xη)2)(1−ρ) (13)

where .

###### Remark 2.

The meaning of the plus/minus sign is typically of critical importance in problems involving linear shear (see for example [34]). However, when simply solving for the interface variables and , an alternative choice in is equivalent to changing the sign of vorticity in the upper layer given by . Furthermore, to ensure that remains a solution to the problem, that is, that we consider traveling wave solutions that bifurcate from flat water with no jump in the tangential velocity, we choose as indicated by (13) where is the speed corresponding to bifurcation from the trivial solution.

The resulting equations to solve for and are given by

 ∫2π0e−ikx(C1ω1+μkS1(σ√ρ(ϵ∂xq2−ϵω2η−c)2+(c2−2ϵη)(1+(ϵμ∂xη)2)(1−ρ)))dx=0, (14)
 ∫2π0e−ikx(C2ω2+μkS2(ϵ(∂xq2−ω2η)−c))dx=0,, (15)

where and are the same as given before This system of two equations can be solved for the unknowns and for an admissible value of the wave-speed . Once a solution and are determined, the is obtained directly by substitution in (13).

### 3.1. Bifurcation from Trivial Flow

Since we ultimately use a numerical continuation scheme in this work, we determine for which wave speeds a non-trivial solution will bifurcation from the zero solution. This is achieved by expanding the free surface , and in powers of of the form

 η(x)=∞∑n=0ϵnηn(x),∂xq2=∞∑n=0ϵn∂xvn,c=∞∑n=0ϵncn.

We assume that solutions will bifurcation from a wave with positive velocity so that and thus, choosing enforces bifurcating from a branch where there is no jump in the tangential velocity. Expanding (14) and (15) in , we find the following linear system for and given by

 (16)

where represents the -th Fourier coefficient of defined by

 Fk{f(x)}=12π∫2π0e−ikxf(x)dx,

and we have introduced the quantity

Since the right-hand side of (16) evaluates to zero, we can then say that the only way for there to be a non-zero solution for and is for the linear operator on the left-hand side of (16) singular for at least one value (henceforth referred to as ). This leads to the equation given by

 c20−~ωTk0c0−Tk0(ρ−1)=0,

where

 Tk=tanh(μk)tanh(δμk)μk(ρtanh(μk)+tanh(δμk)). (17)

Using this notation, we find

 c0=12(−Tk0~ω+√T2k0~ω2+4(ρ−1)Tk0), (18)

where the negative solution is extraneous due to our earlier assumptions of bifurcating from a positive speed .

The corresponding nontrivial solutions are given by

 η0(x)=cos(k0x),v0(x)=μk0c0tanh(μδk0)cos(k0x). (19)

Thus, for a given wave-number , the leading-order approximations to the solution set are given by the expressions in (18) and (19). With the aid of a CAS, we determine higher-order corrections both to the speed whereby

 c≈c0+ϵ2c2+O(ϵ4).

Due to the complicated nature of the resulting expressions though, we obmit the details of computing here. By finding though, we can generate analytic results which provide greater insight into the impact of vorticity on interfacial wave propagation. These results will also help us identify potential regions of difficulty in the numerical results presented later in this paper.

We present our results for computing in Figure 2. As seen in the figure, the response of to the choice of vorticities can be quite sensitive. For example, choosing a weak stratification value of , reflective of oceanic density stratification, , , and choosing the branch, we plot in Figure 2 for and . While both plots exhibit a wide range of scales with respect to the magnitude of , this range is far wider for than . Likewise, the region over which is negative is far larger in the case for . Note, those regions for which are regions in which the wave speed decreases with increasing wave height, which is in contrast to the usual result whereby nonlinear wave speed increases with increasing wave height.

### 3.2. Stagnation Points

Given their connections to high pressure regions, the presence of stagnation points interior to each of the fluid domains impacts both the shape and stability of traveling wave solutions. We proceed to find conditions for the presence of stagnation points interior to each of the fluid domains. In order for there to be a stagnation point inside the fluid domain, there must exist a point inside the domain such that . For simplicity, let represent the location of a stagnation point inside the upper () or lower () fluid domain so that and . At leading order, it is straightforward to show that the relative horizontal profile is given by

 u−c=−ωjz−c0+O(ϵ)

inside of each fluid domain implying that stagnation will occur if there is a inside each fluid domain such that .

For example, in the upper fluid domain, the condition for the existence of a stagnation point is given by

 −ω1z(s)1=c±0→z(s)1=c±0−ω1.

From this perspective, it is clear that there will only be a stagnation point inside the upper fluid domain provided that . Likewise, there will only be a stagnation point inside the lower fluid domain provided that . By enforcing the condition that the stagnation point is inside the appropriate fluid domain, we find the following conditions for the existence of a stagnation point when :

• Upper fluid domain:

• Lower fluid domain:

Figure 4 show the regions in the plane where stagnations points are found in the upper, lower, or both fluid domains. Depending on the configuration of and the presence of the stagnation points for the trivial solution from which we continue greatly impact the shape of the solutions along the corresponding bifurcation branch. These various configurations will be contrasted in the next section with traveling wave solutions.

### 3.3. Constructing Solutions: Numerical Implementation

Using the above formulation, we numerically determine traveling wave solutions by solving for the set for solutions with a fixed spatial period and for fixed wave-height . In the following section, we outline the specific details of the numerical implementation of our pseudo-arclength continuation method111Code is available on GitHub..

To numerically solve (14) & (15), we use a pseudospectral method to solve for the Fourier coefficients of our unknown functions and – differentiation is conducted in Fourier space, while nonlinear operations are computed in physical space. Since we are considering periodic solutions, we numerically represent and by their Fourier series representation with Fourier modes of the form

 η(x)=N∑∑′n=−N^ηneiknx,∂xq2(x)=N∑∑′n=−N^Qneiknx,

where the denotes the summation with as we have both eliminated the average value of as well as the average value of . Enforcing has zero average, we consider zero mean-flow in both the upper and lower fluid. Thus, both and are represented by unknown Fourier coefficients.

Evaluating Equations (14) and (15) for generates equations for the unknowns and for . In order to solve for traveling waves solutions with various amplitudes, we enforce a fixed amplitude for the interface by allowing the wave-speed to vary as a function of the amplitude, thereby introducing both a new equation and new unknown into our system.

Equations (14-15) are solved for iteratively via Newton’s method (though other iterative techniques can also be used). Using a pseudo-arclength continuation, we determine traveling wave solutions for increasingly larger amplitudes.

###### Remark 3.

Due to the exponential nature of the hyperbolic sine and cosine functions, we rewrite the quantities and for as follows

 C1=cosh(μk(ϵη−1))=cosh(μk)(cosh(μkϵη)−tanh(μk)sinh(μkϵη))
 S1=sinh(μk(ϵη−1))=cosh(μk)(sinh(μkϵη)−tanh(μk)cosh(μkϵη))

#### 3.3.1. Visualizing Streamlines and Pressure Contours

For periodic traveling wave solutions, the solutions for the streamlines take the form

 ψj=a0+a1(z+(−1)jδj)−ωj2(z+(−1)jδj)2+∞∑k=−∞eikx^ψ(j)ksinh(k(z+(−1)jδj)).

Choosing the nondimensional function

 φj=e−ikxcosh(μk(z+(−1)jδj)),

where , and , then the analogue of (8) for stationary solutions in a traveling frame is, for , given by

 ∫2π0e−ikx(cosh(μk(ϵη+(−1)jδj))(∂xqj−ωjη−c)+ωksinh(μk(z+(−1)jδj))dx =∫2π0e−ikx∂zψj(x,(−1)jδj)dx (20)

whereas for , we find

 ∫2π0∂zψj(x,(−1)jδj)dx=∫2π0e−ikx(∂xqj−c)dx. (21)

Once we have determined a TWS corresponding to the set , we can then use the above equations to determine the streamfunction within the bulk of each fluid by solving for the coefficient as well as . We also note that the streamfunction can only be determined up to an arbitrary constant within the fluid domain. Here we choose to normalize so that the streamfunctions are zero on the interface. This flexibility allows us to determine the value of .

### 3.4. The Irrotational Case

Before we investigate the effects of shear in each layer, we begin by investigating solutions with no shear (vorticity) in each layer. In Figure 5, we examine the weak stratifications for , and various values of the depth-ratio . For , as we increase the amplitude, both the speed of the traveling wave, and the interface becomes more peaked. However, as increases from to , both the amplitude-speed dependence and qualitative shape of the solution undergo a transformation.

Specifically, for , waves with larger amplitude travel at higher speeds. However, as increases, this relationship changes so that waves of larger amplitude travel at sower speeds than those with lower amplitudes. As continues to increase, the amplitude-speed relationship continues to change until the waves with larger amplitudes again travel at larger speeds. Simultaneously, the qualitative shape of solutions morph from “peaked waves of elevation”, through rounded peaks and troughs, and finally to “peaked waves of depression”. This transition can be accurately predicted by examining how the sign of varies with respect to and for a specific value of as seen in Figure 3.

### 3.5. Including Shear Effects

We now consider the influence of linear shear within each fluid layer. In particular, we focus on the difference between the shape of solutions and the overall bifurcation structure of traveling wave solutions. Specifically, we examine solutions for , and that bifurcate from trivial flows that various configurations of stagnation points within each layer of the fluid.

Noting that we numerically solve for both the free surface and the tangential velocity , we can conjecture a limiting wave-height by imposing that

#### 3.5.1. Bifurcating from trivial flow with no stagnation points in either fluid.

To look at an example in which we bifurcate from a flat surface with a stagnation point in the lower fluid, we choose and . As seen in Figures 6 (a) and (b), the abscence of stagnation points does not appear to bias the interface towards elevation or depression, thus allowing for the development of far more symmetric nonlinear wave profiles than are seen in subsequent sections. We also note that, by referring to Figure 2, the choice of vorticities corresponds to a negative value of . This is corroborated by the speed-amplitude curve in Figure 6 (a), which shows the speed of the wave decreasing with the increase in the amplitude.

#### 3.5.2. Bifurcating from trivial flow with a stagnation point in the upper fluid.

To look at an example in which we bifurcate from a flat surface with a stagnation point in the lower fluid, we choose . As can be seen by examing the interface profiles in Figure 7 (a), the particular case we have chosen corresponds to interfaces of elevation. By examing the impact of the stagnation point on the streamline patterns seen in Figure 7 (b), we can see from where this tendency towards forming interfaces of elevation comes.

#### 3.5.3. Bifurcating from trivial flow with a stagnation point in the bottom fluid.

To look at an example in which we bifurcate from a flat surface with a stagnation point in the lower fluid, we choose . As can be seen by examing the interface profiles in Figure 8 (a), the particular case we have chosen corresponds to interfaces of depression. By examing the impact of the stagnation point on the streamline patterns seen in Figure 8 (b), we can see from where this tendency towards forming interfaces of depression comes.

#### 3.5.4. Bifurcating from trivial flow with a stagnation point in both fluids.

To see an example in which we bifurcate with stagnation points in both the lower and upper fluid regions, we choose and . In this case, we see that the tendency is for waves of depression to form which increase in speed with amplitude; see Figure 9 (a). The role of stagnation becomes increasingly more relevant on the nonlinear profile as its amplitude is increased, and thus has an ever increasing effect on the shape of the wave profile; compare Figures 9 (b) and (c). For lower amplitudes, the presence of stagnation points both above and below the interface allows for the kind of ambiguity between the interface being elevated or depressed for which the case of no stagnation points in the fluid allowed. However, the presence of stagnation points allows for the formation of broader rising portions of heavier fluid, which at larger amplitudes arguably makes the interface one of depression with sharper falling peaks that are shaped by having to move between the stagnation points in the fluid.

## 4. Stability Formulation

We investigate the spectral stability of the traveling wave profiles computed in the previous section with respect to infinitesimal perturbations. First, it is important that we discuss what perturbations we wish to allow. It may appear natural to consider disturbances of the same period as the underlying stationary wave, as is often done in the literature. However, we wish to work with a more general class of disturbances; namely, we choose to include all bounded on the whole real line. It is important to realize that this class is the largest class of perturbations allowed by the physical problem at hand. Indeed, disturbances should be bounded and continuous functions, but there is no physical reason to restrict their spatial dependence to be periodic.

### 4.1. Stability Formulation: Problem Reformulation

In order to investigate the stability of the traveling wave profiles with respect to such perturbations, it is necessary to reformulate the governing equations. Equation (10) is a local statement and does not require modification. However, (11) and (12) are nonlocal and in their current incarnation, apply specifically to waves of period . Thus, we cannot use the equation in its current form for any bounded function on the real line.

Following the method outlined in [13], we reformulate the nonlocal equations on the whole-line via a spatial average value. For a continuous bounded function , i.e. , let represent the spatial average of the function defined by

 ⟨f⟩=limM→∞1M∫M/2−M/2f(x)dx.

It should be noted that the kernel of this operation is quite large. For instance, all functions that approach zero as have zero spatial average. Nevertheless, we may use the spatial average to replace (11) and (12) with the more general nonlocal equations given by

 (22)

and

 ⟨e−ikx((ikϵμ2∂tη−ω2)C2−μk(ϵ(∂xq2−ω2η)−c)S2)⟩=0, (23)

where and are defined as before. These equations are valid for all , as no quantitization condition is imposed by the periodicity of the solutions considered. Further, as solutions of increasingly larger period are considered, the set of values to be considered in (22) and (23) approaches a dense subset of the real line. The equation corresponding to is excluded, as before. Equations (22) and (23) allows us to perturb our traveling wave solution with any bounded perturbation regardless of periodicity.

### 4.2. Stability Formulation: Eigenvalue Problem

Having generalized the dynamical equations to accommodate the perturbations we wish to consider, we briefly discuss the definition of spectral stability. A stationary solution of a nonlinear problem is spectrally stable if there are no exponentially growing modes of the corresponding linearized problem. To determine the spectral stability of the periodic traveling wave solutions, we start by considering a traveling wave solution set , ,). In the same traveling coordinate frame, we add a small perturbation of the form

 η(x−ct,t) = η0(x−ct)+εη1(x−ct)eλt+O(ε2), q1(x−ct,t) = α0(x−ct)+εα1(x−ct)eλt+O(ε2), q2(x−ct,t) = β0(x−ct)+εβ1(x−ct)eλt+O(ε2),

where is a small parameter. The perturbations , , and are moving at the same speed and in the same direction as the original traveling wave solution. Our goal is to determine the time dependence of the perturbation in order to determine how the deviation from the unperturbed solution evolves.

As we are linearlizing about a traveling wave solution, we substitute the above expansions for , and into (10), (22), and (23) keeping only terms. This is rewritten compactly as

 λL1U(x)=L2U(x), (24)

where and are matrices of linear operators and is a vector-valued function with entries . Details are provided in Appendix A.

Since the time dependence of the perturbation depends exponentially on , we can determine information about the stability of the underlying traveling wave by determining all bounded solutions of this generalized eigenvalue problem. If any bounded solutions exist for which the corresponding has a positive real part, the linear approximation of the solution will grow in time and thus the perturbed solution will exponentially diverge from the stationary solution in the linear approximation.

Since the coefficient functions of (24) are periodic in with period , we decompose the perturbations further using Floquet’s Theorem, see for instance [12, 8]. This allows us to further decompose , , and in the form

 η1(x)=eipx¯η1(x),α1(x)=eipx¯α1(x),β1(x)=eipx¯β1(x)

where , , and are periodic with period and (the Floquet exponent) can be restricted to the interval .

Substituting the Floquet decomposition directly into (24) while simultaneously using the Fourier series representation for , and as outlined in [13], we obtain a new equation of the form

 λ~L1(p)^U=~L2(p)^U, (25)

where represents the concatenation of three bi-infinite vectors of the Fourier coefficients for , and .For , are linear operators that now depend on the Floquet exponent . Equation (25) gives a generalized bi-infinite eigenvalue problem for determining the spectrum of the linearized operator about the stationary traveling wave solutions.

Instead of directly solving the eigenvalue problem as stated above, a more stable approach is to reformulate the problem as a quadratic eigenvalue problem for the eigenvalue and the corresponding eigenfunction of the form

 (λ2A2(p)+λA1(p)+A0(p))¯η1=0. (26)

Details regarding (26), the form of the ’s for and are given in Appendix B.

We solve this generalized eigenvalue problem numerically by truncating the Fourier series representation for over to where we have eliminated the zero-mode as we only consider zero-average perturbations to the free-surface. With this truncation, we obtain quadratic equations for unknown Fourier coefficients of . Finally, we note that due to the underlying symmetries in the problem, we may restrict the interval from to and consider both and the conjugate of . Thus, for every , we solve the generalized eigenvalue problem given by

 (λ2A(N)2(p)+λA(N)1(p)+A(N)0(p))^¯η(N)1=0. (27)

where the superscript denotes the projection on to the Fourier modes. This is done via the QZ algorithm [23]. The dimension of the truncation on the Fourier modes is chosen so that both the eigenvalues and eigenvectors converge to 12 digits of accuracy.

## 5. Stability Results

Using the formulation outlined in the previous section, we numerically compute spectra beginning with irrotational interfacial waves. It is well known that the trivial solution is subject to the Kelvin-Helmholtz instability provided that there is a jump in the horizontal velocity at the interface. However, the role that nonlinearities play in the formation of possible Kelvin-Helmholtz instabilities in our configuration is not as well known. In what follows then, we examine the role nonlinearity plays in the formation of instabilities by examining the spectra associated with the non-trivial TWSs presented in the previous section. We begin by considering the spectral stability corresponding to the irrotational problem and then subsequently explore the impact of linear shear.

It is straight-forward to show that the trivial solution will be spectrally stable when there is no jump in the tangential velocity at the interface. However, the spectral stability is not guaranteed for non-trivial solutions that bifurcate from this solution. Throughout this section, we consider two separate classes of instabilities: modulational instabilities (MIs) corresponding to long-wave perturbations (), and high-frequency instabilities (HFIs) corresponding to the Floquet parameter . Note, while computations need only be done for Floquet parameters due to underlying symmetries in the problem, we plot spectra over for illustrative purposes.

For the configuration with , and , Figure 10 shows the maximum growth rate as a function of the Floquet parameter for increasing values of the amplitude of the traveling wave solution . As in the single-layer free-boundary problem, all non-zero, small-amplitude traveling waves appear unstable to perturbations corresponding to HFIs localized around small bands of Floquet parameters . An interesting finding is that for the various depth ratios examined in the irrotational problem with fixed , all solutions are not susceptible to MIs. However, if and are varied appropriately, then even the irrotational problem may become susceptible to both HFI and MI. For example, as seen in Figure 11, when , , traveling wave solutions are susceptible to both HFIs and MIs where magnitude of both instabilities increase as a function of increasing .

Including vorticity has an impact on the strength of HFI as seen in Figure 12, thereby adding a mechanism for the suppression of instability not seen in the stratified shear-free case. Furthermore, the inclusion of vorticity may also impact whether or not a traveling wave is susceptible to MI. This hints at the more complicated role that vorticity plays in understanding and parameterizing regions of instability for the various parameters in the problem. We explore this issue more fully in Figures 13 and 14. Here, we show how different balances between the shear, depth ratio, and either allow or suppress MIs.

For simplicity, in Figure 13, we have chosen