Natural Hazards and Earth System Sciences Large amplitude internal solitary waves over a shelf

Abstract. Dynamics of large amplitude internal waves in two-layers of shallow water is considered. It is demonstrated that in laboratory experiments the subsurface waves of depression over a shelf may be simulated by internal symmetric solitary waves of the mode 2 ("lump-like" waves). The mathematical model describing the propagation and decaying of large internal waves in two-layer fluid is introduced. It is a variant of Choi-Camassa equations with hydrostatic pressure distribution in one of the layers. It is shown that the numerical scheme developed for the Green-Naghdi equations in open channel flows may be applied for the description of large amplitude internal waves over a shelf.


Introduction
Propagation of high amplitude internal waves in a shelf zone is the very important physical mechanism of coastal water ventilation.The internal waves are able to transport cold water, lift up sediments, intensify mixing processes in shallow waters, etc.Such high-energetic mechanisms of the shelf ventilation can effectively intensify the biological and hydrological processes in coastal waters.In particular, they can redistribute the waste waters and influence the water quality in near shore area.Nonlinear internal waves generated by tides as well as by the interaction of flows with topography play an important role in the energy transfer from the large-scale motion to small-scale mixing.Very often the wave fronts take the form of solitary wave trains with the large ratio of the wave amplitude to the upper layer depth (Helfrich and Melville, 2006).Propagating to shore, they transform into the large amplitude internal waves of elevation.The large amplitude internal waves can be identified by "closed streamline" regions in which trapped Correspondence to: V. Liapidevskii (liapid@hydro.nsc.ru)fluid is carried horizontally for long distances for both types of bottom and surface internal waves (Maxworthy, 1980;Tung et al., 1982;Stamp and Jacka, 1995).Scotti and Pineda (2004) have recently observed bore-like structures with trapped cores in the near shore area of the Atlantics.The transition from wave-like motion to the separate moving soliton-like waves ("solibores") containing trapped dense core is the common feature of the run-up process of internal waves.It can be observed in any shelf zone with high internal wave activity as well as in laboratory experiments.
The run up of internal waves in near shore waters is very similar to the run up of surface long waves (tsunami waves), but the process of internal wave breaking and dissipation is not quite understood.In contrast to the energy dissipation mechanism for surface solitary waves in a homogeneous fluid, the energy dissipation in internal waves is closely connected with the entrainment and mixing in stratified shear flows.Recently, the structure of such a flow has been actively studied both theoretically and experimentally.In stratified water reservoirs, internal waves with a trapped core have been found in the bottom layer (Scotti and Pineda, 2004) and the surface layer (Helfrich and Melville, 2006), and at the interface between homogeneous fluid layers of different densities (Stevens et al., 2005).Intrusion flows in the form of symmetric solitary waves at the interface between fluids have been studied using experimental and theoretical methods (Benjamin, 1967;Davis and Acrivos, 1967;Maxworthy, 1980;Kao and Pao, 1980;Tung et al., 1982;Hohji et al., 1995;Stamp and Jacka, 1995;Schmidt and Spigel, 2000).Interest in this class of flows is motivated by their unique ability to transfer mass along high-gradient interlayers in a stratified fluid due to the initial horizontal momentum.Solitary-wave flow is similar to vortex-ring flow (Akhmetov, 2009), with the entrainment of the surrounding fluid and the vertical velocity component in rotational motion being effectively suppressed by stable stratification.The mathematical model, which is an extension of the Green-Naghdi model and includes the turbulent layer, has been developed in (Liapidevskii and Teshukov, 2000) for stratified flows governed by entrainment and mixing processes.Here this model is adapted to unsteady run-up problems for internal waves of depression and elevation.The plan of the paper is as follows.In Sect.2, we describe the laboratory experiments on solitary wave generation in two-layer fluid at rest.In Sect.3, we present a brief derivation of the adequate mathematical models for large amplitude solitary wave dynamics and introduce the numerical algorithm.In Sect.4, numerical calculations of shoaling solitary waves are compared with the experimental data.

Laboratory experiments
Experiments were performed in a laboratory tank (Fig. 1) 140 cm long, 20 cm wide and 35 cm deep.The walls of the test tank were made of Perspex.The test tank was divided by a vertical removable plate in two parts.The geometrical set-up of the experiments is clear from the sketches shown in Fig. 1.A weak solution of sugar in water was used to create the density stratification.In all experiments, the relative difference in density ε = (ρ 1 − ρ 2 )/ρ 1 was less then 0.005 and the viscosity of the sugar solution is comparable to the viscosity of water.In the first series of experiments (experiment I, Fig. 1a), the tank geometry was chosen as symmetric relative to the plane y = H .The density distribution, which was established between the layers at rest, was also symmetric and was well approximated by the formula where ρ 2 = (ρ 0 + ρ 1 )/2, 2δ is the characteristic thickness of the pycnocline, ρ 0 and ρ 1 are asymptotic values of the density in the upper and lower homogeneous layers, respectively.
The length of the compartment with the mixed fluid of density ρ 2 was chosen so that the only one solitary wave at the interface was produced.The wave of amplitude 2A was also symmetric and exceeded the interlayer thickness considerably (A δ).
The aim of experiment I was to simulate the largeamplitude waves of depression spreading on near surface pycnocline.The symmetry of the second mode internal solitary waves generated in experiment I ensured that only the lower part of the flow 0 ≤ y ≤ H was considered.It follows from the flow symmetry that the corresponding Reynolds stress component vanished at y = H .
In experiment II, a solitary wave of elevation was generated in the bottom layer of density ρ 1 and, in contrast to experiment I, was propagating under the action of the bottom friction force (Fig. 1b).It has been shown in Gavrilov and Liapidevskii (2010) that an initially symmetric solitary wave can transform quickly in the highly turbulent gravity flow for thin bottom layers.Nevertheless, the dense bottom layer was chosen rather thick in experiment II so that the pycnocline was above the shelf zone (Fig. 1b).In this case, the wave of elevation kept the form of the solitary wave during its shoaling.
The front propagation speed in both experiments was evaluated from video-records obtained by a digital camera located at 2.5 m from the tank (25 frames per s).The flow visualization is demonstrated in Fig. 2, where a photograph of the bottom solitary wave of elevation (exp.II) is shown against the background of the luminous screen with a grid of inclined lines imposed on it (Ermanyuk and Gavrilov, 2005).In the zones with a high density gradient, a specific distortion of these lines was observed, and in mixing zones, the optical transparency of the fluid was reduced (Ermanyuk and Gavrilov, 2007).In both types of experiments, the fluid of density ρ 2 in the left compartment was slightly coloured with an ink solution for visualizing the fluid core trapped by a solitary wave.The bold and the thin lines in Fig. 2 represent the exact and numerical solutions of model 3, which are discussed in the following Sects.3 and 4. 3 Mathematical models

Model 1
As a basic model for description of nonlinear internal waves, we consider the two-layer shallow water equations derived in Miyata (2000) and Choi and Camassa (1999) for large amplitude, weakly nonhydrostatic long waves.In the Boussinesq approximation (0 < (ρ 1 − ρ 2 )/ρ 1 1), the second mode of internal waves at interfaces generated in the first series of experiments are symmetric with respect to the channel midline y = H for initial stratification shown in Fig. 1a.Therefore, the flow can be considered in the lower part of the tank (0 < y < H ). Let y = h(t,x) and u = u(t,x) be the depth and the mean velocity of the lower homogeneous layer with the density ρ 1 , η = η(t,x) and w = w(t,x) be the depth and the mean velocity of half of the intermediate layer with the density ρ 2 (ρ 2 = (ρ 0 +ρ 1 )/2,ρ 0 < ρ 2 ).
For non-uniform bottom topography z = z(x) we have For the second series of experiments on shoaling solitary internal waves of elevation (Fig. 1b), the notations are the same.In this case ρ 2 = ρ 0 and (2) is fulfilled.
In the Boussinesq approximation, the two-layer shallow water equations over the uneven bottom take the form (Helfrich, 2007) Here g is the gravity acceleration, p is the pressure at interface, β = H /L 1, L is a typical wave length.Note that the parameter β can be put equal to 1, but we keep it in the formulae to have better scaling between horizontal and vertical axes.The friction terms f ± are represented as follows: The friction coefficients c ± , c i at the upper and lower boundaries of flow, as well as at the interface are supposed to be constants.
By virtue of (2), (3) note that In experiments Q ≡ 0 and (3) can be reduced to the system of two equations, say, for the dependent variables h and u.The travelling waves of (2)-(3) describing the large amplitude solitary waves are discussed in the next section.In the following subsections, two models, which are the simplifications of ( 2)-(3), are introduced.Note, that in (3) the terms, which have the order O(β), are kept only.In the case of slowly varying topography (z = z 0 (βx)), we have x have the order O(β 2 ) and they are neglected in models 2 and 3.Moreover, if one of the layers is thin compared to the wave amplitude, the fluid particles in such a wave move almost with the wave velocity and the corresponding material derivative vanishes in the layer.It results in the hydrostatic pressure distribution in the layer.
For symmetric solitary waves realized in the first series of experiments (Fig. 1a), the intermediate layer is supposed to be hydrostatic and the propagation of the large amplitude wave of depression in the lower part of the tank (0 < y < H ) is governed by model 2, which is introduced in Sect.3.2.Solitary waves of elevation, propagating in a thin hydrostatic bottom layer (Fig. 1b), are described by model 3 in Sect.3.3.The effects of the omitted nonhydrostatic components in the simplified models on the structure of travelling waves are considered in Sect.3.4.Further, for models 1-3, we chose the dimensionless variables such that ρ 1 = 1, b = 1, H = 1.

Model 2
The model 1, in which the hydrostatic intermediate layer is supposed and the terms O(β) are kept only, takes in the dimensionless variables form (Gavrilov and Liapidevskii, 2009) In view of ( 2), (4) the system (6) can be written in the form Here We will apply Eqs. ( 7)-( 8) for numerical calculations of non-stationary problems, concerning the propagation of large amplitude solitary waves of depression over a shelf.The variables h and K are considered as the dependent variables, describing time evolution of the flow, and the velocity u = u(t,x) can be restored for given h, h x , K from the linear ODE The realization of such approach, which is analogous to the method of calculations developed in (Le Metayer et al., 2010) for the Green-Naghdi model, will be discussed in the Sect.3.6.

Model 3
The next model describes the propagation of large amplitude waves of elevation in the case of thin initial bottom layer (Fig. 1b).Now the bottom layer is supposed to be hydrostatic and the model 3 can be derived from (3) in a similar way to (6), namely, Analogously, in view of ( 2), (4) the system (10) is reduced to where Note that ( 11), ( 12) can be transformed into ( 6), ( 7) by changing dependent variables h and η, u and w.Therefore, the numerical scheme, developed for model 2, can be easily applied to model 3.In the next section, we show that the solutions of the models 2 and 3 give the effective approximation of solutions of model 1 for large amplitude solitary waves of depression or elevation propagating along thin initially unpertubated layer.

Travelling waves
For the models 1-3, we construct travelling waves, i.e., solutions of Eqs.(3), ( 6), (10), which depend on the variable ξ = x − Dt (D ≡ const).Such solutions exist for nondissipative flows over a flat bottom (c ± = 0, c i = 0, z(x) ≡ 0).Let's consider the solitary waves satisfying for |ξ | → 0 the following conditions (here and below, the prime denotes differentiation with respect to the variable ξ ).For model 1 the solitary waves have been investigated in (Choi and Camassa, 1999), for model 2 (and, consequently, for model 3) they have been constructed in Gavrilov and Liapidevskii (2009).It is shown in Gavrilov and Liapidevskii (2010) that for models 1-3 the velocity distribution and the profiles of solitary waves can be represented in the form and model 1 Here Fr = D/ √ bH = D is the Froude number.It follows from ( 15)-( 17) that the solitary waves exist for where Therefore, the admissible intervals for h are Note also that is the necessary condition of solitary wave existence.
It is worthy to note that the functions G 1 (h) and G 2 (h) have the common limit for h 0 → 1 The same is true for the functions G 1 (h), G 3 (h) when h 0 → 0 For i = 0,1,2,3,4 the equation specifies the profile of symmetric solitary waves for admissible intervals The solution of ( 23) can be found in quadratures Here h m = h + for a wave of depression and h m = h − for a wave of elevation.The solutions (24) for Fr = 0.49 and β = 1 (i = 0,1,2,3,4) are shown in Fig. 3. Since the change of dependent variable h → 1 − h does not vary the function G 1 (h) and transforms the function G 2 (h) into G 3 (h), the solitons of elevation, calculated for h 0 = ζ < h − , and the solitons of depression, calculated for h 0 = 1 − ζ > h + , have the profiles, which are symmetric relative to the line y = 1/2.Therefore, it is sufficient to consider only the waves of elevation.It follows from Fig. 3 that for h 0 = 0.2 the solutions of model 1 (curve 1; Eq. 15) and model 3 (curve 3; Eq. 16) are close to each other, in contrast to the solution of model 2 (curve 2; Eq. 17).Note also that for h 0 < 0.1 the profiles of large amplitude solitary waves practically coincide.Therefore, for numerical calculations of nonstationary problems on nonlinear wave propagation in twolayer flow, we shall use the more simpler models 2 and 3.The limiting case of the solitary wave with h 0 = 0 (curve 0; Eq. 21) is also shown in Fig. 3.It was mentioned above that the wave amplitude is the same as for the waves 1 and 2, but the solution with h 0 = 0 is not soliton since, in view of ( 21), ( 24), the length of the wave is finite in contrast to the waves with h 0 > 0. The solitary waves of depression which are plotted in the upper part of Fig. 3 are analogous to the considered wave of elevation.The curves are represented by ( 24) for i = 1,2,3,4.

Solitary waves in homogeneous fluid
For h H and ρ + = 0, ρ − = 1 the model 2 transforms in the dispersive shallow water equations for open channel flows (Green and Naghdi, 1976;Serre, 1953) Analogously to (7), Eqs. ( 25) can be rewritten as follows (Gavrilyuk and Teshukov, 2001) where We will introduce in the next section the numerical scheme for two-layer flows, which can also be applied for ( 26).This approach is promising in the nonlinear wave shoaling problem due to its insensitivity to the small bottom perturbations.But here we consider only the solitary waves over an even bottom to find the proper value of the scaling coefficient β from the experimental data for homogeneous fluid.
The function G 2 (h) in ( 16) for the case under consideration takes the form (Miles, 1986) In Fig. 4, the dependence h/ h 0 = ϕ(t/ √ gh 0 ), calculated from ( 27) for fixed x 0 = ξ + Dt and for two values of parameter β (β = 0.75 and β = 1) is compared with the experimental data from (Jensen et al., 2003).Experimental solitary waves in (Jensen et al., 2003) (experiment 4, x 0 = 3.37 m, h 0 = 0.2 m) have been generated by moving a wall in an unperturbated layer of water.We can see that the Green-Naghdi equations with β = 1 give the wave profile, which is too wide, and the value β = 0.75 fits better for high amplitude solitary waves.For internal waves considered in experiments I and II, the coefficient β depends on the stratification and two-layer scheme of flow is just the approximation of the real flow.Nevertheless, we could choose the fixed value of the parameter β, namely, β = 0.5, which gives the reasonable agreement between numerical calculations and experimental data for the non-stationary problem considered in next subsection.

Non-stationary problem
In real stratified fluids, solitary waves cannot propagate at a constant velocity due to the dissipation of energy.Friction and entrainment of surrounding fluid lead to a deceleration of wave.Note that energy dissipation in a solitary wave is caused also by mass loss during its propagation.For a two-layer scheme of flow, such effects can be simulated by including the friction terms f ± in the models 1-3.
For numerical treatment of non-stationary problems, we confine ourselves to the models 2 and 3, which are more simpler then model 1, because they are similar in structure to the Green-Naghdi equations for a one-layer fluid (Le Metayer et al., 2010).We apply the numerical scheme developed in (Le Metayer et al., 2010) for ( 6) and (10).As the governing equations, we consider the equivalent forms ( 7) and ( 11) with dependent variables (h, K) for model 2 and (h, M) for model 3 (Gavrilov and Liapidevskii, 2010).The numerical scheme is formally a version of the Godunov's scheme, in which the fluxes through the lateral boundaries of the meshes are calculated by a Riemann solver.For ( 7) and (11) the equilibrium system without dispersion (β = 0) is represented by two-layer shallow water  27) with β = 1, β = 0.5, correspondingly; 3 -experiment in (Jensen et al., 2003).equations The characteristics of ( 28) have the form Equations ( 28) are hyperbolic for |u−w| < 1.The characteristics (29) are used in the Riemann solver construction and for time step definition.The numerical algorithm for ( 7), ( 11) can be realized in a similar way to the constructed one in (Le Metayer et al., 2010).For (7) the time evolution of conservative variables (h, K) is calculated by a Godunov's type scheme, then the functions u and u x are restored from the ordinary differential equation ( 9).This provides the infinite velocity of perturbation propagation for (7) in spite of the fact that ( 28) is hyperbolic for the moderate velocity shear in the layers.A numerical scheme for (11) can be constructed analogously.In all calculations presented in the paper the boundary conditions are u = 0, w = 0 at the side walls of the tank.As the initial conditions, we use the exact solution (24), describing a solitary wave for the corresponding model.The topography y = z(x) is also included, when the problem on the solitary wave shoaling is considered.

Shoaling internal waves
In this section, we consider the results of numerical calculations of the solitary wave propagation over a shelf in the frame of models 2 and 3.The results are compared with corresponding experimental data from laboratory experiments I and II.To specify the model coefficients, the evolution of solitary waves of depression and elevation in a horizontal channel is considered first.

Solitary waves in a horizontal channel
Numerical calculations of the symmetric internal waves of the second mode along the pycnocline has been considered in (Gavrilov andLiapidevskii, 2009, 2010) (the model 2).It has been shown that Eqs. ( 6) represent the form and the decay rate of internal solitary waves in accordance with the laboratory experiments.It was found that the main input in the energy dissipation of waves was given by the interfacial friction term since c + vanished due to the symmetry of flow and the bottom friction could be neglected (c − ≡ 0).In all calculations the value c i = 0.012 was chosen.In the present paper, we again put c i = 0.012 for models 2 and 3.The coefficients c + , c − for model 3 were chosen as c ± = 0.004, but their values didn't influence essentially the numerical solutions of the non-stationary problem for relatively thick initial bottom layers.In models 2 and 3 describing solitary wave propagation in a stratified fluid, we put β = 0.5 for all calculations mentioned in the paper.
The numerical calculation of the solitary wave propagation over a horizontal bottom is shown in Fig. 5 (the model 3).The exact solution (24) with i = 3, h 0 = 0.135, Fr = 0.48, H = 10 cm, was used as initial data.The calculated wave profiles are depicted with the fixed time interval (4 s in dimensional variables corresponding to the experiment II shown in Fig. 2).One of the profiles (thick line in Fig. 5) is shown in Fig. 2 (curve 2) together with the exact solution of (17) (curve 1) and the experimental photo as the background.The bold points in Fig. 5 are the crest positions of the wave of elevation found from the experiment.One can see from Fig. 5 that the mathematical model 3 describes the decay of the solitary wave and its phase position adequately.Note also that the solitary wave lost its initial symmetry during propagation due to friction effects.The shedding rate of mass from the wave can also be found by the numerical calculations.

Solitary waves over a shelf
Let us consider first the transformation of the solitary wave of elevation over a shelf realized in experiment II.The evolution of the solitary wave is shown in Fig. 6.The thick lines represent the numerical solutions of (10).The initial position of the solitary wave in two-layer fluid and the shelf form are shown in Fig. 6d.The wave profiles calculated through the fixed time interval (4 s in dimensional variables of experiment II) are drown in Fig. 6d together with the experimental crest positions of the wave (bold points).The waves depicted by thick lines in Fig. 6d correspond to the waves shown in Fig. 6a-c.We can see from Fig. 6 that in spite of the rather steep shelf slope in the experiment (α = 8 • ), model 3 gives a realistic description of solitary wave transformation over the shelf.
As was mentioned in Sect.2, we consider the second mode solitary wave to simulate a large amplitude subsurface internal wave of depression in a shelf zone.It is possible so far that the waves are symmetric relative to the plane y = H .The transformation of a solitary wave in the second mode in the vertical channel contraction (experiment I) is shown in Fig. 7.The thick lines represent the numerical solution of (6) (model 2), so they are symmetric relative to the plane y = H .The shoaling solitary internal wave, calculated by (6) transformed over the shelf in the "top-table" wave calculated by 2-D model for continuously stratified fluid in (Lamb, 2003).In experiment I (Fig. 7), the wave crest started to flatten over the shelf slope, but due to flow instability the first mode was generated and the wave lost its symmetry.This resulted in the more intensive loss of mass by the wave than predicted by the mathematical model.Nevertheless the main features of the wave evolution was captured by the model.

Conclusions
The amplitude of the internal waves propagating in a twolayer fluid may exceed many times the initial depth of one of the layers without considerable breaking in contrast to the high amplitude surface waves.Therefore, internal waves can transport mass and momentum over a large distance in the coastal waters.The mathematical models 2 and 3 with the hydrostatic pressure distribution in one of the layers show that the formation of large-amplitude solitary waves in a two-layer fluid in a channel of finite depth is due to  dispersion effects in the outer layer.Therefore, for the models considered above, soliton-like solutions have a limit as the thickness of one of the layers tends to zero.Laboratory experiments on shoaling internal waves have demonstrated that the generation and decay of large amplitude internal waves in a two-layer fluid may be described by the shallow water equations of the second order approximation with the terms representing the turbulent friction at the interface.Models 2 and 3 are rather simple variants of the Green-Naghdi equations in open channel flows, which are applied for description of large amplitude internal waves in a shelf zone.It is shown that the numerical scheme developed for the open channel flows may be applied for description of large amplitude internal waves over a shelf.

Fig. 1 .
Fig. 1.A sketch of the lock problem: (a) solitary wave at the interface, experiment I; (b) shoaling solitary wave, experiment II.

Fig. 6 .
Fig. 6.Solitary wave over a shelf: thick lines -numerical calculations by model 3 corresponding to experimental data; bold pointsexperimental positions of the wave crests.

Fig. 7 .
Fig. 7. Transformation of the second mode solitary wave in a vertical channel contraction: thick line -numerical (symmetric) solution by model 2.