Mathematical Modelling of Geophysical Fluid Flow: The Condition for Deep Water Stratification

The deep water stratification can be formulated using various physical and mathematical principles. A common approach is to use the concept of density variation or potential temperature to describe the stratification of the water column. Now for a model formulation for deep water stratification:

The conservation of momentum equation can be used to describe the vertical movement of water masses in response to density differences. This equation can be quite complex and is often solved using a numerical model.

image003.png

Figure 1. Exponential decay in circular motion of deep water

Figure 1 shows the exponential decay of circular trajectory of deep water at the point the vertical amplitude is equal to the horizontal amplitude thereby guaranteeing deep water condition.

image004.png

Figure 2. Geometry of the flow problem in deep water with infinite depth at the thermocline regime

Figure 2 illustrates the zone where the circular orbit dissipates and shows $\eta =\eta \left( x,~y \right)$ representing the progressive surface wave in deep water, with$~\zeta =\zeta \left( x,~y \right)$ defined as the bottom boundary condition.

Consider perturbing the deep water’s surface of amplitude $\eta $($x,y,t$). So $z=\eta $($x,y,t$) is the instantaneous positon of the actual water surface measured from the plane ($z=0$). Where $\eta $ is the surface elevation of the deep water.

3.1 The fluid satisfies Laplacian equation ${{\nabla }^{2}}\varnothing =0$ and continuity equation

$\nabla .\vec{q}=0$                  (1)

There are two boundary condition that will define the regime of deep water stratification and they are free surface and bottom surface. The free surface boundary condition include: Kinematic condition and dynamic condition and the bottom surface condition is at the thermocline regime. Eq. (1) is for momentum conservation which account for the body of deep water and the continuity equation that account for the mass of the deep water.

On the free surface: $pressure={{p}_{atmosphere}}={{p}_{atm}}$

Applying Bernoulli’s equation of motion

$\frac{\partial \phi }{\partial t}+{}^{1}/{}_{2}\left( \phi _{x}^{2}+\phi _{y}^{2} \right)+\frac{P}{\rho }+gy=const=0\left( W.L.G \right)$                   (2)

$\gamma ={{p}_{atm}}~$at $y=\eta \left( x,t \right)$

$\frac{\partial \phi }{\partial t}+{}^{1}/{}_{2}\left( \phi _{x}^{2}+\phi _{y}^{2} \right)+gh=\frac{-{{p}_{atm}}}{\rho }$

Dynamic free surface condition

$y=\eta \left( x,t \right)$

Kinematic condition

$\frac{D}{Dt}\left( y-\eta (x.t \right))=0,~$on $y=\eta \left( x,t \right)$

The equation then simplifies to

$\frac{D}{Dt}=\frac{\partial }{\partial t}+u\frac{\partial }{\partial x}+v\frac{\partial }{\partial y}$

$\frac{D}{Dt}\left( y-\eta  \right),\eta =\eta \left( x,t \right)\frac{\partial \eta }{\partial t}+u\frac{\partial \eta }{\partial x}=v$                         (3)

$y=\eta \left( x,t \right)$ and $u={{\phi }_{x}},v={{\phi }_{y}}$                (4)

$\frac{\partial \phi }{\partial y}=\frac{\partial \eta }{\partial t}+{{\phi }_{x}}\frac{\partial \eta }{\partial x}$

Kinetic free surface condition

$y=\eta \left( x,t \right)$

where, $\phi _{x}^{2},~\phi _{y}^{2},~{{\phi }_{x}},\frac{\partial \eta }{\partial x},\phi \left( x,y,t \right)$ and $\phi \left( x,\eta ,t \right)$ are the two boundary conditions non-linear for 2-dimension.

3.2 Bottom boundary conditions at thermocline regime

In deep water, the thermocline regime as indicated in Figure 2 is associated with the distribution of temperature vertically with depth. So, thermocline is the depth interval where there is a significant change in water temperature, salinity, and other deep water properties. Hence, bottom boundary condition plays very significant role in determining the thermocline regime in a stratified deep water under modified gravity, where, $\eta $ is the free surface elevation and $\phi $ is the surface flux for the flow.

In Figure 2, h is the depth of the water column. Then the substantive derivative or material derivative at the boundary can be expressed as:

$y=h\left( x,t \right),\frac{D}{Dt}\left( y-\eta \left( x,t \right) \right)=0$

$-{{\phi }_{y}}={{h}_{x}}~{{\phi }_{n}}+{{h}_{t}},h\left( x,t \right)=const$

${{\phi }_{y}}=0$ on$\text{ }\!\!~\!\!\text{ }y=-h$, this is the bottom boundary conditions at thermocline regime where gravitational force is different from that on earth giving rise to modified gravity.

The impact of this modified gravity is what our model has sufficiently captured.

Assume that $\eta $ is small, the quantity associated with $\phi $ is equally small

${{\phi }_{x}}/y=\eta ={{\phi }_{x}}/y=0+\eta \frac{\partial {{\phi }_{x}}}{\partial x}/y=0+\frac{{{\eta }^{2}}}{2!}\frac{{{\partial }^{2}}{{\phi }_{x}}}{\partial _{x}^{2}}+\ldots $              (5)

Similarly,

${{\phi }_{y}}/y=\eta ={{\phi }_{y}}/y=0+\eta .\frac{\partial {{\phi }_{y}}}{\partial y}+\frac{{{\eta }^{2}}}{2!}\frac{{{\partial }^{2}}}{\partial _{y}^{2}}\left( {{\phi }_{y}} \right)+\ldots $                    (6)

Substitute for ${{\phi }_{x}},{{\phi }_{y}}~$at$~y=\eta $ as in Eq. (4) in both dynamic and kinematic conditions and neglect the product/higher power terms, which yields:

$\frac{\phi _{x}^{2}}{y}=0,\frac{\phi _{y}^{2}}{y}=0,~\eta .\frac{{{\phi }_{x}}}{y}=0,~{{\eta }^{2}}{{\phi }_{x}}\ldots $              (7)

If the higher power terms are negligible, then the result from the dynamic condition and kinematic condition Eq. (7) become linearized boundary conditions

Dynamic $\rightarrow \phi_t+$ gn $=0$ on $y=0$            (8)

Condition

Kinematic $\rightarrow \eta_t=\phi_y$ on $y=0$          (9)

Now from the above kinematic conditions Eqs. (8) and (9), we obtain

${{\phi }_{tt}}+g{{\eta }_{t}}=0~\text{on }\!\!~\!\!\text{ }y=0$            (10)

${{\phi }_{tt}}+g{{\phi }_{y}}=0~\text{on}~y=0$                 (11)

Free surface boundary condition which is the combination of dynamic condition and kinematic condition.

Eq. (3) necessitate the flow and Eq. (1) which is Laplacian equation is the equation that govern the body of stratified deep water at the thermocline regime.

3.3 Vertical and horizontal amplitude of stratified deep water

At the surface, the vertical displacement $\delta z$ must correspond to a traveling wave

$\delta z=Acos$(k$x$- $\omega $t)

Because vertical displacement of deep water corresponds to the traveling wave due to the principles of wave propagation. In deep water, the wave’s motion is primarily vertical, and this shows that water particles move up and down as this wave passes through. The amount of vertical displacement depends on the wave’s amplitude, which is the maximum distance the water particles move from their resting position.

$\delta z=Acos\left( kx-\omega t \right)$                (12)

$\delta z=-{{A}_{x}}\text{sin}\left( kx-\omega t \right)$             (13)

$\omega =\frac{d}{dt}\left( \delta z \right)={{A}_{z}}\omega \text{sin}\left( kx-\omega t \right)$             (14)

$u=\frac{d}{dt}\left( \delta x \right)={{A}_{x}}\omega \text{cos}\left( kx-\omega t \right)$                     (15)

But the condition that the velocity field is non-divergent and Irrotational arises from conservation of mass in absence of viscosity; hence, Eqs. (16) and (17):

$\frac{du}{dx}+\frac{d\omega }{dz}=0$             (16)

$\frac{d\omega }{dx}-\frac{du}{dz}=0$           (17)

Applying Eq. (14) into Eq. (16) gives

$-{{A}_{z}}k\omega sin~\left( kx-\omega t \right)+\frac{d}{dz}{{A}_{z}}\omega \sin \left( kx-\omega t \right)=0$        (18)

$-k{{A}_{z}}+\frac{d}{dz}{{A}_{z}}=0$        (19)

Using Eq. (15) in Eq. (17), we have

${{A}_{z}}k\omega \cos \left( kx-\omega t \right)-\frac{d}{dz}{{A}_{x}}\omega \cos \left( kx-\omega t \right)=0$           (20)

Differentiating Eq. (19) further gives

$\frac{{{d}^{2}}{{A}_{z}}}{d{{z}^{2}}}-{{k}^{2}}{{A}_{z}}$               (21)

Eq. (21) is a second order differential equation and the solution becomes

${{A}_{z}}=\alpha {{e}^{kz}}+\beta {{e}^{-kz}}$                (22)

where, $\alpha ~$and$~\beta ~$are unknown arbitrary constants and when z$~=0$,

$\Rightarrow {{A}_{z}}=\alpha +\beta =A$                    (23)

When, $z=h$, the bottom condition which is at the thermocline regime where there is no vertical displacement, therefore

$\delta z$($-h$)$=0=\alpha {{e}^{-kh}}+\beta {{e}^{kh}}$             (24)

${{A}_{Z}}=\alpha -\frac{\alpha {{e}^{-kh}}}{{{e}^{kh}}}=\alpha \left( \frac{{{e}^{kh}}-{{e}^{-kh}}}{{{e}^{kh}}} \right)$

Likewise

${{A}_{Z}}=\beta -\frac{\beta {{e}^{kh}}}{{{e}^{-kh}}}=\beta \left( \frac{{{e}^{-kh}}-{{e}^{kh}}}{{{e}^{-kh}}} \right)$            (25)

$\text{or}~\alpha =A~\frac{{{e}^{kh}}}{{{e}^{kh}}-{{e}^{-kh}}},\beta =\frac{-A~{{e}^{-kh}}}{{{e}^{kh}}-{{e}^{-kh}}}$                                (26)

So,

${{A}_{Z}}=\alpha {{e}^{kz}}+\beta {{e}^{kz}}=\frac{A}{{{e}^{kh}}-{{e}^{-kh}}}\left( {{e}^{kh}}{{e}^{kz}}-{{e}^{-kh}}{{e}^{-kz}} \right)=\frac{A{{e}^{k\left( h+z \right)}}-{{e}^{-k\left( h+z \right)}}}{{{e}^{kh}}-{{e}^{-kh}}}$          (27)

${{A}_{Z}}=A\frac{\sin h\left( k\left( h+z \right) \right)}{\sin h\left( kh \right)}$      (28)

Eq. (28) is the vertical amplitude which stands for the vertical movement or oscillation of water particles in deep water. It helps to determine the wave energy and potential impact on coastal structures and marine environment.

Applying Eq. (28) to obtain the horizontal amplitude ${{A}_{x}}$, by substituting solution for ${{A}_{Z}},$ back into Eq. (27), $-KA+\frac{d}{dz}{{A}_{Z}}=0$, gives

$-K{{A}_{x}}+\frac{d}{dz}\left( \frac{A{{e}^{k\left( h+z \right)}}-{{e}^{-k\left( h+z \right)}}}{{{e}^{kh}}-{{e}^{-kh}}} \right)=0,$

$\text{then},~K{{A}_{x}}=AK\left( \frac{A{{e}^{k\left( h+z \right)}}+{{e}^{-k\left( h+z \right)}}}{{{e}^{kh}}-{{e}^{-kh}}} \right)$

So if we take vertical derivative, it yields

${{A}_{x}}=A\frac{\text{cosh}\left( k\left( z+h \right) \right)}{\text{sin}\left( kh \right)}$            (29)

Eq. (29) is the horizontal amplitude of the moving wave which is influenced by the interaction between the waves and stratified water layers. When z becomes negative then $\text{sinh}\left( k\left( z+h \right) \right)$ disappears and the amplitude decreases with depth and when the depth is significantly greater than the wavelength then we have deep water.

Importantly, when Eqs. (27) and (28) are equal then equation of stratified deep water is established at the thermocline regime.

image005.png

Figure 3. Stratification effects in different fluid layers under modified gravity conditions

3.4 Variation of pressure in stratified deep water

In stratified deep water, where multiple layers are stacked atop one another, the pressure within each layer is determined by the hydrostatic approximation. This pressure variation is influenced by several factors, including depth, temperature, salinity, and density gradients. As depicted in Figure 3, these layers possess distinct densities. Consequently, the pressure at any given depth is affected by the density of the layer above it. This results in a layered effect where pressure changes are not uniform across the deep water column.

For the Bosques equation, the hydrostatic approximation is given as

$\frac{dp}{dz}=-gp$                 (30)

The layers in Figure 3 are numbered from top down. The coordinates of the interfaces are denoted by $\eta $, and the layer thickness h_i. So integrating from η₀ to η_z

$\int dp=\mathop{\int }_{{{\eta }_{0}}}^{{{\eta }_{z}}}\left( -g\rho dz \right)\\ {{p}_{1}}=-{{\rho }_{1}}g\left( z-{{\eta }_{0}} \right),~{{p}_{1}}={{\rho }_{1}}g\left( {{\eta }_{0}}-z \right)$       (31)

For the second layer, integrating from$~{{\eta }_{0}}$ to ${{\eta }_{1}}$ and ${{\eta }_{1}}$ to $z$, we have

${{p}_{2}}=-\rho g\mathop{\int }_{{{\eta }_{0}}}^{{{\eta }_{1}}}dz-\rho g\mathop{\int }_{z}^{{{\eta }_{1}}}dz$

${{p}_{2}}={{\rho }_{1}}g\left( {{\eta }_{0}}-{{\eta }_{1}} \right)+{{\rho }_{2}}g\left( {{\eta }_{1}}-z \right)$

$={{\rho }_{1}}g{{\eta }_{0}}-{{\rho }_{1}}g'{{}_{1}}{{\eta }_{1}}+{{\rho }_{2}}g'{{}_{1}}{{\eta }_{1}}-{{\rho }_{2}}gz$

${{p}_{2}}={{\rho }_{1}}g{{\eta }_{0}}+({{\rho }_{2}}-{{\rho }_{1)}}g'{{\eta }_{1}}-{{\rho }_{2}}gz$

At $z=0,$ ${{p}_{2}}gz=0,~$hence

${{p}_{2}}={{\rho }_{1}}g{{\eta }_{0}}+({{\rho }_{2}}-{{\rho }_{1)}}g'{{\eta }_{1}}$                 (32)

From Eq. (32), ${{p}_{1}}={{\rho }_{1}}g{{\eta }_{0}}$ and

${{\rho }_{1}}g'{{}_{1}}{{\eta }_{1}}=g\left( {{\rho }_{2}}-{{\rho }_{1}} \right){{\eta }_{1}}$              (33)

From Eq. (33), we get

$g'{{}_{1}}=g\frac{\left( {{\rho }_{2}}-{{\rho }_{1}} \right)}{{{\rho }_{1}}}=g\frac{\text{ }\!\!\Delta\!\!\text{ }\rho }{{\bar{\rho }}}$                    (34)

Eq. (34) is the modified gravity, where $\frac{\text{ }\!\!\Delta\!\!\text{ }\rho }{{\bar{\rho }}}=\frac{{{\rho }_{i}}_{+1}}{{{\rho }_{i}}}$, $\rho $ is density.

The term involving $z$ is irrelevant for the dynamics, because only the horizontal motion is considered; hence Eqs. (32) and (33) become

${{\begin{matrix}   p  \\ \end{matrix}}_{1}}={{\rho }_{1}}g{{\eta }_{0}},{{p}_{2}}={{\rho }_{1}}g{{\eta }_{0}}+{{\rho }_{1}}g'{{}_{1}}{{\eta }_{1}}$                 (35)

Summing from the top down, we have

${{\eta }_{0}}={{h}_{1}}+{{h}_{2}}+{{\eta }_{b}},~{{\eta }_{1}}={{h}_{2}}+{{\eta }_{b}}$                 (36)

Therefore, putting Eqs. (32) and (33) the pressure in the two layers’ system can be expressed as:

${{\begin{matrix}   p  \\\end{matrix}}_{1}}={{\rho }_{1}}g{{\eta }_{0}}={{\rho }_{1}}g$(${{h}_{1}}+{{h}_{2}}+{{\eta }_{b}}$)               (37)

${{\begin{matrix}   p  \\\end{matrix}}_{2}}={{\rho }_{1}}g{{\eta }_{0}}+{{\rho }_{1}}g'{{}_{1}}{{\eta }_{1}}$$={{\rho }_{1}}g$(${{h}_{1}}+{{h}_{2}}+{{\eta }_{b}}$)$+{{\rho }_{1}}g'{{}_{1}}$(${{h}_{2}}+{{\eta }_{b}}$)                     (38)

Now for the conservation equation of mass, each layer has the same form as the single-layer case and it is expressed as

$\frac{D{{h}_{n}}}{Dt}+{{h}_{n}}\nabla .{{u}_{n}}$                     (39)

Considering the three-layer model, pressure can be given as

${{p}_{1}}={{\rho }_{1}}g{{\eta }_{0}}={{\rho }_{1}}g$(${{h}_{1}}+{{h}_{2}}+{{h}_{3}}+{{\eta }_{b}}$)                (40)

${{p}_{2}}={{\rho }_{1}}g{{\eta }_{0}}+{{\rho }_{1}}{{g'}_{1}}{{\eta }_{1}}={{\rho }_{1}}g$(${{h}_{1}}+{{h}_{2}}+{{h}_{3}}+{{\eta }_{b}}$)+${{\rho }_{1}}{{g'}_{1}}$(${{h}_{2}}+{{h}_{3}}+{{\eta }_{b}}$)                (41)

${{p}_{3}}={{\rho }_{1}}g{{\eta }_{0}}+{{\rho }_{1}}{{g'}_{1}}{{\eta }_{1}}$+${{\rho }_{1}}g'{{}_{2}}{{\eta }_{2}}={{\rho }_{1}}g$(${{h}_{1}}+{{h}_{2}}+{{h}_{3}}+{{\eta }_{b}}$)+${{\rho }_{1}}{{g'}_{1}}$(${{h}_{2}}+{{h}_{3}}+{{\eta }_{b}}$)+${{\rho }_{1}}g'{{}_{2}}$(${{h}_{3}}+{{\eta }_{b}}$)                                 (42)

where, ${{\eta }_{b}}~$is the nth layer of the strata.

Hence,

${{\eta }_{1}}={{\eta }_{b}}+\underset{i=n+1}{\overset{i=N}{\mathop \sum }}\,{{h}_{i}}$            (43)

Therefore, for the nth layer model, the dynamical pressure is given as:

${{p}_{n}}={{\rho }_{i}}\underset{i=0}{\overset{n-1}{\mathop \sum }}\,{{g}_{i}}{{\eta }_{i}}$                (44)

Eq. (39) is a continuity equation in stratified deep water and Eq. (44) is the pressure variation for nth stratification of deep water.

$g'{{}_{i}}=g\frac{\left( {{\rho }_{i+1}}-{{\rho }_{i}} \right)}{{{\rho }_{i}}}=g\frac{\Delta\rho }{{\bar{\rho }}}$          (45)

Eqs. (34) and (45) show the effect of modified gravity on deep water stratification.

It leads to changes in the density and temperature gradients within stratified deep water column.

Taking ${{\rho }_{0}}=0$.

Therefore, incorporating modified gravity into the stratified deep water equations at the thermocline regime become:

$\frac{\partial \left( {{h}_{1}}{{u}_{1}} \right)}{\partial t}+\frac{\partial ({{h}_{1}}u_{1}^{2}+g\frac{\Delta\rho }{{\bar{\rho }}}h_{1}^{2}/2)}{\partial x}+\frac{\partial \left( {{h}_{1}}{{u}_{1}}{{v}_{1}} \right)}{\partial y}=-g\frac{\Delta\rho }{{\bar{\rho }}}{{h}_{1}}\frac{\partial {{h}_{2}}}{\partial x}-g\frac{\Delta\rho }{{\bar{\rho }}}{{h}_{1}}\frac{\partial \left( \xi  \right)}{\partial x}+f{{v}_{1}}$                    (46)

$\frac{\partial \left( {{h}_{1}}{{v}_{1}} \right)}{\partial t}+\frac{\partial \left( {{h}_{1}}{{u}_{1}}{{v}_{1}} \right)}{\partial x}+\frac{\partial ({{h}_{1}}v_{1}^{2}+h_{1}^{2}/2)}{\partial y}=-g\frac{\Delta\rho }{{\bar{\rho }}}{{h}_{1}}\frac{\partial {{h}_{2}}}{\partial y}-g\frac{\Delta\rho }{{\bar{\rho }}}{{h}_{1}}\frac{\partial \left( \xi  \right)}{\partial y}+f{{u}_{1}}$               (47)

$\frac{\partial {{h}_{1}}}{\partial t}+\frac{\partial \left( {{h}_{1}}{{u}_{1}} \right)}{\partial x}+\frac{\partial \left( {{h}_{1}}{{v}_{1}} \right)}{\partial y}=0$             (48)

$\frac{\partial \left( {{h}_{2}}{{u}_{2}} \right)}{\partial t}+\frac{\partial ({{h}_{2}}u_{2}^{2}+g\frac{\Delta\rho }{{\bar{\rho }}}h_{2}^{2}/2+g\frac{\Delta\rho }{{\bar{\rho }}}{{h}_{2}}{{h}_{1}})}{\partial x}+\frac{\partial \left( {{h}_{2}}{{u}_{2}}{{v}_{2}} \right)}{\partial y}=-g\frac{\Delta\rho }{{\bar{\rho }}}~{{h}_{1}}\frac{\partial {{h}_{2}}}{\partial x}-g\frac{\Delta\rho }{{\bar{\rho }}}{{h}_{2}}\frac{\partial \left( \xi  \right)}{\partial x}+f{{v}_{2}}$                 (49)

$\frac{\partial \left( {{h}_{2}}{{v}_{2}} \right)}{\partial t}+\frac{\partial \left( {{h}_{2}}{{u}_{2}}{{v}_{2}} \right)}{\partial x}+\frac{\partial ({{h}_{2}}v_{2}^{2}+g\frac{\Delta\rho }{{\bar{\rho }}}h_{2}^{2}/2+g\frac{\Delta\rho }{{\bar{\rho }}}{{h}_{2}}{{h}_{1}})}{\partial y}=-g\frac{\Delta\rho }{{\bar{\rho }}}~{{h}_{1}}\frac{\partial {{h}_{2}}}{\partial y}-g\frac{\Delta\rho }{{\bar{\rho }}}{{h}_{2}}\frac{\partial \left( \xi  \right)}{\partial y}+f{{u}_{2}}$                (50)

$\frac{\partial {{h}_{2}}}{\partial t}+\frac{\partial \left( {{h}_{2}}{{u}_{2}} \right)}{\partial x}+\frac{\partial \left( {{h}_{2}}{{v}_{2}} \right)}{\partial y}=0$                   (51)

We can now provide solution to the model equation in Eqs. (46)-(51) by considering the characteristic equations in matrix form.

In conservation form,

$\frac{\partial q}{\partial t}+\frac{\partial F}{\partial x}+\frac{\partial G}{\partial y}=\psi $        (52)

where,

$q=\left[ \begin{matrix}   {{q}_{1}}  \\   {{q}_{2}}  \\   {{q}_{3}}  \\ \end{matrix} \right]=\left[ \begin{matrix}   h  \\   hu  \\   hv  \\ \end{matrix} \right],\\ F=\left[ \begin{matrix}  {{q}_{2}}  \\   \frac{q_{2}^{2}}{{{q}_{1}}}+\frac{g\frac{\Delta\rho }{{\bar{\rho }}}q_{1}^{2}}{2}  \\   \frac{{{q}_{2}}{{q}_{3}}}{{{q}_{1}}}  \\ \end{matrix} \right]=\left[ \begin{matrix}   hu  \\   h{{u}^{2}}+\frac{g\frac{\Delta\rho }{{\bar{\rho }}}{{h}^{2}}}{2}  \\   huv  \\ \end{matrix} \right],\\G=\left[ \begin{matrix}   {{q}_{3}}  \\   \frac{{{q}_{2}}{{q}_{3}}}{{{q}_{1}}}  \\   \frac{{{q}_{3}}}{{{q}_{1}}}+\frac{g\frac{\Delta\rho }{{\bar{\rho }}}q_{1}^{2}}{2}  \\\end{matrix} \right]=\left[ \begin{matrix}   hv  \\   huv  \\   h{{v}^{2}}+\frac{g\frac{\Delta\rho }{{\bar{\rho }}}{{h}^{2}}}{2}  \\ \end{matrix} \right],\\ \psi =\left[ \begin{matrix}   0  \\   -g\frac{\Delta\rho }{{\bar{\rho }}}h{{\xi }_{x}}+fv  \\   -g\frac{\Delta\rho }{{\bar{\rho }}}h{{\xi }_{y}}-fu  \\ \end{matrix} \right]$                        (53)

The subscripts $x,y$ denote differentiation with respect to that variable. In vector form, the unknown is $q={{\left[ h;hu;hv \right]}^{T}}$ is the vector of conserved variables.

The right-hand side vector of source terms, ψ, includes the effects of density stratification at the transition zone (thermocline) and the effects of the Coriolis force. To calculate the eigenvalues of the deep water equations, define the Jacobian matrix with differentiated coefficients of $F\left( q \right)={{\left[ {{f}_{1}};{{f}_{2}};{{f}_{3}} \right]}^{T}}$ as $F'\left( q \right)=\frac{\partial {{f}_{1}}}{\partial {{q}_{1}}}$ for $i,j=1,2,3$.

$F'\left( q \right)=\left[ \begin{matrix}   \begin{matrix}   \frac{\partial {{f}_{1}}}{\partial {{q}_{1}}} & \frac{\partial {{f}_{1}}}{\partial {{q}_{2}}} & \frac{\partial {{f}_{1}}}{\partial {{q}_{3}}}  \\   \frac{\partial {{f}_{2}}}{\partial {{q}_{1}}} & \frac{\partial {{f}_{2}}}{\partial {{q}_{2}}} & \frac{\partial {{f}_{2}}}{\partial {{q}_{3}}}  \\   \frac{\partial {{f}_{3}}}{\partial {{q}_{1}}} & \frac{\partial {{f}_{3}}}{\partial {{q}_{2}}} & \frac{\partial {{f}_{3}}}{\partial {{q}_{3}}}  \\\end{matrix}  \\\end{matrix} \right]$

And an equivalent expression for

$G'\left( q \right){{q}_{1}}+F'\left( q \right)\frac{\partial q}{\partial x}+G'\left( q \right)\frac{\partial q}{\partial y}=\psi $              (54)

where,

$F'\left( q \right)=\left[ \begin{matrix}   0 & 1 & 0  \\   g\frac{\Delta\rho }{{\bar{\rho }}}h-{{u}^{2}} & 2u & 0  \\   -uv & v & u  \\\end{matrix} \right]$                 (55)

$G'\left( q \right)=\left[ \begin{matrix}   0 & 0 & 1  \\   -uv & v & u  \\   g\frac{\Delta\rho }{{\bar{\rho }}}h-{{v}^{2}} & 0 & 2v  \\\end{matrix} \right]$ and$\text{ }\!\!~\!\!\text{ }F'\left( q \right).x=\lambda x$

That is

($F'\left( q \right)-\lambda I)x=0$                 (56)

The characteristics determinant:

$F\left( q \right)-\lambda I=\left| g\frac{\Delta\rho }{{\bar{\rho }}}\begin{matrix}   -\lambda  & 1 & 0  \\   h-{{u}^{2}} & 2u-\lambda  & 0  \\   -uv & v & u-\lambda   \\\end{matrix} \right|$                    (57)

Will lead to the characteristic equation:

$F\left( q \right)-\lambda I=\left| \begin{matrix}   -\lambda  & 1 & 0  \\   g\frac{\Delta\rho }{{\bar{\rho }}}h-{{u}^{2}} & 2u-\lambda  & 0  \\   -uv & v & u-\lambda   \\ \end{matrix} \right|=0$              (58)

Solving the 3 by 3 matrix in Eq. (45),

$-\lambda \left\{ \left( 2u-\lambda  \right)\left( u-\lambda  \right) \right\}-1\left\{ \left( g\frac{\Delta\rho }{{\bar{\rho }}}h-{{u}^{2}} \right)\left( u-\lambda  \right) \right\}=0$

$\left( u-\lambda  \right)\left( {{\lambda }^{2}}-2u\lambda -\left( g\frac{\Delta\rho }{{\bar{\rho }}}h-{{u}^{2}} \right) \right)=0$

                  $\left( u-\lambda  \right)=0,\lambda =u$

$\text{or}~\left( {{\lambda }^{2}}-2u\lambda -\left( g\frac{\Delta\rho }{{\bar{\rho }}}h-{{u}^{2}} \right) \right)=0$

Solving using the quadratic formula: $\lambda =\frac{-b\pm \sqrt{{{b}^{2}}-4ac}}{2a}$, where, $a=1,b=-2u$ and$~c=-\left( g\frac{\Delta\rho }{{\bar{\rho }}}h-{{u}^{2}} \right),$ substituting these values in the formula; we obtain:

$\lambda =u\pm \sqrt{g\frac{\Delta\rho }{{\bar{\rho }}}h}$                (59)

This leads us to the eigenvalues,$~$

${{\lambda }_{1}}\left( q \right)=u-\sqrt{g\frac{\Delta\rho }{{\bar{\rho }}}h},~{{\lambda }_{2}}\left( q \right)=u+\sqrt{g\frac{\Delta\rho }{{\bar{\rho }}}h},~{{\lambda }_{3}}\left( q \right)=u$                   (60)

The following eigenvectors are obtained using Maple Software;

For ${{\lambda }_{1}}\left( q \right)$ we have:

${{\lambda }_{1}}\left( q \right)=\left[ \begin{matrix}   \frac{\sqrt{g\frac{\Delta\rho }{{\bar{\rho }}}h}\left( g\frac{\Delta\rho }{{\bar{\rho }}}h-{{u}^{2}} \right)}{uv\left( -u-\sqrt{g\frac{\Delta\rho }{{\bar{\rho }}}h} \right)(u-\sqrt{g\frac{\Delta\rho }{{\bar{\rho }}}h)}}  \\   \frac{\sqrt{g\frac{\Delta\rho }{{\bar{\rho }}}h}\left( g\frac{\Delta\rho }{{\bar{\rho }}}h-{{u}^{2}} \right)}{uv\left(-u-\sqrt{g\frac{\Delta\rho }{{\bar{\rho }}}h}\right)} \\   1  \\\end{matrix} \right]\\{{\lambda }_{2}}\left( q \right)=\left[ \begin{matrix}   \frac{-\sqrt{g\frac{\Delta\rho }{{\bar{\rho }}}h}\left( g\frac{\Delta\rho }{{\bar{\rho }}}h-{{u}^{2}} \right)}{uv\left( -u+\sqrt{g\frac{\Delta\rho }{{\bar{\rho }}}h} \right)(u+\sqrt{g\frac{\Delta\rho }{{\bar{\rho }}}h)}}  \\   \frac{-\sqrt{g\frac{\Delta\rho }{{\bar{\rho }}}h}\left( g\frac{\Delta\rho }{{\bar{\rho }}}h-{{u}^{2}} \right)}{uv\left(-u-\sqrt{g\frac{\Delta\rho }{{\bar{\rho }}}h}\right)}  \\   1  \\\end{matrix} \right]\\{{\lambda }_{3}}\left( q \right)=\left[ \begin{matrix}   0  \\   0  \\   1  \\\end{matrix} \right]$              (61)

Similarly, for$~G'\left( q \right)~$

${{\lambda }_{1}}\left( q \right)=v-\sqrt{g\frac{\Delta\rho }{{\bar{\rho }}}h},~{{\lambda }_{2}}\left( q \right)=v+\sqrt{g\frac{\Delta\rho }{{\bar{\rho }}}h},~{{\lambda }_{3}}\left( q \right)=v$                         (62)

With the corresponding eigenvectors,

${{\chi }_{1}}\left( q \right)=\left[ \begin{matrix}   \frac{1}{v-\sqrt{g\frac{\Delta\rho }{{\bar{\rho }}}h}}  \\   \frac{-ug\frac{\Delta\rho }{{\bar{\rho }}}h+u{{v}^{2}}+uv\left( v-\sqrt{g\frac{\Delta\rho }{{\bar{\rho }}}h} \right)-2u{{v}^{2}}}{\left( g\frac{\Delta\rho }{{\bar{\rho }}}h-{{v}^{2}} \right)\sqrt{g\frac{\Delta\rho }{{\bar{\rho }}}h}}  \\   1  \\ \end{matrix} \right]$         (63)

${{\chi }_{2}}\left( q \right)=\left[ \begin{matrix}   \frac{1}{v+\sqrt{g\frac{\Delta\rho }{{\bar{\rho }}}h}}  \\   \frac{-ug\frac{\Delta\rho }{{\bar{\rho }}}h+u{{v}^{2}}+uv\left( v+\sqrt{g\frac{\Delta\rho }{{\bar{\rho }}}h} \right)-2u{{v}^{2}}}{\left( g\frac{\Delta\rho }{{\bar{\rho }}}h-{{v}^{2}} \right)\sqrt{g\frac{\Delta\rho }{{\bar{\rho }}}h}}  \\   1  \\ \end{matrix} \right]$             (64)

${{\chi }_{3}}\left( q \right)=\left[ \begin{matrix}   0  \\   1  \\   0  \\\end{matrix} \right]$      (65)

Since this system has two spatial dimensions for Eqs. (59) and (60), the waves (solutions) move in horizontal direction establishing the vertical amplitude of the wave and the horizontal amplitude which is given by the unit normal vector $n~$where $n\left( {{n}_{x}};{{n}_{y}} \right).$

$|{{n}_{x}}F'\left( q \right)+{{n}_{y}}G'\left( q \right)-{{\lambda }_{i}}I|=0$ for$~i=1,2,3$           (66)

that yield

${{\lambda }_{1}}=\frac{{{n}_{x}}{{q}_{2}}}{{{q}_{1}}}+\frac{{{n}_{y}}{{q}_{3}}}{{{q}_{1}}}-\sqrt{g\frac{\Delta\rho }{{\bar{\rho }}}{{q}_{1}}},\\{{\lambda }_{2}}=\frac{{{n}_{x}}{{q}_{2}}}{{{q}_{1}}}+\frac{{{n}_{y}}{{q}_{3}}}{{{q}_{1}}},{{\lambda }_{3}}=\frac{{{n}_{x}}{{q}_{2}}}{{{q}_{1}}}+\frac{{{n}_{y}}{{q}_{3}}}{{{q}_{1}}}+\sqrt{g\frac{\Delta\rho }{{\bar{\rho }}}{{q}_{1}}}$            (67)

${{\lambda }_{1}}\left( q \right)=u+v-\sqrt{g\frac{\Delta\rho }{{\bar{\rho }}}h},~{{\lambda }_{2}}\left( q \right)=u+v,~\\{{\lambda }_{3}}\left( q \right)=u+v+\sqrt{g\frac{\Delta\rho }{{\bar{\rho }}}h}$               (68)

where

$c=\sqrt{g\frac{\Delta\rho }{{\bar{\rho }}}h}$           (69)

Eq. (69) provides the speed of stratified deep water at equilibrium position which is when the vertical amplitude is equal to the horizontal amplitude.

The characteristic equations/eigenvalues obtained in Eq. (68) are very useful in deep water analysis as it tells us the behavior of the deep water waves. The eigenvalues help us to predict the growth or decay rates of wave modes in deep water as demonstrated in Figure 1.

The Froude Number is a dimensionless quantity mainly used in fluid dynamics in defining the nature of boundary layer. Fundamentally,roude Number is defined as the ratio of the inertial force to gravitational force acting on a fluid.

The eigenvalues take the form of a convective velocity minus/plus a phase velocity. Dividing convective by phase yields a Froude number for each x- and y- direction, therefore, ${{F}_{r}}~$is Froude number which can be expressed as:

${{F}_{\left( r \right)i}}\equiv \frac{{{U}_{x}}}{\sqrt{\left( g\frac{\Delta\rho }{{\bar{\rho }}} \right)h}}=\frac{{{V}_{y}}}{\sqrt{g\frac{\text{ }\!\!\Delta\!\!\text{ }\rho }{{\bar{\rho }}}}}$                  (70)

Eq. (70) is the Froude number expressed in terms of vertical and horizontal velocities respectively.

Now in x-direction, the Froude number between two stratified layers can be expressed as:

${{F}_{r}}=\frac{u}{\sqrt{g\frac{\left( {{\rho }_{2}}-{{\rho }_{1}} \right)}{{{\rho }_{1}}}h}}=\frac{u}{c}$              (71)

Eq. (71) is the dimensionless quantity which helps to determine the stability and propagation of internal waves which occurs at the interface between the two layers.

$q=\left[ \begin{matrix}   {{h}_{1}}  \\   {{h}_{1}}{{u}_{1}}  \\   {{h}_{1}}{{v}_{1}}  \\   {{h}_{2}}  \\   {{h}_{2}}{{u}_{2}}  \\   {{h}_{2}}{{v}_{2}}  \\\end{matrix} \right]\\S\left( q \right)=\left[ \begin{matrix}   0  \\   -\frac{g\Delta\rho }{{\bar{\rho }}}{{h}_{1}}\frac{\partial {{h}_{2}}}{\partial x}-\frac{g\Delta\rho }{{\bar{\rho }}}{{h}_{1}}\frac{\partial \xi }{\partial x}+f{{v}_{1}}  \\   -\frac{g\Delta\rho }{{\bar{\rho }}}{{h}_{1}}\frac{\partial {{h}_{2}}}{\partial y}-\frac{g\Delta\rho }{{\bar{\rho }}}{{h}_{1}}\frac{\partial \xi }{\partial y}-f{{u}_{1}}  \\   0  \\   -\frac{g\Delta\rho }{{\bar{\rho }}}{{h}_{1}}\frac{\partial {{h}_{2}}}{\partial x}-\frac{g\Delta\rho }{{\bar{\rho }}}{{h}_{2}}\frac{\partial \xi }{\partial x}+f{{v}_{2}}  \\   -\frac{g\Delta\rho }{{\bar{\rho }}}{{h}_{1}}\frac{\partial {{h}_{2}}}{\partial y}-\frac{g\Delta\rho }{{\bar{\rho }}}{{h}_{2}}\frac{\partial \xi }{\partial y}-f{{u}_{2}}  \\ \end{matrix} \right]$

$f\left( q \right)=\left[ \begin{matrix}   {{h}_{1}}{{u}_{1}}  \\   {{h}_{1}}u_{1}^{2}+\frac{1}{2}\frac{g\Delta\rho }{{\bar{\rho }}}h_{1}^{2}  \\   {{h}_{1}}{{u}_{1}}{{v}_{1}}  \\   {{h}_{2}}{{u}_{2}}  \\   {{h}_{2}}u_{2}^{2}+\frac{1}{2}\frac{g\Delta\rho }{{\bar{\rho }}}h_{2}^{2}+\frac{g\Delta\rho }{{\bar{\rho }}}{{h}_{2}}{{h}_{1}}  \\   {{h}_{2}}{{u}_{2}}{{v}_{2}}  \\\end{matrix} \right]$.

$g\left( q \right)=\left[ \begin{matrix}   {{h}_{1}}{{v}_{1}}  \\   {{h}_{1}}{{u}_{1}}{{v}_{1}}  \\   {{h}_{1}}v_{1}^{2}+\frac{1}{2}\frac{g\Delta\rho }{{\bar{\rho }}}h_{1}^{2}  \\   {{h}_{2}}{{v}_{2}}  \\   {{h}_{2}}{{u}_{2}}{{v}_{2}}  \\   {{h}_{2}}v_{2}^{2}+\frac{1}{2}\frac{g\Delta\rho }{{\bar{\rho }}}h_{2}^{2}+\frac{g\Delta\rho }{{\bar{\rho }}}{{h}_{2}}{{h}_{1}}  \\ \end{matrix} \right]$

We can then calculate, the flux Jacobians of $f\left( q \right)$ as:

$f^{\prime}\left( q \right)=\left[ \begin{matrix}   0 & 1 & 0 & 0 & 0 & 0  \\   -u_{1}^{2}+\frac{g\Delta\rho }{{\bar{\rho }}}{{h}_{1}} & 2{{u}_{1}} & 0 & 0 & 0 & 0  \\   -{{u}_{1}}{{v}_{1}} & {{v}_{1}} & {{u}_{1}} & 0 & 0 & 0  \\   0 & 0 & 0 & 0 & 1 & 0  \\   \frac{g\Delta\rho }{{\bar{\rho }}}{{h}_{2}} & 0 & 0 & -u_{2}^{2}+\frac{g\Delta\rho }{{\bar{\rho }}}{{h}_{2}}+\frac{g\Delta\rho }{{\bar{\rho }}}{{h}_{1}} & 2{{u}_{2}} & 0  \\   0 & 0 & 0 & -{{u}_{2}}{{v}_{2}} & {{v}_{2}} & {{u}_{2}}  \\ \end{matrix} \right]$                (72)

$g'\left( q \right)=\left[ \begin{matrix} 0 & 0 & 1 & 0 & 0 & 0  \\   -{{u}_{1}}{{v}_{1}} & {{v}_{1}} & {{u}_{1}} & 0 & 0 & 0  \\-v_{1}^{2}+\frac{g\Delta\rho }{{\bar{\rho }}}{{h}_{1}} & 0 & 2{{v}_{1}} & 0 & 0 & 0  \\   0 & 0 & 0 & 0 & 0 & 0  \\   0 & 0 & 0 & -{{u}_{2}}{{v}_{2}} & {{v}_{2}} & {{u}_{2}}  \\   \frac{g\Delta\rho }{{\bar{\rho }}}{{h}_{2}} & 0 & 0 & -u_{2}^{2}+\frac{g\Delta\rho }{{\bar{\rho }}}{{h}_{2}}+\frac{g\Delta\rho }{{\bar{\rho }}}{{h}_{1}} & 0 & 2{{v}_{2}}  \\ \end{matrix} \right]$             (73)

3.5 Calculating the eigenspace

A key step in solving multi-layer deep water equations is calculating the Eigen space of the system, which consists of two sets: eigenvalues and eigenvectors, one set corresponds to slightly perturb classical deep water gravity waves, while the other set corresponds to internal waves traveling at a much lower speed. By the use of the Maple Software, we have the following eigenvalues and eigenvectors, for $A\left( q \right):$