NUMERICAL ANALYSES ON PROPAGATION OF NONLINEAR INTERNAL WAVES

A set of nonlinear surface/internal-wave equations, which have been derived on the basis of the variational principle without any assumptions concerning wave nonlinearity and dispersion, is applied to compare numerical results with experimental data of surface/internal waves propagating through a shallowor a deep-water region in a tank. Internal waves propagating over a submerged breakwater or a uniformly sloping beach are also simulated. The internal progressive wave shows remarkable shoaling when the interface reaches the critical level, after which physical variables including wave celerity become unstable near the wave-breaking point. In the case of the internal-wave trough reflecting at the vertical wall, the vertical velocities of water particles in the vicinity of the interface are different from that of the moving interface at the wall near the wave breaking, which means that the kinematic boundary condition on the interface of trough has been unsatisfied.


INTRODUCTION
In a lake or an ocean where density stratification is well developed, not only internal long-period waves, e.g.internal seiches and tides, but also internal short-period waves are observed.The sources of the latter include bottom topography and interfacial instability.In coastal zones, internal waves play an important role in nearshore environment including nutrient salts, water temperature, etc., so that it is important to know the characteristics of internal-wave propagation in consideration of both nonlinearity and dispersion of internal waves.In the present study, internal waves in two-layer systems are numerically simulated using a set of surface/internal-wave equations (Kakinuma, 2001).In the derivation process of the equations, no assumption is used for nonlinearity and dispersion of waves, such that the application of this model is expected to be theoretically free from limitations concerning the relative thickness of fluid layers or the frequency band of surface/internal waves.
First the verification of numerical results through the present model is performed in two-layer systems.In shallow-water cases, computational results of interface displacements up to each order on the vertical length scale of motion are compared with the corresponding calculation results through a Boussinesq-type internal-wave model or the existing experimental data.On the other hand, in deepwater cases, experimental data of surface/interface displacements are obtained and compared with the corresponding computational results through the proposed model.
Second the numerical model is applied to internal waves propagating over a submerged breakwater or a uniformly sloping beach, where physical variables including water-particle velocity and internalwave celerity are examined near wave-breaking points.

NONLINEAR EQUATIONS FOR SURFACE/INTERNAL WAVES Multilayer Fluids
Inviscid and incompressible fluids, as shown in Fig. 1, are assumed to be stable in still water.The ilayer thickness in still water is denoted by h i (x).None of the fluids mix even in motion and the density ρ i (ρ 1 < ρ 2 < …< ρ I ) is spatially uniform and temporally constant in each layer.Surface tension and capillary action are neglected.
Fluid motion is assumed to be irrotational, resulting in existence of velocity potential φ i defined as , and z w where ( ) , , i.e., a partial differential operator in the horizontal plane.

Functional of the Variational Problem
The pressure on z = η i,0 , i.e., the lower interface of the i-layer, is written by p i (x,t).In the i-layer, if both the elevation of one interface, z = η i,1-j (x,t) (j = 0 or 1), and the pressure on the other interface, p i-j (x,t), are known, then the unknown variables are the velocity potential φ i (x,z,t) and interface elevation η i,j (x,t) such that the functional for the variational problem in the i-layer, F i , is determined by where ( ) ; g is gravitational acceleration; the plane A which is the orthogonal projection of the object domain on to the x-y plane is assumed to be independent of time.
In comparison with the functional referred to in Luke (1967) for rotational motion, Eq. ( 2) has an additional term of pressure on interfaces, without the terms relating to vorticity.

Vertically Distributed Functions of Velocity Potential
In order to derive a set of equations of a horizontally two-dimensional type, vertical integration is performed analytically.In a manner similar to that in the derivation process of the nonlinear surfacewave model (Isobe, 1995), the velocity potential φ i is expanded into a series in terms of a given set of vertically distributed functions Z i,α multiplied by their weightings f i,α , i.e.,

(
) where the sum rule of product is adopted for subscript α.

Euler-Lagrange Equations
We substitute Eq. (3) into Eq.( 2), after which the functional F i is integrated vertically, see APPENDICES.Then the variational principle is applied to obtain the following Euler-Lagrange equations, i.e., the fully nonlinear equations for surface and internal waves (Kakinuma, 2001) (e = 0 or 1);

Selection of Vertically Distributed Functions in Expanded Velocity Potential
In the present paper, the vertically distributed function Z i,α is determined by

NONLINEAR EQUATIONS FOR SURFACE/INTERNAL WAVES IN A TWO-LAYER SYSTEM
Nonlinear Equations for Internal Waves in a Two-Layer System Covered with a Fixed Horizontal Plate If we consider a two-layer system covered with a fixed horizontal plate, where η 1,1 = 0 and η 2,0 = b(x), and the interface profile is described by z = η (x,t), where η = η 1,0 = η 2,1 , then Eqs. ( 4) and ( 5) are reduced to the following equations.
[Upper layer] (i = 1) In the upper layer, i = 1 and j = 0. Equation ( 6) is substituted into Eqs.( 4) and ( 5), resulting in ( ) In the lower layer, i = 2 and j = 1.Equation ( 6) is substituted into Eqs.( 4) and ( 5), leading to Nonlinear Equations for Surface and Internal Waves in a Two-Layer System If we consider a two-layer system with a free surface where the pressure p 0 is set at zero and the surface profile is described by z = ζ (x,t) = η 1,1 , then Eqs. ( 4) and ( 5) are reduced to the following equations.

NUMERICAL CALCULATION METHOD
Two-layer problems are solved in vertically two-dimensional cases.The set of Eqs. ( 7), ( 8), ( 9), and (10) or the set of Eqs. ( 11), ( 13), ( 14), and ( 16) is rewritten to a set of finite difference equations and the time development is carried out by applying implicit schemes similar to that of Nakayama and Kakinuma (2010).

NUMERICAL SIMULATION OF INTERNAL WAVES REFLECTING IN A TANK BETWEEN TWO FIXED HORIZONTAL PLATES
Shallow-Water Case Horn et al. (2000) performed hydraulic experiments using a tank, where the length L, depth D, and width W were 6.0, 0.29, and 0.3 m, respectively, as shown in Fig. 2(a).Three ultra-sonic wave gauges were set at the positions marked A, B, and C. The tank was filled with a two-layer stratification, where h 1 /D = 0.8, after which it was rotated very slowly around the axis of rotation.At the beginning of experiments, this tilted tank, where the tilt angle was θ, was returned to a horizontal position very quickly, after which internal waves traveled in the two-layer system between two fixed horizontal plates.
In the initial condition of numerical computation, the tank was horizontal and the interface was inclined linearly as shown in Fig. 2(b); the initial velocity potential was assumed to be zero through the computational domain.The grid width ∆x and the time-step interval ∆t were equal to 0.06 m and 0.02 s, respectively.
In Figs. 3 -5, the experimental and calculation results are compared for the time series of interface displacements at the position marked C in Fig. 2(a) in the case where the density ratio ρ 2 /ρ 1 and tilt angle θ were 1.019 and 0.4617°, respectively.Figure 3 shows the experimental interface displacement measured by the wave gauge at position C. Figure 4 shows the corresponding calculation result through a Boussinesq-type model (BT), the fundamental equations of which are described in APPENDICES.
Figure 5 shows the calculation result obtained using the proposed strongly nonlinear model (SN), where the number of vertically distributed functions in expanded velocity potential, N, was equal to 1, 2, 3, 4, or 5.When N = 1, the set of nonlinear internal-wave equations, i.e., Eqs. ( 7), ( 8), (9), and (10), is reduced to a set of nonlinear and non-dispersion internal-wave equations, which shows extreme disintegration around the wave crests as shown in Fig. 5, without dispersion balancing with nonlinearity.
When N = 2, the SN takes into account linear and uniform distributions of u i and w i in the direction of z, respectively, such that the balance between nonlinearity and dispersion is considered, leading to the more accurate result than that when N = 1.
When N = 3, the interface displacement evaluated using the SN is close to that shown in Fig. 4 obtained through the BT, where the parabolic distribution of u i in the direction of z is considered in both the SN and the BT.It should be noted that though the effect due to the linear distribution of w i in   the direction of z is considered in both the SN and the BT, the SN estimates the wave period, or the wave number, more accurately by solving the contribution of each order without perturbation, while the wave period through the BT is longer than that in the experiment.Interface profiles obtained using the SN, where N = 4, and the BT when t = 280 s are shown in Fig. 6, according to which the representative ratio of water depth to wavelength, h 2 /λ, was about 0.06.The SN estimated the wave heights larger and the wavelengths shorter with higher nonlinearity than the BT.

Deep-Water Case
Similar hydraulic experiments on two-layer systems between two fixed horizontal plates were performed using a tank, where the length, depth, and width were 0.6, 0.3, and 0.1 m, respectively.The thickness ratio h 1 /h 2 and density ratio ρ 1 /ρ 2 were equal to 0.25 and 0.802, respectively.
When the tilted tank was returned to a horizontal position, the interface profile was measured and given as the initial interface profile in numerical computation; the initial velocity potential was assumed to be zero through the computational domain.The number of terms in expanded velocity potential, N, was equal to three.Equations ( 7), ( 8), (9), and (10) were solved, where ∆x and ∆t were 0.005 m and 0.00005 s, respectively.
Figure 7 shows interface displacements of both calculation results and experimental data, where the tilt angle θ was 3.8° and the position x of measurement was 0.15 m.An internal seiche was generated.Also in the deep-water case, where the BT based on a perturbation method around the shallow-water condition should not be applied, the calculation results through the SN are in harmony with the corresponding experimental data.

NUMERICAL SIMULATION OF SURFACE AND INTERNAL WAVES REFLECTING IN A TANK
Hydraulic experiments on two-layer systems where the top face, i.e., the upper face of the 1-layer, was a free surface touching the air were performed using a tank of length 1.663 m, height 0.34 m, and width 0.1 m.The density ratio ρ 1 /ρ 2 was equal to 0.802.
In numerical computation, the set of Eqs. ( 11), ( 13), ( 14), and ( 16) was solved, where ∆x and ∆t were equal to 0.01663 m and 0.0001 s, respectively.The initial velocity potential was assumed to be zero through the computational domain.
Figure 8 shows time series of surface and interface displacements obtained through the experiment, as well as the calculation, where the thickness ratio h 1 /h 2 and tilt angle θ were 0.4 and 1.8°, respectively.The position x of measurement was 1.065 m.The calculation results show good accuracy especially in the wave periods.The surface-wave mode was dominant in both the surface and the interface displacements for the thickness of the upper layer was quite thinner than that of the lower layer.
On the other hand, Fig. 9 shows the time series of surface and interface displacements obtained through the experiment and calculation, where the thickness ratio h 1 /h 2 and tilt angle θ were 2.3 and 1.9°, respectively.The position x of measurement was 1.065 m.The predominant modes were different between the surface and the interface displacements.

NUMERICAL SIMULATION OF INTERNAL WAVES OVER A SUBMERGED BREAKWATER Initial Setup Using a Third-Order Theoretical Solution of Internal Solitary Wave
The proposed model was applied to internal waves propagating over a submerged breakwater as shown in Fig. 10, where the breakwater shoulder was located where x = 5.0 m and the lateral boundaries were vertical walls of perfect reflection.The thickness ratio h 1 /h 2 and density ratio ρ 1 /ρ 2 were equal to 0.25 and 0.98, respectively.
A third-order theoretical solution of internal solitary wave in a two-layer system between two fixed horizontal plates derived by Mirie and Pennell (1989) was given in the initial condition of an interface profile and velocity potentials in computation.The initial wave amplitude was equal to 0.02 m.The number of terms in expanded velocity potential, N, was equal to three.The set of Eqs. ( 7), ( 8), (9), and (10) was solved, where ∆x and ∆t were equal to 0.02 m and 0.005 s, respectively.

Wave Shoaling and Disintegration
Numerical results of interface profiles are shown in Fig. 11.The internal-wave trough was inclined backward due to the nonlinearity over the breakwater slope, after which wave disintegration occurred just over the breakwater shoulder.
Figure 12 shows vertical distributions of horizontal velocity u 2 and dynamic pressure p d , i.e., total pressure minus hydrostatic pressure in consideration of the interface displacement, below the internalwave trough when passing the places.At the beginning of wave disintegration, both u 2 and p d showed remarkable curvature in their vertical distributions with wave dispersion.

NUMERICAL SIMULATION OF INTERNAL PROGRESSIVE WAVES OVER A UNIFORMLY SLOPING BEACH
Internal waves propagating over a uniformly sloping beach as shown in Fig. 13(a) were simulated by applying the present numerical model.As mentioned above, a third-order theoretical solution of internal solitary wave was given in the initial condition, where the amplitude of internal wave was 0.02, 0.04, 0.05, 0.07, or 0.08 m.The number of terms in expanded velocity potential, N, was equal to three.The set of Eqs. ( 7), ( 8), ( 9), and (10) was solved, where ∆x and ∆t were equal to 0.02 m and 0.005 s, respectively.In Fig. 14, when the interfaces reached the critical level, the internal-wave profiles began inclined backward, after which the internal wave showed remarkable disintegration.
Figure 15 shows the internal wave height H, i.e., the vertical distance between the first trough and the first crest when the internal-wave trough passed each place.It should be noted that the amplitude of the first crest increased gradually after the internal-wave trough touched the critical level.
Computational results of wave celerity and horizontal velocities of water particles in the vicinity of trough interfaces when the internal-wave troughs passed each place are shown in Fig. 16, where the wave celerity is determined by the horizontal velocity of moving internal-wave trough, while the horizontal velocities of water particles in the vicinity of an interface, i.e., u 1,η and u 2,η in the upper and lower layers, respectively, are ).
According to discrete variables, see Fig. 17, the horizontal velocity of water particle in the vicinity of an interface where  ( ), where the subscript m indicates the grid-point number in the direction of x and ∆z i = η m -z i .
In every case of Fig. 16, there existed a place where the physical variables showed rather discontinuous change; after the internal-wave trough passed this place, the calculation stopped due to numerical divergence.The last profiles of interfaces shown in Fig. 14 are the profiles just before the stops of calculation.According to Fig. 14, the maximum slope of interface was not more than 2.0 near the breaking point in any present case, such that it was not the reason why the calculation stopped that the proposed vertically-integrated model could not treat discontinuous change in an interface profile but that the variatioal principle could not find correct solutions where the kinematic boundary conditions on interfaces were satisfied on the assumption of the expanded velocity potentials.
In addition, the condition where the Kelvin-Helmholtz instability occurs for linear shallow-water waves is described as where k is wave number and ε = (ρ 2 − ρ 1 )/ ρ 1 ; we assume that | U 1 − U 2 | = | u 1,trough − u 2,trough |.In any case of Fig. 16, the Kelvin-Helmholtz instability did not occur, which means that the wave breaking was not caused due to this instability.7), ( 8), ( 9), and (10) was solved, where ∆x and ∆t were equal to 0.02 m and 0.005 s, respectively.
Vertical velocities of water particles in the vicinity of an interface, i.e., w 1,η and w 2,η in the upper and lower layers, respectively, at a vertical wall are Figure 18 shows time variations of interface displacements and vertical velocities of water particles in the vicinity of interfaces at the wall where x = 10.0 m.In the case where the initial wave amplitude a was equal to 0.02 m, w 1,η was in agreement with w 2,η at the vertical wall.
On the other hand, in the case where the initial wave amplitude a was 0.08 m, w 1,η became larger than w 2,η just after the interface displacement showed the minimum value at the wall where x = 10.0 m, which means that in the normal direction of the interface the velocity components of water particles in the vicinity of the interface were different from that of the moving interface at the vertical wall such that the kinematic boundary condition on the interface was unsatisfied to detach the upper layer from the lower one in the wave-breaking process.

CONCLUSIONS
The surface and internal waves in the two-layer systems were simulated by solving numerically the set of nonlinear equations in consideration of both strong nonlinearity and strong dispersion of waves.The calculation results of surface/interface displacements were in harmony with the corresponding experimental data in both the shallow-and the deep-water cases.
The internal waves propagating over the uniformly sloping beach were simulated by applying the present numerical model.When the internal-wave trough reached the critical level, the internal-wave profile began inclined backward, after which the internal wave showed disintegration remarkably, increasing the amplitude of the first crest gradually.The celerity of internal wave, as well as the horizontal velocities of water particles in the vicinity of the interface, showed rather discontinuous change near the wave-breaking point for the variational principle could not find correct solutions where the kinematic boundary condition was satisfied on the interface; then the calculation stopped due to numerical divergence.
In the reflection of internal wave at the vertical wall, the kinematic boundary condition on the interface was unsatisfied in the wave-breaking process, where the initial amplitude of internal solitary wave was large.
In future work the wave-breaking criteria of internal waves will be discussed using more terms in expanded velocity potentials such that we can investigate the details of the wave-breaking mechanisms.

ACKNOWLEDGMENTS
The author would like to extend sincere gratitude to Prof. M. Oikawa, Fukuoka Institute of Technology, and Res.Assoc.H. Tsuiji, Kyushu University, for their beneficial comments on the theoretical solutions of internal solitary waves.

APPENDICES Euler-Lagrange Equations of the Variational Principle
The first variation of the functional F i determined by Eq. ( 2) is ( ) where closed curve C is periphery of plane A; s and n are tangential and outdirected-normal axes to the curve C, respectively, such that if δ F i is equal to zero for any δφ i and δη i,j , the equation of continuity in the i-layer, the kinematic boundary condition on the lower interface of the i-layer, the kinematic boundary condition on the upper interface of the i-layer, and dynamic boundary condition on the unknown interface of the i-layer are satisfied with appropriate initial and lateral-boundary conditions.
We substitute Eq. (3) into Eq.( 2), after which the functional F i is integrated vertically as follows: Then the Euler-Lagrange equations of the variational problem are The numerical schemes with finite difference methods to solve these equations are resemble to those applied in the proposed nonlinear model.
(a) Dimensions of the laboratory tank.(b) Initial condition.The interface is inclined linearly with the angle θ θ θ θ in the horizontal tank.

Figure 3 .
Figure 3.Time series of interface displacement measured using the wave gauge at position C in the hydraulic experiment.

Figure 4 .
Figure 4. Time series of interface displacement obtained through the Boussinesq-type model.

Figure 5 .
Figure 5.Time series of interface displacements obtained through the present model, where N is the number of vertically distributed functions in expanded velocity potential.

Figure 6 .
Figure 6.Interface profiles obtained through the present model (SN), where the number of vertically distributed functions in expanded velocity potential, N, is equal to four, and the Boussinesq-type model (BT), when t = 280 s.

Figure 5
Figure5hardly shows difference in interface displacement obtained through the SN between the cases where N = 4 and 5.Although both the SN and the BT do not include dissipation effects due to friction, which results to wave heights larger than the experimental data, the harmony of results between the SN and the BT indicates the high accuracy of results calculated using the SN in this long-wave condition.Interface profiles obtained using the SN, where N = 4, and the BT when t = 280 s are shown in Fig.6, according to which the representative ratio of water depth to wavelength, h 2 /λ, was about 0.06.The SN estimated the wave heights larger and the wavelengths shorter with higher nonlinearity than the BT.

Figure 10 .
Figure 10.Submerged breakwater and the initial interface profile.

Figure 11 .
Figure 11.Time variation of interface profile.

Figure 12 .
Figure 12.Vertical distributions of horizontal velocity u2 and dynamic pressure pd below the internal-wave trough when the trough passed the places.

Figure 13 .
Figure 13.Computational domains including a uniformly sloping beach.The slope of the beach is 1/50.

Figure 14 .
Figure 14.Time variations of Interface profiles.The broken line shows the critical level obtained through the KdV theory.

Figure 15 .
Figure 15.Wave height when the first internal-wave trough passed each place.

Figure 16 .
Figure 16.Horizontal velocities of water particles in the vicinity of the interface and wave celerity when the first internal-wave trough passed each place.

Figure 17 .
Figure 17.Interface profile.The filled circles denote water particles in the vicinity of the interface, where m is a grid-point number in the direction of x.

Figure 18 .
Figure 18.Time series of interface displacements and vertical velocities of water particles in the vicinity of interfaces at the vertical wall where x = 10.0 m.
. (4) and (5) are derived as the nonlinear surface/internal wave equations.Boussinesq-Type Equations for Internal WavesWe also utilize the numerical model based on the Boussinesq-type equations, i.e., φ 2 are velocity potentials in the upper and lower layers, respectively.