Kevin G. Lamb

Department of Applied Mathematics, University of Waterloo
Waterloo, Ontario, Canada         N2L 3G1

Two formulations of shallow-water weakly nonlinear theory for internal waves are presented. The structure of solitary waves predicted by these two theories are compared with fully nonlinear solitary waves.


In this talk two different formulations of shallow-water, weakly-nonlinear theory will be considered and the predicted structure of solitary waves will be compared with fully nonlinear solitary waves. The latter are computed numerically using an iterative method (1). The two formulations of the theory are obtained by using different sets of independent variables. For the Eulerian formulation the independent variables are (x,z,t) where x and z are the horizontal and vertical Cartesian coordinates and t is time (2,3). For the Euler-Lagrange formulation of the theory z is replaced by the Lagrangian coordinate $y=z-\eta(x,z,t)$ where $\eta(x,z,t)$ is the vertical displacement of the fluid particle at (x,z,t) from its rest position (4).

Eulerian Formulation

Making the incompressible and Boussinesq approximations, neglecting rotation and viscosity, the equations of motion in Cartesian coordinates are

\begin{align}\nabla^2\psi_t - b_x &=\psi_x\nabla^2\psi_z-\psi_z\nabla^2\psi_x , \tag{1a}\\b_t+N^2(z)\psi_x &= \psi_x b_z-\psi_z b_x, \tag{1b}\end{align}


where$\psi$ is the stream function for the velocity $(u,w) =(\psi_z,-\psi_x)$$b = g\rho'/\rho_o$ where g=9.81 ms-1is the constant gravitational acceleration, $\rho_o$ is a constant reference density and $\rho'(x,z,t)$ is the density perturbation due to the wave, and $\nabla^2$ is the Laplacian operator. The buoyancy frequency N(z) is given by

\begin{displaymath}N^2(z) = -g \frac{d\bar{\rho}}{dz} \geq 0, \tag{2}\end{displaymath} (2)

where $\rho_o\bar\rho$ is the undisturbed density field. For simplicity it is assumed that there is no background flow. The fluid is assumed to lie between two horizontal, rigid planes at z=0,H so that the boundary conditions are $b=\psi_x=0$ at z=0,H.

Under the assumption of small, long waves an asymptotic expansion of the above equations leads to a weakly-nonlinear solution of the form

\begin{displaymath}\begin{split}\begin{pmatrix}\psi \\ [-2mm] _{b}\\ \overline{......{0,2} \end{pmatrix} \biggr\} + \cdots \, .\end{split} \tag{3}\end{displaymath} (3)

Here B = B(x,t) and the $\phi$ and E with the various indices are functions of z$\epsilon$ is a parameter whose value we take to be one. Nondimensional equations would have the same form, in which case $\epsilon$ appears as a small parameter measuring the strength of the nonlinearity and dispersion. Here it is included to indicate the relative ordering of the various terms. The leading-order terms in (3) describe infinitesimal, long waves. These, together with the $O(\epsilon)$ terms, comprise first-order theory. Second-order theory includes the$O(\epsilon^2)$ terms. $\phi(z)$ and the linear long-wave propagation speed c are given by the eigenvalue problem

\begin{displaymath}\mathcal{L} \phi \equiv \left( \frac{d^{2}}{dz^{2}} +\frac{N^{2}}{c^{2}} \right) \phi = 0 \tag{4}\end{displaymath} (4)

with boundary conditions $\phi(0) = \phi(H) = 0$. Only mode-1 waves are considered so c is the largest eigenvalue of (4) with $\phi$ the corresponding eigenfunction. The Ei,jq can be expressed in terms of the $\phi^{i,j}_q$which are determined by equations of the form

\begin{displaymath}\mathcal{L} \phi_{a}^{i,j} = - r_{ijq} \frac{2N^{2}}{c^{3}} \phi+ S_{q}^{i,j} \tag{5}\end{displaymath} (5)

where Si,jq is a known function of the lower-order vertical structure functions. The general solution includes an undetermined constant multiple of $\phi$ (see 3). The constants rijq are determined by the solvability conditions $2r_{ijq}\int_0^H N^2\phi^2\,dz = c^3\int_0^H\phiS^{i,j}_q\,dz$. The first-order constants r10 and r01are given by

\begin{displaymath}Ir_{10} = - \frac{3}{2} \int_{0}^{H} \phi^{\prime 3} dz, \quadIr_{01} = - c \int_{0}^{H} \phi^{2} dz. \tag{6}\end{displaymath} (6)


\begin{displaymath}I = 2 \int_{0}^{H} \phi^{\prime 2} dz. \tag{7}\end{displaymath} (7)

The second order constants are

\begin{displaymath}\begin{split}Ir_{20} &= - \int_{0}^{H} \left\{ \tfrac{8}{3}{r_{01}}{c} \phi^{\prime 3} \right\} dz.\end{split} \tag{8}\end{displaymath} (8)

The theory also yields an evolution equation for B(x,t), namely

\begin{displaymath}\begin{split}B_{t} = - c B_{x} + \epsilon ( 2 r_{10} c B \, ......} B_{xx} + r_{02} B_{xxxxx}) + \cdots\, .\end{split} \tag{9}\end{displaymath} (9)

To $O(\epsilon)$ this is the KdV equation.

Euler-Lagrange Formulation

To use the Lagrangian coordinate $y=z-\eta(x,z,t)$ as the independent variable instead of z we must be able to solve z = z(x,y,t) for y. This requires that $\eta_z<1$ and in particular rules out waves with overturning isopycnals. The governing equations (1a)-(1b) are replaced by

\begin{align}z_{t} & = - \Psi_{x} , \tag{10a} \\\Theta_{t} &= - \frac{\Psi_{y}}{z_{y}} \Theta_{x} - N^{2}(y) \frac{z_{x}}{z_{y}} , \tag{10b}\end{align}


where$\Psi$ and $\Theta$ are the streamfunction $\psi$ and vorticity $\nabla^2\psi$ expressed as functions of (x,y,t).

Solutions of (10a)-(10b) have the form

\begin{displaymath}\begin{split}\begin{pmatrix}\zeta \\ \Psi \\ \Theta \end{pma......0,2}\end{pmatrix} \Biggr\} + \cdots \, .\end{split} \tag{11}\end{displaymath} (11)

where $\zeta =\eta(x,z(x,y,t),t)$ is the isopycnal displacement written as a function of yG = G(x,t), and Z, P and T with the various indices are functions of y. G(x,t) satisfies (9) with identical values for the coefficients. Z and $\phi$ satisfy the same eigenvalue problem so, if they are both scaled to have a maximum value of 1, they are identical. In that case we can take B=G. The remaining vertical structure functions can be expressed in terms of the $\phi^{i,j}_q(y)$

Comparisons with fully-nonlinear waves

The theoretical approximations are now compared with fully nonlinear waves. Two different density profiles are considered in a fluid of depth H=100 m. The first,

\begin{displaymath}\bar{\rho}_{1} (z) = 1 + 0.01 \left(\frac{1-e^{\frac{z-H}{15}}}{1-e^{\frac{-H}{15}}} \right) ,\tag{12}\end{displaymath} (12)

has an exponentially increasing stratification. The second,

\begin{displaymath}\bar{\rho}_{2} (z) = 1.005 - 0.005 \tanh \left( \frac{z -80}{5} \right), \tag{13}\end{displaymath} (13)

has a well defined pycnocline 20 m below the surface. The bottom to surface density difference is 0.01 in both cases. For both stratifications the wave half-widths are approximately 100 m for the largest waves. Comparisons of the isopycnal displacement along the vertical line through the centre of a number of waves are shown in figures 1 and 2. Solid curves are profiles from the fully nonlinear waves. Dashed curves are the theoretical profiles. In the Euler-Lagrange formulation the leading order isopycnal displacement is $\zeta= GZ(y)$. For the Eulerian formulation the isopycnal displacement is obtained by expanding $\bar\rho(z) +\rho'(x,z,t) = \bar\rho(z-\eta(x,z,t))$ which gives $\eta=B\phi(z)$ at leading-order. Since $Z(y)=\phi(y)$ both expression have the identical form however there is a crucial difference between the two. In the Eulerian expression $\eta$is the vertical displacement of an isopycnal as a function of the displaced isopycnal height. In the Euler-Lagrange expression the vertical displacement of an isopycnal is given as a function of the rest height of the isopycnal. In figure 1 it can be seen that for density 1 the Eulerian description is much better than the Euler-Lagrange description while figure 2 shows that the opposite is true for density 2. Using first-order theory the theoretical isopycnal profiles obtained using the Eulerian formulation are almost indistinguishable from the fully-nonlinear wave profiles. This is remarkable as the largest wave considered is very close to breaking (i.e., $\max{\eta_z}$ is close to 1). Use of first-order theory for density 2 also reduces the error significantly but the agreement is not as good. For density 2 the largest wave is far from breaking.

In figures 3 and 4 vertical profiles of the horizontal velocity through the centre of the waves are shown. Solid curves are the profiles for the fully-nonlinear waves. Dashed profiles are the theoretical results. In figure 3 leading- and first-order Eulerian approximations are shown for density 1. In figure 4 leading- and first-order Euler Lagrange approximations are shown for density 2. For the Eulerian formulation the first-order velocity is given by

\begin{displaymath}u = B \phi' + \epsilon (B^{2} c^{2} \phi_z^{1,0} + B_{xx} c\phi_z^{0,1}) . \tag{14}\end{displaymath} (14)

For the Euler-Lagrange formulation u is given, at first-order, by

\begin{displaymath}u = \frac{Bc \phi' + \epsilon (B^{2} P_y^{1,0} + B_{xx} P_y^{0,1})}{1+ \epsilon B \phi'} . \tag{15}\end{displaymath} (15)

For leading-order theory the $O(\epsilon)$ terms are dropped. For the first-order expressions for u the value of Bxx is needed in the centre of the wave. For this the well-known solitary wave solution $B=b_o\hbox{sech}^2(\theta)$ of the KdV equation is used to express Bxx(0,0) in terms of bo. Here $\theta=(x-Vt)/\lambda = 0$ at the wave centre. For the vertical profiles the second-order coefficients r20, r11a, r11b and r02 are not required, however for horizontal profiles they are required to calculate $O(\epsilon)$corrections to $\lambda$ and to the expression for B.

For density 1 the theoretical estimate is very good for even the largest wave. For density 2 the magnitude of the velocities is underestimated for the larger waves well below the pycnocline. The theory does correctly predict a minimum in uimmediately below the pycnocline, something that is not predicted by leading-order theory. The Eulerian formulation incorrectly predicts a subsurface maximum in u for the largest 2 waves (not shown).

These results show that the best formulation of shallow-water weakly-nonlinear theory for solitary internal waves is stratification dependent. Comparisons for more complicated stratifications taken from oceanographic observations will be shown.


(1) Turkington, B., Eydeland, A., and S. Wang, Stud. Appl. Math., 85, 93-127, 1991.

(2) Benney, D. J., J. Math. Phys., 45, 52-63, 1966.

(3) Lamb, K. G., and L. Yan, J. Phys. Oceanogr., 26, 2712-2734, 1996.

(4) Gear, J. A., and R. Grimshaw, Phys. Fluids. 26, 14-29, 1983.

  • Figures ...

  • nextupprevious
    Next:Figures ...