Ripples from a Splashing Drop

Ripples from a Splashing Drop

When Jordan Spieth failed to make the green on the 12th at Augusta in 2016 and saw his ball plunge into the lake, his dismay was perhaps tempered by the beautifully symmetric ripples or waves emanating from the point where the ball entered the water. An example of such a splash can be seen in Figure 1.

Figure 1: A splashing drop displaying the concentric wave motion

This is a classical free surface fluid mechanics problem involving the Navier–Stokes equations together with the conservation of mass. These may be solved by a variety of methods. We choose a finite difference method employing marker and cell (MAC method) on a staggered grid, which we have been developing over many years. A review of the MAC method may be found in [1]. Figure 2 displays the numerical solution of a splashing drop at three different times using the MAC method (albeit, in a confined region with wave reflection). Note that in Figure 1, a series of concentric crests and troughs are observed whose amplitudes are not constant, nor are the distances between the crests.

Figure 2: Numerical computation at three different times

It is, however, well known but nonetheless somewhat surprising that substantial progress may be made analytically using asymptotic arguments. In the remainder of this note we shall show how using Hankel transforms, the method of stationary phase and asymptotic analysis, can lead to a simple analytic expression which encapsulates the main features of the flow and indeed compares tolerably with the computation of Figure 2.

An asymptotic expression

In cylindrical polar coordinates with z-axis vertically upwards, the linearised water wave equations are
Ripples-from-a-Splashing-Drop-maths 1where \zeta denotes the elevation (above z = 0) of the water and \phi is the velocity potential (i.e. \bm{u} = -\nabla \phi). See, e.g. Coulson [2].

Note that from Jordan Spieth’s point of view, time is initiated, not when his ball hits the water, but at the time when the resultant splash occurs.

Define the Hankel transforms

    \begin{align*} \tilde\phi (k, z, t) &= \int^\infty_0 r\, J_0 (kr) \phi (r, z, t) \mathrm{d} r\\ \tilde\zeta (k, t) &= \int^\infty_0 r\, J_0 (kr) \zeta (r, t) \mathrm{d} r. \end{align*}

Taking Hankel transforms of (1) yields

    \[ \frac{\partial^2 \tilde \phi}{\partial z^2} - k^2 \tilde \phi = 0 \]

or, using (4),

(6)   \begin{equation*} \tilde \phi = A(k, t) \mathrm{e}^{|k|z}. \end{equation*}

Equations (2) and (3) provide

(7)   \begin{equation*} \frac{\partial^2 \tilde \phi}{\partial t^2} + g\frac{\partial \tilde \phi}{\partial z} = 0 \quad \mbox{at}\quad z = 0, \end{equation*}


    \[ \frac{\partial^2 A}{\partial t^2} + g|k| A = 0 \]

and so

    \[ A(k, t) = a(k) \mathrm{e}^{\mathrm{i} wt} + b(k) \mathrm{e}^{-\mathrm{i}wt}, \ w^2 = g|k|. \]

Substituting back into (6) then gives

    \[ \tilde\phi (k, z, t) = \left[ a(k) \mathrm{e}^{\mathrm{i}wt} + b(k) \mathrm{e}^{-\mathrm{i}wt} \right]\mathrm{e}^{|k|z}. \]


    \[ \tilde\phi = 0 \quad \mbox{when} \quad t = 0 \]

and so

    \[ \tilde \phi (k, z, t) = a(k) (\mathrm{e}^{\mathrm{i}wt} - \mathrm{e}^{-\mathrm{i}wt}) \mathrm{e}^{|k|z}. \]

Moreover, substitution into the Hankel transform of (2) results in

(8)   \begin{equation*} \tilde\zeta (k, t) = \frac{\mathrm{i}w}{g}\, a(k) (\mathrm{e}^{\mathrm{i}wt} + \mathrm{e}^{-\mathrm{i}wt}). \end{equation*}


    \[ \tilde\zeta (k, t) = \tilde \zeta_0 (k) \]

when t = 0 and so

    \[ a(k) = -\frac{\mathrm{i}g \tilde\zeta_0 (k)}{2w}. \]

Thus, from (8), we have

    \[ \tilde\zeta (k, t) = \frac{1}{2} (\mathrm{e}^{\mathrm{i}wt} + \mathrm{e}^{-\mathrm{i}wt}) \tilde\zeta_0 (k). \]

Taking the inverse Hankel transform results in

    \[ \zeta (r, t) = \int^\infty_0 kJ_0 (kr) \cos wt\, \tilde \zeta_0 (k) \mathrm{d} k \]


    \[ \tilde\zeta_0 (k) = \int^\infty_0 r\, J_0 (kr) \zeta_0 (r) \mathrm{d} r. \]

Let us define the initial disturbance to be

(9)   \begin{equation*} \zeta_0 (r) = \begin{cases} a, \quad r < \rho,\\ 0, \quad r > \rho, \end{cases} \end{equation*}

at time t = 0.


    \[ \tilde\zeta_0 (k) = a\int^\rho_0 r\, J_0 (kr) \mathrm{d} r = \frac{a\rho}{k}\, J_1 (k\rho). \]


(10)   \begin{equation*} \zeta (r, t) = a\rho \int^\infty_0 J_0 (kr) J_1 (k\rho)\cos wt \mathrm{d} k. \end{equation*}

The integral form of the Bessel function is

    \[ J_0 (kr) = \frac{2}{\pi} \int^{\pi/2}_0 \cos (kr \cos \theta) \mathrm{d}\theta. \]

We may therefore express (10) as follows:

(11)   \begin{equation*} \zeta (r, t) = \operatorname{Re} \bigg\{ \frac{a\rho}{\pi}\int^\infty_0 \int^{\pi/2}_0 J_1 (k\rho) \left(\mathrm{e}^{\mathrm{i}(wt + kr \cos \theta)} + \mathrm{e}^{\mathrm{i}(-wt + kr \cos \theta)}\right) \bigg\} \mathrm{d}\theta \mathrm{d} k \end{equation*}

where \operatorname{Re} denotes the real part of the double integral.

Interchange the order of integration and define

(12)   \begin{equation*} k(\theta) = \int^\infty_0 J_1 (k\rho) \left[\mathrm{e}^{\mathrm{i}(wt + kr \cos \theta)} + \mathrm{e}^{\mathrm{i}(-wt + kr \cos \theta)}\right] \mathrm{d} k. \end{equation*}

With the change of variables \tau = t\sqrt{g\rho}/2r and p = \sqrt{k\rho}, then (12) becomes

(13)   \begin{equation*} k(\theta) =\frac{2}{\rho} \int^\infty_0 J_1 (p^2) \left[\mathrm{e}^{\mathrm{i}(p^2 \cos \theta + 2p \tau)r/\rho} + \mathrm{e}^{\mathrm{i}(p^2 \cos \theta - 2p\tau)r/\rho}\right] p \mathrm{d} p. \end{equation*}

Here we assume {r}/{\rho} \gg 1. In this case, the trigonometric functions will vary extremely rapidly from -1 to +1 with the positive and negative components cancelling each other out. This will happen except where the function f(p) = p^2 \cos \theta \pm 2p \tau is stationary in the interval of integration. We observe that this will not be the case for f(p) = p^2 \cos \theta + 2p \tau since the stationary point is outside the range of integration. Thus the leading order term in (13) is the second (and consequently the first may be neglected). Furthermore

    \[ f' (p) = 2p \cos \theta - 2\tau \]

implying that the stationary point p_0 is

    \[ p_0 = \tau/\cos \theta. \]

Excluding \theta in the neighbourhood of \pi/2, we may write

    \[ f(p) = [(p - p_0)^2 - p^2_0] \cos \theta. \]

Then the integral (13) may be approximated by

    \begin{align*} k(\theta) &\sim \frac{2}{\rho} \int^{p_0 +\epsilon}_{p_0-\epsilon} J_1 (p^2) \\ & \qquad \exp \left[\mathrm{i}\left(\cos \theta \, (p-p_0)^2 \frac{r}{\rho} - \cos \theta \, p^2_0 \frac{r}{\rho} \right)\right] p \mathrm{d} p \\ &= \frac{2}{\rho} \exp \left(-\mathrm{i}\cos \theta \, p^2_0 \frac{r}{\rho}\right) \\ & \qquad \int^{p_0 +\epsilon}_{p_0-\epsilon} J_1 (p^2) \exp \left[\mathrm{i}\cos \theta \, (p-p_0)^2 \frac{r}{\rho}\right] p \mathrm{d} p, \end{align*}

with the integral outside the range [p_0 - \epsilon, p_0 + \epsilon] cancelling to zero.

The factors p and J_1 (p^2) are slowly varying compared with the rapidly fluctuating cosine and sine functions and so may be replaced by p_0 and J_1 (p^2_0), respectively.


    \begin{multline*} k(\theta) \sim \frac{2}{\rho} \exp \left[-\mathrm{i}\cos \theta\, p^2_0\, \frac{r}{\rho}\right] p_0 J_1 (p^2_0) \\ \int^{p_0 +\epsilon}_{p_0-\epsilon} \exp \left[\mathrm{i}\cos \theta \, (p-p_0)^2 \frac{r}{\rho}\right] \mathrm{d} p. \end{multline*}

We now make the following further change of variables:

    \[ s = q(p-p_0), \quad q^2 = \frac{r}{\rho} \cos \theta \]

so that

    \[ k(\theta) \sim \frac{2}{\rho} \exp \left[-\mathrm{i}\cos \theta \, p^2_0 \frac{r}{\rho}\right] \frac{p_0}{q}\, J_1 (p^2_0) \int^{q\epsilon}_{-q\epsilon} \mathrm{e}^{\mathrm{i}s^2} \mathrm{d} s. \]

The limits in the integral can be considered infinite if r/\rho \to \infty, implying q \to \infty for a small, but fixed \epsilon assuming \cos \theta \ne 0 (that is, excluding \theta in the neighbourhood of \pi/2).


(14)   \begin{equation*} \int^\infty_{-\infty} \mathrm{e}^{\pm \mathrm{i}s^2} \mathrm{d} s = \sqrt{\pi} \mathrm{e}^{\pm \mathrm{i}\pi/4} \quad \mbox{(see [3])} \end{equation*}

so that we may write

    \[ k(\theta) \sim \frac{2\sqrt{\pi}}{\rho} \exp \left[-\mathrm{i}\left(\cos \theta \, p^2_0 \,\frac{r}{\rho} - \frac{\pi}{4}\right)\right] \frac{p_0}{q}\, J_1 (p^2_0). \]


    \[ p_0 = \dfrac{t\sqrt{g\rho}}{2r \cos \theta}, \]

so p_0 \ll 1 when {\rho}/{r} \ll 1 and consequently J_1 (p^2_0) \sim {p^2_0}/{2}, and we can write

(15)   \begin{equation*} k(\theta) \sim \frac{\sqrt{\pi}}{\rho q}\, p^3_0 \exp \left[-\mathrm{i}\left(\cos \theta \, p^2_0 \frac{r}{\rho} - \frac{\pi}{4}\right)\right]. \end{equation*}

Substituting back into equation (11) results in

(16)   \begin{equation*} \zeta (r, t) \sim \operatorname{Re} \left\{\frac{a}{\sqrt{\pi}} \int^{\pi/2}_0 \frac{p^3_0}{q} \exp \left[-\mathrm{i}\left(\cos \theta \, p^2_0 \frac{r}{\rho} - \frac{\pi}{4}\right)\right] \mathrm{d}\theta\right\}. \end{equation*}

On introducing the volume displaced by the initial pulse V_0 = \pi \rho^2 a, equation (16) becomes

(17)   \begin{multline*} \zeta(r, t) \sim \\ \operatorname{Re} \left\{\frac{V_0}{\pi^{3/2} \rho^2} \int^{\pi/2}_0 \frac{p^3_0}{q} \exp \left[-\mathrm{i}\left(\cos\theta \, p^2_0 \frac{r}{\rho} - \frac{\pi}{4}\right)\right] \mathrm{d}\theta\right\}. \end{multline*}

We shall consider \rho \to 0 such that V_0 remains finite (and constant). Introducing a final change of variables u = gt^2/4r, (17) becomes

(18)   \begin{equation*} \zeta (r, t) \sim\operatorname{Re} \left\{\frac{V_0}{r^2} \left(\frac{u}{\pi}\right)^{3/2} \int^{\pi/2}_0 \frac{1}{\cos^{7/2} \theta} \mathrm{e}^{-\mathrm{i}(u/\cos \theta - \pi/4)} \mathrm{d}\theta \right\}. \end{equation*}

We now consider the asymptotic limit as u\to \infty, i.e. t^2/r \gg 1, or t^2 \gg r \gg \rho.

Once again we note that for large u (and \cos \theta \ne 0) the integrand oscillates extremely rapidly and so defining

    \[ g(\theta) = u/\cos \theta - \frac{\pi}{4}, \]

we observe that g'(\theta) = 0 implies \theta = \theta_0 = 0. Further g''(\theta_0) = 1 and so

    \[ g(\theta) \simeq g(\theta_0) + \frac{(\theta - \theta_0)^2}{2!} g'' (\theta_0) = u\left( 1 + \frac{\theta^2}{2}\right) - \frac{\pi}{4}. \]

We may therefore replace (18) by

    \[ \zeta (r, t) \sim \operatorname{Re} \left\{\frac{V_0}{r^2} \left(\frac{u}{\pi}\right)^{3/2} \mathrm{e}^{-\mathrm{i}\pi/4}\int^\epsilon_{-\epsilon} \mathrm{e}^{-\mathrm{i}u(1+\theta^2/2)} \mathrm{d}\theta \right\} \]

and, since the integral is essentially annihilated (cancels itself out) for \theta \in (-\infty, -\epsilon) \cap (\epsilon, \infty) we may write

    \[ \zeta (r, t) \sim \operatorname{Re} \left\{\frac{V_0}{r^2} \left(\frac{u}{\pi}\right)^{3/2} \mathrm{e}^{-\mathrm{i}\pi/4} \int^\infty_{-\infty} \mathrm{e}^{-\mathrm{i}u(1+\theta^2/2)} \mathrm{d}\theta\right\}. \]

Note that (\cos \theta)^{-7/2} \simeq 1 for \theta \in [-\epsilon, \epsilon] and is slowly varying for \theta \in [0, \pi/2).


    \begin{align*} \zeta (r, t) &\sim \frac{V_0}{r^2} \left(\frac{u}{\pi}\right)^{3/2} \operatorname{Re} \left\{\mathrm{e}^{-\mathrm{i}(u+\pi/4)} \int^\infty_{-\infty} \mathrm{e}^{-\mathrm{i}{u\theta^2}/{2}} \mathrm{d}\theta\right\}\\ & \sim \frac{V_0}{r^2} \left(\frac{u}{\pi}\right)^{3/2} \operatorname{Re} \left\{\mathrm{e}^{-\mathrm{i}u}\sqrt{\frac{2\pi}{u}} \mathrm{e}^{-\mathrm{i}\pi/2}\right\} \end{align*}


    \[ \int^\infty_{-\infty} \mathrm{e}^{-\mathrm{i}u\theta^2/2} \mathrm{d}\theta = \sqrt{\frac{2\pi}{u}} \mathrm{e}^{-\mathrm{i}\pi/4}. \]

Thus, upon taking the real part, we obtain

(19)   \begin{equation*} \zeta (r, t) \sim \frac{\sqrt{2} V_0}{r^2} \left(\frac{u}{\pi}\right) \cos \left(u + \frac{\pi}{2}\right) \end{equation*}

where, recall, u = gt^2/4r.

Figure 3: Concentric waves (expression (19)) resulting from a disturbance at r = 0

The asymptotic expression (19) for the resulting concentric waves are displayed in Figure 3. It should be stressed that this is only a valid approximation when the fluid velocity is small (i.e. |\boldsymbol{u}|^2 can be neglected) and when t^2 \gg r and \rho \ll r.

We note that \zeta becomes infinite at r = 0 (a direct result of letting \rho \to 0 for a fixed V_0); however, in Figure 3, a smoothing function has been added for exposition. The amplitudes of the crests decrease according to u/r^2 or r^{-3}. Furthermore, it is easily seen that the crests follow each other at a distance \Delta r \simeq {8\pi r^2}/{gt^2}. Hence the wavelength is no longer a constant, but increases at constant t with r^2 and decreases at constant r with t^2.

For further reading on the subject and more complicated wave patterns, for example, the wake behind a ship, the reader is referred to [4] and [5].

Sean McKee CMath FIMA
University of Strathclyde

Murilo F. Tomé
University of São Paulo


  1. McKee, S., Tomé, M.F., Ferreira, V.G., Cuminato, J.A., Castelo, A. and Sousa, F.S. (2008) The MAC method, Comput. Fluids, vol. 37, pp. 907930.
  2. Coulson, C.A. (1961) Waves, Oliver and Boyd, Edinburgh.
  3. Spiegel, M.R. (1964) Theory and Problems of Complex Variables, Schaum Publishing Co, NY.
  4. Whitham, G.B. (1999) Linear and Nonlinear Waves, Wiley, New Jersey.
  5. Johnson, R.S. (1997) A Modern Introduction to the Mathematical Theory of Water Waves, Cambridge University Press.

Reproduced from Mathematics Today, April 2018

Download the article, Ripples from a Splashing Drop (pdf)

Image credit: Water Drop. Crown, close. © Evan Spiler |



Leave a Reply

Your email address will not be published. Required fields are marked *