61
61

A Numerical Calculation for Internal Waves over a Slope or a Mound
Taro KAKINUMA

Graduate Schooı of Science and Engineering, Kagoshima University
1‐21‐40 Korimoto, Kagoshima, Kagoshima 890‐0065, JAPANtaro@oce.kagoshima‐u.ac.jp
1. Introduction

Intemal waves show interface deformation when they reach a continental shelf after travelling in the deep

ocean, as studied by e.g. Ostrovsky& Stepanyants (ı989). Grimshaw et al. (2004) numericalıy simulated thetransformation of internal solitary waves using the EKdV eq., obtained by Lee & Beardsley (1974) byextending theKdV

eq. However, the applicability of the nonlinear wave equations, derived using

perturbation methods for two‐ıayer systems, is restricted owing to the perturbation assumptions, such that,for example, the Benjamin‐Ono (BO) equation (Benjamin, 1967; Ono, 1975), as well as the deeper version

derived by Choi & Camassa (1999), can be adopted for internal waves in offshore deep regions, while theKdV equation, and the shallower version by Choi & Camassa (1999), for those in shaılower water over acontinental shelf. In the derivation process for the deeper version by Choi & Camassa (1999), O(hı)<<

O(h_{2}) is assumed, where h_{1} and h_{2} are the thickness ofthe upper and lower layers in stiıl water, respectively,

whiıe for the shaııower version by Choi & Camassa (1999), O(alh)=1 , O(hıll) <<1 , and O(h_{2}ll)<<1 areassumed, where alh, and l, are the representative values for the ratio between wave height and layer
thickness, and wavelength, respectively. Conversely, the set of equations based on a variational principle(Kakinuma, 2001), without any assumption concerning wave nonlinearity or dispersion, is expected to be
applicable to nonlinear waves propagating from deep to shallow water, ifthe vertical distribution ofvelocity

potential is described appropriately.
In the present study, the set of nonlinear equations based on the variational principle, is solved

numerically using the implicit scheme developed by Nakayama & Kakinuma (2010), to simulate internal
waves in two‐layer systems over topography, in consideration of both the strong nonıinearity, and the strong
dispersion, of internal waves.

2. Nonlinear equations for surface/internal waves and numerical method
The illustration in Fig. 1 is our schematic for the multi‐ıayer system of fluids, represented as i(i=1, 2, \cdot\cdot, 1)
from top to bottom, respectively, showing only inviscid and incompressible motion. The thickness of the i‐layer in still water is denoted by h_{l}(x) . None of the fluids mix even in motion, and the density in the i‐layer,\rho_{l} , is spatiaıly uniform and temporalıy constant in each layer.

62y\ap roxL_{\vec{\bulet}

Fig. 1. A schematic for a multi‐layer fluid system.

In the i ‐layer, if both the dispıacement of one interface, z=\eta_{l,1-j}(x,t) ( j=0 or 1), and the pressure on
the other interface, p_{-j}(x,t) , are known, then the unknown variables are the velocity potential \phi, (x,z,t) andthe interface displacement \eta_{i_{j}}(x,t) , such that the definition of the functional for the variational problem in

the i‐ıayer, F_{l} , is as folıows (Kakinuma, 2001), by adjusting the functional introduced by Luke (ı965):

F_{l}[ \phi_{l},\eta_{lj}]=\int_{t_{0} ^{t_{1} \i nt_{A}L_{l}[\phi_{l}, \eta_{lj}]dAdt ,L_{l}[\phi_{l},\eta_{i j},]=\int_{\eta_{l}^{\eta,1}0[\frac{\partial\phi_{l} }{\partialt}+\frac{1}{2}(\nabla\phi_{l})^{2}+\frac{\imath}{2} (\frac{\partial\phi_{l}{\partialz})^{2}+gz+\frac{p_{l-j}+P_{l}{\rho_{l}]dz

(1)

,

(2)

where \eta_{l},0(x,t) and \eta_{l},1(x,t) are the displacements of the lower and upper interfaces of the i ‐layer,
respectively; p_{l},0(x,t) and p_{l},1(x,t) are the pressures at the ıower and upper interfaces of the i‐ıayer,
respectively;

P_{l}= \sum_{k={\imath} ^{l-1}(\rho_{l}-\rho_{k})gh_{k};\nabla is a partial differential operator in the horizontal plane, i.e.,\nabla=(\partial/\partial x,\partial/\partial y);(\nabla\phi,)^{2}\equiv|\nabla \phi_{i}|^{2} . The gravitational acceleration g is 9.8 m/s^{2} . The plane A , which is theorthogonal projection ofthe object domain on to the x‐y plane, is assumed to be independent oftime.

Fluid motion is assumed to be irrotational, resuıting in the existence of velocity potentiaı \phi_{I} , which is
expanded into a power series of vertical position z with weightings f_{l,\alpha_{I} , in a manner similar to that for thederivation of surface wave equations by Isobe (1995), as\phi_{l}(x,z t)=\sum_{\alpha_{l}=0}^{N_{l}-1}[f_{l,\alpha_{\ovalbox{\t \smal REJECT} (x,t)\cdot z^{\alpha_{t} ]\equiv f_{l,\alpha_{I} \cdot z^{\alpha_{I} ,(3)

where N_{l} is the number ofterms for an expanded velocity potential in the i‐layer.
If we consider a two‐ıayer fluid system, over a seabed of a fixed profile, i.e., z=\eta_{2,0}=b(x) , then the
Euler‐Lagrange equations for the variational problem are reduced to the nonlinear surface/internal wave

equations for two‐layer fluid motion as63[Upper layer]\zeta^{\alpha_{1}\frac{\partial\zeta}{\partialt}-\eta^{\alpha_{1} \frac{\partial\eta}{\partialt}+\frac{\imath}{\alpha_{1}+\beta_{1}+1} \nabla[(\zeta^{\alpha_{1}+\beta_{1}+1}-\eta^{\alpha_{\imath}+\beta_{1}+1}) \nablaf_{1.\beta_{1}]-\frac{\alpha_{1}\beta_{1}{\alpha_{1}+\beta_{1}-1} (\zeta^{\alpha_{1}+\beta_{1}-1}-\eta^{\alpha_{\imath}+\beta_{1}-1}) f_{1\beta_{1}=0\zeta^{\beta_{1}\frac{\partialf_{1\beta_{1} {\partialt}+\frac{1}{2}\zeta^{ \beta_{1}+r_{1}\nablaf_{\imath},\beta_{1}\nablaf_{\imath},\gam a_{1}+ \frac{1}{2}\beta_{\imath}\gam a_{\imath}\zeta^{\beta_{1}+\gam a_{1}-2 f_{\imath},\beta_{1}f_{\imath}r_{1}+g\zeta=0 ,\eta^{\beta_{1}\frac{\partialf_{\imath},\beta}{\partialt}+\frac{\imath} {2}\eta^{\beta_{\imath}+\gam a_{1}\nablaf_{\imath},\beta}\nablaf_{1, \gam a_{1}+\frac{1}2\beta_{\imath}\gam a_{1}\eta^{\beta_{1}+\gam a_{1}-2f_ {1\beta_{\imath} f_{1\gam a_{1}+g\eta+\frac{p}\rho_{1}=0

,

(4)

(5)

,

(6)

[Lower layer]\eta^{\alpha_{2} \frac{\partial\eta}{\partialt}+\frac{ \imath} {\alpha_{2}+ \beta_{2}+1}\nabla[(\eta^{\alpha_{2}+\beta_{2}+1}-b^{\alpha_{2}+\beta_{2}+1}) \nablaf_{2,\beta_{2} ]-\frac{\alpha_{2}\beta_{2} {\alpha_{2}+\beta_{2}-1}(\eta^ {\alpha_{2}+\beta_{2}-1}-b^{\alpha_{2}+\beta_{2}-1})f_{2\beta_{2} =0\eta^{\beta_{2} \frac{\partialf_{2\beta_{2} {\partialt}+\frac{1}{2} \eta^{\beta_{2}+\gam a_{2} \nablaf_{2,\beta_{2} \nablaf_{2,\gam a_{2} + \frac{1}{2}\beta_{2}\gam a_{2}\eta^{\beta_{2}+r_{2}-2}f_{2,\beta_{2} f_{2, \gam a_{2} +g\eta+\frac{1}{\rho_{2} [p+(\rho_{2}-\rho_{ \imath} ) ghı ]=0 ,

,

(7)

(8)

where the surface, and the interface, profiıes are described by z=\eta_{1,1}=\zeta(x,t), and z =\eta ı,0 =\eta_{2,1}=\eta(x,t) ,respectively; the pressure at the interface is denoted by p=p_{1,0}=p_{2,1} . The sum rule ofproduct is adopted forsubscripts \alpha_{l}, \beta_{l} , and \gamma_{l} , such that, for example, \alphaı in the first term on the ıeft‐hand side of Eq. (4) denotes
the power of \zeta

After eliminating p from Eqs. (6) and (8), we obtain\eta^{\beta_{2}\frac{\partialf_{2\beta_{2} {\partialt}+\frac{1}{2} \eta^{\beta_{2}+\gam a_{2}\nablaf_{2,\beta_{2}\nablaf_{2,\gam a_{2}+ \frac{1}{2}\beta_{2}\gam a_{2}\eta^{\beta_{2}+\gam a_{2}-2}f_{2,\beta_{2}f_{2}, r_{2}+(1-\frac{\rho_{1}{\rho_{2})g(\eta+h_{1})-\frac{\rho_{\imath} {\rho_{2} (\eta^{\beta_{\imath} \frac{\parti l f_{1\beta_{1} {\parti l }+\frac{1}2\eta^{\beta_{1}+\gam a_{1}\nabl f_{\imath},\beta_{\imath} \nabl f_{1,\gam a_{1}+\frac{1}2\beta_{1}\gam a_ {\imath}\eta^{\beta_{1}+\gam a_{\imath}-2f_{\imath} \beta{\imath} f_{\imath}\gam a_{\imath} )=0.

(9)

In this paper, the top face of the upper layer is assumed to touch a fixed horizontal plate at z=0 , suchthat \zeta=0 without any surface wave. In order to represent both the nonlinearity and dispersion of internaı
waves propagating from deep to shallow, or shallow to deep, regions, the numbers ofterms for the expanded
velocity potentials introduce in Eq. (3), are N_{1}=3 and N_{2}=5 for the upper and lower layers, respectively,
referrin g to the accuracy of the numerical results obtamined by Yamashita & Kakinuma (20ı5), using thesame equations for the internaı solitary waves.
The fundamental equations mentioned above, are solved using the fmite‐difference implicit scheme

developed by Nakayama & Kakinuma (2010). The incident wave is a soıitary‐wave solution, obtained forthe same equations, using the numerical method presented by Yamashita and Kakinuma (2015). In thefollowing discussion, the fluid density ratio \rho_{1}/\rho_{2} is 0.98; in the offshore deep region, the thickness of theupper and ıower ıayers in stilı water, is hı
depth.

=

0.01h and h_{2}=0.99h , respectiveıy, where h is the total water

64\bigwedg_{\tilde{\aleph}

10 20 30 40

0xlh

Fig. 2. The target domain with a slope in Case 1‐A, where h_{2}lh_{1}=99.0 in the deep area I, while h_{2}lh_{1}=4.0 in the
shalıow area III. The broken line shows the critical level z_{c} , defmed by Eq. (10).

\leq_{\grave{N}^{-}
‐

0

10

20

30

40

\lambda/hFig. 3. The target domain with a mound in Case 2‐A, where h2lhı =99.0 in the deep area I, while h_{2}lh_{1}=4.0 in theshallow area III. The broken line shows the critical level z_{c} , defined by Eq. (10).
Tab.1. The still‐water thickness ratios between the lower and upper layers, h_{2}lh_{1} , in the shallow area III.

3. Conditions for numerical computation
Figures 2 and 3 show the target domains for computation, where the totaı water depth is uniform in the areas
I, III, and V, while the slope gradient is 0.ı and -0.1 in the areas II and IV, respectively. The still‐water

thickness ratios between the lower and upper layers, h_{2}lh_{1} , in the shallow area III, are shown in Tab. 1. The
broken ıines in the figures show the elevation of criticaı leveı z_{c} , defmed by Funakoshi & Oikawa (1986) asz_{c}=b/({\imath}+\sqrt{\rho_{2}/\rho_{ \imath} }) .

(10)

65
-0.01

‐0.0105<_{\tilde{P}‐0.011—: Present model----‐0.01150

: BO

1

2

xlh

3

4

5

Fig. 4. The interface profile ofthe incident internal solitary wave, given in the deep area I, in comparison with thecorresponding BO solution.

In the deep area I, the interface exists above the critical level, such that downward convex intemal
waves are stable, while in the shalıow area III for Cases 1‐C and 1‐D, as well as Cases 2‐C and 2‐D, theinterface exists below the critical level, therefore upward convex internal waves are stable.

The grid width, and time step, for computation are \Delta x=0.01h , and \Delta t=0.01\Delta xlC_{\dot{{\imath}},0,I} , respectively,
where C_{\dot{ \imath} ,0,I} is the phase velocity of linear intemal waves in the area I, on the assumption that the internal

waves propagate in shallow water. The phase velocity of linear internal waves in shallow water, C_{i,0} , isevaluated by

C_{i,0}=\sqrt{(\rho_{2}-\rho_{1})gh_{1}h_{2}/(\rho_{2}h_{1}+\rho_{1}h_{2})}(11)Shown in Fig. 4 is the interface profile of the incident minternaı solitary wave, given in the deep area I
for all the present cases, in comparison with the corresponding BO solution, where the ratio between wave

height and upper layer thickness in still water, alhı, is 0.1.

4. Internal waves propagating from deep to shallow regions
Figure 5 shows the interface profile at each time in Case 1‐O, where h_{2}lh_{1}=8.0 in the shallow area III. The
intemaı wave shows disintegration over the continentaı shelf in the area III after passing the sıope, resulting
in two downward convex intemal waves.

In Fig. 6, the interface profile at tC_{i,0,I} lh=34.9 in Case 1‐A, where h_{2}lh_{1}=4.0 in the shallow area III, iscompared with the theoretical solutions for the corresponding intemal solitary wave. According to the figure
for Case 1‐A, although the BO‐type downward convex internal wave is stable in the deep area I, the KdV‐

type downward convex wave appears in the shaııow area III.
Figure 7 shows the interface profile at tC_{i,0,I} lh=30.5 in Case ı‐B, where h_{1}=h_{2} in the shallow area III.

In this case, the intemal wave shows instabiıity, ıeading to a long wave train.
Shown in Fig. 8 is the time variation of the interface profiıe in Case 1‐D, where h_{2}lh_{1}=0.4 in theshallow area III. Partial reflection of the internal wave occurs at the boundary between the areas II and III,

66
.9
.5

.2

9^{g}\underline{x}\dot{e.}5r^{-}\sim
.1

.7
<-

\sim\wedge
+

0. 1

4

\overline{\approx} 0 0\vee\vee-0.1
slopc

xlh

Fig. 5. The time variation ofthe interface profile in Case 1‐O, where h_{2}lh_{1}=8.0 in the shallow area III.

0.02
0\leq\wedge^{-0.02}-<--0.04+-0.06\approx_{-0.08}
-0.1
-0.12
35

36

X/h

37

38

Fig. 6. The interface profile at tC_{\dot{{\imath}},0,I}lh=34.9 in Case 1‐A, where h_{2}lh_{1}=4.0 in the shallow area III, in comparison
with the BO and KdV solutions for the corresponding intemal solitary wave.
0.04—: Present model<\sim0leq\smilewdge+^{-02}\aprox-04 60.\cdot02arow^{\neg}cup^{-}c(/t\backslh\backslh/^{r}|_C/^{} /
-0.08

26

27

28

29

X/h

Fig. 7. The interface profile at tC_{\dot{{\imath}},0,I}lh=30.5 in Case 1‐B, where hı=

h2 in the shalıow area III.

owing to the change in slope gradient. The total length of the internal wave train increases as the waves
propagate.

Depicted in Fig. 9 is the interface profile at tC_{\dot{{\imath}},0,I}lh=34.9 in Case 1‐C, where h_{2}lh_{1}=0.5 in the shallow

area III, in comparison with the KdV solution for the corresponding internal solitary wave. According to the67\neg, |_{\wedge}^{34.9}5.29\simeq-\approx.5\dot{o.}\prime\sim-.1.7\overline{\leq}\sim\wedge
+

0. 1

4

\leq-

0

0

\vee\approx-0.1
slope

xlh

Fig. 8. The time variation ofthe interface profile in Case 1‐D, where h_{2}lh_{1}=0.4 in the shallow area III.
0.

\sim<\wedge-

O.

<-

+-0.

\vee\approx_{-0}.
‐O.

26.5

27

27.5

28

X/h
Fig. 9. The interface profile at tC_{i,0,I}lh=34.9 in Case 1‐C, where h_{2}lh_{1}=0.5 in the shallow area III, in
comparison with the KdV solution for the corresponding internal solitary wave.

figure, the BO‐type downward convex intemal wave in the deep area has changed to a KdV ‐type upward
convex intemal wave in the shallow area, for the interface elevation exceeds the critical level over the slope.

The reversal of convex direction, i.e., convex flipping, of internal waves were shown in the numerical results
obtamined by Helfrich et aı. (1984), as welı as the field data from the observation by Orr & Mignerey (2003).
Figure 10 shows the relative energy, i.e., the ratios ofkinetic and potential energy to the total energy E_{i},in Case ı‐C, where h2lhı=0.5 in the shallow area III. For instance, E_{i,K,I+II} is the kinetic energy of internal

waves in both the deep area I, and the area II with the slope, while E_{\dot{ \imath} ,P,I I} is the potentiaı energy of intemal
waves in the shallow area III. In Case 1‐C, the potential energy becomes larger than the kinetic energy forthe intemal waves in the shallow area III, aıthough the kinetic energy is slightıy ıarger than the potentiaı
energy for the incident intemal solitary wave in the deep area I.

5. Internal waves propagating over a moundIntemal waves passing through a region over a mound, as shown in Fig. 3, have also been numericallysimulated: the lowest elevation of the interfaces at each location is shown in Fig. ıl. In Cases 1‐A and 2‐A,

68
0Q^{-()}q_{A}^-0}\backslah\backslash q^{-0}\ltimes\overline{i}A00tC_{i.0.1}/hFig. 10. The relative kinetic and potential energy, i.e., E_{\dot{ \imath} ,K}lE_{i} and E_{i,P}lE_{i} , respectiveıy, in Case ı‐C, where h2lhı=0.5 in the shallow area III; E_{i} is the total energy of internal waves. For example, E_{i,K,1+1I} denotes the
kinetic energy of internaı waves in both the areas I and II.
t^{1}

-00

-00

\overline{\leq}\backsla h\overline{*}\wedge^{-}\sim-0+\approx\xi-00
-0.

‐O. 1

‐ () . 1,

xlh

Fig. 11. The distributions ofthe lowest elevation \eta_{m\dot{ \imath} n} ofthe interfaces in Cases 1‐A, 1‐B, and 1‐C, where h_{2}lh_{1}=4.0,1.0 , and 0.5, respectiveıy, in the shallow area III shown in Fig. ı, as welı as Cases 2‐A, 2‐B, and 2‐
C , where h_{2}lh_{1}=4.0,1.0 , and 0.5, respectively, in the shallow area III shown in Fig. 2.
where the interface exists above the critical level also in the shallow area III, the ıowest elevation of wave‐

trough interface in Case 2‐A is higher than that in Case 1‐A at xlh>25.0 , which means that the wave heightof the first internal wave decreases after its passing the shallow area III. Conversely, the lowest elevation of
the interfaces in Cases 2‐B and 2‐C is lower than that in Cases 1‐B and 1‐C, respectively, at xlh>25.0.
6. Conclusions

The intemal waves propagating from the deep to shalıow, or the shallow to deep, areas in the two‐layer fluid
systems, were numerically simulated by solving the set of nonlinear equations, based on the variational

69
principle in consideration of both the strong nonlinearity and strong dispersion of internal waves. When the

BO‐type downward convex internaı waves propagated into the shalıow areas, the convex direction ofinterfaces in the shalıow areas, depended on the elevational relationship between the interface, and thecritical level, in the shallow area. In the cases where the elevation of the still‐water interface was below, orequaı to, the criticaı ıeveı in the shaılow area, the disintegration of intemaı waves occurred remarkabıy,leading to a long wave train. When the intemal wave travelled in the deeper area after passing the shallow

area, the wave height decreased when the interface eıevation is higher than the critical leveı also in the
shaııow area.

Acknowledgments
This work was supported by both JSPS Grant‐in‐Aid for Scientific Research (B) Grant Number JP17H02856,
and (C) Grant Number JP17K06585.

References

Benjamin, T. B.: Intemal waves ofpermanent form in fluids ofgreat depth, J. Fluid Mech., Vol. 29, pp. 559‐592, 1967.

Choi, W. and Camassa, R.: Fully nonlinear internal waves in a two‐fluid system, J. Fluid Mech., Vol. 396, pp. 1‐36,
1999.

Funakoshi, M. and Oikawa, M.: Long internal waves oflarge amplitude in a two‐layer fluid, J. Phys. Soc. Jpn, vol. 55,
pp. 128‐144, 1986.

Grimshaw, R., Pelinovsky, E., Talipova, T. and Kurkin, A.: Simulation ofthe transformation of internal solitary waves
on oceanic shelves, J. Phys. Oceanogr., Vol. 34, pp. 2774‐2791, 2004.
Helfrich, K. R., Melville, W. K. and Miles, J. W.: On interfacial solitary waves over slowly varying topography, J.
Fluid Mech., Vol. 149, pp. 305‐317, 1984.
Isobe, M.: Time‐dependent mild‐slope equations for random waves, In Coastal Eng. ı994, pp. 285‐299, 1995.

Kakinuma, T.: A set of fulıy nonlinear equations for surface and internal gravity waves, Proc. 5th Int. Conf. on
Computer Modelling of Seas and Coastal Regions, Wessex Insti. Tech. Press, pp. 225‐234, 200ı.
Lee, C.‐Y. and Beardsley, R. C.: The generation of long nonlinear intemal waves in a weakly stratified shear flow, J.
Geophys. Res., Vol. 79, No. 3, pp. 453‐462, 1974.
Luke, J. C.: A variational principle for a fluid with a free surface, J. Fluid Mech., Vol. 27, pp. 395 ‐ 397 , 1967.
Nakayama, K. and Kakinuma, T.: Intemal waves in a two‐layer system using fully nonlinear internal‐wave equations,Int. J. Numer. Meth. Fluids, Vol. 62, pp. 574‐590, 2010.
Ono, H.: Algebraic solitary waves in stratified fluids, J. Phys. Soc. Jpn., Vol. 39, pp. 1082‐1091, 1975.
Orr, M. H. and Mignerey, P. C.: Nonlinear intemal waves in the South China Sea: Observation of the conversion ofdepression internal waves to elevation internal waves, J. Geophys. Res., Vol. 108, No. C3, pp. 9‐1‐9‐16, 2003.
Ostrovsky, L. A. and Stepanyants, Yu. A.: Do internal solitons exist in the ocean?, Rev. Geophys, Vol. 27, pp. 293 ‐ 310,
1989.

Yamashita, K. and Kakinuma, T.: Properties of surface and intemal solitary waves, In Coastal Eng. 2014, waves. 45, 15
pages, 20ı5.