Weakly Nonlinear Wave Theory for Periodic Waves (Stokes Expansion)

From WikiWaves
Jump to: navigation, search


The solution for Stokes waves is valid in deep or intermediate water depth. It is assumed that the wave steepness is much smaller than one.

[math] \varepsilon = a k \ll 1 \, [/math]

[math] k h = O (1) \, [/math]

where [math] k \, [/math] is the wavenumber and [math] h \, [/math] is the water depth which is assumed constant.

  • Nondimensional Variables

[math] X = x k, \quad Z = z k, \quad Y = y k, \, [/math]

[math] \bar{t} = t \sigma, \quad \bar{\eta} = \eta / a, \, [/math]

[math] \Phi = \frac{\sigma\phi}{a g}, \quad \bar{C} = \frac{C}{a g}, \quad D = \frac{g k}{\sigma^2}, \, [/math]

where [math] x, \ z, \ t, \ h, \ \phi \, [/math] and [math] C \, [/math] are dimensional variables and [math] X, \ Z, \ \bar{t}, \ \bar{\eta}\, \ \Phi \, [/math] and [math] \bar{C} \, [/math] are corresponding nondimensional variables.

Nondimensional Governing Equation & Boundary Conditions

[math] \frac{\partial^2\theta}{\partial X^2} + \frac{\partial^2\phi}{\partial Y^2} + \frac{\partial^2\theta}{\partial Z^2} = 0 \, [/math] [math] - k h \lt Z \lt \varepsilon \bar{\eta} \, [/math] [math] (3.1.1) \, [/math]
[math] \frac{\partial\Phi}{\partial Z} = 0 \, [/math] at [math] Z = - k h \, [/math] [math] (3.1.2) \, [/math]
[math] \frac{\partial\bar{\eta}}{\partial\bar{t}} + \varepsilon D \nabla_h \Phi \cdot \nabla_h \bar{\eta} = D \frac{\partial\Phi}{\partial Z} [/math] at [math] Z = \varepsilon \bar{\eta} \, [/math] [math] (3.1.3) \, [/math]
[math] \frac{\partial\Phi}{\partial\bar{t}} + \frac{1}{2} \varepsilon D [ \nabla \Phi ]^2 = \bar{C} \, [/math] at [math] Z = \varepsilon \bar{\eta} \, [/math] [math] (3.1.4) \, [/math]

where [math] \nabla \, [/math] and [math] \nabla_h \, [/math] stand for gradient and horizontal gradient, respectively.

Perturbation (Stokes Expansion)

Assuming the wave train is weakly nonlinear [math] ( \varepsilon = a k \ll 1 ) \, [/math], its potential and elevation can be perturbed in the order of [math] \varepsilon \, [/math].

[math] \Phi = \Phi^{(1)} + \varepsilon\Phi^{(2)} + \varepsilon^2\Phi^{(3)} + \ldots + \varepsilon^{j-1}\Phi^{(j)} + \ldots [/math]

[math] \bar{\eta} = {\bar{\eta}}^{(1)} + \varepsilon {\bar{\eta}}^{(2)} + \varepsilon^2 {\bar{\eta}}^{(3)} + \ldots + \varepsilon^{j-1} {\bar{\eta}}^{(j)} + \dots [/math]

[math] \bar{C} = {\bar{C}}^{(1)} + \varepsilon {\bar{C}}^{(2)} + \varepsilon^2 {\bar{C}}^{(3)} + \ldots + \varepsilon^{j-1} {\bar{C}}^{(j)} + \dots [/math]

Hierarchy Equations

Using the Taylor expansion, the free-surface boundary conditions (Equations(3.1.3) and (3.1.4) are expanded at the still water level [math] (Z=0) \, [/math]. Then we substitute perturbation forms of potential and elevation into the Laplace Equation, bottom and free-surface boundary conditions. The equations are sorted and grouped according to the order in wave steepness [math] \varepsilon^{(j)} \, [/math]. The governing equations for [math] j-th \, [/math] order solutions is given by:

[math] \nabla^2\Phi^{(j)}=0 \quad -h k \leq Z \leq 0 \qquad (3.1.5) [/math]

[math] \bar{\eta}^{(j)} + \frac{\partial{\Phi}^{(j)}}{\bar{\partial{t}}} = P^{(j)} \left[ \Phi^{(j-1)}, {\bar{\eta}}^{(j-1)} \right] \ \mbox{at} \quad Z = 0 \qquad (3.1.6) [/math]

[math] \frac{{\bar{\partial{\eta}}}^{(j)}}{\bar{\partial{t}}} - D \frac{{\partial{\Phi}}^{(j)}}{\partial{Z}} = Q^{(j)} \left[ \Phi^{(j-1)}, {\bar{\eta}}^{(j-1)} \right] \quad \mbox{at} \quad Z = 0 \quad (3.1.7) [/math]

[math] \frac{\partial\Phi^{(j)}}{\partial Z} = 0 \qquad \mbox{at} \quad Z = - k h \qquad (3.1.8) [/math]

where the [math] P^{(j)} \, [/math] and [math] Q^{(j)} \, [/math] can be derived in terms of the solutions for the potential and elevation of order [math] (j-1) \, [/math] or lower. Therefore, the above hierarchy equations must be solved sequentially from lower to higher order until the required accuracy is reached. To derive the third-order solution for a Stokes wave train, it is adequate to truncate the equations at [math] j = 3 \, [/math].

Up to [math] j = 3, \ P^{(j)} \, [/math] and [math] Q^{(j)} \, [/math] are given below.

[math] P^{(1)} = {\bar{C}}^{(1)} \quad \mbox{and} \quad Q^{(1)} = 0 [/math]

[math] P^{(2)} = - {\bar{\eta}}^{(1)} \frac{\partial^2\Phi^{(1)}}{\partial{Z}\partial{\bar{t}}} - \frac{D}{2} \left| \nabla\Phi^{(1)} \right|^2 + {\bar{C}}^{(2)} [/math]

[math] Q^{(2)} = D {\bar{\eta}}^{(1)} \frac{\partial^2\Phi^{(1)}}{\partial{Z}^2} - D \nabla_h \Phi^{(1)} \cdot \nabla_h {\bar{\eta}}^{(1)} [/math]

[math] P^{(3)} = - D \nabla\Phi^{(1)} \cdot \nabla\Phi^{(2)} - {\bar{\eta}}^{(1)} \frac{\partial}{\partial{Z}} \left[ \frac{D}{2} \left| \nabla\Phi^{(1)} \right|^2 + \frac{\partial\Phi^{(2)}}{\partial\bar{t}} \right] [/math]

[math] - {\bar{\eta}}^{(2)} \frac{\partial^2\Phi^{(1)}}{\partial{Z}\partial{\bar{t}}} - \frac{1}{2} \left( {\bar{\eta}}^{(1)} \right)^2 \frac{\partial^3\Phi^{(1)}}{\partial{Z}^2 \partial{\bar{t}}} + {\bar{C}}^{(3)} [/math]

[math] Q^{(3)} = D {\bar{\eta}}^{(2)} \frac{\partial^2\Phi^{(1)}}{\partial{Z}^2} + \frac{D}{2} \left( {\bar{\eta}}^{(1)} \right)^2 \frac{\partial^3\Phi^{(1)}}{\partial{Z}^3} [/math]

[math] +D {\bar{\eta}}^{(1)} \frac{\partial}{\partial{Z}} \left( \frac{\partial\Phi^2}{\partial{Z}} - \nabla_h \Phi^{(1)} \cdot \nabla_h {\bar{\eta}}^{(1)} \right) [/math]
  • Solving the non-dimensional Equations from lower order [math] (j=1) \, [/math] to higher order [math] (j=3) \, [/math] for the non-dimensional solutions ( wave advances in the x-direction ).

[math] \Phi^{(1)} = \frac{ \cosh ( k h + Z ) }{\cosh k h} \sin ( X - \bar{t} ), [/math]

[math] {\bar{\eta}}^{(1)} = \cos ( X - \bar{t} ), \quad {\bar{C}}^{(1)} = 0 [/math]

[math] \Phi^{(2)} = \frac{3}{8} \frac{\cosh (2kh+2Z)}{\sinh^3kh\cosh kh} \sin(2X-2\bar{t}) [/math]

[math] {\bar{\eta}}^{(2)} = \frac{\alpha}{4} ( 3\alpha^2 - 1 ) \cos ( 2X - 2 \bar{t} ) \, [/math]

[math] {\bar{C}}^{(2)} = \frac{1}{2\sinh 2 kh} \, [/math]

[math] \Phi^{(3)} = \frac{1}{64} (\alpha^2-1) (\alpha^2+3) (9\alpha^2-13) \cdot \frac{\cosh(3kh+3Z)}{\cosh kh} \sin (3X - 3\bar{t}) [/math]

[math] {\bar{\eta}}^{(3)} = \frac{3}{8} (\alpha^4 - 3\alpha^2 + 3) \cos (X-\bar{t}) + \frac{3}{64} (8\alpha^6 + (\alpha^2-1)^2) \cos (3X-3\bar{t}) [/math]

[math] {\bar{C}}^{(3)} = 0 \, [/math]

[math] D^{-1} = \alpha^{-1} \left\{ 1 + \varepsilon^2 \left[ \frac{9}{8} \left( \alpha^2 -1 \right)^2 + \alpha^2 \right] \right\} [/math]

where [math] \alpha = \coth k h \, [/math]

  • The non-dimensional solutions are then transferred back to the dimensional form.


[math] \phi^{(1)} = A \frac{\cosh [k(z+h)]}{\cosh (kh)} \sin\theta [/math]
[math] \eta^{(1)} = a \cos\theta \, [/math]


[math] \theta = kx - \sigma t + \beta \, [/math] and [math] a = A \frac{\sigma}{g} \, [/math]


[math] \phi^{(2)} = \frac{3akA}{8\alpha} (\alpha^2-1)^2 \cosh [ 2k (z+h)] \sin 2 \theta [/math]
[math] \eta^2 = \frac{1}{4} \alpha (3\alpha^2-1) a^2 k \cos 2 \theta [/math]

Bernoulli Constant: [math] C_o = \frac{1}{4\alpha} a^2 kg (\alpha^2-1) [/math]


[math] \phi^{(3)} = \frac{1}{64} (\alpha^2-1)(\alpha^2+3)(9\alpha^2-13) \cdot \frac{\cosh[3k(z+h)]}{\cosh3kh} a^2 k^2 A \sin 3 \theta [/math]
[math] \eta^{(3)} = - \frac{3}{8} ( \alpha^4 - 3 \alpha^2 + 3 ) a^3 k^2 \cos \theta + \frac{3}{64} ( 8 \alpha^6 + ( \alpha^2-1 )^2 ) a^3 k^2 \cos 3 \theta [/math]

Nonlinear Dispersion Relation:

[math] \sigma^2 = g k \tanh (kh) \left\{ 1 + k^2 a^2 \left[ \frac{9}{8} \left( \alpha^2 - 1 \right)^2 + \alpha^2 \right] \right\} [/math]


For the fast convergence of the perturbed coefficient, [math] \varepsilon \, [/math], must be much smaller than unity, which is consistent with weakly nonlinear assumption. However, when the ratio of depth to wave length is small, the Stokes perturbation may not be valid.

Convergence rate:

[math] R_\phi = \frac{\left| \phi^{(2)} \right|_{mag}}{\left| \phi^{(1)} \right|_{mag}}, \, [/math]

[math] R_\phi \, [/math] is the ratio of the potential magnitude of second-order to that of first order solution at [math] Z = 0 \, [/math].

[math] R_\phi = \frac{ 3 \varepsilon ( \alpha^2 - 1 )^2 }{ 8 a } \, [/math]

For fast convergence, [math] R_\phi \, [/math] should be [math] \ll 1 \, [/math]. This is true when [math] k h \sim O (1) \, [/math]. When [math] k h \ll 1 \, [/math], we have: [math] \alpha \sim ( kh )^{-1} \, [/math], hence [math] R_\phi \sim O \left( \frac{3}{8} \varepsilon ( kh )^{-3} \right) \, [/math]

[math] R_\phi \, [/math] may be much greater than unity [math] U_r \, [/math] sell number [math] U_r = \frac{a}{h} \frac{1}{(kh)^2} = \frac{\varepsilon}{(kh)^3} \, [/math]

For [math] R_\phi \ll 1 [/math], then [math] U_r \ll \frac{8}{3} \, [/math].

A few striking features of a nonlinear wave train can be described for the above equation:

  • The crests are steeper and troughs are flatter; (see applet (Nonlinear Wave Surface)).
  • Phase velocity increases with the increase in wave steepness.
  • Non-closed trajectories of particles movement. (see applet (N-Trajectory)).
  • Nonlinear wave characteristics (up to 2nd order).

Wave advancing in the x-direction

Particle velocity [math] \vec{V} = \vec{i} u + \vec{k} w \, [/math]

[math] u = \frac{akg}{\sigma} \frac{\cosh[k(z+h)]}{\cosh kh} \cos\theta + \frac{3}{4} \frac{a^2 k^2 g}{\sigma} \frac{\cosh[2k(z+h)]}{\sinh^3 kh \cosh kh} \cos 2 \theta [/math]

[math] w = \frac{akg}{\sigma} \frac{\sinh[k(z+h)]}{\cosh kh} \sin\theta + \frac{3}{4} \frac{a^2 k^2 g}{\sigma} \frac{\sinh[2k(z+h)]}{\sinh^3 kh \cosh kh} \sin 2 \theta [/math]

Acceleration [math] \vec{a} = \vec{i} a_x + \vec{k} a_z \, [/math]

[math] a_x = \frac{\partial u}{\partial t} + {\vec{V}}^{(1)} \cdot \nabla u^{(1)} = akg \frac{\cosh[k(z+h)]}{\cosh kh} \sin\theta + a^2 k^2 g \left[ \frac{3}{2} \frac{\cosh [2k(z+h)]}{\sinh^3kh \cosh kh} - \frac{1}{\sinh 2kh} \right] \sin 2 \theta [/math]

[math] a_z = \frac{\partial w}{\partial t} + {\vec{V}}^{(1)} \cdot \nabla w^{(1)} = -akg \frac{\sinh[k(z+h)]}{\cosh kh} \cos\theta + a^2 k^2 g \left[ \frac{3}{2} \frac{\sinh [2k(z+h)]}{\sinh^3kh \cosh kh} - \frac{1}{\sinh 2kh} \right] \cos 2 \theta [/math]

Particle Trajectory

Denoting the mean position of a particle by [math] (x,z) \, [/math], and its instantaneous displacement from the mean position by [math] (\zeta, \zeta) \, [/math], the Lagrangian velocities of the particle are hence [math] u ( x + \zeta, z + \zeta ) \, [/math] and [math] w ( x + \zeta, z + \zeta ) \, [/math], they are related to the Eulurian velocities through a Taylor Expansion:

[math] u (x+\zeta, z+\zeta) = u^{(1)}(x,z) + u^{(2)}(x,z) + \frac{\partial u^{(1)}}{\partial x} \zeta^{(1)} + \frac{\partial u^{(1)}}{\partial z} \zeta^{(1)} + O ( \varepsilon^2 ) u^{(1)} [/math]
[math] w (x+\zeta, z+\zeta) = w^{(1)}(x,z) + w^{(2)}(x,z) + \frac{\partial w^{(1)}}{\partial x} \zeta^{(1)} + \frac{\partial w^{(1)}}{\partial z} \zeta^{(1)} + O (\varepsilon^2) w^{(1)} [/math]

where [math] u^{(1)}(x,z),u^{(2)}(x,z),w^{(1)}(x,z) \, [/math] and [math] w^{(2)}(x,z) \, [/math] are first- and second-order horizontal and vertical velocities.

[math] ( \zeta, \xi) \, [/math] are calculated by integrating the related Lagrangian velocities. [math] \zeta (t) = \int_0^t u (x+\zeta,z+\xi,\tau)\mathrm{d}t + \zeta_0; \ \xi(t) = \int_0^t w(x+\zeta,z+\xi,\tau)\mathrm{d}\tau + \xi_0 [/math]

We intend to compute [math] (\zeta, \xi) \, [/math] up to second order in wave steepness

[math] \zeta(t) = \zeta^{(1)} (x,z,t) + \zeta^{(2)} (x,z,t) + {\overline{\zeta}}^{(2)} (x,z) + O (\varepsilon^2) \zeta^{(1)} [/math]

[math] \xi(t) = \xi^{(1)} (x,z,t) + \xi^{(2)} (x,z,t) + O ( \varepsilon^2) \xi^{(1)} [/math]

where superscripts stand for orders and overbar denotes a secular term. At lending-order, the solution is the same as that in LWT,

[math] \zeta^{(1)} = \int_0^t u^{(1)} (x+\zeta, z+\xi, \tau) \mathrm{d}\tau + \zeta_0^{(1)} = - a \frac{\cosh [ k(z+h)]}{\sinh kh} \sin\theta [/math]

[math] \xi^{(1)} = \int_0^t w^{(1)} (x+\zeta, z+\xi, \tau) \mathrm{d}\tau + \xi_0 = a \frac{\sinh [ k(z+h)]}{\sinh kh} \cos\theta [/math]

[math] \frac{\partial{u}^{(1)}}{\partial{x}} \zeta^{(1)} + \frac{\partial{u}^{(1)}}{\partial{z}} \xi^{(1)} = \frac{a^2k^2g}{\sigma} \left[ \frac{\cosh^2 [k(z+h)]}{\sinh kh \cosh kh} \sin^2 \theta + \frac{\sinh^2 [k(z+h)]}{\sinh kh \cosh kh} \cos^2\theta \right] [/math]
[math] = \frac{a^2k^2\sigma}{\sigma} \left\{ - \frac{1}{\sinh 2kh} \cos 2 \theta + \frac{ \cosh [ 2k(z+h)]}{\sinh 2kh} \right\} [/math]

[math] \frac{\partial w^{(1)}}{\partial{x}} \zeta^{(1)} + \frac{\partial{w}^{(1)}}{\partial{z}} \xi^{(1)} = 0 \, [/math]

[math] u^{(2)} (x,z) + \frac{\partial{u}^{(1)}}{\partial{x}} \zeta^{(1)} + \frac{\partial{u}^{(1)}}{\partial{z}} \xi^{(1)} = \frac{a^2k^2g}{\sigma} \left[ \frac{3}{4} \frac{\cosh [2k(z+h)]}{\sinh^3 kh \cosh kh} - \frac{1}{\sinh 2 kh} \right] \cos 2 \theta [/math]
[math] + \frac{a^2 k^2 g}{\sigma} \frac{\cosh [2 k(z+h)]}{\sinh 2 kh}, [/math]

[math] w^{(2)} (x,z) + \frac{ \partial{w}^{(1)}}{\partial x} \zeta^{(1)} + \frac{ \partial{w}^{(1)}}{\partial z} \xi^{(1)} = \frac{3}{4} \frac{a^2 k^2 g}{\sigma} \frac{ \sinh [ 2 k (z+h)]}{\sinh^3 kh \cosh kh} \sin 2 \theta [/math]

The leading-order trajectory of a particle is an ellipse of the center at [math] (x,z) \, [/math] and a major-axis [math] a \frac{\cosh [ k(z+h)]}{\sinh kh} [/math] and minor-axis [math] a \frac{\sinh[k(z+h)]}{\sinh kh} [/math].

[math] \frac{ \left(\zeta^{(1)} \right)^2}{\left[a\frac{\cosh[k(z+h)]}{\sinh kh} \right]^2} + \frac{\left(\xi^{(1)}\right)^2}{\left[a\frac{\sinh[k(z+h)]}{\sinh kh} \right]^2} = 1 [/math]

The second-order solutions for the displacement are calculated by integrating the related second-order lagrangian velocities.

[math] \zeta^{(2)} = a^2 k \left[ - \frac{3}{8} \frac{\cosh [ 2k (z+h)]}{\sinh^4 kh} + \frac{1}{4\sinh^2 kh} \right] \sin 2\theta [/math]

[math] {\bar{\zeta}}^{(2)} = a^2 k \sigma \frac{ \cosh [ 2k (z+h)]}{2\sinh^2 kh} t [/math]

[math] \xi^{(2)} = \frac{3}{8} a^2 k \frac{ \sinh [ 2k (z+h) ] }{ \sinh^4 kh } \cos 2 \theta [/math]

The secular term [math] \left( \bar{\zeta^{(2)}} \right) \, [/math] in the horizontal displacement indicates the particles will continuously move in the wave direction. Hence, the trajectory of a particle is no longer an ellipse. Because the horizontal mean position of a particle is not fixed at [math] x \, [/math] but change with time, we re-define the horizontal mean position by

[math] x' = x + a^2 k \sigma \frac{ \cosh[2k(z+h)]}{2\sinh^2kh} t [/math] and [math] \theta' = k x' - \sigma t \, [/math].

Correspondingly, the displacement with respect to the instantaneous mean position [math] (x', z) \, [/math] is given by,

<center>[math] \zeta^{(1)} = - a \frac{\cosh[k(z+h)]}{\sinh kh} \sin\theta', \quad \xi^{(1)} = a \frac{\sinh [k (z+h)]}{\sinh kh} \cos \theta', [/math]

[math] \zeta^{(2)} = a^2 k \left[ - \frac{3}{8} \frac{\cosh [ 2k (z+h)]}{\sinh^4 kh} + \frac{1}{4\sinh^2 kh} \right] \sin 2 \theta', [/math]

[math] \xi^{(2)} = \frac{3}{8} a^2 k \frac{\sinh[2k(z+h)]}{\sinh^4kh} \cos 2 \theta' [/math]

The trajectory of a particle based on the solution is plotted in Applet (N-trajectory). The time average Lagrangian velocity of a particle is equal to the derivative of the secular term [math] \left( \bar{\zeta}^{(2)}\right) \, [/math] with respect to time.

[math] \bar{u}_l = a^2 k \sigma \frac{ \cosh [2k(z+h)]}{2\sinh^2 kh} [/math]

The integral of the average Lagrangian velocity with respect to water depth renders the average mass flux induced by a periodic wave train over a unit width

Mass flux [math] = \int_{-h}^0 \rho \bar{u}_l \mathrm{d}z = \frac{1}{2} \rho a^2 \sigma \alpha = \frac{1}{2} \rho a^2 kg/\sigma = E/C_p, [/math]

which is consistent with the result derived using Eulurian approach.

Dynamic Pressure

Using the Bernoulli equation, dynamic pressure head induced by a periodic wave train can be calculated up to second-order,

[math] \frac{p}{\rho g} = - \frac{1}{g} \left[ \frac{\partial\phi^{(1)}}{\partial t} + \frac{ \partial\phi^{(1)} }{\partial t} \right] - \frac{1}{2g} \left[ \left| \nabla \phi^{(1)} \right|^2 \right] + C_0, \, [/math]

[math] \frac{p^{(1)}}{\rho g} = a \frac{\cosh [ k (z+h)]}{\cosh (kh)} \cos \theta, \, [/math]

[math] \frac{p^{(2)}}{\rho g} = \frac{3a^2 k}{4\alpha} \left( \alpha^2 -1 \right)^2 \cosh [2k(z+h)] \cos 2 \theta - \frac{a^2k}{4} \alpha ( \alpha^2 -1 ) \cos 2 \theta, \, [/math]

[math] \frac{ {\bar{p}}^{(2)}}{\rho g} = \frac{\alpha}{4} a^2 k ( \alpha^2 - 1 ) \left[ 1 - \cosh [ 2 k (z+h)] \right]. \, [/math]

Radiation Stress

Radiation stress: defined as the time average of excess quasi momentum flux due to the presence of a periodic wave train.

Up to second order, a wave train advancing in the x-axis, [math] S = \begin{bmatrix} S_{xx} & S_{xy} \\ S_{yx} & S_{yy} \end{bmatrix} [/math]

[math] S_{xx} = \overline{\int_{-h}^{\eta} \left( p + \rho u^2 \right) \mathrm{d}z} - \int_{-h}^0 p_0 \mathrm{d}z = \int_{-h}^0 \rho \overline{u^2} \mathrm{d}z + \int_{-h}^0 \overline{ \left( p - p_0 \right) } \mathrm{d}z + \overline{ \int_0^\eta p \mathrm{d}z } [/math]

Noticing that [math] \overline{ \left( p + \rho w^2 \right) } = - \rho g z = p+0, \, [/math]

[math] S_{xx} = \int_{-h}^0 \rho \overline{ \left( u^{(1)2} - w^{(1)2} \right) } \mathrm{d}z + \overline{ \int_0^\eta p^{(1)} \mathrm{d}z } = \frac{ \rho g a^2 k h }{ \sinh 2 k h } + \frac{1}{4} \rho g a^2 [/math]

[math] S_{yy} = \overline{ \int_{-h}^\eta \left( p + \rho v^2 \right) \mathrm{d}z } - \int_{-h}^0 p_0 \mathrm{d}z = \int_{-h}^0 \overline{ (p - p_0) } \mathrm{d}z + \overline{ \int_0^\eta p \mathrm{d}z } = \frac{ \rho g a^2 kh}{ 2 \sinh 2 k h } [/math]
[math] S_{xy} = S_{yx} = 0 \, [/math]

[math] S = E \begin{bmatrix} \frac{2kh}{\sinh 2kh} + \frac{1}{2} & 0 \\ 0 \frac{kh}{\sinh 2kh} \end{bmatrix} [/math], where [math] E = \frac{1}{2} \rho g a^2 \, [/math], the energy density.

In deep water In shallow water
[math] S= E \begin{bmatrix} 1/2 & 0 \\ 0 & 0 \end{bmatrix} [/math] [math] S= E \begin{bmatrix} 3/2 & 0 \\ 0 & 1/2 \end{bmatrix} [/math]

In the case of a wave train having an angle, [math] \theta \, [/math] with respect to the x-axis,

[math] S = E \begin{bmatrix} \left( \frac{kh}{\sinh 2kh} + \frac{1}{2} \right) \frac{(3+\cos 2\theta)}{2} - \frac{1}{2} & \left( 1 + \frac{2kh}{\sinh 2kh} \right) \frac{\sin 2\theta}{4} \\ \left(1 + \frac{2kh}{\sinh 2kh} \right) \frac{\sin 2 \theta}{4} & \left( \frac{kh}{\sinh 2kh} + \frac{1}{2} \right) \frac{(3-\cos 2 \theta)}{2} - \frac{1}{2} \end{bmatrix} [/math]

This article is based on the lecture notes of Dr. Jun Zhang and the original article can be found here