Spherical waves

Spherically symmetric three-dimensional waves propagate in the radial direction \(r\) only so that \(u = u(r,t)\). The fully three-dimensional wave equation

\[\frac{\partial^2u}{\partial t^2}=\nabla\cdot (c^2\nabla u) + f\]

then reduces to the spherically symmetric wave equation

\[\begin{split}\frac{\partial^2u}{\partial t^2}=\frac{1}{r^2}\frac{\partial}{\partial r} \left(c^2(r)r^2\frac{\partial u}{\partial t}\right) + f(r,t),\quad r\in (0,R),\ t>0 {\thinspace .}\end{split}\]

One can easily show that the function \(v(r,t) = ru(r,t)\) fulfills a standard wave equation in Cartesian coordinates if \(c\) is constant. To this end, insert \(u=v/r\) in

\[\frac{1}{r^2}\frac{\partial}{\partial r} \left(c^2(r)r^2\frac{\partial u}{\partial t}\right)\]

to obtain

\[r\left(\frac{d c^2}{dr}\frac{\partial v}{\partial r} + c^2\frac{\partial^2 v}{\partial r^2}\right) - \frac{d c^2}{dr}v {\thinspace .}\]

The two terms in the parenthesis can be combined to

\[r\frac{\partial}{\partial r}\left( c^2\frac{\partial v}{\partial r}\right),\]

which is recognized as the variable-coefficient Laplace operator in one Cartesian coordinate. The spherically symmetric wave equation in terms of \(v(r,t)\) now becomes

\[\begin{split}\frac{\partial^2v}{\partial t^2}= \frac{\partial}{\partial r} \left(c^2(r)\frac{\partial v}{\partial r}\right) -\frac{1}{r}\frac{d c^2}{dr}v + rf(r,t),\quad r\in (0,R),\ t>0 {\thinspace .}\end{split}\]

In the case of constant wave velocity \(c\), this equation reduces to the wave equation in a single Cartesian coordinate called \(r\):

(1)\[\begin{split} \frac{\partial^2v}{\partial t^2}= c^2 \frac{\partial^2 v}{\partial r^2} + rf(r,t),\quad r\in (0,R),\ t>0 {\thinspace .}\end{split}\]

That is, any program for solving the one-dimensional wave equation in a Cartesian coordinate system can be used to solve (1), provided the source term is multiplied by the coordinate, and that we divide the Cartesian mesh solution by \(r\) to get the spherically symmetric solution. Moreover, if \(r=0\) is included in the domain, spherical symmetry demands that \(\partial u/\partial r=0\) at \(r=0\), which means that

\[\frac{\partial u}{\partial r} = \frac{1}{r^2}\left( r\frac{\partial v}{\partial r} - v\right) = 0,\quad r=0,\]

implying \(v(0,t)=0\) as a necessary condition. For practical applications, we exclude \(r=0\) from the domain and assume that some boundary condition is assigned at \(r=\epsilon\), for some \(\epsilon >0\).

The linear shallow water equations

The next example considers water waves whose wavelengths are much lager than the depth and whose wave amplitudes are small. This class of waves may be generated by catastrophic geophysical events, such as earthquakes at the sea bottom, landslides moving into water, or underwater slides (or a combination, as earthquakes frequently release avalanches of masses). For example, a subsea earthquake will normally have an extension of many kilometers but lift the water only a few meters. The wave length will have a size dictated by the earthquake area, which is much lager than the water depth, and compared to this wave length, an amplitude of a few meters is very small. The water is essentially a thin film, and mathematically we can average the problem in the vertical direction and approximate the 3D wave phenomenon by 2D PDEs. Instead of a moving water domain in three space dimensions, we get a horizontal 2D domain with an unknown function for the surface elevation and the water depth as a variable coefficient in the PDEs.

Let \(\eta(x,y,t)\) be the elevation of the water surface, \(H(x,y)\) the water depth corresponding to a flat surface (\(\eta =0\)), \(u(x,y,t)\) and \(v(x,y,t)\) the depth-averaged horizontal velocities of the water. Mass and momentum balance of the water volume give rise to the PDEs involving these quantities:

(2)\[ \eta_t = - (Hu)_x - (Hv)_x\]

(3)\[ u_t = -g\eta_x,\]

(4)\[ v_t = -g\eta_y,\]

where \(g\) is the acceleration of gravity. Equation (2) corresponds to mass balance while the other two are derived from momentum balance (Newton’s second law).

The initial conditions associated with (2)-(4) are \(\eta\), \(u\), and \(v\) prescribed at \(t=0\). A common condition is to have some water elevation \(\eta =I(x,y)\) and assume that the surface is at rest: \(u=v=0\). A subsea earthquake usually means a sufficiently rapid motion of the bottom and the water volume to say that the bottom deformation is mirrored at the water surface as an initial lift \(I(x,y)\) and that \(u=v=0\).

Boundary conditions may be \(\eta\) prescribed for incoming, known waves, or zero normal velocity at reflecting boundaries (steep mountains, for instance): \(un_x + vn_y =0\), where \((n_x,n_y)\) is the outward unit normal to the boundary. More sophisticated boundary conditions are needed when waves run up at the shore, and at open boundaries where we want the waves to leave the computational domain undisturbed.

Equations (2), (3), and (4) can be transformed to a standard, linear wave equation. First, multiply (3) and (4) by \(H\), differentiate (3)) with respect to \(x\) and (4) with respect to \(y\). Second, differentiate (2) with respect to \(t\) and use that \((Hu)_{xt}=(Hu_t)_x\) and \((Hv)_{yt}=(Hv_t)_y\) when \(H\) is independent of \(t\). Third, eliminate \((Hu_t)_x\) and \((Hv_t)_y\) with the aid of the other two differentiated equations. These manipulations results in a standard, linear wave equation for \(\eta\):

(5)\[ \eta_{tt} = (gH\eta_x)_x + (gH\eta_y)_y = \nabla\cdot (gH\nabla\eta)\]\[ {\thinspace .}\]

In the case we have an initial non-flat water surface at rest, the initial conditions become \(\eta =I(x,y)\) and \(\eta_t=0\). The latter follows from (2) if \(u=v=0\), or simply from the fact that the vertical velocity of the surface is \(\eta_t\), which is zero for a surface at rest.

The system (2)-(4) can be extended to handle a time-varying bottom topography, which is relevant for modeling long waves generated by underwater slides. In such cases the water depth function \(H\) is also a function of \(t\), due to the moving slide, and one must add a time-derivative term \(H_t\) to the left-hand side of (2). A moving bottom is best described by introducing \(z=H_0\) as the still-water level, \(z=B(x,y,t)\) as the time- and space-varying bottom topography, so that \(H=H_0-B(x,y,t)\). In the elimination of \(u\) and \(v\) one may assume that the dependence of \(H\) on \(t\) can be neglected in the terms \((Hu)_{xt}\) and \((Hv)_{yt}\). We then end up with a source term in (5), because of the moving (accelerating) bottom:

(6)\[ \eta_{tt} = \nabla\cdot(gH\nabla\eta) + B_{tt}\]\[ {\thinspace .}\]

The reduction of (6) to 1D, for long waves in a straight channel, or for approximately plane waves in the ocean, is trivial by assuming no change in \(y\) direction (\(\partial/\partial y=0\)):

(7)\[ \eta_t = (gH\eta_x)_x + B_{tt}\]\[ {\thinspace .}\]

Wind drag on the surface

Surface waves are influenced by the drag of the wind, and if the wind velocity some meters above the surface is \((U,V)\), the wind drag gives contributions \(C_V\sqrt{U^2+V^2}U\) and \(C_V\sqrt{U^2+V^2}V\) to (3) and (4), respectively, on the right-hand sides.

Bottom drag

The waves will experience a drag from the bottom, often roughly modeled by a term similar to the wind drag: \(C_B\sqrt{u^2+v^2}u\) on the right-hand side of (3) and \(C_B\sqrt{u^2+v^2}v\) on the right-hand side of (4). Note that in this case the PDEs (3) and (4) become nonlinear and the elimination of \(u\) and \(v\) to arrive at a 2nd-order wave equation for \(\eta\) is not possible anymore.

Effect of the Earth’s rotation

Long geophysical waves will often be affected by the rotation of the Earth because of the Coriolis force. This force gives rise to a term \(fv\) on the right-hand side of (3) and \(-fu\) on the right-hand side of (4). Also in this case one cannot eliminate \(u\) and \(v\) to work with a single equation for \(\eta\). The Coriolis parameter is \(f=2\Omega\sin\phi\), where \(\Omega\) is the angular velocity of the earth and \(\phi\) is the latitude.