## Abstract

We introduce a set of coupled equations for multi-layer water waves that removes the ill-posedness of the multi-layer Green–Naghdi (MGN) equations in the presence of shear. The new well-posed equations are Hamiltonian and in the absence of imposed background shear, they retain the same travelling wave solutions as MGN. We call the new model the square root depth () equations from the modified form of their kinetic energy of vertical motion. Our numerical results show how the equations model the effects of multi-layer wave propagation and interaction, with and without shear.

## 1. Introduction

The propagation and interactions of internal gravity waves on the ocean thermocline may be observed in many areas of strong tidal flow, including the Gibraltar Strait and the Luzon Strait. These waves are strongly nonlinear and may even be seen from the Space Shuttle (Liu *et al.* 1998) with their crests moving in great arcs hundreds of kilometres in length and traversing sea basins thousands of kilometres across. The MGN model—the multi-layer extension of the well-known Green–Naghdi (GN) equations (Green & Naghdi 1976; Choi & Camassa 1999)—has been used with some success to model the short-term behaviour of these waves (Jo & Choi 2002). Nevertheless, MGN and its rigid-lid version, the Choi–Camassa equation (Choi–Camassa 1996, hereafter CC), have both been shown by Liska & Wendroff (1997) to be *ill-posed* in the presence of background shear. That is, background shear causes the linear growth rate of a perturbation to increase without bound as a function of wave number. Needless to say, this ill-posedness has made numerical modelling of MGN problematic. In particular, ill-posedness prevents convergence of the numerical solution, since the energy cascades to smaller scales and builds up at the highest resolved wave number. Grid refinement only makes the problem worse. Regularization by keeping higher order expansion terms is possible (Barros & Choi 2009), but such methods tend to destroy the Hamiltonian property of the system and thus may degrade its travelling wave structure. Nevertheless, if one is to consider the wave generation problem, one must consider the effects of both topography and shear.

Thus, the MGN equations must be modified to make them well-posed. We shall require that the new system (i) is both linearly well-posed and Hamiltonian; (ii) preserves the MGN linear dispersion relation for fluid at rest; (iii) has the same travelling wave solutions as MGN in the absence of imposed background shear.

The equation we introduce here satisfies these requirements. The remainder of the paper is organized as follows. In §2, we derive the model in the same Euler–Poincaré variational framework as for MGN (Percival *et al.* 2008). In §3, we compare the linear dispersion analysis of and MGN, and thus show that linear ill-posedness has been removed. In §4, we show that the and the shear-free MGN equations possess the same travelling wave solutions. In §5, we present numerical results that compare the wave propagation and interaction properties of the and GN solutions for a single layer. Finally, in §6, we present numerical results for the equations that model the effects of two-layer wave propagation and interaction, both with and without a background shear.

## 2. The governing equations

In this section, we derive the equations by approximating the kinetic energy of vertical motion in Hamilton’s principle for a multi-layer ideal fluid. Such a system consists of *N* homogeneous fluid layers of densities *ρ*_{i},*i*∈[1,…,*N*], where *i*=1 is the top layer and *i*=*N* is the bottom layer. Thus, for stable stratification, *ρ*_{i+1}>*ρ*_{i}. The *i*-th layer has a *horizontal* velocity, *u*_{i} and thickness, *D*_{i}; the interface between the *i*-th and *i*−1-th layer is at depth , for a prescribed bathymetry, *b*(*x*,*y*). We assume columnar motion within each layer (horizontal velocity independent of vertical coordinate, *z*). Incompressibility then implies that vertical velocity is linear in *z*. Under this ansatz and after a vertical integration, the Lagrangian for Euler’s fluid equations with a free surface becomes,2.1In this Lagrangian, the *w*_{i} should be viewed as functions of layer velocities, layer thicknesses and their partial derivatives. The fluid momentum equations are obtained from Euler–Poincaré theory (Holm *et al.* 1998) asThe system is closed by the layer continuity equations,2.2

To obtain the MGN equations, one sets *w*_{i} in equation (2.1) equal to the vertical component of the fluid velocity, represented in terms of the material time derivative as2.3In each layer, we define a vertical *material* coordinate, 0≤*s*_{i}≤1, by *s*_{i}=(*h*_{i}−*z*)/*D*_{i}. Then, by using *ds*_{i}/*dt*_{i}+*w*_{i}∂*s*_{i}/∂*z*=0 with *d*/*dt*_{i}=(∂/∂*t*+*u*_{i}⋅∇), we note2.4

To obtain the new equations, we will approximate the physics of the underlying system contained in the Lagrangian for the MGN equations. We thus replace *w*_{i} in equation (2.1) with a different linear approximation of the vertical motions within each layer, as follows2.5This choice introduces a set of *N* new length scales, *d*_{i}, the far-field fluid thicknesses or equivalently the local thickness of the unperturbed fluid at rest. The new measure of vertical motion in equation (2.5) is related to the natural MGN measure for non-breaking waves in equation (2.4) through two complementary approximations: (i) neglecting the effect of the advective term on the layer coordinate, and (ii) premultiplying by the ratio *d*_{i}/*D*_{i}. Each approximation may be justified individually in the limit of small amplitude and shallow water flow. It can be shown, either though an asymptotic analysis of the MGN equation set, or directly by a scaling argument in the Lagrangian (2.1) that the system approaches the linear shallow water equations in the combined limit of small aspect ratio and long waves. Small amplitude wave-like disturbances can thus have arbitrarily small fluid velocities, but must be supercritical with respect to the linear baroclinic wave speed and are hence of sufficiently small Mach number to justify neglect of advection. Merely dropping advective terms in *u*_{i}⋅∇*s*_{i} in the vertical measure of velocity in the *i*-th layer (i.e. setting *W*_{i}=*D*_{i}(∂*s*_{i}/∂*t*)) gives, on taking variations, an equation set from the Boussinesq long wave family, which of itself is linearly well-posed. However, the travelling wave solutions of this Boussinesq equation set are not those of the MGN equations. This discrepancy between solutions is corrected through the additional rescaling of amplitude. The wave amplitude rescaling, *d*_{i}/*D*_{i}, is obviously of order unity for small amplitude waves.

It will turn out that the composition of these two approximations actually has considerably more validity for wave-like solutions than either of them does separately. We show in §4 that for shear free travelling waves, the quantities defined in equations (2.4) and (2.5) are the same, even for arbitrarily large amplitude waves. Our new measure of vertical motion (2.5) may thus be viewed as a linearization of the full velocity (2.4) around a superposition of shear free travelling wave solutions.

The final quantity in equation (2.5) also equals the vertical fluid velocity, expressed in the *convective* representation (Holm *et al.* 1986). The spatial and convective representations of fluid dynamics are the analogues, respectively, of the body and spatial representations of rigid body dynamics. This analogy arises because the configuration spaces for fluid dynamics and for rigid bodies are both Lie groups. However, further discussion of this analogy is beyond our present scope, so we refer to e.g. Holm *et al.* (1986, 1998).

On taking variations, the final equations arising from Hamilton’s principle for the new Lagrangian may be written compactly in terms of a shallow water equation, plus additional nonlinear dispersive terms that represent a non-hydrostatic pressure gradient,2.6

In comparison, the MGN equations have a similar form, but the non-hydrostatic forces no longer correspond to the gradient of a pressure. Instead, the MGN equations take the form2.7

For both the MGN and equation sets, there exist equivalent continuum systems defined in isopycnal coordinates, for which the relevant vertical material coordinate is given in terms of the physical density, *ρ*, and its vertical stratification, ∂*ρ*/∂*z*, rather than a layer thickness. These continuum equation sets are obtainable from the layer equations in the limit of an infinite number of layers of infinitesimal thickness, in which the summations given above become integrals, while the relationship between height and density is assumed to remain invertible.

Having arisen from Hamilton’s principle, the and MGN equations may both be written in Hamiltonian form. In addition, since their Lagrangians are both invariant under particle relabelling, the two sets of equations each possesses the corresponding materially conserved quantities,The quantity *q*_{i} is the potential vorticity in the *i*-th layer. For the equations,just as for the unmodified shallow water equations. Hence, equation (2.6) admits potential flow solutions for which ∇×*u*_{i}=0. In contrast, the MGN potential vorticity contains higher derivatives of *u*_{i}.

When restricted to a single layer, the Lagrangian for the equation (2.6) isand the single-layer motion equation becomesThe single-layer Lagrangian contains a term that coincides with the Fisher–Rao metric in probability theory (see Brody & Hughston 1998) for an excellent review of the subject). This feature suggested the name ‘ equations’ for the new system.

## 3. Linear dispersion analysis

This section shows that the equation in (2.6) possesses the same linear dispersion relation at rest as MGN, but that unlike MGN they remain linearly well-posed when background shear is present. As discussed in Liska & Wendroff (1997), linear well-posedness requires that the phase speed remain bounded as for all background shear profiles. In what follows, we specialize to the case *N*=2, and impose a rigid lid constraint, *h*_{1}=0, through an additional barotropic pressure term, although the analysis generalizes to an arbitrary number of layers and to the free surface case. Linearizing equations (2.6) and (2.2) for far-field depths, *d*_{1}, *d*_{2} and background velocities, *U*_{1}, *U*_{2} produces the following dispersion relationHere, the phase speed is *λ*=*ω*/*k* for frequency *ω* and wave number *k*. Meanwhile, under the same assumptions, the equivalent CC equations linearize toWhen the background state is at rest, *U*_{1}=*U*_{2}=0, both sets of equations exhibit the same dispersion relation,in which the phase speed of the high wave number modes converges towards zero. This is not surprising, since when linearized at rest, the two measures of vertical motion coincide in equations (2.3) and (2.5), with *w*_{i}=*W*_{i}.

In the presence of background shear, *U*_{1}≠*U*_{2}, the solution for the phase speed *λ* of the rigid lid equations isand for the CC equations isHere the coefficients *A*_{0}, *B*_{0}, *C*_{0}, etc., are functions of the mean depth ratio, *d*_{1}/*d*_{2} and of the shear (*U*_{2}−*U*_{1}). Liska & Wendroff (1997) showed that the discriminant is strictly negative in the CC case, thereby ensuring the existence of a root with growth rate, ℑ(*λk*), which is positive and scales linearly with *k*. Thus, the CC equations are linearly ill-posed. In contrast, the magnitude of growth rates in the case are bounded from above, independently of wave number *k*. Thus, the equations are linearly well-posed. Indeed, provided a Richardson-number condition is satisfied, that3.1then the linearized system remains stable to perturbations at *any* wave number. The dimensionless quantity *Ri* expresses the ratio of the contributions of potential and kinetic energy in the dispersion relation, in terms of the bulk shear flow, *ΔU*/*Δz*=|*U*_{2}−*U*_{1}|/(*d*_{1}+*d*_{2}). Upon taking the standard oceanographic definition,for the Brunt–Väisälä frequency,we see that the bulk density for the equations is not a simple arithmetic mean; instead, it is given by the geometric quantity

## 4. Travelling wave solutions

At this point, we have shown that the equations satisfy requirements (i) and (ii) in §1. We will now show that the and MGN equations also satisfy requirement (iii); that is, they share the same travelling wave solutions for any number of layers in the absence of imposed background shear. We take a travelling wave ansatz for a choice of wave speed, *c*,Integrating the continuity equation (2.2) in a frame moving at the mean barotropic velocity implies for both the MGN and systems that4.1

For such solutions we see from equations (2.3) and (2.5) that *w*_{i} and *W*_{i} coincide,4.2To obtain travelling wave solutions, one varies the Lagrangian subject to the relation (4.1). Using relation (4.2) between the vertical velocities implies that the travelling wave Lagrangians for the and MGN equations coincide. Consequently, the travelling wave solutions of the two equation sets must also coincide. This finishes the derivation of the new equation set and the verification of the three requirements set out in §1. In the presence of background shear, however (*w*_{i}−*W*_{i})≠0, and the coincidence of the travelling waves of the two models no longer holds.

## 5. Numerical experiments (single layer)

We now present some results of numerical computation with the equations for a single layer. The equations are discretized using standard finite volume techniques, without any filtering. In this section, we compare results for the one-layer and single-layer GN equations. These results show that the two systems have the same qualitative behaviour.

### (a) Single-layer lock-release

We first consider a ‘lock-release’ experiment, in which a perturbation to a flat equilibrium interface is held initially at rest, then allowed to flow freely in a domain with periodic boundary conditions. The physical parameters and initial conditions are taken asandWe integrated this initial condition forward in time for both the GN and systems. Figure 1 shows a snapshot of a train of travelling waves emerging from the lock release in each equation set. In both cases, the first wave is the largest and the wave profiles of the two systems track each other closely. This may have been expected, because the two systems share the same travelling wave solutions and possess identical linear dispersion relations for a quiescent background.

### (b) Single layer interaction of solitary waves

We now consider the interaction between pairs of solitary travelling wave solutions in the two systems. The and GN systems have the same travelling wave solutions, but their PDEs differ; so one may expect to see differences in their wave interaction properties. For both systems, the travelling wave takes a sech^{2} form,and the velocity profile is given by equation (4.1). Figure 2 shows a snapshot of the results for the case of two equal amplitude, oppositely directed travelling waves after suffering a symmetric head-on collision. Both the and GN descriptions show a near-elastic collision followed by small-amplitude radiation. The travelling wave behaviour is essentially identical. The main difference lies in the details of the slower linear waves that are radiated and left behind during the nonlinear wave interactions. Similar results were obtained in the case of asymmetric collisions.

## 6. Numerical experiments (two-layer rigid lid)

We now consider numerical solutions of the two-layer rigid lid equations in one spatial dimension. Unlike the result of Jo & Choi (2002) for the CC equations, there is no value of resolution at which the code for numerically integrating the equations becomes unstable. This might have been hoped, because the equations are linearly well-posed.

### (a) Two-layer rigid lid lock-release experiments

We next present the result of a lock-release experiment analogous to that of the previous section, but for the two-layer rigid lid equations. Once again an initial perturbation to the layer height fields is introduced, this time on the interface between the lighter top layer and heavier bottom layer, while the total mean velocity across both both layers is initially held fixed at zero. In this case, we choose *d*_{1}=0.1, *d*_{2}=0.9, *g*=1 andThe initial configuration is taken as a depression of the fluid interface, because both CC and travelling waves extend from the thinner layer into the thicker one. Since our motivating interest was in modelling waves on the ocean thermocline, we will concentrate here on the situation in which the top layer is thinner than the bottom layer. We show results for the following velocity profiles:

— No shear:

*u*_{1}(*x*,0)=0,*u*_{2}(*x*,0)=0 and*x*∈[−100,100].— Background shear flow:

*u*_{1}(*x*,0)=0.01,*u*_{2}=−*D*_{1}*u*_{1}/*D*_{2}and*x*∈[−100,100].

The Richardson number for the case with shear, calculated using equation (3.1), is *Ri*=5.02, and the flow is still linearly stable for all wave numbers.

The result of the integration is shown in figure 3. Unlike the single-layer case, the linear long wave speed in the two-layer case is close to the wave speed of the resulting travelling waves. In figure 4 we compare the leading short wave disturbances of one lobe in the shear-free lock-release experiment to the calculated solitary travelling wave solutions in the same parameter regime, as defined by the maximum amplitude of the disturbance. For the three largest waves the differences are indistinguishable, while the smaller waves, which are still interacting with the long wave signal are wider than the equivalent true travelling wave. A similar result holds for the asymmetric crests generated in the presence of shear and for the respective velocity profiles. The leading order solution is thus well approximated as a train of independent solitary travelling waves propagating in order of amplitude.

### (b) Interaction of two-layer solitary waves with rigid lid

Finally, we consider wave interaction experiments in the two-layer case. Here we initialize with two numerical solutions of the travelling wave problem with oppositely directed velocities. Results for these two experiments are similar to those in the single layer case and also similar to those in Jo & Choi (2002) in the case without shear. Figure 5 shows our numerical solution following the wave collision in the case of equal magnitude wave speed with and without a background shear. We again see nearly elastic collisions, after which the waves approximately regain both their shape and amplitude. In the case with shear, the wave profiles are no longer symmetric. Defining upstream and downstream direction in terms of the background flow in the thicker layer, we see that the downstream gains amplitude and narrows compared with the shear-free wave profile at the same wave speed, while the wave travelling upstream loses amplitude and widens. Nevertheless, the waves are seen to nearly conserve amplitude, and hence momentum through the collision.

## 7. Summary

We have introduced a system of Hamiltonian multi-layer water wave equations that is linearly well-posed in the presence of background shear but possesses the same travelling wave solutions as the MGN equations, whose ill-posedness had previously caused difficulties in their numerical integration. The new system also has been found numerically to generate fast-moving trains of large-amplitude coherent waves that exhibit ballistic, nearly elastic, nonlinear scattering behaviour amidst a background of slow, small-amplitude, weakly interacting or linear wave radiation. Future work will explore wave–wave and wave–topography interactions in the presence of vertical shear between the layers, as well as two-dimensional multi-layer wave interactions and wave generation by flow over topography by using the new well-posed equation set.

## Acknowledgements

We are grateful to W. Choi, F. Dias and R. Grimshaw for fruitful discussions on this topic. J.R.P. was supported by the ONR grant no. N00014-05-1-0703. D.D.H. was partially supported by the Royal Society of London Wolfson Award.

## Footnotes

- Received March 1, 2010.
- Accepted April 28, 2010.

- © 2010 The Royal Society