2. Model

We examine the response of the ocean to an atmospheric perturbation (cyclone or anticyclone) as this crosses the coast. We model this process using the linear rotating shallow-water equations forced by a skew-gradient stress. In dimensionless form, the governing equations are

(2.1a)$$\begin{gather} \varepsilon u_t - v ={-} h_x - \varepsilon \varPhi_y, \end{gather}$$

(2.1b)$$\begin{gather}\varepsilon v_t + u ={-} h_y + \varepsilon \varPhi_x, \end{gather}$$

(2.1c)$$\begin{gather}\varepsilon \lambda^2 h_t + u_x + v_y = 0. \end{gather}$$

Here $\varepsilon = U/(\,fL)$ is the Rossby number, defined in terms of the inertial frequency $f>0$ and the velocity $U$ and length scale $L$ of the perturbation, and $\lambda ^2 = L^2/L_D^2 = f^2 L^2/(gH)$ is an inverse Burger number (e.g. Vallis Reference Vallis2017, § 5.1). The domain is the half-plane $x \ge 0$. With the non-dimensionalisation chosen, the gravity wave speed is $(\varepsilon \lambda )^{-1}$. The boundary condition

(2.2)\begin{equation} u=0 \quad \text{at} \ x = 0 \end{equation}

applies along the coast.

The forcing is proportional to $\varepsilon$ – so the ocean response is primarily geostrophic – and defined by the $O(1)$ scalar function $\varPhi$ taken as a localised function of $x \mp t$ and $y$. The upper sign corresponds to a perturbation that travels westwards over the ocean for $t<0$, makes landfall at $t=0$ and continues overland for $t>0$; the lower sign corresponds to a perturbation travelling eastwards overland for $t<0$, crossing the coast at $t=0$ and continuing over the ocean for $t>0$. Figure 1 illustrates the case of the westward-travelling perturbation. To fix ideas, and for the numerical simulations below, we take $\varPhi$ to be the Gaussian

(2.3)\begin{equation} \varPhi(x \mp t,y) = (2{\rm \pi})^{{-}1} \exp[-\tfrac{1}{2}((x\mp t)^2+y^2)]. \end{equation}

This corresponds to an anticyclonic wind stress, which raises the sea surface in its wake. We consider the evolution from a state of rest as $t \to -\infty$. This eliminates any wave generation by adjustment and implies that $u, v, h \to 0$ as $y \to \pm \infty$.

Figure 1. Schematic of the model. An atmospheric perturbation indicated by the white circle travels westwards over the ocean, leaving a geostrophically balanced flow in its wake (the height field $h$ is shown by the colour scale), before crossing the coast and continuing over the land (shown in brown). The case of an eastward-travelling perturbation starting overland is also considered.

We can eliminate $u$ and $v$ from (2.1) to obtain the single equation

(2.4)\begin{equation} (\nabla^2 - \lambda^2) h_t - \varepsilon^2 \lambda^2 h_{ttt} = \nabla^2 \varPhi, \end{equation}

with $\nabla ^2 = \partial _x^2 + \partial _y^2$ being the Laplacian. The associated boundary condition is found by combining (2.1a), (2.1b) and (2.2) to obtain

(2.5)\begin{equation} h_y + \varepsilon h_{xt}= \varepsilon \varPhi_x - \varepsilon^2 \varPhi_{yt} \quad \text{at} \ x = 0. \end{equation}

Note that (2.4) and (2.5) are consistent with the conservation of the total mass

(2.6)\begin{equation} m(t) = \iint_{x \ge 0} h \, \mathrm{d}\kern0.7pt x \,\mathrm{d} y, \end{equation}

which follows from (2.1c) and (2.2): integrating (2.4) and taking (2.5) into account gives

(2.7)\begin{equation} m_t + \varepsilon^2 m_{ttt} = \lambda^{{-}2} \int_{x=0} (\varPhi_x - h_{xt} ) \, \mathrm{d} y = \lambda^{{-}2} \int_{x=0} (\varepsilon^{{-}1} h_y + \varepsilon \varPhi_{yt} ) \, \mathrm{d} y = 0 , \end{equation}

since $h, \varPhi _t \to 0$ as $y \to \pm \infty$, with $m_t = 0$ the only solution once initial conditions are inferred from (2.1).

3. Small-Rossby-number asymptotics

We now solve (2.4) and (2.5) in the small-Rossby-number regime $\varepsilon \ll 1$ with $\lambda =O(1)$ corresponding to standard QG scaling. Expanding

(3.1)\begin{equation} h = {{h}}^{(0)} + \varepsilon {{h}}^{(1)} + \cdots \end{equation}

and introducing this into (2.4) and (2.5) gives, at leading order, the linearised QG equation

(3.2)\begin{equation} (\nabla^2 - \lambda^2) {{h}}^{(0)}_t = \nabla^2 \varPhi \end{equation}

with boundary condition

(3.3)\begin{equation} {{h}}^{(0)}_y = 0 \quad \text{at} \ x = 0. \end{equation}

Since $h \to 0$ as $y \to \pm \infty$, this implies that

(3.4)\begin{equation} {{h}}^{(0)} = 0 \quad \text{at} \ x = 0. \end{equation}

As Dorofeyev & Larichev (Reference Dorofeyev and Larichev1992), Reznik & Grimshaw (Reference Reznik and Grimshaw2002) and Reznik & Sutyrin (Reference Reznik and Sutyrin2005) point out, this boundary condition is inconsistent with the conservation of the total mass associated with ${{h}}^{(0)}$. Indeed integrating (3.2) gives

(3.5)\begin{equation} {{m}}^{(0)}_t= \lambda^{{-}2} \int_{x=0} ( \varPhi_x - {{h}}^{(0)}_{xt} ) \, \mathrm{d} y, \end{equation}

which is not constrained to vanish by (3.4). This is typical of QG dynamics in large domains, with boundary lengths that are $O(\varepsilon ^{-1})$ or larger compared with the deformation radius. In contrast, for $O(1)$ boundary lengths, as often assumed implicitly, (3.3) implies only that ${{h}}^{(0)} = C(t)$ on the boundary, with a function $C(t)$ determined to ensure total mass conservation. See Reznik & Sutyrin (Reference Reznik and Sutyrin2005) for a detailed discussion of the two regimes. Note that the non-vanishing of (3.5) also implies a failure of QG to represent correctly the evolution of the circulation along the coast.

Mass conservation is resolved by considering the correction $\varepsilon {{h}}^{(1)}$ to ${{h}}^{(0)}$. At $O(\varepsilon )$, (2.4) and (2.5) give

(3.6a,b)\begin{equation} (\nabla^2 - \lambda^2) {{h}}^{(1)}_t = 0 \quad \text{with}\ {{h}}^{(1)}_y = \varPhi_x - {{h}}^{(0)}_{xt} \quad \text{at}\ x=0. \end{equation}

These equations can be solved, but not in a way that ensures ${{h}}^{(1)} \to 0$ as $y \to \pm \infty$, since (3.6b) implies the jump

(3.7)\begin{equation} [{{h}}^{(1)}(0,y,t)]_{y \to -\infty}^{y \to \infty} = \int_{x=0} ( \varPhi_x - {{h}}^{(0)}_{xt} ) \, \mathrm{d} y = \lambda^2 {{m}}^{(0)}_t. \end{equation}

The difficulty arises because (3.6b) is not uniformly valid as $y \to \pm \infty$. The solution ${{h}}^{(1)}$ to (3.6) should be regarded as an inner solution, to be matched to an outer solution, ${{H}}^{(1)}(x,Y,t)$ say, valid in the region $Y=\varepsilon y = O(1)$. This outer solution is found from (2.4) and (2.5) to satisfy

(3.8a,b)\begin{equation} {{H}}^{(1)}_{xxt} - \lambda^2 {{H}}^{(1)}_t = 0 \quad \text{with}\ {{H}}^{(1)}_Y + {{H}}^{(1)}_{xt} = 0 \quad \text{at} \ x= 0, \end{equation}

since $\varPhi$ is exponentially small for $Y=O(1)$, together with matching conditions as $Y \to 0^\pm$. It takes the form

(3.9)\begin{equation} {{H}}^{(1)}(x,Y,t) = K(Y,t) \, \mathrm{e}^{- \lambda x}, \end{equation}

as required by (3.8a) and the condition ${{H}}^{(1)} \to 0$ as $x \to \infty$, and corresponds to a Kelvin wave. From (3.8b) the amplitude $K(Y,t)$ satisfies $K_Y - \lambda K_t = 0$, with a jump across $Y=0$ that matches (3.7). This can be rewritten compactly as

(3.10)\begin{equation} K_Y - \lambda K_t = \lambda^2 {{m}}^{(0)}_t \delta(Y), \end{equation}

with $\delta (Y)$ the Dirac distribution, and solved to obtain

(3.11)\begin{equation} K(Y,t) ={-} \lambda^2 {{m}}^{(0)}_t(t + \lambda Y) \, \varTheta({-}Y), \end{equation}

where $\varTheta (Y)$ denotes the Heaviside step function. This represents a long Kelvin wave, propagating with the coast to its right (here southwards) from the point where the atmospheric perturbation makes contact with the coast. The large Kelvin-wave speed $-(\varepsilon \lambda )^{-1}$ ensures that, although the wave has a small $O(\varepsilon )$ amplitude, it carries an $O(1)$ mass, thus balancing the change of the mass ${{m}}^{(0)}$ associated with the leading-order QG solution. We verify that ${{m}}^{(0)} + \varepsilon {{m}}^{(1)} = \mathrm {const.}$ by computing

(3.12)\begin{equation} \varepsilon {{m}}^{(1)}_t = \varepsilon \iint_{x \ge 0} {{h}}^{(1)}_t \, \mathrm{d}\kern0.7pt x \,\mathrm{d} y = \iint_{x \ge 0} {{H}}^{(1)}_t \, \mathrm{d}\kern0.7pt x\,\mathrm{d} Y = \lambda^{{-}1} \int K_t \, \mathrm{d} Y ={-} {{m}}^{(0)}_t, \end{equation}

where the last equality follows from integrating (3.10) in $Y$.

4. Gaussian perturbation

We now obtain explicit predictions and compare them with the results of numerical simulations of the linear shallow-water equations (2.1) in the case of a Gaussian perturbation (2.3). We solve the QG equation (3.2) for ${{h}}^{(0)}$ with boundary condition (3.4) by writing

(4.1a)\begin{equation} \varPhi(x\mp t,y) = \iint\hat\varPhi(k,l)\exp[\mathrm{i}(k(x\mp t) + ly)]\,\mathrm{d} k\,\mathrm{d} l \end{equation}

and

(4.1b)\begin{equation} {{h}}^{(0)}_t(x,y,t) = \iint\hat g(x,k,l)\exp( \mathrm{i}({\mp} k t + ly) ) \, \mathrm{d} k \,\mathrm{d} l, \end{equation}

where $\hat \varPhi (k,l)=(2{\rm \pi} )^{-2} \exp ({-(k^2+l^2)/2})$ and the function $\hat g$ is to be determined. Introducing (4.1) into (3.2) yields

(4.2)\begin{equation} \hat g_{xx} - (l^2 + \lambda^2) \hat g ={-}(k^2 + l^2) \hat \varPhi(k,l) \mathrm{e}^{\mathrm{i} k x}. \end{equation}

Solving and imposing that $\hat g(x=0,k,l)=0$ as required by (3.4) leads to

(4.3)\begin{equation} \hat g(x,k,l) = \frac{k^2+l^2}{k^2 + l^2 + \lambda^2} \left(\mathrm{e}^{\mathrm{i} k x} - \exp\left({-\sqrt{l^2+\lambda^2} x}\right) \right) \hat{\varPhi}(k,l). \end{equation}

This gives a Fourier representation for ${{h}}^{(0)}_t$ and, by integration in time from $t \to -\infty$, ${{h}}^{(0)}$.

We focus on the Kelvin-wave response, which, according to (3.11), is determined by the mass rate of change ${{m}}^{(0)}_t$. We compute ${{m}}^{(0)}_t$ from its definition (3.5), using (4.1), (4.3) and $\int \mathrm {e}^{\textrm {i}ly} \, \mathrm {d} y = 2 {\rm \pi}\delta (l)$, to find

(4.4)\begin{equation} {{m}}^{(0)}_t = 2 {\rm \pi}\int \frac{\mathrm{i} \lambda k - k^2}{\lambda(k^2 + \lambda^2)} \hat \varPhi(k,0)\, \mathrm{e}^{{\mp} \mathrm{i} k t} \, \mathrm{d} k. \end{equation}

Introducing this into (3.11) gives an explicit form for the Kelvin wave generated by the atmospheric perturbation. In particular, at the coast, for $y<0$ and away from the inner region $y=O(1)$,

(4.5)\begin{equation} h(0,y,t) \sim \varepsilon K(\varepsilon y,t) ={-} 2 {\rm \pi}\varepsilon \int \frac{\mathrm{i} \lambda^2 k - \lambda k^2}{k^2 + \lambda^2} \exp({\mp \mathrm{i} k (t + \varepsilon \lambda y)}) \hat \varPhi(k,0) \, \mathrm{d} k + O(\varepsilon^2), \end{equation}

since ${{h}}^{(0)}=0$ for $x=0$.

Figure 2 shows $K(0,t) = - \lambda ^2 {{m}}^{(0)}_t$ as a function of $t$ for $\lambda =0.5$, 1 and 2 and for an atmospheric perturbation travelling westwards, from the open ocean towards the coast (i.e. with the factor $\mathrm {e}^{\mathrm {i} k t}$ in (4.4) and (4.5)). This function approximates the evolution of the height at the location of the landfall as described by the outer solution ${{H}}^{(1)}$; it can be interpreted as the amplitude of a wavemaker generating the Kelvin-wave response to the atmospheric perturbation. Note the asymmetry of the evolution about $t=0$, with a depression phase ($t<0$) that is weaker and lasts longer than the elevation phase ($t>0$). This asymmetry is constrained by the vanishing of the integral of $K(0,t)$ for $t \in (-\infty ,\infty )$, which reflects the fact that the change in the QG mass ${{m}}^{(0)}$ is purely transient with the forcing chosen. The asymmetry increases as $\lambda$ decreases. The case of an atmospheric perturbation travelling eastwards, from overland towards the ocean, is deduced by reversing the sign of $t$.

Figure 2. Evolution of the height at $(x,y)=0$, where the atmospheric perturbation crosses the coast, according to the outer solution (4.5) for $\lambda =0.5$ (orange), $1$ (red) and $2$ (blue).

We illustrate and verify the above predictions by carrying out numerical simulations of the linear shallow-water equations (2.1). The numerical model used discretises the fields $(u,v,h)$ on a staggered grid in the $x$-direction, with $u$ represented at grid points $j \Delta x$ for $j=1,2,\ldots$ and $v$ and $h$ at grid points $(\,j-1/2)\Delta x$. A Fourier expansion is used in the $y$-direction with an assumption of periodicity. The $x$-derivatives are approximated by central differences; the $y$-derivatives are computed spectrally. The domain size is $10 \times 160$, with the $y=0$ axis, corresponding to the path of the atmospheric perturbation, placed at a distance 40 from the northern boundary of the domain. These choices ensure that the distance from the path of the perturbation to the southern boundary is large enough for the long Kelvin wave to propagate south along the coast without being affected by the $y$-periodicity of the domain. The boundary condition $u_x=0$ is imposed at the eastern boundary $x=10$. The physical parameters for the simulations reported are $\lambda =1$ and $\varepsilon =0.1$; the numerical parameters are 128 grid points in $x$, 256 Fourier modes in $y$, and a time step $\Delta t = 0.1 \varepsilon \Delta x$ with $\Delta x = 10/128$.

Figure 3 shows successive snapshots of the height field $h$ for a westward-travelling perturbation as this makes landfall. At $t=-4$, when the centre of the perturbation is a distance 4 away from the coast (recall that the speed of the perturbation is $1$ in dimensionless units), the response of the ocean is well balanced, QG to a good approximation and, in particular, symmetric about the $x$-axis as the QG solution (4.1)–(4.3) indicates. At $t=-2$, there is a weak signal of negative $h$ along the coast and south of the atmospheric perturbation, breaking the $\pm y$ symmetry. This is the signature of the $O(\varepsilon )$ Kelvin wave generated by the mass imbalance associated with the QG response. This initial depression wave propagates south, with its maximum reaching $y=-20$ at $t=0$. It is followed by the propagation of an elevation visible for $t=2$ and $4$.

Figure 3. Height field $h$ for a Gaussian atmospheric perturbation propagating westwards from the open ocean towards the coast, shown at times $t=-4$, $-2$, $0$ (landfall time), $2$ and $4$ from (a) to (e). The parameters are $\lambda =1$ and $\varepsilon =0.1$. The location of the perturbation is indicated in the first three panels by the black circle with radius $1$ corresponding to the length scale of the Gaussian.

A clearer depiction of the Kelvin-wave propagation is given in figure 4, which shows the height field along the coast, where the QG contribution ${{h}}^{(0)}$ vanishes. The figure confirms the validity of the matched-asymptotics prediction (4.5), which approximates the numerical solution with remarkable accuracy except, of course, in the ‘inner’ region $y=O(1)$ where (3.6) applies. Thus the simple picture of a Kelvin wave generated by a wavemaker with time dependence as shown in figure 2 is well justified.

Figure 4. Height field along the coast $x=0$ at times $t=-4$ (blue), $-2$ (red), $0$ (landfall time, orange), $2$ (purple) and $4$ (green) for the simulation in figure 3, showing the generation near $y=0$ and southward propagation of a Kelvin wave. The predictions based on the matched-asymptotics result (4.5) are shown by the black curves closely matching the coloured curves except near $y=0$.

Figures 5 and 6 are the analogues of figures 3 and 4 for an eastward-travelling atmospheric perturbation coming from overland to reach the coast at $t=0$. In this case, the Kelvin-wave generation, visible at $t=-1$ in figure 5, precedes the bulk of the QG response, which appears only around $t=3$, when the perturbation is well away from the coast. As expected from (4.4) and (4.5), the Kelvin wave leads to an initial elevation of the sea surface, followed by a depression, reversing the evolution compared with the westward-travelling case. Figure 6 shows that the Kelvin-wave amplitude is again predicted with high accuracy by the asymptotic formula (4.5).

Figure 5. Same as figure 3 but for an atmospheric perturbation travelling eastwards, from overland to the open ocean, crossing the coast at $t=0$. The height field $h$ is shown at times $t=-1$, $1$, $3$, $5$ and $7$ from (a) to (e). The location of the perturbation is indicated by the black circles.

Figure 6. Height field along the coast $x=0$ at times $t=-1$ (blue), $1$ (red), $3$ (landfall time, orange), $5$ (purple) and $7$ (green) for the simulation in figure 5. The predictions based on the matched-asymptotics result (4.5) are shown by the black curves.