Energy Spectrum of Linear Internal Wave Field in the Vicinity of Continental Slope

Ranis N. Ibragimov¹, Austin Biondi¹, Nathan Arndt², Maria Castillo², Guang Lin³, Vesselin Vatchev⁴, Sheng Zhang³

¹Wenatchee Valley College, Wenatchee, WA, USA.
²University of Wisconsin-Parkside, Kenosha, WI, USA.
³Purdue University, West Lafayette, IN, USA.
⁴University of Texas Rio Grande, Edinburg, TX, USA.

DOI: 10.4236/jamp.2020.810169                

Abstract

The purpose of the research was to investigate two-dimensional modeling of efficiency of mixing, resulting from the reflection of a linear internal wave field (IWF) off a continental slope. Efficiency of deep ocean mixing was associated with the energy balance of the radiating IWF into an interior of the ocean in the vicinity of a sloping bottom topography. Since waves are generated not only at the fundamental frequency but also at all of its harmonics ω_n = nω less than buoyancy frequency N and greater than Coriolis frequency f, our analysis includes, in general, an infinite number of discrete internal wave modes n satisfying the dispersion relationship for internal waves. However, since we are interested only in the radiating part of the field, the mode numbers are limited. Due to multiple singularities of order two caused by resonance in the vicinity of critical slope, the energy is visualized in L_δ -norm with δ > 2. Research results include the visualization of the effects of the continental slope and the Earth’s rotation on resulting energy in the vicinity of the slope.

Keywords

Internal Wave Field, Effects of Rotation, Continental Slope, Energy Radiation

1. Introduction

Internal waves occur in density stratified fluids in the presence of a gravitational field. They arise as a result of perturbations which force the stratified fluid to move vertically (i.e. against gravity). Internal wave field (IWF) has periods between the inertial period $2\pi /f$ and the local buoyancy period $2\pi /N$. In this notation f is the inertial frequency, which depends on the rotation rate of the earth (angular velocity $\Omega$ ) and the latitude $\varphi$ :

$f=2\Omega \sin \left(\varphi \right)$ (1)

and N is the (depth-dependent) buoyancy frequency, defined by:

$N\left(z\right)={\left(\frac{-g}{{\rho }_0}\frac{\partial \rho }{\partial z}\right)}^{\frac{1}{2}},$ (2)

where g is the gravitational acceleration and ${\rho }_0$ is the mean density of the medium and $\Omega =2\pi \text{rad}/\text{day}\approx 0.73\times {10}^{-4}{\text{s}}^{-1}$ is the rate of the earth’s rotation. In the mid/high latitude upper ocean, N is typically one or two orders of magnitude larger than f ( $N~{10}^{-3}{\text{s}}^{-1}$ ; $f~{10}^{-4}{\text{s}}^{-1}$ ). The quantity N measures the degree of density stratification of a fluid with average potential density $\rho \left(z\right)$.

Linear IWF is characterized by an amplitude A, a wavenumber $\stackrel{\to }{k}$, and frequency $\omega$, which are related by a dispersion relation $\omega =\omega \left(\stackrel{\to }{k}\right)$.

The seafloor is one of the critical controls on the ocean’s general circulation. Near bottom topography such as a continental slope, and internal waves, sometimes reaching storm proportions, causes that a strong mixing precesses playing an important role in redistributing sediments and organic matter in deep water, leading to the formation of vast sediment drifts [1] [2] [3] [4]. Mixing efficiency in the ocean due to interaction of IWF near the critical slope of a continental slope makes the IWF generation more effective than that due to supercritical or subcritical slopes [4] [5]. The enhancement of vertical mixing over regions of rough topography has important implications for the abyssal stratification and circulation of the ocean [6] [7] [8].

The subject of mixing processes due to breaking of internal waves in the ocean and atmosphere has received much attention in recent years. In particular, the problem of determining the sources of mechanical energy for mixing in the oceans and atmosphere is an important problem for local micrometeorology [3] [9] [10] [11] [12] and in controlling large-scale ocean circulation [13] [14] [15] [16]. Recent measurements in regions of continental slope characterized by ridges and valleys running up and down the slope by [17] and [18] suggest that 40% - 50% of the required energy available for mixing is produced by tide-topography interactions. This has recently been confirmed in our previous analytical studies [14] [19], where the process of internal wave generation by the interaction of an oscillatory three-dimensional background flow over a uniform slope has been considered (in the model considered in [19], the upper radiation condition has been used). Observations by Polzin et al. in [20] and [21] suggest that diapycnal mixing is strongly enhanced over rough bottom topography which can be associated with the conversion of barotropic tidal energy into baroclinic tidal energy and the subsequent instability and breaking of internal waves in the ocean.

The length of the world’s continental slopes in the oceans sums up to 300,000 km. About 65% of all continental slopes descend in deep sea regions. The rest of the slopes located in the shallow regions of the oceans (see Figure 1). It has been observed that the “Pacific Continental Slope” is steeper than the “Atlantic Continental Slope” whereas the continental slope gradients are flattest in the Indian Ocean [18]. The continental slopes occupy almost 8.8% of the world’s surface.

As has been pointed out in [22], increased habitat heterogeneity in canyons is responsible for enhancing benthic biodiversity and creating biomass hotspots. So canyons in the U.S. East Coast were a high priority for federal and state agencies research programs.

The primary focus of this paper is to investigate two-dimensional efficiency of mixing, resulting from the reflection of a linear IWF off a continental slope. The effects of rotation on a linear IWF waves are investigated numerically and using an analytical non-hydrostatic Boussinesq and rigid lid approximations. However, because of the inclusion of high-frequency waves, we cannot make the hydrostatic approximation. Efficiency of deep ocean mixing is associated with the energy balance of the radiating IWF into an interior of the ocean in the vicinity of a sloping bottom topography. Since waves are generated not only at the fundamental frequency but also at all of its harmonics ${\omega }_n=n\omega$ spanning the range of possible frequencies $\omega \in \left[f\mathrm{,}N\right]$, our analysis includes, in general, an infinite number of discrete internal wave modes n satisfying the dispersion relationship for internal waves. However, since we are interested only in the radiating part of the IWF, the mode numbers are limited. In order to separate the spectrum irregularities caused by waves of fundamental frequency of mode $n=1$ from resonance due to higher harmonic waves ( $n>1$ ) in the vicinity of critical slope, the energy is visualized in $L_{\delta }$ -norm with $\delta >2$.

Figure 1. Canyons on the edge of the Continental Shelf show former river channels of Pliocene/Pleistocene age. Source: National Oceanic and Atmospheric Administration (NOAA), Deep-Water Mid-Atlantic Canyons Exploration 2011.

2. Two-Dimensional Boussinesq Model and Linearization

In the domain of continuous motion (i.e. in the domain where the velocity $\stackrel{\to }{v}=\left(u\mathrm{,}v\mathrm{,}w\right)$, pressure p, density $\rho$ and their first derivatives are continuous), integral conservation laws are equivalent (see e.g. [23] ) to the following system of differential equations [6] [24] [25]:

${\stackrel{\to }{v}}_t+\stackrel{\to }{v}\cdot \nabla \stackrel{\to }{v}=-\frac{1}{\rho }\nabla p-2\Omega \times \stackrel{\to }{v}-g\stackrel{\to }{k}\mathrm{,}$ (3)

$\nabla \cdot \stackrel{\to }{v}=\mathrm{0,}$ (4)

${\rho }_t+\stackrel{\to }{v}\cdot \nabla \rho =\mathrm{0,}$ (5)

which are the momentum balance equation, the continuity equation and conservation of mass.

In Equations (3)-(5), $\nabla$ is the gradient operator and $\stackrel{\to }{k}$ denotes a unit vector in the vertical direction. By convention, the xy-plane is tangent to the surface of the ocean with $\left(x\mathrm{,}y\mathrm{,}z\right)$ pointing (east, north, up). The term $2\Omega \times \stackrel{\to }{v}$ is referred to as the Coriolis acceleration. Note that the momentum balance Equation (3) is a vector equation of motion valid for any coordinate system rotating with the earth.

In most geophysical systems, the fluid density varies, but not greatly, around a mean value (more specifically, the density variations in the ocean do not exceed 3% - 4%). This observation is the basis for the Boussinesq approximation, which replaces the exact density by a constant representative value in many terms of the equations of motion (see also [25] ).

Without approximation, we may divide the exact ocean density into a constant part ${\rho }_0$, a part $\bar{\rho }\left(z\right)$ that varies only with depth, and a perturbation ${\rho }^{\prime }\left(x\mathrm{,}y\mathrm{,}z\mathrm{,}t\right)$ which represents the density variation caused by the waves. We thus write

$\rho ={\rho }_0+\hat{\rho }\left(x\mathrm{,}y\mathrm{,}z\mathrm{,}t\right)\mathrm{,}$

where $\hat{\rho }=\bar{\rho }\left(z\right)+{\rho }^{\prime }$. Since typical oceanic sizes are

${\rho }_0~1.04\text{gm}\cdot {\text{cm}}^{-3}\mathrm{,}\left|\bar{\rho }\right|~0.03\text{gm}\cdot {\text{cm}}^{-3}\mathrm{,}\left|{\rho }^{\prime }\right|~0.003\text{gm}\cdot {\text{cm}}^{-3}\mathrm{,}$

we can assume $\left|{\rho }^{\prime }\right|\mathrm{,}\left|\bar{\rho }\right|\ll {\rho }_0$ and the Boussinesq approximation is quite accurate. Similarly,

$p\left(x\mathrm{,}y\mathrm{,}z\mathrm{,}t\right)=\bar{p}\left(z\right)+p^{\prime }\left(x\mathrm{,}y\mathrm{,}z\mathrm{,}t\right)\mathrm{,}$

where $\bar{p}\left(z\right)=P_0-{\rho }_{0}gz-g{\displaystyle \int }\bar{\rho }\text{d}z$ and we require ${\rho }_0+\bar{\rho }$ and $\bar{p}$ to be consistent with the state of rest, i.e.

$\frac{\text{d}\bar{p}}{\text{d}z}=-\left({\rho }_0+\bar{\rho }\right)g\mathrm{.}$

With these definitions, since

$-\nabla p-\rho g\stackrel{\to }{k}\equiv -\nabla \left(-{\rho }_{0}gz-\bar{\rho }gz+p^{\prime }\right)-{\rho }_{0}g\stackrel{\to }{k}-\hat{\rho }g\stackrel{\to }{k}=-\nabla p^{\prime }-{\rho }^{\prime }g\stackrel{\to }{k}$

we can write the momentum equation as

${\rho }_0\left(1+\frac{\bar{\rho }}{{\rho }_0}+\frac{{\rho }^{\prime }}{{\rho }_0}\right)\left[\frac{D\stackrel{\to }{v}}{Dt}+2\Omega \times \stackrel{\to }{v}\right]=-\nabla p^{\prime }-{\rho }^{\prime }g\stackrel{\to }{k}.$

Thus, to a very good approximation,

$\frac{D\stackrel{\to }{v}}{Dt}+2\Omega \times \stackrel{\to }{v}=-\nabla \left(\frac{p^{\prime }}{{\rho }_0}\right)-\frac{{\rho }^{\prime }g}{{\rho }_0}\stackrel{\to }{k}\mathrm{,}$ (6)

where $D=\frac{\partial }{\partial t}+\stackrel{\to }{u}\cdot \nabla$ is the material derivative.

The approximation (6) is called the Boussinesq approximation [25].

In summary, the Boussinesq approximation, rooted on the assumption that the density does not depart much from a mean value, has allowed us to replace the exact density $\rho$ by its reference value ${\rho }_0$ everywhere, except in front of the gravitational acceleration and in the energy equation, which has become an equation governing density variations.

In two dimensions, the incompressibility condition is taken into account by introducing a stream function $\psi$ which is related to the two velocity components u and w by $u={\psi }_z,w=-{\psi }_x$. We write Equations (3)-(5) as the system of three equations for $\psi \mathrm{,}\rho$ and v as follows:

${\nabla }^2{\psi }_t-g{\rho }_x-f{v}_z=J\left(\psi \mathrm{,}{\nabla }^2\psi \right)\mathrm{,}$ (7)

$v_t+f{\psi }_z=J\left(\psi \mathrm{,}v\right)\mathrm{,}$ (8)

${\rho }_t+\frac{N^2}{g}{\psi }_x=J\left(\psi \mathrm{,}\rho \right)\mathrm{,}$ (9)

where the nonlinear advective terms are written by means of the Jacobian operator $J\left(a\mathrm{,}b\right)=a_x{b}_z-a_z{b}_x$ and subscripts denote partial differentiation. We look for solutions

$\left(\psi ,v,\rho \right)=\stackrel{\infty }{\underset{i=1}{{\displaystyle \sum }}}{\epsilon }^i\left({\psi }_i,v_i,{\rho }_i\right),$ (10)

where $\epsilon$ characterizes the size of the perturbation due to the waves. We then have

$\begin{array}{l}\left\{{\nabla }^2\frac{\partial {\psi }_1}{\partial t}-g\frac{\partial {\rho }_1}{\partial x}-f\frac{\partial v_1}{\partial z}\right\}+\epsilon \left\{{\nabla }^2\frac{\partial {\psi }_2}{\partial t}-g\frac{\partial {\rho }_2}{\partial x}-f\frac{\partial v_2}{\partial z}\right\}+0\left({\epsilon }^2\right)\\ =\epsilon J\left({\psi }_1,{\nabla }^2{\psi }_1\right)+0\left({\epsilon }^2\right),\end{array}$ (11)

$\left\{\frac{\partial v_1}{\partial t}+f\frac{\partial {\psi }_1}{\partial z}\right\}+\epsilon \left\{\frac{\partial v_2}{\partial t}+f\frac{\partial {\psi }_2}{\partial z}\right\}+0\left({\epsilon }^2\right)=\epsilon J\left({\psi }_1,v_1\right)+0\left({\epsilon }^2\right),$ (12)

$\left\{\frac{\partial {\rho }_1}{\partial t}+\frac{N^2}{g}\frac{\partial {\psi }_1}{\partial x}\right\}+\epsilon \left\{\frac{\partial {\rho }_2}{\partial t}+\frac{N^2}{g}\frac{\partial {\psi }_2}{\partial x}\right\}+0\left({\epsilon }^2\right)=\epsilon J\left({\psi }_1,{\rho }_1\right)+0\left({\epsilon }^2\right),$ (13)

In the linear approximation, obtained by setting $\epsilon =0$, we can derive a single equation for ${\psi }_1$ which we will call $\psi$,

${\nabla }^2{\psi }_{tt}+N^2{\psi }_{xx}+f^2{\psi }_{zz}=0.$ (14)

If N is a constant we can look for plane wave solutions of the form

$\psi =A\exp \left\{i\left(kx+mz-\omega t+\phi \right)\right\},$ (15)

where A is a constant and $\phi$ is the phase constant. Substitution of representation (15) into (14) leads to the dispersion relation of internal waves

${\omega }^2=\frac{N^2{k}^2+f^2{m}^2}{k^2+m^2},$ (16)

where k and m are wavenumbers of the wave vector $\stackrel{\to }{k}=\left(k,m\right)$. Dispersion relation (16) relates the frequency $\omega$ to the parameters $N\mathrm{,}f\mathrm{,}k$ and m. Unlike for surface waves, the wave vector of internal plane waves is perpendicular to the group velocity and, thus, to the direction of the energy propagation. Hereafter it is assumed that N is a constant value, which is set to be 10⁻³ s⁻¹. As follows from the relation (16), the range of possible frequencies is limited,

$f\le \omega \le N\mathrm{.}$ (17)

This means that internal waves with frequencies greater than buoyancy frequency or with frequencies less than inertial frequency cannot exist. Internal waves of inertial frequency are horizontal rotations of water parcels; internal waves at the buoyancy frequency (when m/k is small, $\omega \approx N$ and, hence, rotation effects are negligible) are purely vertical oscillations of water parcels. In both these limiting cases the group velocity is zero and there is no wave propagation. At intermediate frequencies, the parcels move in oblique elliptical orbits. It is worth noting that f goes to zero as the equator is approached and therefore the rotational effect becomes smaller and this has an important effect on the internal wave structure and dynamics.

2.1. Vertical Modes and Rigid Lid Approximation

With the boundary conditions,

$\psi \left(0\right)=\psi \left(H\right)=0$ (18)

appropriate for a flat bottom at $z=0$ and rigid lid at $z=H$, where H is the (range-independent) depth of the water column, we can look for a separable solution of the form

$\psi \left(x,z,t\right)=\Psi \left(z\right){\text{e}}^{i\left(kx-\omega t+\phi \right)}.$ (19)

To a good approximation, a rigid lid approximation holds for IWF which produces very small vertical displacement at the free surface. Because of the small density variations ( $\Delta \rho /\rho ~{10}^{-3}$ ), the vertical displacement involved in internal wave motions can normally reach several meters in amplitude. In some isolated cases, the displacement can be several tens of meters, e.g. internal solitary waves generated by tidal flow over topography. For example, [26] presents fascinating observations of large solitary waves in the Adaman Sea.

Substituting this into Equation (14) we obtain a differential equation for the vertical structure of internal waves:

${\Psi }^{″}+k^2\left(\frac{N^2-{\omega }^2}{{\omega }^2-f^2}\right)\Psi =\mathrm{0,}0\le z\le H\mathrm{.}$ (20)

Solving this equation using the boundary conditions (18) provides the modes, the eigenfunctions of the system. These modes are oscillating only in the region of the water column where the internal wave frequency $\omega$ is between f and N. Outside this region the solutions can’t satisfy both boundary conditions and they are evanescent (i.e. exponentially decreasing). Each frequency $\omega$ has its own set of modes ${\Psi }_{\omega \mathrm{,}n}\left(z\right)$, where n indicates the order of the mode. The general solution of equation (20) can be expressed as a linear combination of these modes:

$\Psi \left(x\mathrm{,}z\mathrm{,}t\right)=\stackrel{N}{\underset{f}{{\displaystyle \int }}}\underset{n}{{\displaystyle \sum }}{A}_n\left(\omega \mathrm{,}n\right){\Psi }_{\omega \mathrm{,}n}\left(z\right){\text{e}}^{i\left(k_{\omega \mathrm{,}n}x-\omega t+\phi \right)}\text{d}\omega \mathrm{,}$ (21)

where $A_n\left(\omega \mathrm{,}n\right)$ is the amplitude of the n-th mode and $k_{\omega \mathrm{,}n}$ the corresponding wave number.

Equation (20) can also be used to obtain solutions for a wave propagating horizontally with a large enough wavelength to feel a bottom.

When the boundary conditions (18) are used one can readily verify that

$\Psi =Csin\left(mz\right)\mathrm{,}C=\text{const}\mathrm{.}$ (22)

is a solution of (20) and thus there exist an infinite number of discrete internal wave modes satisfying the relationship (16), i.e.

${\omega }_n^2=\frac{k^2{N}^2+m_n^2{f}^2}{m_n^2+k^2},m_{n}H=n\pi ,n=1,2,\cdots$ (23)

Thus in our model (14) we can set

$\psi =A\frac{\omega }{k}sin\left(kx-\omega t+\phi \right)sin\left(mz\right)\mathrm{,}$ (24)

where A is an amplitude of the isopycnal displacement and $\phi$ is a constant phase.

Since N is assumed to be a constant and the coefficients of Equation (14) are independent of spatial variables, we can superimpose solutions of the form (24) to get more general solutions to the $0\left(\epsilon \right)$ problem, i.e.

${\psi }_1=\underset{i}{{\displaystyle \sum }}{A}_i\frac{{\omega }_i}{k_i}\sin \left(k_{i}x-{\omega }_{i}t+{\phi }_i\right)\sin \left(m_{i}z\right),$ (25)

${\rho }_1=\frac{N^2}{g}\underset{i}{{\displaystyle \sum }}{A}_i\sin \left(k_{i}x-{\omega }_{i}t+{\phi }_i\right)\sin \left(m_{i}z\right),$ (26)

$v_1=-f\underset{i}{{\displaystyle \sum }}{A}_i\frac{m_i}{k_i}\cos \left(k_{i}x-{\omega }_{i}t+{\phi }_i\right)\cos \left(m_{i}z\right).$ (27)

2.2. Energy of IWF

To get the formula expressing the total energy of a IWF, we find it useful to rewrite the linearized governing equations combining the Boussinesq approximation (6) in the component form, i.e.

$\frac{\partial u}{\partial t}-fv=-\frac{1}{{\rho }_0}\frac{\partial p}{\partial x},$ (28)

$\frac{\partial v}{\partial t}+fu=-\frac{1}{{\rho }_0}\frac{\partial p}{\partial y},$ (29)

$\frac{\partial w}{\partial t}=-\frac{1}{{\rho }_0}\frac{\partial p}{\partial z}-\frac{{\rho }^{\prime }g}{{\rho }_0},$ (30)

$\frac{\partial \rho }{\partial t}-\frac{N^2{\rho }_0}{g}w=0,$ (31)

$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0,$ (32)

Following [24] we obtain the energy equation of linear internal waves by multiplying (28) by u, (29) by v, (30) by w, and (31) by $g^2{\rho }^{\prime }/{\rho }_0{N}^2$, then adding the results. With the use of the continuity equation this gives

$\frac{\partial }{\partial t}\left\{\frac{1}{2}{\rho }_0\left(u^2+v^2+w^2\right)+\frac{1}{2}\frac{g^2}{{\rho }_0{N}^2}{{\rho }^{\prime }}^2\right\}+\frac{\partial }{\partial x}\left(pu\right)+\frac{\partial }{\partial y}\left(pv\right)+\frac{\partial }{\partial z}\left(pw\right)=0.$ (33)

From (33) we define the averaged linear energy density E by

$E=\frac{1}{2}{\rho }_0\overline{\left(u^2+v^2+w^2\right)}+\frac{1}{2}\frac{g^2}{{\rho }_0{N}^2}\overline{{{\rho }^{\prime }}^2}\mathrm{,}$ (34)

where an overbar is used to define the mean over a wavelength. Substituting solutions (25)-(27) into (34) we obtain the averaged linear energy density for each wave in the form

$E_i=\frac{1}{8}{A}_i^2\left[{\omega }_i^2\left(1+\frac{m_i^2}{k_i^2}\right)+\frac{N^2}{{\rho }_0}+f^2\frac{m_i^2}{k_i^2}\right],$ (35)

where i is the wave number in our model.

Using dispersion relation (23) we can simplify (35) as

$E_{ni}=\frac{1}{4}{A}_i^2\left(N^2+f^2\frac{m_{ni}^2}{k_i^2}\right).$ (36)

2.3. Reflection of Waves

We approximate a continental slope by a uniform slope $z=\alpha x$. The problem consists of determining the reflected wave off the slope if the components of the impinging wave on the slope are known. For completeness, in this section we summarize that part of the Thorpe’s solution [27] to the posed problem which is essential for our analysis.

Let us associate the incident wave with the wave propagating in the $\left(x\mathrm{,}z\right)$ -plane with the group velocity vector directed down to the right. The incident wave is given by

$\left(u_I\mathrm{,}{w}_I\right)=a_I{\text{e}}^{i{\varphi }_I}\mathrm{,}$ (37)

where the subscript I is used to indicate the component of the incident wave, $\left(u\mathrm{,}w\right)$ are the horizontal and vertical components of the velocity vector, a is a constant amplitude and the phase ${\varphi }_I$ of the incident wave is given by

${\varphi }_I=k_{I}x+m_{I}z-{\omega }_{I}t.$ (38)

Since the normal velocity vanishes on the slope, the boundary condition is

$w_I+w_R=\alpha \left(u_I+u_R\right)\text{on}z=\alpha x\mathrm{,}$ (39)

where the subscript R is used to denote the component of the reflected wave. Since the frequencies of the incident and reflected waves must be the same, ${\omega }_I={\omega }_R=\omega$, the components of the reflected wave are determined as

$k_R=\frac{k_I\left(1+{\eta }_I\right)}{\left(1-{\eta }_I\right)},m_R=-\frac{m_I\left(1+{\eta }_I\right)}{\left(1-{\eta }_I\right)},a_R=-\frac{a_I\left(1+{\eta }_I\right)}{\left(1-{\eta }_I\right)},$ (40)

where ${\eta }_I=\alpha m_I/k_I$ is the aspect ratio.

In general, if the incident wave has a unit amplitude, the reflected wave does not. It is only true ( $\left|a_R\right|=1$ ) for horizontal boundary ( $\alpha =0$ ) or for the limiting cases $m_I/k_I=0$ ( $\omega =N$ ) and $m_I/k_I\to \infty$ ( $\omega =f$ ) which follows from dispersion relation (23). As one can see from (40), the incenedent and reflected IWFs are the same if $\alpha =0$.

3. Visualization of IWF Reflected Wave off the Slope

We study the IWF reflected wave off the slope if the incident wave field is known. The energy (36) of each wave is calculated as follows:

First, according to (22), we define $m_{in}$ if wave i by formula

$m_{in}=\frac{n\pi }{H}\mathrm{,}$ (41)

where $n=1,2,3,\cdots ,n_{\max }$ is a mode number and H is an average ocean depth. Hereafter we set set $n_{\max }=50$ and $H=5\text{km}$.

Next, we subdivide the interval $\left[f\mathrm{,}N\right]$ into frequency bands $\left(f+\left(i-1\right)\Delta \omega \mathrm{,}f+i\Delta \omega \right)$, where the wave number i is the integer $i\in \left[\mathrm{1,}M\right]$ and M is the number of frequency bands. We set $M=100$. For each i we choose ${\omega }_{in}$ such that each harmonic wave i with the frequency ${\omega }_{i,n}=n{\omega }_i$ is within the range of possible frequencies. Thus,

$n=1:\omega \in \left[f,N\right],$ (42)

$n=2:\omega \in \left[2f,N/2\right],$ (43)

$\cdots$ (44)

$n=50:\omega \in \left[50f,N/50\right].$

We choose the wave distribution (42)-(44) for the following reason: Consider, for example the first wave ( $i=1$ ) of mode 2 ( $n=2$ ). From one hand,

${\omega }_{1,2}=2{\omega }_1=2f.$ (45)

From the other hand, the frequency ${\omega }_{\mathrm{1,2}}$ must satisfy the inequality (17), i.e.

$f\le 2{\omega }_{\mathrm{1,2}}\le N\text{or}\frac{f}{2}\le {\omega }_{\mathrm{1,2}}\le \frac{N}{2}\mathrm{.}$ (46)

Since $2f>f/2$, the choice $\omega \in \left[2f\mathrm{,}N/2\right]$ for $n=2$ will satisfies both requirements (45)-(46) and the inequality

$\frac{f}{2}<2f\le {\omega }_{1,2}\le \frac{N}{2}$

is valid.

Finally, we find $k_i$ from the dispersion relation (23) as follows:

$k_i=+m_{in}\sqrt{\frac{{\omega }_{in}^2-f^2}{N^2-{\omega }_{in}^2}}$ (47)

We next define the auxiliary function

$S\left({\omega }_{i,n}\right)=\frac{N^2-n^2{\omega }_{i,n}^2}{n^2{\omega }_{i,n}^2-f^2}$ (48)

and rewrite the energy (36)

$E\left({\omega }_{i,n}\right)=\frac{1}{4}{\left(\frac{1+\alpha \sqrt{S}}{1-\alpha \sqrt{S}}\right)}^2{A}^2\left(N^2+f^{2}S\right),$ (49)

with the corresponding domain for each n. For example

$\begin{array}{l}n=1:S\left({\omega }_{i,1}\right)=\frac{N^2-1^2{\omega }_{i,1}^2}{1^2{\omega }_{i,1}^2-f^2},\\ E\left({\omega }_{i,1}\right)=\frac{1}{4}{\left(\frac{1+\alpha \sqrt{S}}{1-\alpha \sqrt{S}}\right)}^2{A}^2\left(N^2+f^{2}S\right),\omega \in \left[f,N\right],\end{array}$ (50)

$\begin{array}{l}n=50:S\left({\omega }_{i,50}\right)=\frac{N^2-{50}^2{\omega }_{i,1}^2}{{50}^2{\omega }_{i,1}^2-f^2},\\ E\left({\omega }_{i,50}\right)=\frac{1}{4}{\left(\frac{1+\alpha \sqrt{S}}{1-\alpha \sqrt{S}}\right)}^2{A}^2\left(N^2+f^{2}S\right),\omega \in \left[50f,\frac{N}{50}\right].\end{array}$ (51)

Without loss of generality we set $A_i=1$. We note that the function E has a singular point

$\xi \left(n\right)=\sqrt{\frac{f^2+{\alpha }^2{N}^2}{n^2\left(1+{\alpha }^2\right)}}\mathrm{.}$ (52)

Indeed we can rewrite E as

$E\left(\omega \mathrm{,}n\right)={\left(\omega -\sqrt{\frac{f^2+{\alpha }^2{N}^2}{n^2\left(1+{\alpha }^2\right)}}\right)}^{-2}G\left(\omega \mathrm{,}n\right)$ (53)

for the positive function

$G\left(\omega ,n\right)=\frac{1}{4}{\left(1+\alpha \sqrt{S}\right)}^4\left(N^2+f^{2}S\right)$ (54)

so that one can see that for any n the function $E\left(\omega \mathrm{,}n\right)$ has a singularity of order 2 at the point $\xi \left(n\right)$. Since we are interested in the total energy for $\omega \in \left[f\mathrm{,}N\right]$ for a fixed $\alpha$ and any integer

$n\le \frac{1}{1+\alpha }\sqrt{1+{\alpha }^2{\left(\frac{N}{f}\right)}^2}$ (55)

the corresponding energy $E\left({\omega }_{i\mathrm{,}n}\right)$ will have a singularity. Thus the domain of exiting internal wave filed is very limited, as shown in Figure 2. This figure shows the restriction (55) of mode numbers n for increasing values of latitude $\varphi$. In particular, as panel (d) shows, at high latitude $\varphi ={70}^{\circ }$, only waves of fundamental frequency ( $n=1$ ) contribute to the energy mixing in the vicinity of continental slope.

Any averaging of the energy involving $E\left({\omega }_{i\mathrm{,}n}\right)$ for a range of n defined by (55) and involving integration with respect to $\omega$ has to be adjust for the singularities. We suggest to consider a parameter $0<\delta <\frac{1}{2}$ and define the contribution of each frequency band by

$E\left({\omega }_{i\mathrm{,}n}\mathrm{<mo>\; ;\; </mo>}\delta \right)=\frac{1}{\Delta \omega }{\displaystyle {\int }_{f+\left(i-1\right)\Delta }^{f+i\Delta }}{E}^{\delta }\left({\omega }_{i\mathrm{,}n}\right)\text{d}\omega \mathrm{.}$ (56)

For the range of $\delta$ all integrals above are convergent and the quantities $E\left({\omega }_{i\mathrm{,}n}\mathrm{<mo>\; ;\; </mo>}\delta \right)$ are finite positive constants. To account for $\delta$ we define the total weigthed energy through $L_{\delta }$ vector norm,

$E^{\text{tot}}\left({\omega }_{i\mathrm{,}n}\mathrm{<mo>\; ;\; </mo>}\delta \right)=\frac{1}{R_i}{\left({\displaystyle \underset{j=1}{\overset{R_j}{\sum }}}E\left({\omega }_{i\mathrm{,}n}\mathrm{<mo>\; ;\; </mo>}\delta \right)\right)}^{\frac{1}{\delta }}\mathrm{,}$ (57)

where $R_i$ is the discrete number of waves in the frequency band i.

Figure 2. The existence domain of continuous IWF for different values of latitude $\varphi$.

3.1. Fundamental Frequency (n = 1)

First, we visualize the dominant contribution to energy density (53) of IWF as a function of frequency when $n=1$. In this case, since higher harmonic waves don’t contribute, we set $\delta =1$. Figure 3 shows the energy density for fixed values of slope $\alpha$ and various latitudes $\varphi =0^{\circ },\varphi ={20}^{\circ },\varphi ={40}^{\circ }$ and $\varphi ={70}^{\circ }$. We can observe the critical slopes ${\alpha }_{\ast }$ at different latitudes. For example, the critical slope at the equator (when latitude $\varphi =0^{\circ }$ ) is ${\alpha }_{\ast }^{\left(0\right)}\approx {0.075}^{\circ }$, i.e. energy is decreasing monotonically for the values of $\alpha <{\alpha }_{\ast }^{\left(0\right)}$. Panel (a) shows the energy density for subcritical slope $\alpha ={0.07}^{\circ }<{\alpha }_{\ast }^{\left(0\right)}$ for $\varphi =0^{\circ }$. Panel (b) shows the energy density for supercritical slope $\alpha ={0.1}^{\circ }>{\alpha }_{\ast }^{\left(0\right)}$ for $\varphi =0^{\circ }$. The energy of low frequency waves is increasing and it’s decreasing for higher frequency waves. However, the value of slope $\alpha ={0.1}^{\circ }$ is still subcritical for higher values of latitude. So we can find numerically the critical slope for latitude $\varphi ={20}^{\circ }$ is ${\alpha }_{\ast }^{\left(20\right)}\approx {0.92}^{\circ }$. Panel (c) shows the energy density for subcritical slope $\alpha ={0.07}^{\circ }<{\alpha }_{\ast }^{\left(20\right)}$ for latitude $\varphi ={20}^{\circ }$. Similarly, the panel (d) shows the energy density for supercritical slope $\alpha ={1.5}^{\circ }>{\alpha }_{\ast }^{\left(20\right)}$ for latitude $\varphi ={20}^{\circ }$. To better visualize the effects of rotation on the values of critical slope on IWF of fundamental frequency, we plot the energy density for fixed values of latitude $\varphi ={40}^{\circ }$ in the vicinity of a critical slope for latitude $\varphi ={40}^{\circ }$, which is ${\alpha }_{\ast }^{\left(40\right)}\approx {1.62}^{\circ }$. Figure 4 shows the energy density for two cases: 1) subcritical slopes $\alpha =1^{\circ }<{\alpha }_{\ast }^{\left(40\right)}$ and $\alpha ={1.61}^{\circ }<{\alpha }_{\ast }^{\left(40\right)}$ and 2) supercritical slopes $\alpha =5^{\circ }>{\alpha }_{\ast }^{\left(40\right)}$ and $\alpha ={10}^{\circ }>{\alpha }_{\ast }^{\left(40\right)}$.

We can can also observe that the energy is a monotonically increasing function of frequency for subcritical values of a continental slope $\alpha$ and the energy has singularities for supercritical values of a continental slope $\alpha$.

Figure 3. Energy density for fixed values of slope $\alpha$ and various latitudes $\varphi =0^{\circ },\varphi ={20}^{\circ },\varphi ={40}^{\circ }$ and $\varphi ={70}^{\circ }$. Panel (a) shows the energy density for subcritical slope $<{\alpha }_{\ast }^{\left(0\right)}$. Panel (b) shows the energy density for supercritical slope $>{\alpha }_{\ast }^{\left(0\right)}$. Panels (c) and (d) show the energy distribution for increasing values of the continental slope $\alpha$.

Figure 4. Illustration of critical slope. Energy density for fixed values of latitude $\varphi ={40}^{\circ }$ in the vicinity of a critical slope slope for latitude $\varphi ={40}^{\circ }$, which is ${\alpha }_{\ast }^{\left(40\right)}\approx {1.62}^{\circ }$.

3.2. Harmonic Frequency (n ≥ 1)

The contribution of higher harmonic waves into energetics of the reflected IWF off a slope is analyzed with $\delta =3$ and visualized in Figure 5. The panels (a) and (b) show the energy density for fixed values of subcritical slope $\alpha$ and various latitudes $\varphi =0^{\circ },\varphi ={20}^{\circ },\varphi ={40}^{\circ }$ and $\varphi ={70}^{\circ }$. The panels (c) and (d) show the energy density for fixed values of latitudes $\varphi ={10}^{\circ }$ and various values of subcritical and supercritical slopes. In particular, the panel (d) shows the dominant contribution to the fundamental vertical mode $n=1$ as a main peak at higher frequency bands and secondary contribution of higher harmonic waves as smaller peaks at lower frequency bands.

Namely, in the top two plots of Figure 5, the effect of latitude is displayed for a fixed slope. As latitude approaches 90˚, the energy increases more rapidly as frequency approaches f. In the bottom two plots, a clear difference between the effects of subcritical and supercritical angles on the behavior of the graph can be observed. In the bottom left plot where the angle is subcritical, the energy decreases monotonically. In the bottom right plot where the angle is supercritical, a large spike in energy can be observed. This situation is visualized in more details using a log scale as shown in Figure 6. This figure is used to illustrate the log of energy density for fixed values of latitude and various values of supercritical slope $\alpha$.

In particular, the upper panel shows the secondary contribution of higher harmonic waves as smaller amplitude spikes of energy at lower frequency bands.

In order to visualize the effects of rotation on harmonic IWF in the vicinity of a supercritical slope $\alpha$, we plot the log of energy density for the fixed value of the supercritical slope $\alpha =5^{\circ }$ and various values of latitude $\varphi ={10}^{\circ }$ and $\varphi ={10}^{\circ }$, as shown in Figure7. We can observe from this Figurethat the main spike in energy is moved to the higher frequency domain if latitude is increased. So we can conclude that the increasing both latitude and the slope transfers the peak of the energy density of harmonic IWF to a higher frequency region.

As follows from (55), for small values of the continental slope ( $\alpha <{\alpha }_{\ast }^{\left(0\right)}$ ), the radiating IWF exists for limited mode numbers n, i.e. if n satisfies the following asymptotic inequality:

$n\le \left(1+\alpha \right)\left[1-\frac{1}{2}{\left(\frac{N}{f}\right)}^2{\alpha }^2+\frac{3}{8}{\left(\frac{N}{f}\right)}^4{\alpha }^4+O\left({\alpha }^6\right)\right]\mathrm{.}$ (58)

Figure 5. The energy density for fixed values of subcritical slope $\alpha$ and various latittudes $\varphi =0^{\circ },\varphi ={20}^{\circ },\varphi ={40}^{\circ }$ and $\varphi ={70}^{\circ }$.

Figure 6. The log of energy density for fixed values of latitude and various values of supercritical slope $\alpha$.

Figure 7. The log of energy density for the fixed value of the supercritical slope $\alpha$ and various values of latitude $\varphi$.

Overall, the energy density of fundamental frequency $n=1$ increases with $\alpha$, reaching a maximum at ${\alpha }_{\ast }^{\left(0\right)}$, the cutoff frequency for the lowest harmonic. It then drops significantly, due to the loss of waves of the fundamental frequency, and then increases as $\alpha$ continues to increase reaching a second maximum (singular frequency) at which point it drops again as waves of second harmonic frequency are lost. The process continues until all the harmonics are lost.

As an additional; illustration to our modeling, the total energy density in $\left(\omega \mathrm{,}n\right)$ and in $\left(k\mathrm{,}n\right)$ -spaces is shown in Figure 8 and Figure 9. Figure 8 shows the log of energy distribution scale in $\left(\omega \mathrm{,}n\right)$ -space for latitude latitude $\varphi ={20}^{\circ }$, for which the critical slope ${\alpha }_{\ast }^{\left(20\right)}\approx {0.92}^{\circ }$. The upper panel shows the energy density for subcritical slope $\alpha =0^{\circ }<{\alpha }_{\ast }^{\left(20\right)}$ and $\alpha ={0.7}^{\circ }<{\alpha }_{\ast }^{\left(20\right)}$. The lower panel shows the energy density for supercritical slope $\alpha =1^{\circ }>{\alpha }_{\ast }^{\left(20\right)}$ and $\alpha =5^{\circ }>{\alpha }_{\ast }^{\left(20\right)}$. We can see that the most energetic part of the spectrum is contained at lower mode numbers and frequencies of $\left(\omega \mathrm{,}n\right)$ -space. Namely, the effects of increasing slope, both approaching and exceeding the critical angle, can be observed in Figure 8 in the harmonic case. After the slope exceeds the critical angle, a large spike in energy occurs at lower mode numbers and frequencies.

Thus the effects of increasing slope, both approaching and exceeding the critical angle, can be observed in Figure 8 in the harmonic case. After the slope exceeds the critical angle, a large spike in energy occurs at lower mode numbers frequencies.

Similarly, Figure 9 shows the log of energy distribution in $\left(k\mathrm{,}n\right)$ -space for latitude latitude $\varphi ={20}^{\circ }$, for which the critical slope ${\alpha }_{\ast }^{\left(20\right)}\approx {0.92}^{\circ }$. The upper panel shows the energy density for subcritical slope $\alpha =0^{\circ }<{\alpha }_{\ast }^{\left(20\right)}$ and $\alpha ={0.7}^{\circ }<{\alpha }_{\ast }^{\left(20\right)}$. The lower panel shows the energy density for supercritical slope $\alpha =1^{\circ }>{\alpha }_{\ast }^{\left(20\right)}$ and $\alpha =5^{\circ }>{\alpha }_{\ast }^{\left(20\right)}$.

We can see that the most energetic part of the spectrum is contained at the higher horizontal wavenumber—higher mode number region of $\left(k\mathrm{,}n\right)$ -space.

Figure 8. The log of energy distrinution in log scale in $\left(\omega \mathrm{,}n\right)$ -space for latitude latitude $\varphi ={20}^{\circ }$.

Figure 9. The log of energy distrinution in log scale in $\left(k\mathrm{,}n\right)$ -space for latitude latitude $\varphi ={20}^{\circ }$.

4. Conclusions

The results of studies in [16] [17] [18] and [20] show that the ocean mixing is most efficient near bottom irregularities. One particular region of ocean mixing occurs near Mid-Atlantic Ridge (see e.g. [18] or [28] ). Our recent numerical results and analytical studies in [19] of internal tide generation over a continental slope confirm the results in [29], showing that the turbulent dissipation is increased due to nonlinear effects and mixing near Mid-Atlantic Ridge. The problem of reflection of internal waves from a plane boundary in different contexts was studied in [27] [28] and [30] [31] [32]. Recent laboratory observations of reflected internal waves on sloping boundaries at a moderately large Reynolds number were reported in [33]. The particular question of reflecting internal wave beams was studied in [34] and [35].

We have developed the numerical algorithm designed to analyze the effects of rotation on the energy distribution of linear IWF in the vicinity of the critical slope. Since waves are generated not only at the fundamental frequency but also at all of its harmonics ${\omega }_n=n\omega$ spanning the range of possible frequencies $\omega \in \left[f\mathrm{,}N\right]$, our analysis includes, in general, an infinite number of discrete internal wave modes n satisfying the dispersion relationship for internal waves. Our algorithm allows to analyze separately the contribution of fundamental frequency and higher harmonic waves on the linear energy distribution. Higher harmonic waves contribute to multiple singularities of the function G defined by (54). So we have introduced the $L_{\beta }$ -norm in order to smooth out the spectral irregularities of the reflected IWF spectrum.

To conclude, our numerical algorithm gave us insight into the effects of continental slope and latitude on energy of a linear reflected internal wave field. In this model, our bathymetry is given by a simple slope, which is used to approximate continental slope. Also, the developed algorithm could be used to expand the existing model to incorporate small amplitude, gentle slope bell-shape topography, which would yield a model much more representative of the ocean floor, and therefore the internal wave fields contributing to mixing in the ocean. This analysis will be reported in the forthcoming studies published elsewhere.

Acknowledgements

We would like to thank the Dean of STEM & Social Sciences at Wenatchee Valley College, Holly Bringman for providing us with all the necessary facilities for the research, great support, and encouragement. This publication is supported by the National Science Foundation (DMS-1555072, DMS-1736364, CMMI-1634832, and CMMI-1560834), and Brookhaven National Laboratory Subcontract 382247.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Hillister, C.D. and McCave, I.N. (1984) Sedimentation under Deep-Sea Storms. Nature, 309, 220-225. https://doi.org/10.1038/309220a0
[2]	Garrett, C.J.R. and Munk, W.H. (1972) Space-Time Scales of Internal Waves. Geophysical Fluid Dynamics, 3,225-264. https://doi.org/10.1080/03091927208236082
[3]	McComas, C.H. (1977) Equilibrium Mechanisms within the Internal Wave Field. Journal of Physical Oceanography, 7, 836-845. https://doi.org/10.1175/1520-0485(1977)007<0836:EMWTOI>2.0.CO;2
[4]	Müller, P. (1977) Spectral Features of the Energy Transfer between Internal Waves and a Larger-Scale Shear Flow. Dynamics of Atmospheres and Oceans, 2, 49-72. https://doi.org/10.1016/0377-0265(77)90015-X
[5]	Zhang, H.P., King, B. and Swinney, H.L. (2008) Resonant Generation of Internal Waves on a Model Continental Slope. Physical Review Letters, 100, Article ID: 244504. https://doi.org/10.1103/PhysRevLett.100.244504
[6]	Bell, T.H. (1975) Lee Waves in Stratified Flows with Simple Harmonic Time Dependence. Journal of Fluid Mechanics, 67, 705-722. https://doi.org/10.1017/S0022112075000560
[7]	Baines, P.G. and Fang, X.H. (1985) Internal Tide Generation at a Continental Shelf/Slope Junction: A Comparison between Theory and a Laboratory Experiment. Dynamics of Atmospheres and Oceans, 9, 297-314. https://doi.org/10.1016/0377-0265(85)90025-9
[8]	Morozov, E.G. (1995) Semidiurnal Internal Wave Global Field. Deep Sea Research Part I: Oceanographic Research Papers, 42, 135-148. https://doi.org/10.1016/0967-0637(95)92886-C
[9]	Shutts, G. (1995) Gravity Wave Drag Parameterization over Complex Terrain: The Effect of Critical Level Absorption in Directional Wind-Shear. Quarterly Journal of the Royal Meteorological Society, 121, 1005-1021. https://doi.org/10.1002/qj.49712152504
[10]	McComas C.H. and Bretherton, F.P. (1977) Resonant Interaction of Oceanic Internal Waves. Journal of Geophysical Research, 82, 1397-1412. https://doi.org/10.1029/JC082i009p01397
[11]	Furue, R. (2003) Energy Transfer within the Small-Scale Oceanic Internal Wave Spectrum. Journal of Physical Oceanography, 33, 267-282. https://doi.org/10.1175/1520-0485(2003)033<0267:ETWTSS>2.0.CO;2
[12]	Cummins, P.F. (1991) The Deep Water Stratification of Ocean General Circulation Models. Atmosphere-Ocean, 29, 563-575. https://doi.org/10.1080/07055900.1991.9649417
[13]	Khatiwala, S. (2003) Generation of Internal Tides in an Ocean of Finite Depth: Analytical and Numerical Calculations. Deep Sea Research Part I: Oceanographic Research Papers, 50, 3-21. https://doi.org/10.1016/S0967-0637(02)00132-2
[14]	Ibragimov, R.N. (2007) Oscillatory Nature and Dissipation of the Internal Wave Energy Spectrum in the Deep Ocean. The European Physical Journal of Applied Physics, 40, 315-334. https://doi.org/10.1051/epjap:2007164
[15]	Munk, W. and Wunsch, C. (1998) Abyssal Recipes II: Energetics of Tidal and Wind Mixing. Deep Sea Research Part I: Oceanographic Research Papers, 45, 1977-2010. https://doi.org/10.1016/S0967-0637(98)00070-3
[16]	Legg, S. and Adcroft, A. (2003) Internal Wave Breaking at Concave and Convex Continental Slopes. Journal of Physical Oceanography, 33, 2224-2246. https://doi.org/10.1175/1520-0485(2003)033<2224:IWBACA>2.0.CO;2
[17]	Wunsch, C. and Ferrari, R. (2004) Vertical Mixing, Energy, and the General Circulation of the Oceans. Annual Review of Fluid Mechanics, 36, 281-314. https://doi.org/10.1146/annurev.fluid.36.050802.122121
[18]	Legg, S. (2004) Internal Tides Generated on a Corrugated Continental Slope. Part II: Along-Slope Barotropic Forcing. Journal of Physical Oceanography, 34, 1824-1838. https://doi.org/10.1175/1520-0485(2004)034<1824:ITGOAC>2.0.CO;2
[19]	Ibragimov, R.N. (2008) Generation of Internal Tides by an Oscillating Background Flow along a Corrugated Slope. Physica Scripta, 78, Article ID: 065801. https://doi.org/10.1088/0031-8949/78/06/065801
[20]	Polzin, K.L., Toole, J.M., Ledwell, J.R. and Schmitt, R.W. (1997) Spatial Variability of Turbulent Mixing in the Abyssal Ocean. Science, 276, 93-96. https://doi.org/10.1126/science.276.5309.93
[21]	Llewellyn Smith, S.G. and Young, W.R. (2002) Conversion of the Barotropic Tide. Journal of Physical Oceanography, 32, 1554-1566. https://doi.org/10.1175/1520-0485(2002)032<1554:COTBT>2.0.CO;2
[22]	Obelcz, J., et al. (2014) Geomorphic Characterization of Four Shelf-Sourced Submarine Canyons along the U.S. Mid-Atlantic Continental Margin. Deep Sea Research Part II: Topical Studies in Oceanography, 104, 106-119. https://doi.org/10.1016/j.dsr2.2013.09.013
[23]	Stoker, J.J. (1957) Water Waves, the Mathematical Theory with Applications. Interscience Publishers, New York.
[24]	Gill, A.E. (1983) Atmosphere-Ocean Dynamics. Academic Press, New York.
[25]	Kundu, P.K. (1990) Fluid Mechanics. Academic Press, New York.
[26]	Osborne, A.R. and Burch, T.L. (1980) Internal Solitons in the Andaman Sea. Science, 208, 451-460. https://doi.org/10.1126/science.208.4443.451
[27]	Thorpe, S.A. (1997) On the Interactions of Internal Waves Reflecting from Slopes. Journal of Physical Oceanography, 27, 2072-2078. https://doi.org/10.1175/1520-0485(1997)027<2072:OTIOIW>2.0.CO;2
[28]	Thorpe, S.A. (2001) On the Reflection of Internal Wave Groups from Sloping Topography. Journal of Physical Oceanography, 31, 3121-3126. https://doi.org/10.1175/1520-0485(2001)031<3121:OTROIW>2.0.CO;2
[29]	Lien, R.C. and Gregg, M.C. (2001) Observations of Turbulence in a Tidal Beam and Across a Costal Ridge. Journal of Geophysical Research: Oceans, 106, 4575-4591. https://doi.org/10.1029/2000JC000351
[30]	Eriksen, C.C. (1982) Observation of Internal Wave Reflection from Sloping Bottoms. Journal of Geophysical Research: Oceans, 87, 525-538. https://doi.org/10.1029/JC087iC01p00525
[31]	Ivey, G.N. and Nokes, R.I. (1989) Vertical Mixing due to Breaking of Critical Waves on Sloping Boundaries. Journal of Fluid Mechanics, 204, 479-500. https://doi.org/10.1017/S0022112089001849
[32]	Slinn, D.N. and Riley, J.J. (1998) Turbulent Dynamics of a Critically Reflecting Internal Gravity Wave. Theoretical and Computational Fluid Dynamics, 11, 281-303. https://doi.org/10.1007/s001620050094
[33]	Gostiaux, L. and Dauxois, T. (2006) Quantitative Laboratory Observations of Internal Wave Reflection on Ascending Slopes. Physics of Fluids, 18, Article ID: 056602. https://doi.org/10.1063/1.2197528
[34]	Gerkema, T., Staquet, C. and Bouruet-Aubertot, P. (2006) Non-Linear Effects in Internal-Tide Beams, and Mixing. Ocean Modelling, 12, 302-318. https://doi.org/10.1016/j.ocemod.2005.06.001
[35]	Tabaei, A., Akylas, T.R. and Lamb, K.G. (2005) Non-Linear Effects in Reflecting and Colliding Internal Wave Beams. Journal of Fluid Mechanics, 526, 217-243. https://doi.org/10.1017/S0022112004002769