Mildslope equation
In fluid dynamics, the mildslope equation describes the combined effects of diffraction and refraction for water waves propagating over bathymetry and due to lateral boundaries—like breakwaters and coastlines. It is an approximate model, deriving its name from being originally developed for wave propagation over mild slopes of the sea floor. The mildslope equation is often used in coastal engineering to compute the wavefield changes near harbours and coasts.
The mildslope equation models the propagation and transformation of water waves, as they travel through waters of varying depth and interact with lateral boundaries such as cliffs, beaches, seawalls and breakwaters. As a result, it describes the variations in wave amplitude, or equivalently wave height. From the wave amplitude, the amplitude of the flow velocity oscillations underneath the water surface can also be computed. These quantities—wave amplitude and flowvelocity amplitude—may subsequently be used to determine the wave effects on coastal and offshore structures, ships and other floating objects, sediment transport and resulting geomorphology changes of the sea bed and coastline, mean flow fields and mass transfer of dissolved and floating materials. Most often, the mildslope equation is solved by computer using methods from numerical analysis.
A first form of the mildslope equation was developed by Eckart in 1952, and an improved version—the mildslope equation in its classical formulation—has been derived independently by Juri Berkhoff in 1972.^{[1]}^{[2]}^{[3]} Thereafter, many modified and extended forms have been proposed, to include the effects of, for instance: wave–current interaction, wave nonlinearity, steeper seabed slopes, bed friction and wave breaking. Also parabolic approximations to the mildslope equation are often used, in order to reduce the computational cost.
In case of a constant depth, the mildslope equation reduces to the Helmholtz equation for wave diffraction.
Contents
Formulation for monochromatic wave motion
For monochromatic waves according to linear theory—with the free surface elevation given as \zeta(x,y,t)=\Re\left\{\eta(x,y)\,\text{e}^{i\omega t}\right\} and the waves propagating on a fluid layer of mean water depth h(x,y)—the mildslope equation is:^{[4]}
 \nabla\cdot\left( c_p\, c_g\, \nabla \eta \right)\, +\, k^2\, c_p\, c_g\, \eta\, =\, 0,
where:
 \eta(x,y) is the complexvalued amplitude of the freesurface elevation \zeta(x,y,t);
 (x,y) is the horizontal position;
 \omega is the angular frequency of the monochromatic wave motion;
 i is the imaginary unit;
 \Re\{\cdot\} means taking the real part of the quantity between braces;
 \nabla is the horizontal gradient operator;
 \nabla\cdot is the divergence operator;
 k is the wavenumber;
 c_p is the phase speed of the waves and
 c_g is the group speed of the waves.
The phase and group speed depend on the dispersion relation, and are derived from Airy wave theory as:^{[5]}
 \begin{align} \omega^2 &=\, g\, k\, \tanh\, (kh), \\ c_p &=\, \frac{\omega}{k} \quad \text{and} \\ c_g &=\, \frac12\, c_p\, \left[ 1\, +\, kh\, \frac{1  \tanh^2 (kh)}{\tanh\, (kh)} \right] \end{align}
where
 g is Earth's gravity and
 \tanh is the hyperbolic tangent.
For a given angular frequency \omega, the wavenumber k has to be solved from the dispersion equation, which relates these two quantities to the water depth h.
Transformation to an inhomogeneous Helmholtz equation
Through the transformation
 \psi\, =\, \eta\, \sqrt{c_p\, c_g},
the mild slope equation can be cast in the form of an inhomogeneous Helmholtz equation:^{[4]}^{[6]}
 \Delta\psi\, +\, k_c^2\, \psi\, =\, 0 \qquad \text{with} \qquad k_c^2\, =\, k^2\, \, \frac{\Delta\left(\sqrt{c_p\,c_g}\right)}{\sqrt{c_p\,c_g}},
where \Delta is the Laplace operator.
Propagating waves
In spatially coherent fields of propagating waves, it is useful to split the complex amplitude \eta(x,y) in its amplitude and phase, both real valued:^{[7]}
 \eta(x,y)\, =\, a(x,y)\, \text{e}^{i\, \theta(x,y)},
where
 a=\eta\, is the amplitude or absolute value of \eta\, and
 \theta=\arg\{\eta\}\, is the wave phase, which is the argument of \eta.\,
This transforms the mildslope equation in the following set of equations (apart from locations for which \nabla\theta is singular):^{[7]}
 \begin{align} \frac{\partial\kappa_y}{\partial{x}}\, \, \frac{\partial\kappa_x}{\partial{y}}\, =\, 0 \qquad &\text{ with } \kappa_x\, =\, \frac{\partial\theta}{\partial{x}} \text{ and } \kappa_y\, =\, \frac{\partial\theta}{\partial{y}}, \\ \kappa^2\, =\, k^2\, +\, \frac{\nabla\cdot\left( c_p\, c_g\, \nabla a \right)}{c_p\, c_g\, a} \qquad &\text{ with } \kappa\, =\, \sqrt{\kappa_x^2 \, +\, \kappa_y^2} \quad \text{ and} \\ \nabla \cdot \left( \boldsymbol{v}_g\, E \right)\, =\, 0 \qquad &\text{ with } E\, =\, \frac12\, \rho\, g\, a^2 \quad \text{and} \quad \boldsymbol{v}_g\, =\, c_g\, \frac{\boldsymbol{\kappa}}{k}, \end{align}
where
 E is the average waveenergy density per unit horizontal area (the sum of the kinetic and potential energy densities),
 \boldsymbol{\kappa} is the effective wavenumber vector, with components (\kappa_x,\kappa_y),\,
 \boldsymbol{v}_g is the effective group velocity vector,
 \rho is the fluid density, and
 g is the acceleration by the Earth's gravity.
The last equation shows that wave energy is conserved in the mildslope equation, and that the wave energy E is transported in the \boldsymbol{\kappa}direction normal to the wave crests (in this case of pure wave motion without mean currents).^{[7]} The effective group speed \boldsymbol{v}_g is different from the group speed c_g.
The first equation states that the effective wavenumber \boldsymbol{\kappa} is irrotational, a direct consequence of the fact it is the derivative of the wave phase \theta, a scalar field. The second equation is the eikonal equation. It shows the effects of diffraction on the effective wavenumber: only for moreorless progressive waves, with \nabla\cdot(c_p\, c_g\, \nabla a)\ll k^2\, c_p\, c_g\, a, the splitting into amplitude a and phase \theta leads to consistentvarying and meaningful fields of a and \boldsymbol{\kappa}. Otherwise, κ^{2} can even become negative. When the diffraction effects are totally neglected, the effective wavenumber κ is equal to k, and the geometric optics approximation for wave refraction can be used.^{[7]}
When \eta=a\mbox{ }\exp(i\theta) is used in the mildslope equation, the result is, apart from a factor \exp(i\theta):
 c_p\,c_g\, \left( \Delta a\, +\, 2i\, \nabla a \cdot \nabla\theta\, \, a\, \nabla\theta \cdot \nabla\theta\, +\, i\, a\, \Delta\theta \right)\, +\, \nabla \left( c_p\, c_g \right) \cdot \left( \nabla a\, +\, i\, a\, \nabla\theta \right)\, +\, k^2\, c_p\, c_g\, a\, =\, 0.
Now both the real part and the imaginary part of this equation have to be equal to zero:
 \begin{align} & c_p\,c_g\, \Delta a\, \, c_p\, c_g\, a\, \nabla\theta \cdot \nabla\theta\, +\, \nabla \left( c_p\, c_g \right) \cdot \nabla a\, +\, k^2\, c_p\, c_g\, a\, =\, 0 \quad \text{and} \\ & 2\, c_p\,c_g\, \nabla a \cdot \nabla\theta\, +\, c_p\, c_g\, a\, \Delta\theta\, +\, \nabla \left( c_p\, c_g \right) \cdot \left( a\, \nabla\theta \right)\, =\, 0. \end{align}
The effective wavenumber vector \boldsymbol{\kappa} is defined as the gradient of the wave phase:
 \boldsymbol{\kappa}\, =\, \nabla\theta and its vector length is \kappa=\boldsymbol{\kappa}.
Note that \boldsymbol{\kappa} is an irrotational field, since the curl of the gradient is zero:
 \nabla\, \times\, \boldsymbol{\kappa}\, =\, 0.
Now the real and imaginary parts of the transformed mildslope equation become, first multiplying the imaginary part by a:
 \begin{align} &\kappa^2\, =\, k^2\, +\, \frac{\nabla(c_p\, c_g)}{c_p\, c_g} \cdot \frac{\nabla a}{a}\, +\, \frac{\Delta a}{a} \quad \text{and} \\ &c_p\, c_g\, \nabla\left(a^2\right) \cdot \boldsymbol{\kappa}\, +\, c_p\, c_g\, \nabla\cdot\boldsymbol{\kappa}\, +\, a^2\, \boldsymbol{\kappa} \cdot \nabla \left( c_p\, c_g \right)\, =\, 0. \end{align}
The first equation directly leads to the eikonal equation above for \kappa\,, while the second gives:
 \nabla \cdot \left( \boldsymbol{\kappa}\, c_p\, c_g\, a^2 \right)\, =\, 0,
which—by noting that c_p=\omega/k in which the angular frequency \omega is a constant for timeharmonic motion—leads to the waveenergy conservation equation.
Derivation of the mildslope equation
The mildslope equation can be derived by the use of several methods. Here, we will use a variational approach.^{[4]}^{[8]} The fluid is assumed to be inviscid and incompressible, and the flow is assumed to be irrotational. These assumptions are valid ones for surface gravity waves, since the effects of vorticity and viscosity are only significant in the Stokes boundary layers (for the oscillatory part of the flow). Because the flow is irrotational, the wave motion can be described using potential flow theory.
Luke's variational principle
Luke's Lagrangian formulation gives a variational formulation for nonlinear surface gravity waves.^{[9]} For the case of a horizontally unbounded domain with a constant density \rho, a free fluid surface at z=\zeta(x,y,t) and a fixed sea bed at z=h(x,y), Luke's variational principle \delta\mathcal{L}=0 uses the Lagrangian
 \mathcal{L} = \int_{t_0}^{t_1} \iint L\; \text{d}x\; \text{d}y\; \text{d}t,
where L is the horizontal Lagrangian density, given by:
 L = \rho\, \left\{ \int_{h(x,y)}^{\zeta(x,y,t)} \left[ \frac{\partial\Phi}{\partial t} +\, \frac{1}{2} \left( \left( \frac{\partial\Phi}{\partial x} \right)^2 + \left( \frac{\partial\Phi}{\partial y} \right)^2 + \left( \frac{\partial\Phi}{\partial z} \right)^2 \right) \right]\; \text{d}z\; +\, \frac{1}{2}\, g\, (\zeta^2\, \, h^2) \right\},
where \Phi(x,y,z,t) is the velocity potential, with the flow velocity components being \partial\Phi/\partial{x}, \partial\Phi/\partial{y} and \partial\Phi/\partial{z} in the x, y and z directions, respectively. Luke's Lagrangian formulation can also be recast into a Hamiltonian formulation in terms of the surface elevation and velocity potential at the free surface.^{[10]} Taking the variations of \mathcal{L}(\Phi,\zeta) with respect to the potential \Phi(x,y,z,t) and surface elevation \zeta(x,y,t) leads to the Laplace equation for \Phi in the fluid interior, as well as all the boundary conditions both on the free surface z=\zeta(x,y,t) as at the bed at z=h(x,y).
Linear wave theory
In case of linear wave theory, the vertical integral in the Lagrangian density L is split into a part from the bed z=h to the mean surface at z=0, and a second part from z=0 to the free surface z=\zeta. Using a Taylor series expansion for the second integral around the mean freesurface elevation z=0, and only retaining quadratic terms in \Phi and \zeta, the Lagrangian density L_0 for linear wave motion becomes
 L_0 = \rho\, \left\{ \zeta\, \left[ \frac{\partial\Phi}{\partial t} \right]_{z=0}\, +\, \int_{h}^0 \frac12 \left[ \left( \frac{\partial\Phi}{\partial x} \right)^2 + \left( \frac{\partial\Phi}{\partial y} \right)^2 + \left( \frac{\partial\Phi}{\partial z} \right)^2 \right]\; \text{d}z\; +\, \frac{1}{2}\, g\, \zeta^2\, \right\}.
The term \partial\Phi/\partial{t} in the vertical integral is dropped since it has become dynamically uninteresting: it gives a zero contribution to the Euler–Lagrange equations, with the upper integration limit now fixed. The same is true for the neglected bottom term proportional to h^2 in the potential energy.
The waves propagate in the horizontal (x,y) plane, while the structure of the potential \Phi is not wavelike in the vertical zdirection. This suggests the use of the following assumption on the form of the potential \Phi:
 \Phi(x,y,z,t) = f(z;x,y)\, \varphi(x,y,t) with normalisation f(0;x,y)=1 at the mean freesurface elevation z=0.
Here \varphi(x,y,t) is the velocity potential at the mean freesurface level z=0. Next, the mildslope assumption is made, in that the vertical shape function f changes slowly in the (x,y)plane, and horizontal derivatives of f can be neglected in the flow velocity. So:
 \begin{pmatrix} \displaystyle \frac{\partial\Phi}{\partial{x}} \\[2ex] \displaystyle \frac{\partial\Phi}{\partial{y}} \\[2ex] \displaystyle \frac{\partial\Phi}{\partial{z}} \end{pmatrix}\, \approx\, \begin{pmatrix} \displaystyle f\, \frac{\partial\varphi}{\partial{x}} \\[2ex] \displaystyle f\, \frac{\partial\varphi}{\partial{y}} \\[2ex] \displaystyle \frac{\partial{f}}{\partial{z}}\, \varphi \end{pmatrix}.
As a result:
 L_0 = \rho\, \left\{ \zeta\, \frac{\partial\varphi}{\partial t}\, +\, \frac12\, F\, \left[ \left( \frac{\partial\varphi}{\partial{x}} \right)^2\, +\, \left( \frac{\partial\varphi}{\partial{y}} \right)^2 \right]\, +\, \frac12\, G\, \varphi^2\, +\, \frac12\, g\, \zeta^2\, \right\}, with F\, =\, \int_{h}^0 f^2\; \text{d}z and G\, =\, \int_{h}^0 \left(\frac{\text{d}f}{\text{d}z}\right)^2\; \text{d}z.
The Euler–Lagrange equations for this Lagrangian density L_0 are, with \xi(x,y,t) representing either \varphi or \zeta:
 \frac{\partial{L_0}}{\partial\xi}  \frac{\partial}{\partial{t}}\left( \frac{\partial{L_0}}{\partial(\partial\xi/\partial{t})} \right)  \frac{\partial}{\partial{x}}\left( \frac{\partial{L_0}}{\partial(\partial\xi/\partial{x})} \right)  \frac{\partial}{\partial{y}}\left( \frac{\partial{L_0}}{\partial(\partial\xi/\partial{y})} \right) = 0.
Now \xi is first taken equal to \varphi and then to \zeta. As a result, the evolution equations for the wave motion become:^{[4]}
 \begin{align} \frac{\partial\zeta}{\partial t}\, &+ \nabla \cdot \left( F\, \nabla\varphi \right)\, \, G\, \varphi\, =\, 0 \quad \text{and} \\ \frac{\partial\varphi}{\partial t}\, &+\, g\, \zeta\, =\, 0, \end{align}
with ∇ the horizontal gradient operator: ∇ ≡(∂/∂x ∂/∂y)^{T} where T denotes the transpose.
The next step is to choose the shape function f and to determine F and G.
Vertical shape function from Airy wave theory
Since the objective is the description of waves over mildly sloping beds, the shape function f(z) is chosen according to Airy wave theory. This is the linear theory of waves propagating in constant depth h. The form of the shape function is:^{[4]}
 f = \frac{\cosh\, \bigl( k\, (z+h) \bigr)}{\cosh\, (k h)},
with k(x,y) now in general not a constant, but chosen to vary with x and y according to the local depth h(x,y) and the linear dispersion relation:^{[4]}
 \omega_0^2\, =\, g\, k\, \tanh\, (k h).
Here \omega_0 a constant angular frequency, chosen in accordance with the characteristics of the wave field under study. Consequently, the integrals F and G become:^{[4]}
 \begin{align} F &= \int_h^0 f^2\; \text{d}z = \frac{1}{g}\, c_p\, c_g \quad \text{and} \\ G &= \int_h^0 \left( \frac{\partial{f}}{\partial{z}} \right)^2\; \text{d}z = \frac{1}{g} \left( \omega_0^2\, \, k^2\, c_p\, c_g \right). \end{align}
The following timedependent equations give the evolution of the freesurface elevation \zeta(x,y,t) and freesurface potential \phi(x,y,t):^{[4]}
 \begin{align} g\, \frac{\partial\zeta}{\partial{t}} &+ \nabla\cdot\left( c_p\, c_g\, \nabla \varphi \right) + \left( k^2\, c_p\, c_g\, \, \omega_0^2 \right)\, \varphi = 0, \\ \frac{\partial\varphi}{\partial{t}} &+ g \zeta = 0, \quad \text{with} \quad \omega_0^2\, =\, g\, k\, \tanh\, (kh). \end{align}
From the two evolution equations, one of the variables \varphi or \zeta can be eliminated, to obtain the timedependent form of the mildslope equation:^{[4]}
 \frac{\partial^2\zeta}{\partial{t^2}} + \nabla\cdot\left( c_p\, c_g\, \nabla \zeta \right) + \left( k^2\, c_p\, c_g\, \, \omega_0^2 \right)\, \zeta = 0,
and the corresponding equation for the freesurface potential is identical, with \zeta replaced by \varphi. The timedependent mildslope equation can be used to model waves in a narrow band of frequencies around \omega_0.
Monochromatic waves
Consider monochromatic waves with complex amplitude \eta(x,y) and angular frequency \omega:
 \zeta(x,y,t)\, =\, \Re\left\{ \eta(x,y)\; \text{e}^{i\, \omega\, t} \right\},
with \omega and \omega_0 chosen equal to each other, \omega=\omega_0. Using this in the timedependent form of the mildslope equation, recovers the classical mildslope equation for timeharmonic wave motion:^{[4]}
 \nabla \cdot \left( c_p\, c_g\, \nabla \eta \right)\, +\, k^2\, c_p\, c_g\, \eta\, =\, 0.
Applicability and validity of the mildslope equation
The standard mild slope equation, without extra terms for bed slope and bed curvature, provides accurate results for the wave field over bed slopes ranging from 0 to about 1/3.^{[11]} However, some subtle aspects, like the amplitude of reflected waves, can be completely wrong, even for slopes going to zero. This mathematical curiosity has little practical importance in general since this reflection becomes vanishingly small for small bottom slopes.
Notes
 ^
 ^ Berkhoff, J. C. W. (1972), "Computation of combined refraction–diffraction", Proceedings 13th International Conference on Coastal Engineering, Vancouver, pp. 471–490
 ^ Berkhoff, J. C. W. (1976), Mathematical models for simple harmonic linear water wave models; wave refraction and diffraction (PhD. Thesis), Delft University of Technology
 ^ ^{a} ^{b} ^{c} ^{d} ^{e} ^{f} ^{g} ^{h} ^{i} ^{j} See Dingemans (1997), pp. 248–256 & 378–379.
 ^ See Dingemans (1997), p. 49.
 ^ See Mei (1994), pp. 86–89.
 ^ ^{a} ^{b} ^{c} ^{d} See Dingemans (1997), pp. 259–262.
 ^ Booij, N. (1981), Gravity waves on water with nonuniform depth and current (PhD. Thesis), Delft University of Technology
 ^ Luke, J. C. (1967), "A variational principle for a fluid with a free surface", Journal of Fluid Mechanics 27 (2): 395–397,
 ^
 ^ Booij, N. (1983), "A note on the accuracy of the mildslope equation", Coastal Engineering 7 (1): 191–203,
References
 Dingemans, M. W. (1997), Water wave propagation over uneven bottoms, Advanced Series on Ocean Engineering 13, World Scientific, Singapore, , 2 Parts, 967 pages.
 Liu, P. L.F. (1990), "Wave transformation", in B. Le Méhauté and D. M. Hanes, Ocean Engineering Science, The Sea 9A, Wiley Interscience, pp. 27–63,
 , 740 pages.
 Porter, D.; Chamberlain, P. G. (1997), "Linear wave scattering by twodimensional topography", in J. N. Hunt, Gravity waves in water of finite depth, Advances in Fluid Mechanics 10, Computational Mechanics Publications, pp. 13–53,
 Porter, D. (2003), "The mildslope equations", Journal of Fluid Mechanics 494: 51–63,
