2 Model equations

The model considered in this work applies to stationary, incompressible flows where viscosity and surface tension effects are ignored. Our problem is readily converted to non-dimensional form using the profile depth $h$ and surface current velocity $U_{s}$ , as follows:

(2.1a-d ) $$\begin{eqnarray}(x,y,z)\mapsto (x,y,z)h,\quad \boldsymbol{k}\mapsto \boldsymbol{k}h^{-1},\quad \hat{\boldsymbol{u}}_{t}\mapsto \hat{\boldsymbol{u}}_{t}U_{s},\quad \hat{p}_{t}\mapsto \hat{p}_{t}\unicode[STIX]{x1D70C}U_{s}^{2},\end{eqnarray}$$

$\boldsymbol{k}=(k_{x},k_{y})$ is the wavenumber in the surface plane, $\hat{\boldsymbol{u}}_{t}$ the fluid velocity and $\hat{p}$ the pressure. A solution will be constructed in Fourier space by considering a single mode. Intuitively, one may in this frame regard the problem configuration as that of water flowing over a sinusoidal bed, as sketched in figure 1. The shear flow profile is represented by the power law family (1.1). The mean vertical position of the bed itself is placed at some elevation $z=\unicode[STIX]{x1D6FF}>0$ above the shear profile origin. Thus, the two parameters $q$ and $\unicode[STIX]{x1D6FF}$ govern the current profile felt by the bed, and the Froude number $Fr$ its strength. The bathymetry undulations redirects the current energy to perturb the velocity field and otherwise flat free surface. With the prescribed non-dimensionalisation, the elevation of the surface, $\hat{\unicode[STIX]{x1D701}}_{s}(\boldsymbol{r})$ , is oriented about $z=1$ , just as the bed elevation, $\hat{\unicode[STIX]{x1D701}}_{b}(\boldsymbol{r})$ , is oriented about $z=\unicode[STIX]{x1D6FF}$ . Symbols $\hat{\unicode[STIX]{x1D702}}$ are used when referring the elevation relative to the unperturbed surface height.

Figure 1.

Sketch of the problem set-up.

The problem to be solved consists of the stationary incompressible Euler equations, along with an impermeability condition at the bed and dynamic and kinematic boundary conditions at the surface;

(2.2a ) $$\begin{eqnarray}\displaystyle & \displaystyle \left.\begin{array}{@{}c@{}}(\hat{\boldsymbol{u}}_{t}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}})\hat{\boldsymbol{u}}_{t}+\hat{\unicode[STIX]{x1D735}}\hat{p}_{t}=-Fr^{-2}\boldsymbol{e}_{z}\\ \hat{\unicode[STIX]{x1D735}}\boldsymbol{\cdot }\hat{\boldsymbol{u}}_{t}=0\end{array}\right\}\!;\quad \hat{\unicode[STIX]{x1D701}}_{b}\leqslant z\leqslant \hat{\unicode[STIX]{x1D701}}_{s}, & \displaystyle\end{eqnarray}$$

(2.2b ) $$\begin{eqnarray}\displaystyle & \displaystyle \left.\begin{array}{@{}c@{}}\hat{\boldsymbol{u}}_{t}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}\hat{\unicode[STIX]{x1D701}}_{s}={\hat{w}}\\ \hat{p}_{t}=0\end{array}\right\}\!;\quad z=\hat{\unicode[STIX]{x1D701}}_{s}, & \displaystyle\end{eqnarray}$$

(2.2c ) $$\begin{eqnarray}\displaystyle & \hat{\boldsymbol{u}}_{t}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}\hat{\unicode[STIX]{x1D701}}_{b}={\hat{w}};\quad z=\hat{\unicode[STIX]{x1D701}}_{b}, & \displaystyle\end{eqnarray}$$

$Fr=U_{s}/\sqrt{gh}$

$\hat{\unicode[STIX]{x1D735}}=(\unicode[STIX]{x2202}_{x},\unicode[STIX]{x2202}_{y},\unicode[STIX]{x2202}_{z})$

3 Linearised solution

Begin by separating the perturbed part of the internal flow from the total,

(3.1a ) $$\begin{eqnarray}\displaystyle & \hat{\boldsymbol{u}}_{t}(\boldsymbol{r},z)=\boldsymbol{U}(z)+\hat{\boldsymbol{u}}(\boldsymbol{r},z), & \displaystyle\end{eqnarray}$$

(3.1b ) $$\begin{eqnarray}\displaystyle & \hat{p}_{t}(\boldsymbol{r},z)=(1-z)Fr^{-2}+\hat{p}(\boldsymbol{r},z), & \displaystyle\end{eqnarray}$$

(3.1c ) $$\begin{eqnarray}\displaystyle & \hat{\unicode[STIX]{x1D701}}_{b}(\boldsymbol{r})=\unicode[STIX]{x1D6FF}+\hat{\unicode[STIX]{x1D702}}_{b}(\boldsymbol{r}), & \displaystyle\end{eqnarray}$$

(3.1d ) $$\begin{eqnarray}\displaystyle & \hat{\unicode[STIX]{x1D701}}_{s}(\boldsymbol{r})=1+\hat{\unicode[STIX]{x1D702}}_{s}(\boldsymbol{r}), & \displaystyle\end{eqnarray}$$

(3.2a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D702}_{b}(\boldsymbol{k})={\mathcal{F}}_{\boldsymbol{k}}\hat{\unicode[STIX]{x1D702}}_{b}(\boldsymbol{r}),\quad \unicode[STIX]{x1D702}_{b}(\mathbf{0})=0,\end{eqnarray}$$

where ${\mathcal{F}}_{\boldsymbol{k}}$ denotes a Fourier transform onto a wave vector plane  $\boldsymbol{k}$ . Real-space flow fields become

(3.3) $$\begin{eqnarray}(\hat{\boldsymbol{u}},\hat{p},\hat{\unicode[STIX]{x1D702}}_{s})={\mathcal{F}}_{\boldsymbol{k}}^{-1}(\boldsymbol{u},p,\unicode[STIX]{x1D702}_{s}),\end{eqnarray}$$

by virtue of the lower boundary condition (2.2c ).

3.1 Internal flow

After linearisation, the Euler equations read

(3.4a ) $$\begin{eqnarray}\displaystyle & \text{i}k_{x}U\boldsymbol{u}+\boldsymbol{U}^{\prime }w+\unicode[STIX]{x1D735}p=O[\ldots ], & \displaystyle\end{eqnarray}$$

(3.4b ) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D735}=(\text{i}k_{x},\text{i}k_{y},\unicode[STIX]{x2202}_{z})$ . The error estimate is $O[(k^{3}+k^{-1})w^{2}]$ for the $x$ -component and $O[(k^{2}+1)w^{2}]$ for the $y$ - and $z$ -components, assuming $k_{x}\sim k_{y}\sim k$ . Eliminating  $u$ , $v$ and $p$ yields the Rayleigh equation

(3.5) $$\begin{eqnarray}w^{\prime \prime }-(k^{2}+U^{\prime \prime }/U)w=O[(k^{3}+k^{-1})w^{2}],\end{eqnarray}$$

$k=(k_{x}^{2}+k_{y}^{2})^{1/2}$ . Upon substituting $w(z)=\sqrt{z}W(z)$ (Phillips et al. Reference Phillips, Wu and Lumley1996), the linear part of (3.5) reduces to the Bessel equation

(3.6) $$\begin{eqnarray}W^{\prime \prime }+\frac{W^{\prime }}{z}-\left(k^{2}+\frac{(1-2q)^{2}}{4z^{2}}\right)W=0\end{eqnarray}$$

whose homogeneous solutions are well known; two such exist which are superposed with the coefficients $c^{+}$ and  $c^{-}$ , to be determined from the boundary conditions. The horizontal velocity components are found by integrating (3.4a ). We write the solution for a generic flow variable $\unicode[STIX]{x1D719}$ in the form

(3.7) $$\begin{eqnarray}\unicode[STIX]{x1D719}=\unicode[STIX]{x1D716}\mathop{\sum }_{\pm }\unicode[STIX]{x1D719}^{\pm }c^{\pm },\end{eqnarray}$$

with

(3.8a-c ) $$\begin{eqnarray}w^{\pm }=\text{i}\sqrt{z}I_{\pm (q-1/2)}(kz),\quad p^{\pm }=\frac{k_{x}}{k}U(z)\sqrt{z}I_{\pm (q+1/2)}(kz),\quad \boldsymbol{u}_{h}^{\pm }=\frac{\text{i}\boldsymbol{e}_{x}U^{\prime }w^{\pm }-\boldsymbol{k}p^{\pm }}{k_{x}U},\end{eqnarray}$$

$I_{\unicode[STIX]{x1D6FC}}$ being the modified Bessel function of the first kind of order $\unicode[STIX]{x1D6FC}$ , $\boldsymbol{u}_{h}=(u,v)$ the horizontal velocity vector and $\boldsymbol{e}_{x}=(1,0)$ is the horizontal unit vector. The function $\unicode[STIX]{x1D716}$ represents the intensity of the internal flow perturbations, to be related to the intensity of the flow disturbance from the bathymetry. This we have separated from the weight coefficients $c^{\pm }$ which are related to the unperturbed flow alone. The summation over ‘ $\pm$ ’ means a sum is taken over ‘ $+$ ’ and ‘ $-$ ’, being either signs or labels as context requires. (This form of (3.8), written in terms of the modified Bessel function, becomes linearly dependent at $q=1/2$ , but is adopted here for symmetry rather than resorting to the second kind of modified Bessel function. Limit values $q\rightarrow 1/2$ are evaluated in place of $q=1/2$ .) Note that a stream function $\hat{u} =\unicode[STIX]{x2202}_{z}\hat{\unicode[STIX]{x1D713}}$ , ${\hat{w}}=-\unicode[STIX]{x2202}_{x}\hat{\unicode[STIX]{x1D713}}$ may easily be constructed in the case of two-dimensional flow. One finds

(3.9) $$\begin{eqnarray}\hat{\unicode[STIX]{x1D713}}(x,z)=\frac{z^{q+1}}{q+1}+{\mathcal{F}}_{\boldsymbol{k}}^{-1}(\text{i}w/k_{x});\quad [k_{y}\equiv 0].\end{eqnarray}$$

3.2 Boundary condition

The lower kinematic boundary conditions (2.2c ), linearised about $z=\unicode[STIX]{x1D6FF}$ , reads

(3.10) $$\begin{eqnarray}w(\unicode[STIX]{x1D6FF})=\unicode[STIX]{x1D716}+O[(k+k^{3})\unicode[STIX]{x1D702}_{b}^{2}]\end{eqnarray}$$

in Fourier space, where

(3.11) $$\begin{eqnarray}\unicode[STIX]{x1D716}=\text{i}k_{x}U(\unicode[STIX]{x1D6FF})\unicode[STIX]{x1D702}_{b}.\end{eqnarray}$$

Linearised about $z=1$ , the upper boundary conditions (2.2b ) yield

(3.12a ) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D702}_{s}-Fr^{2}p(1)=O[Fr^{2}w(1)^{2}], & \displaystyle\end{eqnarray}$$

(3.12b ) $$\begin{eqnarray}\displaystyle & \text{i}k_{x}\unicode[STIX]{x1D702}_{s}-w(1)=O[(k+k^{-1})w(1)^{2}]. & \displaystyle\end{eqnarray}$$

Equation (3.4) has here been used and $U$ derivatives assumed order unity to simplify error estimates. It is tempting to relate these estimates back to the amplitudes of the bed, but, as will be seen, these variables may not scale in a linear manner.

Either of the surface conditions (3.12) yields $\unicode[STIX]{x1D702}_{s}$ directly. The remaining coefficients $c^{+}$ and $c^{-}$ are readily determined by the other surface condition and the kinematic condition at the bed, equation (3.10). In terms of the functions in (3.8),

(3.13) $$\begin{eqnarray}c^{\pm }=\left[w^{\pm }(\unicode[STIX]{x1D6FF})-\frac{w^{\pm }(1)-\text{i}k_{x}Fr^{2}p^{\pm }(1)}{w^{\mp }(1)-\text{i}k_{x}Fr^{2}p^{\mp }(1)}w^{\mp }(\unicode[STIX]{x1D6FF})\right]^{-1}.\end{eqnarray}$$

Note that $c^{\pm }$ approaches finite values both as $Fr^{2}\rightarrow 0$ and $Fr^{2}\rightarrow \infty$ .

Figure 2.

$(k_{x}^{2}/k)Fr_{crit}^{2}$ and $Fr_{crit}(k_{y}=0)$ as functions of $q$ and wavenumber for two values of  $\unicode[STIX]{x1D6FF}$ .

3.3 Critical Froude number

From (3.13) it is seen that there may exist particular parameter values for which both field coefficients diverge. We shall later examine the impact of such singular cases in a spectral bathymetry representation. At present though, a single mode is considered. The critical Froude number $Fr_{crit}$ , at which the coefficient denominator is zero, is

(3.14) $$\begin{eqnarray}\frac{k_{x}^{2}}{k}Fr_{crit}^{2}(\boldsymbol{k})=\frac{I_{q-1/2}(k)I_{-q+1/2}(k\unicode[STIX]{x1D6FF})-I_{-q+1/2}(k)I_{q-1/2}(k\unicode[STIX]{x1D6FF})}{I_{q+1/2}(k)I_{-q+1/2}(k\unicode[STIX]{x1D6FF})-I_{-q-1/2}(k)I_{q-1/2}(k\unicode[STIX]{x1D6FF})},\end{eqnarray}$$

which is found to be always real positive for $q>0$ , $0<\unicode[STIX]{x1D6FF}<1$ . Plots of $Fr_{crit}^{2}$ are shown in figure 2. For uniform current profiles, $q=0$ , one recovers the well known result (Thomson 1886) $(k_{x}Fr_{crit})^{2}/k=\tanh [k(1-\unicode[STIX]{x1D6FF})]$ and $\{\coth [k(1-\unicode[STIX]{x1D6FF})]-{k^{-1}\}}^{-1}$ for linear ones ( $q=1$ ). All profiles $0\leqslant q\leqslant 1$ have the short wave asymptote $(k_{x}Fr_{crit})^{2}/k\rightarrow 1$ as $k\rightarrow \infty$ . From figure 2 it is further seen that $Fr_{crit}$ diminishes monotonically with increasing wavenumber. The largest critical Froude number, graphed in figure 3, is

(3.15) $$\begin{eqnarray}Fr_{crit,max}=\lim _{k\rightarrow 0}Fr_{crit}=\frac{k}{|k_{x}|}\sqrt{\frac{1-\unicode[STIX]{x1D6FF}^{1-2q}}{1-2q}}\end{eqnarray}$$

and marks the distinction between subcritical and supercritical flows over bathymetries represented with an infinite Fourier spectrum, as shall be expounded upon later. Also, $(k_{x}Fr_{crit})^{2}/k\rightarrow I_{-(q-1/2)}(k)/I_{-(q+1/2)}(k)$ in the limit $\unicode[STIX]{x1D6FF}\rightarrow 0$ if $q<1/2$ , agreeing with the result of Lighthill (Reference Lighthill1953) for a standing wave over a flat bed with a $q=1/7$ profile. The critical Froude number for the $q=1/2$ profile is obtained by considering the limit, yielding

(3.16) $$\begin{eqnarray}\left.\frac{k_{x}^{2}}{k}Fr_{crit}^{2}\right|_{q\rightarrow 1/2}=\frac{I_{0}(k)K_{0}(k\unicode[STIX]{x1D6FF})-K_{0}(k)I_{0}(k\unicode[STIX]{x1D6FF})}{I_{1}(k)K_{0}(k\unicode[STIX]{x1D6FF})+K_{1}(k)I_{0}(k\unicode[STIX]{x1D6FF})},\end{eqnarray}$$

where $K_{\unicode[STIX]{x1D6FC}}$ is the modified Bessel function of the second kind. In the limit $k_{x}=k\rightarrow 0$ and $q=0$ one retrieves the expected $Fr_{crit}^{2}=1-\unicode[STIX]{x1D6FF}$ , or just 1 if the Froude number is instead defined based on the actual depth $h-\unicode[STIX]{x1D6FF}$ in physical units. (It is also common to define the Froude number based on the mean velocity for which $Fr_{m}=Fr(1-\unicode[STIX]{x1D6FF}^{q+1})/[(1-\unicode[STIX]{x1D6FF})(1+q)]$ .)

Singular states are a common artefact of resonances when linearisation has been performed. With the allowance of transient dynamics, these singularities are often associated with algebraic wave growth as current energy is transferred to the perturbed flow field without feedback (Benney Reference Benney1961; Craik Reference Craik1970; Akselsen & Ellingsen Reference Akselsen and Ellingsen2019). The asymptotic growth of critical flow over a bump or under a pressure patch goes in linear theory as the cube root of time (Akylas Reference Akylas1984; Cole Reference Cole1985), but a linear steady state exists if the forcing is compact and the flow unrestricted in the spanwise direction (Katsis & Akylas Reference Katsis and Akylas1987). In the setting of a monochromatic bathymetry with a critical Froude number, resonance can be interpreted as the phase velocity being zero relative to the bed, denying dispersion of the bathymetry induced disturbances. This issue can, however, be resolved and a steady state reached when considering the full nonlinear problem. Classical Stokes wave theory provides a conceptual indication of how; as the level of nonlinearity increases with growing wave amplitude so does the Stokes wave phase velocity relative to that of a linear wave. The shift in phase velocity detunes the Stokes wave such that it cannot remain stationary relative to the bed. Resonance is thus broken. Mei (Reference Mei1969) was the first to consider the Stokes wave solution near criticality. He discovered triple-valued steady state solutions. In order to deduce which solutions are likely occur in nature, Sammarco, Mei & Trulsen (Reference Sammarco, Mei and Trulsen1994) examined the transient problem and its initial stability and long-time evolution.

4 Asymptotic behaviour

Consider next the composite solution as a spectrum of modes in a Fourier integral. It was shown in § 3.3 that the linear theory will, in the frame of a monochromatic bathymetry, break down for some particular combinations of parameters as a form of resonance is encountered. Excluding these parameter combinations, and their immediate neighbourhoods, the derived solution is valid. When, on the other hand, unevenness in bed topography is of limited spatial extent the Fourier integral must cover all values of  $\boldsymbol{k}$ . If the flow is subcritical (there exists wavenumbers $\boldsymbol{k}$ so that $Fr_{crit}(\boldsymbol{k})>Fr$ ) then a pair of singular points will be encountered during the integration. The problem as it stands is indeterminate until a radiation condition is applied. Following a procedure due to Rayleigh (1883) (see also Lamb (Reference Lamb1932, p. 399)) the momentum equations are furnished with an artificial frictional force, $-\unicode[STIX]{x1D707}(\hat{\boldsymbol{u}}-\boldsymbol{U})$ . The only requirement for this force is that it always dampens perturbations and that it is vanishingly small, i.e., that $\unicode[STIX]{x1D707}$ be infinitesimally small yet positive everywhere. This provides a direction of time and ensures that causality is obeyed.

Without loss of generality, let this force vary in $z$ according to $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{1}k_{x}U(z)$ where $\unicode[STIX]{x1D707}_{1}=0^{+}\operatorname{sgn}(k_{x})$ . The resulting modified Euler equations are

(4.1) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\text{i}k_{x}U(1-\text{i}\unicode[STIX]{x1D707}_{1})\boldsymbol{u}+\boldsymbol{U}^{\prime }w+\unicode[STIX]{x1D735}p=0,\\ \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0,\end{array}\right\}\end{eqnarray}$$

which, to linear order in $\unicode[STIX]{x1D707}_{1}$ , alters the solution (3.7)–(3.8) as follows:

(4.2a,b ) $$\begin{eqnarray}w^{\pm }=\text{i}z^{(1/2)+(\text{i}/2)q\unicode[STIX]{x1D707}_{1}}I_{\pm (\tilde{q}-1/2)}(kz),\quad \displaystyle p^{\pm }=\frac{k_{x}}{k}(1-\text{i}\unicode[STIX]{x1D707}_{1})z^{1/2+\tilde{q}}I_{\pm (\tilde{q}+1/2)}(kz),\end{eqnarray}$$

(4.3) $$\begin{eqnarray}\boldsymbol{u}_{h}=\frac{\text{i}\boldsymbol{e}_{x}U^{\prime }w-\boldsymbol{k}p}{(1-\text{i}\unicode[STIX]{x1D707}_{1})k_{x}U},\end{eqnarray}$$

where $\tilde{q}=q(1+\text{i}\unicode[STIX]{x1D707}_{1}/2)$ . Now rewrite the $c^{\pm }$ -coefficients by inserting (4.2) into (3.13),

(4.4) $$\begin{eqnarray}c^{\pm }=\frac{F^{\pm }(\boldsymbol{k})}{H(\boldsymbol{k})},\end{eqnarray}$$

with

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle F^{\pm }(\boldsymbol{k})=\pm \left[I_{\mp (\tilde{q}-1/2)}(k)-Fr^{2}\frac{k_{x}^{2}}{k}(1-\text{i}\unicode[STIX]{x1D707}_{1})I_{\mp (\tilde{q}+1/2)}(k)\right], & \displaystyle\end{eqnarray}$$

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle H(\boldsymbol{k})=\text{i}\unicode[STIX]{x1D6FF}^{(1/2)+(\text{i}/2)q\unicode[STIX]{x1D707}_{1}}\mathop{\sum }_{\pm }I_{\pm (\tilde{q}-1/2)}(k\unicode[STIX]{x1D6FF})F^{\pm }(\boldsymbol{k}). & \displaystyle\end{eqnarray}$$

The denominator $H$ is shared by all flow variables and is a dispersion relation for surface waves atop the flowing water. Its roots may be seen as representing waves propagating upstream with phase velocity exactly equal to surface velocity, so that the flow is stationary.

Consider now a two-dimensional flow ( $k_{y}=0$ ). A localised unevenness in the bathymetry generates, for any flow variable, a perturbation of the surface of the form

(4.7) $$\begin{eqnarray}\int _{-\infty }^{\infty }\frac{\text{d}k_{x}}{2\unicode[STIX]{x03C0}}\frac{G(k_{x})}{H(k_{x})}\text{e}^{\text{i}k_{x}x}.\end{eqnarray}$$

For example, if computing the surface elevation $\hat{\unicode[STIX]{x1D702}}_{s}$ from (3.12a ) then

(4.8) $$\begin{eqnarray}G=-\unicode[STIX]{x1D716}\,Fr^{2}\frac{k_{x}}{k^{2}}(1-\text{i}\unicode[STIX]{x1D707}_{1})\frac{2}{\unicode[STIX]{x03C0}}\cos (\unicode[STIX]{x03C0}\tilde{q})\end{eqnarray}$$

– note that Lord Kelvin’s solution for a uniform current (Thomson 1886, p. 520; Lamb Reference Lamb1932, p. 409),

(4.9) $$\begin{eqnarray}\unicode[STIX]{x1D702}_{s}=\unicode[STIX]{x1D702}_{b}\left\{\cosh [k(1-\unicode[STIX]{x1D6FF})]-\frac{k}{Fr^{2}k_{x}^{2}}\sinh [k(1-\unicode[STIX]{x1D6FF})]\right\}^{-1},\end{eqnarray}$$

is recovered when putting $q=0$ . The integral (4.7) may be evaluated using well known techniques from complex analysis. The denominator $H(k_{x})$ has conjugate roots along the imaginary axis whose contributions make up the near field. Two roots $k_{x}=\pm \unicode[STIX]{x1D6FC}+\text{i}0^{+}$ also appear if the flow is subcritical. These are both shifted slightly away from the real axis into the upper imaginary plane by the artificial friction force. Subcritical flow here means that there exists a wavenumber such that $Fr<Fr_{crit}(\boldsymbol{k})$ . As seen before, this is equivalent to checking whether $Fr<Fr_{crit,max}$ , equation (3.15). Above this Froude number current transport is too great to permit disruptions of the surface to propagate upstream or remain stationary.

Although it is possible to express the integral (4.7) considering these singular points only (Lamb Reference Lamb1932, p. 410), a form more suitable for computation is obtained by adding and subtracting the pole terms (leading term of the integrand’s Laurent series expansion) inside the integral whenever the flow is subcritical and using Cauchy’s integral theorem to write (4.7) in the form

(4.10) $$\begin{eqnarray}\int _{-\infty }^{\infty }\frac{\text{d}k_{x}}{2\unicode[STIX]{x03C0}}\left[\frac{G(k_{x})}{H(k_{x})}-\mathop{\sum }_{\pm }\frac{G(\pm \unicode[STIX]{x1D6FC})}{H^{\prime }(\pm \unicode[STIX]{x1D6FC})(k_{x}\mp \unicode[STIX]{x1D6FC})}\right]\text{e}^{\text{i}k_{k}x}+\text{i}\mathop{\sum }_{\pm }\frac{G(\pm \unicode[STIX]{x1D6FC})}{H^{\prime }(\pm \unicode[STIX]{x1D6FC})}\text{e}^{\pm \text{i}\unicode[STIX]{x1D6FC}x}\unicode[STIX]{x1D6E9}(x)\end{eqnarray}$$

when $Fr_{crit,max}>Fr$ . A prime denotes differentiation with respect to $k_{x}$ and $\unicode[STIX]{x1D6E9}$ is the Heaviside function. The poles that hinder numerical convergence are thus eliminated from the integrand and their effect, a wave train emitted downstream, appears explicitly in the solution. If the bathymetry obstruction is localised and centred near $x=0$ , we may regard the integral term of (4.7) as a ‘near field’ and the oscillating term as a ‘far field’. No alteration to (4.7) is needed when $Fr_{crit,max}<Fr$ (supercritical) since this flow consists of a near field only, the real axis being free of poles.

It is worth noting that when the flow is subcritical the shape of the bed deformation is of consequence only for the amplitude of the emitted wave train while the current profile alone dictates the wavelength. In a dynamic sense the downstream far field wave will have exactly the wavelength of an upstream wave propagating towards the left at the surface velocity so as to be stationary in the ‘laboratory’ frame. The dispersion relation of such a wave is a functional of the shear profile $U(z)$ and is, for arbitrary $U(z)$ , readily calculated numerically; cf., e.g., Li & Ellingsen (Reference Li and Ellingsen2019). Figure 4, which is essentially an inversion of figure 2, shows the wavenumber of the wave train for given parameters  $q$ , $\unicode[STIX]{x1D6FF}$ and  $Fr$ . Supercritical flow states, for which wave trains cease to appear, are excluded from the figure.

In two-dimensional problems, equation (4.10) is readily computed numerically with the derivative of $H(k_{x})$ evaluated with discrete differentiation. An example is shown in figure 5. The extended method for obtaining accuracy in adherence to conditions at the lower boundary, presented in § 5, is here employed.

In three dimensions the field separation procedure which yielded (4.10) becomes more involved as roots of $G(\boldsymbol{k})$ form a continuous path in $\boldsymbol{k}$ -space along which an integration must be performed. Instead, we opt in three-dimensional problems for enforcing the radiation condition by allocating a small but finite value of the artificial friction coefficient $\unicode[STIX]{x1D707}_{1}$ such that the original Fourier integral converges. This is a tried and trusted method which is equivalent to introducing a low level of artificial viscosity, see, e.g., Moisy & Rabaud (Reference Moisy and Rabaud2014). The artificial friction thus introduced must be made strong enough for the far field to die out before re-entering the periodically bounded numerical domain, yet weak enough not to mask any salient features pertaining to the inviscid flow or noticeably affect the far field wavelength, which obtains a correction of order  $\unicode[STIX]{x1D707}_{1}^{2}$ . The strategy can add to the computational load since a large computational domain may be required. Many more sophisticated methods are available such as non-reflecting boundary conditions, yet given the modest overall computational effort in our examples the simplest approach is adequate.

A final observation is warranted. The limit $\unicode[STIX]{x1D6FF}\rightarrow 0$ of the denominator $H$ in (4.6) is finite if $q<1/2$ and infinite if $q>1/2$ . This means that infinitesimally small bathymetry perturbations centred around the current stagnation point $z=0$ are felt within the flow only if $\text{d}U/\text{d}z>1/\sqrt{z}$ at $z=0$ . Profiles developing slower than this are not felt by a bed of infinitesimal amplitude. The feature can be seen in (3.15) and figure 3 where critical Froude numbers are infinite for $q>1/2$ as $\unicode[STIX]{x1D6FF}\rightarrow 0$ .

Figure 3.

Critical Froude number $Fr_{crit,max}$ , equation (3.15), for the existence of a far field wave train.

Figure 4.

Wavenumber $\unicode[STIX]{x1D6FC}$ of far field wave train as a function of  $q$ , $\unicode[STIX]{x1D6FF}$ and  $Fr$ .

Figure 5.

Uniform flow over a cylinder of radius 0.2. The surface, bed, streamlines and dynamic pressure field $\hat{p}+Fr^{2}z$ are shown. $Fr=0.5$ . The extended lower boundary conditions of § 5 low are applied, which does not affect the features of the wave train.

5 A lower boundary condition of $n$ th-order accuracy

Recall that the bathymetry $\hat{\unicode[STIX]{x1D702}}_{b}$ is prescribed and notice that the lower boundary condition (2.2c ) therefore contains no products of unknowns. This is convenient for describing the lower boundary condition precisely, without requiring that the bathymetry perturbations be very small.

Many classical works have employed coordinate transformation techniques to two-dimensional problems (Benjamin Reference Benjamin1959; Stokes Reference Stokes1880; Forbes & Schwartz Reference Forbes and Schwartz1982), but these are not so readily extendible to three-dimensional non-potential flow. A perturbation approach is instead viable, and the evaluation of (2.2c ) at $z=\hat{\unicode[STIX]{x1D701}}_{b}$ is achieved by Taylor expansion about $z=\unicode[STIX]{x1D6FF}$ . This generates products between variable derivatives and powers of  $\hat{\unicode[STIX]{x1D702}}_{b}$ , and so a perturbation solution

(5.1) $$\begin{eqnarray}\hat{\unicode[STIX]{x1D719}}=\mathop{\sum }_{n=1}^{\infty }\hat{\unicode[STIX]{x1D719}}_{n},\quad \text{where}~\hat{\unicode[STIX]{x1D719}}_{n}=O(\hat{\unicode[STIX]{x1D702}}_{b}^{n});\quad \unicode[STIX]{x1D719}\in \{\boldsymbol{u},p,\unicode[STIX]{x1D702}_{s}\}\end{eqnarray}$$

is postulated. Inserting this series into the bottom boundary condition (2.2c ) and using Taylor expansion about $z=\unicode[STIX]{x1D6FF}$ , one obtains

(5.2) $$\begin{eqnarray}\mathop{\sum }_{m=0}^{\infty }\frac{\hat{\unicode[STIX]{x1D702}}_{b}^{m}}{m!}\mathop{\left[U^{(m)}\unicode[STIX]{x2202}_{x}\hat{\unicode[STIX]{x1D702}}_{b}+\mathop{\sum }_{n=1}^{\infty }(\hat{\boldsymbol{u}}_{n}^{(m)}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}\hat{\unicode[STIX]{x1D702}}_{b}-{\hat{w}}_{n}^{(m)})\right]}\nolimits_{z=\unicode[STIX]{x1D6FF}}=0,\end{eqnarray}$$

where parenthesised superscripts indicate order of derivative in  $z$ . Sorting into consecutive orders of $\hat{\unicode[STIX]{x1D702}}_{b}$ and imposing that equality holds for each order individually, one finds the lower boundary condition

(5.3) $$\begin{eqnarray}{\hat{w}}_{n}(\unicode[STIX]{x1D6FF})=\hat{\unicode[STIX]{x1D716}}_{n},\end{eqnarray}$$

where

(5.4) $$\begin{eqnarray}\hat{\unicode[STIX]{x1D716}}_{n}=\hat{\unicode[STIX]{x1D702}}_{b}^{n}\mathop{\left[\frac{U^{(n-1)}}{(n-1)!}\frac{\unicode[STIX]{x2202}_{x}\hat{\unicode[STIX]{x1D702}}_{b}}{\hat{\unicode[STIX]{x1D702}}_{b}}+\mathop{\sum }_{m=1}^{n-1}\frac{1}{m!}\left(m\frac{\hat{\boldsymbol{u}}_{h,n-m}^{(m-1)}}{\hat{\unicode[STIX]{x1D702}}_{b}^{n-m}}\boldsymbol{\cdot }\frac{\hat{\unicode[STIX]{x1D735}}\hat{\unicode[STIX]{x1D702}}_{b}}{\hat{\unicode[STIX]{x1D702}}_{b}}-\frac{{\hat{w}}_{n-m}^{(m)}}{\hat{\unicode[STIX]{x1D702}}_{b}^{n-m}}\right)\right]}\nolimits_{z=\unicode[STIX]{x1D6FF}}.\end{eqnarray}$$

For bathymetries consisting of many Fourier modes (typically an obstruction localised in space) the most computationally efficient procedure is to go back and forth between physical and Fourier space using fast Fourier transforms. That is, having computed $\boldsymbol{u}_{n-1}$ in Fourier space, the physical field is readily obtained via the transformation $\hat{\boldsymbol{u}}_{n-1}^{(m)}={\mathcal{F}}\boldsymbol{u}_{n-1}^{(m)}$ . After evaluating (5.4), the result is transformed back to Fourier space, $\unicode[STIX]{x1D716}_{n}={\mathcal{F}}^{-1}\hat{\unicode[STIX]{x1D716}}_{n}$ . The only required adjustment to the linear solution presented in § 3 is to replace in (3.7) $\unicode[STIX]{x1D716}$ with increasing orders of  $\unicode[STIX]{x1D716}_{n}$ . The horizontal derivatives of the velocity field are readily available via the Fourier transformation, or simply using finite differences. The $z$ -derivatives are computed in Fourier space – we have constructed a building block function $f_{a,b,m}\equiv \unicode[STIX]{x2202}_{z}^{m}[z^{b-1/2}I_{a}(kz)]$ evaluated using binomial statements,

(5.5) $$\begin{eqnarray}f_{a,b,m}(z)=\mathop{\sum }_{j=0}^{m}\left(\frac{k}{2}\right)^{j}\binom{m}{m-j}\frac{\unicode[STIX]{x1D6E4}\left(b+{\textstyle \frac{1}{2}}\right)z^{b+j-m-1/2}}{\unicode[STIX]{x1D6E4}\left(b+j-m+{\textstyle \frac{1}{2}}\right)}\mathop{\sum }_{i=0}^{j}\binom{j}{i}I_{a+j-2i}(kz),\end{eqnarray}$$

whence

(5.6a ) $$\begin{eqnarray}\displaystyle & \displaystyle w_{n}^{(m)}(z)=\text{i}\unicode[STIX]{x1D716}_{n}\mathop{\sum }_{\pm }f_{\pm (q-1/2),1,m}(z)c^{\pm }, & \displaystyle\end{eqnarray}$$

(5.6b ) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}_{h,n}^{(m)}(z)=-\unicode[STIX]{x1D716}_{n}\mathop{\sum }_{\pm }\left[\frac{q\boldsymbol{e}_{x}}{k_{x}}f_{\pm (q-1/2),0,m}(z)+\frac{\boldsymbol{k}}{k}f_{\pm (q+1/2),1,m}(z)\right]c^{\pm }, & \displaystyle\end{eqnarray}$$

(5.6c ) $$\begin{eqnarray}\displaystyle & \displaystyle p_{n}^{(m)}(z)=\unicode[STIX]{x1D716}_{n}\frac{k_{x}}{k}\mathop{\sum }_{\pm }f_{\pm (q+1/2),q+1,m}(z)c^{\pm }, & \displaystyle\end{eqnarray}$$

(5.6d ) $$\begin{eqnarray}\displaystyle & \displaystyle U^{(m)}(z)=z^{q-m}\mathop{\prod }_{j=0}^{m-1}(q-j)=z^{q-m}\frac{\unicode[STIX]{x1D6E4}(q+1)}{\unicode[STIX]{x1D6E4}(q-m+1)}, & \displaystyle\end{eqnarray}$$

$m=0,1,2,\ldots .$

The convergence the Taylor expansion of the above expressions will by Taylor’s theorem be limited to the disc within which the represented function is analytical. There is a singularity at $z=0$ if the current profile is curved ( $q$ different from 0 or 1) which imposes the restriction $\unicode[STIX]{x1D702}_{b}\leqslant \unicode[STIX]{x1D6FF}$ . Even so, since the first few terms of the Taylor series typically decrease in magnitude before growing again, good approximations are often obtained by truncating the Taylor series one term before the smallest in the manner recommended by Bender & Orszag (Reference Bender and Orszag1991, § 3.5).

We can still compute cases where fully curved current profiles encounter large-amplitude obstructions by instead adopting an iterative approach. One then evaluates an error to the boundary condition (2.2c ) in real space after having computed $\hat{\boldsymbol{u}}(\hat{\unicode[STIX]{x1D701}}_{b})={\mathcal{F}}^{-1}\{\boldsymbol{u}(\hat{\unicode[STIX]{x1D701}}_{b})\}$ at all points $(x,y)$ . This error is then Fourier transformed and fed back into the amplitudes as an augmentation of  $\unicode[STIX]{x1D716}$ ; $\unicode[STIX]{x1D716}_{new}=\unicode[STIX]{x1D716}_{old}+error$ . Convergence properties are usually good, presumably because $\unicode[STIX]{x1D716}=w(\unicode[STIX]{x1D6FF})\sim \text{i}k_{x}U(\unicode[STIX]{x1D6FF})$ dominates over $\boldsymbol{u}_{h}(\unicode[STIX]{x1D6FF})\ast \text{i}\boldsymbol{k}\unicode[STIX]{x1D702}_{b}$ in the cases tested herein.

The products in (5.4) constitute convolutions in Fourier space. Going back-and-forth between physical and Fourier space circumvents computing the multi-dimensional mode interactions which appear within the convolution. If instead individual mode interactions are considered, the solution takes the form

(5.7) $$\begin{eqnarray}\hat{\boldsymbol{u}}_{n}={\mathcal{F}}_{\boldsymbol{k}_{1}}^{-1}{\mathcal{F}}_{\boldsymbol{k}_{2}}^{-1}\ldots {\mathcal{F}}_{\boldsymbol{ k}_{n}}^{-1}\boldsymbol{u}_{n},\end{eqnarray}$$

where $\boldsymbol{u}_{n}$ is a function of the interaction wave vector

(5.8) $$\begin{eqnarray}\unicode[STIX]{x1D73F}_{n}=\mathop{\sum }_{m=1}^{n}k_{n},\end{eqnarray}$$

(replacing $\boldsymbol{k}$ in (5.5)–(5.6)) except for in $\unicode[STIX]{x1D716}_{n}$ where

(5.9) $$\begin{eqnarray}\unicode[STIX]{x1D716}_{n}=\mathop{\left[\text{i}k_{nx}\frac{U^{(n-1)}}{(n-1)!}\mathop{\prod }_{j=1}^{n}\unicode[STIX]{x1D702}_{b}(\boldsymbol{k}_{j})+\mathop{\sum }_{m=1}^{n-1}\frac{1}{m!}(m\text{i}\boldsymbol{k}_{n}\boldsymbol{\cdot }\boldsymbol{u}_{h,n-m}^{(m-1)}-w_{n-m}^{(m)})\mathop{\prod }_{j=n-m+1}^{n}\unicode[STIX]{x1D702}_{b}(\boldsymbol{k}_{j})\right]}\nolimits_{z=\unicode[STIX]{x1D6FF}}.\end{eqnarray}$$

The lower-order velocity components $\boldsymbol{u}_{m-n}(z;\unicode[STIX]{x1D73F}_{m-n})$ are in turn functions of $\unicode[STIX]{x1D716}_{m-n}$ and so on. Higher harmonics are therefore present even when the bathymetry is a sinusoidal shape provided its amplitude is finite. We will later demonstrate the appearance of higher harmonics in a flow over a sinusoidal bed of finite amplitude and the possibility of higher-order resonance.

Although the nonlinear lower boundary condition greatly extends the range of well approximated flow scenarios, the solution will still be restricted to linear order within the flow field and at the surface. However, there are two cases, $q=0$ and two-dimensional flow with $q=1$ , where the undulated velocity field $\hat{\boldsymbol{u}}$ is inherently irrotational. Potential theory is permissible in these two cases which are consequently the two scenarios where the most theoretical headway has been made. Irrotationality will in the Rayleigh equation (3.5) manifest by its inhomogeneous right-hand term, assumed small in linear theory, becoming identically zero so that the Rayleigh equation reflects the Laplace equation of potential theory (see e.g., Akselsen & Ellingsen Reference Akselsen and Ellingsen2019). In the first case, when also the unperturbed current is uniform and irrotational ( $q=0$ ), flow irrotationality is intuitive since inviscid perturbations are incapable of generating new vorticity. The second case is when the unperturbed current is linear and perturbations are everywhere aligned with the current (two-dimensional flows with $q=1$ ). This latter case becomes rotational, however, once the flow has non-trivial spanwise dependence, in particular if the bathymetry is three-dimensional or not parallel to the current because current vortex lines are made to twist if they pass through an undulated velocity field (Ellingsen Reference Ellingsen2016). Except for in these special cases, the linearised internal flow puts a limit on the steepness $k\unicode[STIX]{x1D702}_{b0}$ of the bathymetry, meaning that perturbations should be shallow ( $k\ll 1$ ) when current profiles are curved. Surface boundary conditions are also approximate to linear order so that perturbations in the region $z=1$ ought to be made small.

Perturbation strategies for resolving the weakly nonlinear dynamics in the entirety of the flows are straightforward in principle (see, e.g., Akselsen & Ellingsen (Reference Akselsen and Ellingsen2019) for an example), but they quickly become unwieldy.

8 Summary

We have presented a perturbation solution for a free-surface shear flow following a $z^{q}$ velocity profile flowing over an arbitrarily but moderately sloped bed. For various positive values of $q$ this analytical form can model bottom boundary layers, linearly sheared flows and near-surface shear layers. Radiation conditions are incorporated to generate the appropriate asymptotic behaviour.

Bathymetries of finite amplitude are found to introduce higher-order harmonics into the solution. For particular combinations of flow profile and bathymetry shapes these may resonate with the surface wave modes, producing nonlinear contributions to the downstream surface pattern of higher-order which may become dominant. Resonances occurring in flows over a sinusoidal bed manifest, in two-dimensional flows, as endless wave trains trailing an obstruction. The expressions for criticality have been derived and a number of observations concerning the relationship between current vorticity, flow criticality, solution nonlinearity and the shape of the surface have been made related to the sinusoidal-bed configuration. Recirculation is observed to occur within deep depressions in the bathymetry when shear is strong. Further, three-dimensional bed–current vorticity interaction phenomena have been demonstrated, revealing for the case of flow at oblique angles to a sinusoidal bed both helical curving and spanwise migration of streamlines. Finally, examples of flows encountering localised obstructions of various shapes were presented in both two and three dimensions.