# Attached and separated rotating flow over a finite height ridge

S. Frei Department of Mathematics & Statistics, University of Konstanz,
Universitätsstr. 10, 78457 Konstanz, Germany<EMAIL_ADDRESS>E.
Burman  E. Johnson Department of Mathematics, University College London, UK,
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

This paper discusses the effect of rotation on the boundary layer in high
Reynolds number flow over a ridge using a numerical method based on stabilised
finite elements that captures steady solutions up to Reynolds number of order
$10^{6}$. The results are validated against boundary layer computations in
shallow flows and for deep flows against experimental observations reported in
Machicoane et al. (Phys. Rev. Fluids, 2018). In all cases considered the
boundary layer remains attached, even at large Reynolds numbers, provided the
Rossby number of the flow is sufficiently small. At any fixed Rossby number
the flow detaches at sufficiently high Reynolds number to form a steady
recirculating region in the lee of the ridge. At even higher Reynolds numbers
no steady flow is found. This disappearance of steady solutions closely
reproduces the transition to unsteadiness seen in the laboratory.

##  1 Introduction

Mason and Sykes Mason and Sykes (1979) observe that topography with horizontal
length scales of the order of a few kilometres is often sufficiently steep to
induce flow separation in the atmosphere and that such flows might be expected
to produce significant changes in the momentum exchange between the atmosphere
and the surface. In addition to exerting a form drag on the fluid, separated
flow over topography leads to enhanced mixing in the lee of the topography,
affecting the lee-wave field and may make a significant contribution to the
overall vertical mixing in the region Richards  _et al._ (1992). More
recently, Machicoane  _et al._ Machicoane  _et al._ (2018) present theoretical
and experimental results for flow past a horizontal cylinder in horizontal
motion relative to a frame rotating rapidly about a vertical axis. Their work
discusses in detail the inertial lee-wave field generated by the cylinder,
extending significantly the results of Johnson (1982), and also includes a
discussion of experimental results for the onset, with increasing translation
speed of the cylinder, of separation and subsequent unsteadiness in the flow.
They note that rotation strongly suppresses separation and suggest, from the
experimental results, that with increasing flow speed the flow changes
directly from being steady and attached to being unsteady and detached as the
product of the Rossby and Reynolds numbers ($Ro=U/2\Omega L$, $Re=\nu/UL$ for
flow speed $U$, rotation rate $\Omega$, kinematic viscosity $\nu$ and typical
horizontal scale $L$) for the flow passes through $275\pm 25$ . This contrasts
with non-rotating flow where separation occurs at $Re\approx 3$ but wake
instability at $Re\approx 47$, Williamson (1996); Machicoane  _et al._ (2018).
The purpose of the present work is to examine this lee separation in rotating
flows in greater detail.

To obtain the critical Reynolds numbers for this nonlinear flow over a wide
range of Rossby numbers, a stable numerical method is required, which is
robust independently of the Reynolds number. For this purpose, we employ a
stabilised equal-order finite element method using Continuous Interior Penalty
Stabilisation, going back to Douglas & Dupont Douglas and Dupont (1976).
Concerning fluid equations, it has first been used by Burman, Hansbo &
Fernández for the Stokes, Darcy and Oseen equations Burman and Hansbo (2006);
Burman  _et al._ (2006) and later for more complex problems, such as flows at
high Reynolds numbers Tong  _et al._ (2022), including turbulent flows Moura
_et al._ (2022), as well as fluid-structure interactions Fernández (2013);
Burman  _et al._ (2020).

The structure of this article is as follows. Section 2 introduces the geometry
and the governing equations. Section 3 introduces a streamfunction formulation
for finite height ridges and obtains the outer and inner, boundary-layer
solution at large Reynolds number for shallow flow.The numerical method for
the boundary layer solution is given in Appendix A. Section 4 introduces the
numerical method, showing that for attached flow the results for increasing
Reynolds number converge towards the boundary layer solutions, but also
capturing steady separated flows at larger Rossby numbers. Section 5 extends
the results to deep flow allowing inertial waves to appear in the outer flow
but showing the same qualitative behaviour of separation and unsteadiness as
shallow flows, albeit at smaller Reynolds numbers. Section 6 investigates deep
rotating flow over a horizontal cylinder, the geometry of Machicoane  _et al._
Machicoane  _et al._ (2018), and shows that the numerical simulations capture
much of the experimental observations. Section 7 gives brief conclusions.

##  2 Governing equations

Consider a fluid of constant density $\rho$ and viscosity $\nu$ in solid body
rotation about a vertical axis at angular velocity $\Omega$. Introduce
Cartesian axes $Ox^{*}y^{*}z^{*}$ fixed in the fluid with $Oz^{*}$ vertical.
Let the fluid lie above a horizontal boundary perturbed by a ridge of height
$h$ and horizontal scale $L$, given by $z^{*}=h\;b(x^{*}/L)$, where the
function $b$ gives the ridge profile and has a maximum value unity. Let the
fluid be bounded above at height $z^{*}=2D$ by a horizontal boundary perturbed
downwards by a complementary ridge so that the flow geometry is symmetric
about the line $z^{*}=D$ and the motion has aspect ratio $\delta=D/L$ (Figure
1).

2*DhxzU2*DyL

Figure 1: Flow over a ridge. We assume a constant inflow $U$ in horizontal
direction, while the geometry rotates around the vertical $z$-axis.

Now suppose the fluid is set in steady motion, relative to the rotating frame,
so that at large distance, outside viscous layers on the boundaries, the flow
is uniform with constant speed $U$ in the positive $x^{*}$ direction. The
Coriolis force associated with this uniform flow is balanced by a constant
pressure gradient of magnitude $2\rho\Omega U$ in the negative $y^{*}$
direction. In the absence of the ridge ($h=0$), and far upstream of the ridge
when $h>0$, the flow is the usual Ekman layer flow with a boundary layers of
thickness $(\nu/\Omega)^{1/2}$ on the upper and lower boundaries. The flow is
independent of the coordinate $y^{*}$ and so non-dimensionalising lengths on
$L$, velocities on $U$ and the pressure on $\rho U^{2}$ gives the two-
dimensional governing equations

$\displaystyle uu_{x}+wu_{z}-\text{Ro}^{-1}v$
$\displaystyle=-p_{x}+\text{Re}^{-1}\nabla^{2}u,$ (1a) $\displaystyle
uv_{x}+wv_{z}+\text{Ro}^{-1}(u-1)$ $\displaystyle=\text{Re}^{-1}\nabla^{2}v,$
(1b) $\displaystyle uw_{x}+ww_{z}$
$\displaystyle=-p_{z}+\text{Re}^{-1}\nabla^{2}w,$ (1c) $\displaystyle
u_{x}+w_{z}$ $\displaystyle=0,$ (1d)

with the constant in (1b) giving the constant far-field pressure gradient in
the $y$-direction. On the top and bottom boundaries the flow satisfies the
viscous boundary conditions that the velocity vanishes., i.e. $u=v=w=0$. The
two non-dimensional parameters appearing are the Reynolds number
$\text{Re}=UL/\nu$, measuring the ratio of inertial to viscous forces, and the
Rossby number $\text{Ro}=U/2\Omega L$, measuring the ratio of inertial to
Coriolis forces. Since the nonlinear inertial terms are always significant in
the flows discussed below, Re and Ro are taken, together with the aspect ratio
$\delta$, as the fundamental parameters. The Ekman number $E=\nu/2\Omega
D^{2}=\text{Ro}/\delta^{2}\text{Re}$ appears naturally with $E^{1/2}$
measuring the ratio of the Ekman layer thickness to the fluid depth and is
treated as derived from the fundamental parameters.

Integrating (1d) across the fluid depth gives the dimensional result

$\int u^{*}\;\text{d}z^{*}=2UD,$ (2)

showing the mass flux per unit width in the $y$ direction is constant and
equal to its upstream value. The dimensional form of (1b), outside boundary
layers, can be written

$(u^{*}v^{*})_{x^{*}}+(w^{*}v^{*})_{z^{*}}=-2\Omega(u^{*}-U),$ (3)

Integrating (3) over the flow domain, using the divergence theorem, the
vanishing of normal velocity on the upper and lower boundaries, the uniform
conditions at large distance and (2), gives

$2UD(V_{+}-V_{-})=-2\Omega\int(2UD-UH^{*}(x^{*}))\;\text{d}x^{*}=-4\Omega UA,$
(4)

where $H^{*}(x^{*})$ is the local depth, $A$ is the cross-sectional area of
the ridge and $V_{-}$ and $V_{+}$ are the values of the cross-flow speed far
upstream and far downstream. The generation of negative relative vertical
vorticity as fluid columns are compressed moving over the ridge thus generates
a clockwise change in flow direction determined by the difference

$V_{+}-V_{-}=-2\Omega A/D,$ (5)

which becomes progressively less important the deeper the flow. Since the flow
geometry is independent of $y$ any value of $V_{-}$ can be enforced by an
appropriate choice of the upstream pressure gradient $p_{x}$.

Sufficiently far upstream of the ridge the velocity field is uniform
horizontally and an Ekman layer solution satisfies the full nonlinear system
(1) for arbitrary $E$ and Ro, giving the solution

$\displaystyle
u^{*}_{0}+\text{i}v^{*}_{0}=(U+\text{i}V_{-})\\{1-\frac{\cosh[(1+\text{i})(2E\delta^{2})^{-1/2}(z-\delta)]}{\cosh[(1+\text{i})(2E)^{-1/2}]}\\}.$
(6)

The non-dimensional shear stress on the lower boundary is thus

$\frac{\text{d}\phantom{z}}{\text{d}z}(u_{0}+\text{i}v_{0})\big{|}_{z=0}=(2E\delta^{2})^{-1/2}(1+\text{i})(1+\text{i}V_{-}/U)\tanh[(1+\text{i})(2E)^{-1/2}],$
(7)

rotated by an angle $\pi/4$ from the outer velocity. For sufficiently small
Ekman number the $x$-component of shear stress is negative if $V_{-}$ exceeds
$U$, which it will for any ridge of non-zero cross-sectional area if the
rotation is sufficiently rapid and the interior flow is symmetric about the
ridge. The equations are no longer parabolic in the direction of increasing
$x$. This difficulty is removed here by taking advantage of the $y$
independence of the flow to set $V_{-}=0$ far upstream (so $p_{x}\to 0$ as
$x\to-\infty$) so the upstream shear stress is always in the $x$ direction.
For small Ekman numbers the $x$-component of wall-shear far downstream of the
ridge is then, from (5),

$(2E\delta^{2})^{-1/2}(U+2\Omega A/D)$ (8)

and so increases due to the clockwise rotation of the exterior flow. This
rotation also affects the structure of the boundary layer and is discussed in
greater detail in §3. The boundary conditions on the flow can thus be written

$\displaystyle(u,v,w)$ $\displaystyle=(0,0,0),\qquad\text{on }z=\alpha\;b(x),$
(9) $\displaystyle(u_{z},v_{z},w)$ $\displaystyle=(0,0,0),\qquad\text{on
}z=\delta,$ (10) $\displaystyle(u,v,w)$
$\displaystyle\to(u_{0},v_{0},0)\big{|}_{V_{-}=0},\qquad\text{as
}x\to-\infty,$ (11)

where the non-dimensional parameter, $\alpha=h/L$, gives the average slope of
the ridge and $\delta=D/L$ gives the channel aspect ratio. The flow conditions
far downstream, outside the Ekman layers, follow from (5) as

$v\to-2\Omega A/UD,\qquad p_{x}\to-2\Omega^{2}AL/DU^{2},\qquad\text{as
}x\to\infty.$ (12)

##  3 Boundary layer analysis, $\text{Re}\gg 1$

The continuity equation (1d) allows the introduction of a streamfunction
$\psi$ defined through

$u=\psi_{z},\qquad w=-\psi_{x}.$ (13)

Then the momentum equations become

$\displaystyle-\frac{\partial(\psi,u)}{\partial(x,z)}-\text{Ro}^{-1}v$
$\displaystyle=-p_{x}+\text{Re}^{-1}\nabla^{2}u,$ (14a)
$\displaystyle-\frac{\partial(\psi,v)}{\partial(x,z)}+\text{Ro}^{-1}(u-1)$
$\displaystyle=\text{Re}^{-1}\nabla^{2}v,$ (14b)
$\displaystyle-\frac{\partial(\psi,w)}{\partial(x,z)}\phantom{+\text{Ro}^{-1}u}$
$\displaystyle=-p_{z}+\text{Re}^{-1}\nabla^{2}w,$ (14c)

where $\frac{\partial(f,g)}{\partial(x,y)}:=f_{x}g_{y}-f_{y}g_{x}$. Cross-
differentiation of (14a) and (14c) gives with (14b)

$\displaystyle-\frac{\partial(\psi,\eta)}{\partial(x,z)}$
$\displaystyle=\text{Ro}^{-1}v_{z}+\text{Re}^{-1}\nabla^{2}\eta,$ (15a)
$\displaystyle-\frac{\partial(\psi,v)}{\partial(x,z)}$
$\displaystyle=\text{Ro}^{-1}(\psi_{z}+1)+\text{Re}^{-1}\nabla^{2}v,$ (15b)

where $\eta=u_{z}-w_{x}=\psi_{xx}+\psi_{zz}=\nabla^{2}\psi$ is the $y$
component of vorticity. The boundary conditions on $\psi$ and $v$ are

$\displaystyle(\psi,\psi_{z},v)$ $\displaystyle=(0,0,0),\qquad\text{on
}z=\alpha\;b(x),$ (16) $\displaystyle(\psi,\psi_{zz},v)$
$\displaystyle=(\delta,0,0),\qquad\text{on }z=\delta,$ (17)
$\displaystyle(\psi,v)$ $\displaystyle\to(\int
u_{0}\text{d}z,v_{0}),\qquad\text{as }x\to-\infty.$ (18)

###  3.1 The outer flow at large Reynolds number

In the limit, $Re\to\infty$ with all variables of order unity, (15),(16),(18)
reduce to

$\displaystyle-\frac{\partial(\psi,\eta)}{\partial(x,z)}$
$\displaystyle=\text{Ro}^{-1}\frac{\partial(x,v)}{\partial(x,z)},$ (19a)
$\displaystyle-\frac{\partial(\psi,v)}{\partial(x,z)}$
$\displaystyle=\text{Ro}^{-1}\frac{\partial(x,z-\psi)}{\partial(x,z)},$ (19b)

giving an inviscid outer flow, subject to the inviscid boundary conditions

$\displaystyle\psi$ $\displaystyle=0,\qquad\text{on }z=\alpha\;b(x),$ (20)
$\displaystyle\psi$ $\displaystyle=\delta,\qquad\text{on }z=\delta,$ (21)
$\displaystyle(\psi,v)$ $\displaystyle\to(z,0),\qquad\text{as }x\to-\infty.$
(22)

For unidirectional motion (where $\psi_{z}=u\neq 0$ throughout the flow
domain) it is convenient to map the flow domain to a rectangle by taking
independent variables $(x,\psi)$ and dependent variable $Z(x,\psi)$.
Multiplying equations (19) by the Jacobian
$|\partial(x,z)/\partial(x,\psi)|=Z_{\psi}$ gives

$\displaystyle-\frac{\partial(\psi,\eta)}{\partial(x,\psi)}$
$\displaystyle=\text{Ro}^{-1}\frac{\partial(x,v)}{\partial(x,\psi)},$ (23a)
$\displaystyle-\frac{\partial(\psi,v)}{\partial(x,\psi)}$
$\displaystyle=\text{Ro}^{-1}\frac{\partial(x,z-\psi)}{\partial(x,\psi)},$
(23b)

i.e.

$\displaystyle\eta_{x}$ $\displaystyle=\text{Ro}^{-1}v_{\psi},$ (24a)
$\displaystyle v_{x}$ $\displaystyle=-\text{Ro}^{-1}Z_{\psi}+\text{Ro}^{-1},$
(24b)

which can be combined to give

$Z_{\psi\psi}+\text{Ro}^{2}\eta_{xx}=0.$ (25)

Although the system (24) and equation (25) appear linear, in $(x,\psi)$
coordinates

$\displaystyle\eta$
$\displaystyle=\nabla^{2}\psi=(1/2)[(1+Z_{x}^{2})/Z_{\psi}^{2}]_{\psi}+(Z_{x}/Z_{\psi})_{x},$
(26)

closing the problem for $Z(x,\psi)$ in $(x,\psi)$ space but rendering the
field equation nonlinear.

The far-field boundary condition (22) becomes

$Z(x,\psi)\to\psi,\qquad\text{as }|x|\to\infty.$ (27)

The boundary condition on the ridge becomes simply a condition on the
coordinate line $\psi=0$, giving a problem in a rectangular domain with

$\displaystyle Z(x,0)$ $\displaystyle=\alpha\;b(x),$ (28a) $\displaystyle
Z(x,\delta)$ $\displaystyle=\delta,$ (28b) $\displaystyle Z(x,\psi)$
$\displaystyle\to\psi,\qquad\text{as }|x|\to\infty.$ (28c)

These equations are equivalent to those derived for stratified, rotating flow
by Jacobs Jacobs (1964), although the change in transverse velocity noted here
is not discussed there.

In $(x,\psi)$ coordinates the Cartesian velocity components become

$u(x,\psi)=1/Z_{\psi},\qquad w(x,\psi)=Z_{x}/Z_{\psi},$ (29)

giving the outer flow velocity component along a streamline as

$U_{\text{out}}(x,\psi)=(u^{2}+w^{2})^{1/2}=(1+Z_{x}^{2})^{1/2}/Z_{\psi}.$
(30)

###  3.2 The inner, boundary-layer solution

For Rossby numbers of order unity, the outer solution matches the no-slip
condition on the ridge in a viscous layer of thickness $\text{Re}^{-1/2}$.
Thus introduce the coordinate $\zeta=\text{Re}^{1/2}\psi$ orthogonal to the
surface $\psi=0$ and the scaled velocity component normal to the surface
$W=\text{Re}^{1/2}u_{\text{n}}$ for normal velocity component $u_{\text{n}}$
which vanishes as $\zeta\to\infty$. Then the boundary layer is governed by the
equations

$\displaystyle UU_{\xi}+WU_{\zeta}-\text{Ro}^{-1}V+P_{\xi}$
$\displaystyle=U_{\zeta\zeta},$ (31a) $\displaystyle
UV_{\xi}+WV_{\zeta}+\text{Ro}^{-1}(1+U)$ $\displaystyle=V_{\zeta\zeta},$ (31b)
$\displaystyle P_{\zeta}$ $\displaystyle=0,$ (31c) $\displaystyle
U_{\xi}+W_{\zeta}$ $\displaystyle=0,$ (31d)

where $\xi$ is the distance along the surface and $U$ the speed in the $\xi$
direction. Choosing the origins of $\xi$ and $x$ to coincide gives

$\xi(x)=\int_{0}^{x}\sqrt{1+\alpha^{2}b^{{}^{\prime}2}(x^{\prime})}\;\text{d}x^{\prime}.$
(32)

The boundary conditions on system (31) are

$\displaystyle[U,V,W]$ $\displaystyle=[0,0,0],$
$\displaystyle\text{on}\quad\zeta=0,$ (33a) $\displaystyle[U,V]$
$\displaystyle\to[U_{e}(\xi),V_{e}(\xi)],$
$\displaystyle\text{as}\quad\zeta\to\infty,$ (33b) $\displaystyle[U,V]$
$\displaystyle=[U_{0}(\zeta),V_{0}(\zeta)],$
$\displaystyle\text{at}\quad\xi=\xi_{0},$ (33c) $\displaystyle P$
$\displaystyle\to P_{e}(\xi),$ $\displaystyle\text{as}\quad\zeta\to\infty,$
(33d)

where, from (30), the exterior velocity components are given by
$U_{e}(\xi)=U_{\text{out}}(x,0)$, $V_{e}(\xi)=v(x,0)$, $P_{e}(\xi)=p(x,0)$ for
$x=x(\xi)$, obtained by inverting (32) and $U_{0}(\zeta)$, $V_{0}(\zeta)$ are
the usual Ekman layer velocities, valid for arbitrary Re sufficiently far
upstream of the ridge, at some $-\xi_{0}\gg 1$, so, from (6),

$U_{0}(\zeta)+\text{i}V_{0}(\zeta)=U\\{1-\exp[-(1+\text{i})(2\text{Ro})^{-1/2}\zeta]\\}.$
(34)

For sufficiently low ridges, system (31) is parabolic and so can be
numerically integrated by stepping downstream from the known upstream profile
at $\xi=\xi_{0}$.

###  3.3 The shallow flow limit, $\delta\ll 1$

In general the solution of (25) and (26) can be found only numerically.
However for shallow flows where $\delta\ll 1$, (25) reduces to
$Z_{\psi\psi}=0$ with solution

$Z(x,\psi)=\alpha b(x)(1-\psi/\delta)+\psi,$ (35)

giving the depth-independent horizontal velocity components

$u(x,\psi)=[1-\Delta b(x)]^{-1},\quad
v(x,\psi)=-(\Delta/\text{Ro})\int_{-\infty}^{x}b(x^{\prime})\text{d}x^{\prime},$
(36)

depending solely on the ridge shape and the fractional height
$\Delta=\alpha/\delta$. The velocity along the lower boundary is thus, from
(30),

$\displaystyle U_{e}(\xi)$
$\displaystyle=u(x,0)\sqrt{1+\alpha^{2}b^{{}^{\prime}2}}$
$\displaystyle=\sqrt{1+\alpha^{2}b^{{}^{\prime}2}}/[1-\Delta b(x)],\qquad
x=x(\xi).$ (37)

The shallow flow solution (36) provides precisely the outer flow required for
the numerical solution of the boundary layer equations (31) as formulated in
Appendix A.

. (a) (b) .

Figure 2: The shear stress (red) in the $x$-direction at the wall and the
boundary-layer displacement thickness (blue) as functions of position $x$ for
a Gaussian ridge of maximum fractional depth $\Delta=0.5$. (a)
$\text{Ro}=0.5.$ (b) $\text{Ro}=1.1.$ Figure 3: Profiles of $U$, the
$x$-component of velocity, as a function of the scaled normal coordinate,
$\zeta$ for a Gaussian ridge of maximum fractional depth $\Delta=0.5$ and
$\text{Ro}=1.1$, as in figure 2(b). The red curve gives the profile at
$x=0.63$ corresponding to the minimum in the wall shear and the blue curve
gives the profile at $x=4.23$ corresponding to the minimum in the displacement
thickness.

Two quantities are of particular interest: the scaled boundary-layer
displacement thickness

$d(x)=\int_{0}^{\infty}[1-U(x,\zeta)/U_{e}(\xi(x))]\;\text{d}\zeta,$ (38)

and the the magnitude of the shear stress at the wall where the shear stress
is given by

$\bm{\tau}(x)=U_{\zeta}(x(\xi,0))\bm{\hat{\xi}}+V_{\zeta}(x(\xi,0))\bm{\hat{y}}.$
(39)

Figure 2 shows the first component of ${\tau}(x)$ and $d(x)$ for a Gaussian
ridge $b(x)=\exp(-x^{2}/2)$ of fractional height $\Delta=0.5$. Negative
$\tau_{1}(x)$ corresponds to recirculation in the planes $y=$ constant. In
figure 2(a) rotation is reasonably strong ($\text{Ro}=0.5$). From (12) the
flow direction thus rotates by $\tan^{-1}[(\Delta/\text{Ro})\int
b(x)\text{d}x]=\tan^{-1}\sqrt{2/\pi}$ clockwise on passing over the ridge and
the component of the velocity in the $x$-direction increases, as does the
shear stress at the wall as in (8). The adverse pressure gradient in the lee
of the ridge causes a small local decrease in wall stress. The rotation of the
external velocity also leads to an internal maximum of the downstream velocity
which exceeds the free-stream velocity (as in figure 3 for $\text{Ro}=1.1$).
This increases the volume flux in the boundary layer compared with that in the
free stream and leads to a negative displacement thickness downstream of the
ridge. The adjustment to the new displacement thickness is through a viscously
damped oscillation with wavelength determined by a balance between the
advection and Coriolis terms in (31a), (31b) and so is of order Ro, as can be
seen by comparing figures 2(a) with 2(b) (for $\text{Ro}=1.1$) where the
wavelength is approximately 2.2 times longer. In figure 2(b) nonlinear
advection is more important making the adverse pressure gradient in the lee of
the ridge larger and causing the wall shear stress to fall to a minimum of
approximately 0.04. For larger Ro the wall shear stress becomes negative
corresponding to a region of reversed flow near the wall. The thickening of
the boundary layer near the point of minimum wall shear, even for
unidirectional flow, is visible in the rapid increase of the displacement
thickness there. Figure 3 shows profiles of $U$ as a function of the scaled
normal coordinate, $\zeta$ for the solution of figure 2(b). The red curve
gives the profile at $x=0.63$ corresponding to the minimum in the wall shear
and so is almost vertical at the wall. The blue curve gives the profile at
$x=4.23$ corresponding to the minimum in the displacement thickness, showing
the jet-like flow in the interior of the layer which increases the mass flux
above the corresponding mean flow flux, thinning the layer and giving the
negative displacement thickness.

##  4 Numerical simulation of the full equations for shallow flow

We present here numerical simulations of the full nonlinear equations at high
$Re$ and small channel aspect ratio $\delta$, comparing them to the boundary
layer integrations (for the limit $Re\to\infty,\delta\to 0$) in Section 3. The
geometry is fixed independently of $\delta$ by the vertical scaling
$\hat{z}=z^{*}/D=z/\delta$. The lower boundary is taken as the Gaussian ridge
$\hat{z}=\Delta\exp(-x^{2}/2)$, with $\Delta=0.5$, as in Figs. 2 and 3. For
the computations, we cut the domain at $x\pm 10$ and define

$\displaystyle\Omega^{s}:=\\{(x,z)\in\mathbb{R}^{2},\,-10<x<10,\,b(x)<\hat{z}<1\\}.$

We introduce a sub-tesselation ${\cal T}_{h}$ of $\Omega^{s}$ into triangles
$T\in{\cal T}_{h}$ and use the following finite element spaces of continuous
piecewise linear finite elements on ${\cal T}_{h}$

$\displaystyle{\cal V}_{h}$ $\displaystyle:=\\{v\in
C(\bar{\Omega}^{s})^{3},\,v_{i}|_{T}\in P_{1}(T)\,\forall T\in{\cal
T}_{h},\,\,v_{i}=0\text{ on }\Gamma_{i}^{d},\,i=1,2,3\\},$ $\displaystyle{\cal
L}_{h}$ $\displaystyle:=\\{p\in C(\bar{\Omega}^{s}),\,p|_{T}\in
P_{1}(T)\,\forall T\in{\cal T}_{h}\\}.$

The Dirichlet part of the boundary $\Gamma_{i}^{d}$ differs slightly between
the components. For $i=1,2$ it is given by the lower boundary, where no-slip
conditions (9) are imposed and the left (upstream) boundary, where the inflow
profile (34) is imposed for $U=1$. $\Gamma_{3}^{d}$ includes additionally the
upper (symmetry) boundary, where the symmetry condition (10) is imposed. On
the right downstream boundary the do-nothing conditions

$\displaystyle Re^{-1}\partial_{n}u-pn=0,\quad\partial_{n}v=\partial_{n}w=0$

are used, where $n$ denotes the exterior normal vector on $\partial\Omega$.

In order to abbreviate notation, we write
$U:=(u,v,w),\bar{U}:=(u,w),\bar{\nabla}:=(\partial_{x},\partial_{\hat{z}})$
and $\bar{\nabla}_{\delta}:=(\partial_{x},\delta^{-1}\partial_{\hat{z}})$. By
$U^{d}=(u^{d},v^{d},w^{d})$ we denote an extension of the Dirichlet data to
$\Omega^{s}$, i.e. the inflow profile (6) for $u$ and $v$ on the upstream
boundary and the zero function for $w$ on the remaining Dirichlet parts.

The discrete variational formulation corresponding to (1a)-(1c) reads: Find
$U\in U^{d}+{\cal V}_{h},p\in{\cal L}_{h}$, such that

$\displaystyle\begin{aligned}
(\bar{U}&\cdot\bar{\nabla}U,\phi)_{\Omega^{s}}+Re^{-1}(\bar{\nabla}_{\delta}U,\bar{\nabla}_{\delta}\phi)_{\Omega^{s}}-Ro^{-1}(v,\phi_{1})_{\Omega^{s}}\\\
&+Ro^{-1}(u-1,\phi_{2})_{\Omega^{s}}-(p,\bar{\nabla}\cdot\phi)_{\Omega^{s}}+(\bar{\nabla}\cdot\bar{U},\eta)_{\Omega^{s}}\\\
&\qquad+s_{u}(U,\phi)+s_{p}(p,\eta)=0\qquad\forall\phi\in{\cal
V}_{h},\eta\in{\cal L}_{h}.\end{aligned}$ (40)

Due the lack of inf-sup-stability of the discrete spaces and to handle
convection-dominated regimes for high $Re$, we add the weakly consistent
Continuous Interior Penalty stabilisation terms defined by

$\displaystyle s_{u}(U,\phi)$ $\displaystyle=\gamma_{u}\sum_{e\in{\cal
E}_{h}}|e|^{2}(\llbracket\nabla
U\rrbracket,\llbracket\nabla\phi\rrbracket)_{e},$ $\displaystyle
s_{p}(p,\eta)$ $\displaystyle=\gamma_{p}\sum_{e\in{\cal
E}_{h}}\max\big{\\{}|e|^{3}\frac{\delta}{\nu},\frac{|e|^{2}}{\|\bar{U}\|}\big{\\}}(\llbracket\bar{\nabla}p\rrbracket,\llbracket\bar{\nabla}\eta\rrbracket)_{e},$

where $\llbracket\cdot\rrbracket$ denotes the jump operator over interior
edges $e\in{\cal E}_{h}$ and $\gamma_{u},\gamma_{p}$ are positive parameters
Burman  _et al._ (2006).

All the simulation results have been obtained with the finite element library
FreeFem++ Hecht (2013). For the visualisation, we use the visualisation tool
ParaView Ayachit (2015). We use a triangulation with $192\,682$ vertices
($379\,024$ triangles), which is highly refined along the lower boundary, in
order to resolve the boundary layer for high $Re$. To solve the non-linear
system of equations, we use a fixed-point iteration based on the Oseen
linearisation of the Navier-Stokes equations: For $l=1,\dots$ find $U^{l}\in
U^{d}+{\cal V}_{h},p^{l}\in{\cal L}_{h}$, such that

$\displaystyle\big{(}\bar{\beta}^{l}$
$\displaystyle\cdot\bar{\nabla}U^{l},\phi\big{)}_{\Omega^{s}}+Re^{-1}(\bar{\nabla}_{\delta}U^{l},\bar{\nabla}_{\delta}\phi)_{\Omega^{s}}-Ro^{-1}(v^{l},\phi_{1})_{\Omega^{s}}$
$\displaystyle+Ro^{-1}(u^{l}-1,\phi_{2})_{\Omega^{s}}-(p^{l},\bar{\nabla}\cdot\phi)_{\Omega^{s}}+(\bar{\nabla}\cdot\bar{U}^{l},\eta)_{\Omega^{s}}$
$\displaystyle\qquad+s_{u}(U^{l},\phi)+s_{p}(p^{l},\eta)=0\qquad\forall\phi\in{\cal
V}_{h},\eta\in{\cal L}_{h},$

where $\beta^{l}:=\frac{1}{2}(U^{l-1}+U^{l-2})$ for $l\geq 2$ and
$\beta_{1}=U^{0}$ for a suitable starting value $U^{0}$. Although its rate of
convergence is in general worse compared to a standard Newton iteration, we
expect this fixed-point iteration to be more robust, and for this reason
better suitable to give an indication on the existence or non-existence of
stationary solutions.

###  4.1 The attached boundary-layer regime, $\text{Ro}<1.2$

Figure 4 compares the non-dimensional wall shear stress,

$\tau=\partial_{n}u/(E^{1/2}\delta),$

for $\delta=0.1$ and $\text{Ro}\in\\{0.5,\,1.1\\}$ at increasing Reynolds
numbers $\text{Re}=4000,\,16000,\,64000\text{ and }256000$ with the boundary
layer solutions of Figure 2, corresponding to the limit $\delta\to 0$,
$\text{Re}\to\infty$.

Figure 4:  The wall shear stress for $\delta=0.1$ and $Ro=0.5$ (left) and
$Ro=1.1$ (right) for different $Re$, compared to the boundary layer
calculations for $\delta\to 0,Re\to\infty$.

We observe convergence for $Re\to\infty$, although high Reynolds numbers are
needed to get close to the results obtained from the boundary layer
computations. Notice that we do not expect to fully reach the latter curve, as
$\delta$ is fixed to 0.1. The solution stays attached to the boundary for all
tested Reynolds numbers, i.e. up to $Re=1.024\cdot 10^{6}$. The finite element
method appears to give reasonable results up to very large Reynolds numbers.

###  4.2 Separated, steady, recirculating flow, $Ro\geq 1.2$

For $Ro\geq 1.2$ (31a)-(31b) can not be used to compute the inner boundary-
layer solution, as the solution detaches in the lee of the ridge (for
$\delta\to 0,Re\to\infty$) and the system is no longer parabolic. The finite-
element solution has no such restrictions. Figure 5 shows numerical results
for fixed $\delta=0.1$ and $Ro=1.2$ and $Ro=1.5$, at various values of Re. For
$Ro=1.2$ the boundary layer detaches at $Re=1.024\cdot 10^{6}$, i.e. the wall
shear stress takes negative values in the lee of the ridge, at $x\approx 0.8$.
For $Ro=1.5$, this happens earlier, by $Re\geq 64\,000$. For both Rossby
numbers, the method gives stationary states up to (at least) $Re=1.024\cdot
10^{6}$.

Figure 5:  The wall shear stress for $\delta=0.1$ and $Ro=1.2$ (left) and
$Ro=1.5$ (right) for various $Re$.

Figure 6 (top) shows isolines in the $x-z$ plane of the streamfunction for
$Re=256\,000$ and $Ro=1.2$, $Ro=1.5$ and $Ro=2$. Since the $y$-component of
velocity is non-zero these are vertical sections of stream surfaces, referred
to subsequently, for brevity, as streamlines. While for $Ro=1.2$ the flow is
still attached to the boundary, the flow detaches for $Ro=1.5$ and $Ro=2$ and
a separation bubble becomes clearly visible. Figure 6 (bottom) shows vertical
profiles through the eddy of the horizontal velocity $u$. The negative values
show the reverse flow for $Ro=1.5$ and $Ro=2$.

Figure 6:  Left: Streamlines in the lee of the ridge at $Re=256000$ for $Ro=$
1.2, 1.5 and 2. The flow detaches for $Ro=1.5$ and $Ro=2$ at $x\approx 0.75$
and $x\approx 0.6$, respectively. Right: Profiles normal to the lower boundary
through the center of the eddy of the horizontal velocity $u$ for these cases.

Figure 7 shows the critical Reynolds numbers (red solid curve) for boundary
layer separation as a function of Rossby number. Separation occurs at all
points above the curve while at points below the curve the flow remains
attached. For the ridge geometry here, the flow remains attached for Reynolds
numbers, $Re\leq 5\,600$, at all rotation rates, with separation depending
solely on the Reynolds number for $\,Ro\geq 6$. At smaller Rossby numbers,
$Ro<1.2$, when rotational effects become significant, the flow remained
attached to the boundary for all tested configurations, i.e. up to at least
$Re=1.024\cdot 10^{6}$. With decreasing $Ro$ the curve approaches the vertical
dashed line showing the upper bound at $Ro=1.1$ for attached flow in the
boundary layer analysis of §3. At any fixed $Ro>1.1$ there exists a critical
Reynolds number, higher than the Reynolds number for separation, above which
the fixed-point algorithm is unable to find a stationary solution. This gives
an indication that above some value of $Re$ no stable stationary solution
exists. The values for the disappearance of stationary solutions are shown by
the blue dotted curve in Figure 7 which shows the same qualitative behaviour
as the separation Reynolds number curve. For Rossby numbers greater than 1.2
there is thus a considerable range of Reynolds numbers where a steady
recirculating separation bubble is present in the lee of the ridge.

Figure 7: Shallow flow, $\delta=0.1$. Red solid curve: The smallest value of
$Re$, at fixed $Ro$, for which the flow separates. For all pairs $(Ro,Re)$
above the line a separation bubble occurs in the lee of the ridge. Blue dotted
curve: The largest value of $Re$, at fixed $Ro$, for which a stationary
solution was obtained by the fixed point algorithm. The vertical dashed line
shows the upper bound at $Ro=1.1$ for attached flow in the boundary layer
analysis of §3.

##  5 Deep flow over a ridge

The shallow domain of §4 suppresses inertial waves. To allow inertial waves
and assess the effect of the absence of an upper rigid boundary we cut the
computational domain at $z=100$ and impose the boundary conditions
$(u,v,w)=(1,0,0)$ at the left inflow boundary and on the top boundary. The
ridge profile is taken as $z=\alpha\exp(-x^{2}/2)$, with $\alpha=1/2$, giving
ridge slopes that are 10 times those of §4. This is reflected below in
separation at much smaller values of the Reynolds number. The variational
formulation is given by (40) for $\delta=1$, where the space ${\cal V}_{h}$ is
modified to accommodate for the Dirichlet boundary conditions on the upper
boundary for all variables.

Figure 8 shows the vertical velocity component $w$ for $Ro=0.1$ at $Re=400$
and $10\,000$. The flow patterns are very similar. This accords with the
analysis in Johnson (1982); Machicoane  _et al._ (2018) showing that the
leading order viscous effect in deep flow at small Rossby number is the
introduction of a downstream decay of the lee wave wake with a scale of order
$Re^{-1}$ and a downstream displacement of the phase lines at large $Re$ of
order $(RoRe)^{-1}$.

Figure 8:  The vertical velocity component $w$ for deep flow with Ro=0.1 at
Re=400 (left) and Re=10 000 (right). The flow field is almost unaltered when
the Reynolds changes by a factor of 25. The mesh is highly refined around the
ridge and the figure shows approximately half the simulation domain.

Figure 9, as for figure 7, shows at fixed Rossby number the largest values of
Reynolds number which the flow remained attached (red curve) and at which the
fixed point iteration converged (blue curve). The qualitative behaviour of the
curves is very similar but the critical Reynolds numbers are an order of
magnitude smaller in the deep flow. As noted above, this is a consequence of
the ridge slope’s being an order of magnitude larger here. At fixed Rossby
number, greater than 1.5, there again exists a considerable range of Reynolds
numbers where a steady recirculating separation bubble is present in the lee
of the ridge.

Figure 9:  As figure 7 but for deep flow. Red solid curve: The smallest value
of $Re$, at fixed $Ro$, for which the flow separates. Blue dotted curve: The
largest value of $Re$, at fixed $Ro$, for which a stationary solution was
obtained by the fixed point algorithm.

Figure 10 shows the vertical velocity $w$ and streamlines in the lee of the
ridge for deep flow and weak rotation with $Ro=3$ at $Re=200,1000$. The flow
is attached for $Re=200$ but has separated by $Re=1000$ with a steady
recirculating separation bubble visible.

\begin{overpic}[width=195.12767pt]{w_deep_Ro3_Re200_v2.png}\hfil\put(0.0,5.0){(a)}
\end{overpic}

\begin{overpic}[width=195.12767pt]{w_deep_Ro3_Re1000_v2.png}
\put(0.0,5.0){(b)} \end{overpic}

Figure 10:  The vertical velocity $w$ and streamlines in the lee of the ridge
for deep flow and weak rotation, $\text{Ro}=3$. (a) $\text{Re}=200$, attached
flow, (b): $\text{Re}=1000$, separated flow with a stable separation bubble.

##  6 Flow over a horizontal cylinder

To capture more closely the geometry of the experiments reported in Machicoane
_et al._ (2018), this section considers flow over a horizontal right circular
cylinder. The variational formulation is again given by (40) with $\delta=1$.
We use an ellipsoidal simulation domain with semi-axes of size 60 (vertically)
and 150 (horizontally) from which a circle of radius $0.5$ (representing an
infinitely long cylinder in transversal direction), is extracted. We set
$U^{d}=(1,0,0)$ on the outer boundary of the ellipsoid and $U^{d}=(0,0,0)$ on
the the cylinder.

Figure 11 shows a pair of streamlines for each of (a,b) $Ro=0.5$ and (c,d)
$Ro=2$ at Reynolds numbers $Re$ chosen so the flow is attached to the cylinder
at the lower $Re$ (a,c) and detached at the higher (b,d). The lower $Ro$ flow
remains attached at significantly higher $Re$ but in both cases the steep
obstacle gradient at the rear of the cylinder, and the concomitant strong
adverse pressure gradient, induces separation at far lower Reynolds numbers
than in flow over the smooth ridge. The inertial wave wake alternatively
compresses and accelerates and then expands and decelerates a jet-like flow
along the axis in the lee of the cylinder. Figure 12 shows vertical profiles
through the jet of the horizontal velocity $u$ at downstream stations
$x=1,2,3$ and $4$ for the flow in Fig. 11(b). The jet is narrow with almost
twice the free-stream speed at $x=2$ but wide with a velocity deficit giving a
speed of two-thirds the free stream speed at $x=4$.

\begin{overpic}[width=173.44534pt]{U_Ro05_Re200_v2.png} \put(0.0,1.0){(a)}
\end{overpic}

\begin{overpic}[width=173.44534pt]{U_Ro05_Re600_v2.png} \put(0.0,1.0){(b)}
\end{overpic}

\begin{overpic}[width=173.44534pt]{U_Ro2_Re15_v2.png} \put(0.0,1.0){(c)}
\end{overpic}

\begin{overpic}[width=173.44534pt]{U_Ro2_Re50_v2.png} \put(0.0,1.0){(d)}
\end{overpic}

Figure 11:  Streamlines for $(\text{Ro},\text{Re})=$ (a) (0.5,200), (b)
(0.5,600), (c) (2,15) and (d) (2,50). The colouring shows the magnitude of the
velocity $(u,w)$. The figure shows only a small part of the simulation domain.
\begin{overpic}[width=173.44534pt]{profiles_jet.png} \put(-3.0,50.0){u}
\put(50.0,-2.0){y} \end{overpic} Figure 12:  Vertical profiles through the lee
jet of the horizontal velocity $u$ at downstream stations $x=1,2,3$ and $4$
for the flow in Fig. 11(b) where $\text{Ro}=0.5$ and $\text{Re}=600$.

As figures 7 and 9, figure 13 shows for fixed Rossby number the lowest $Re$
for which separation occurs and the highest $Re$ for which the fixed point
iteration converged. The curves behave qualitatively very similarly to those
of Figures 7 and 9: as $Ro\to 0$ the critical Reynolds number for convergence,
$Re_{\text{crit}}$ increases rapidly, showing the existence of steady
solutions at high $Re$, and for weakly rotating flow, $Ro\to\infty$,
$Re_{\text{crit}}$ asymptotes to the constant value ($Re_{\text{crit}}\approx
160$ here) of non-rotating flow. For Reynolds numbers of order 200 and larger
the range of values of Rossby numbers for which steady separated solutions
exist is small and becomes smaller with increasing Reynolds number, almost
vanishing by $Re=1000$. These Reynolds numbers include the parameter regime
investigated in Machicoane  _et al._ (2018) and supports their observation
that separation and unsteadiness appear simultaneously. Figure 14, shows the
$Re_{\text{crit}}$ curve over a smaller range of Reynolds numbers with
superposed the experimental observations of Machicoane  _et al._ (2018)
(Figure 7) that are closest to the transition between steady and unsteady
solutions. The qualitative form of the numerical curve appears to capture the
experimentally observed separation between steady and unsteady flow. Two
experimental points for unsteady flow lie just below the numerical curve. The
deviation is within the range of experimental fluctuations and numerical
error. In the numerical solution, numerical diffusion could lead to
convergence to stationary states on the discrete level that are not stationary
states on the continuous level so the critical curve in the continuous limit
would lie slightly below the plotted curve. In the experiments even small
disturbances to the flow might be of sufficient magnitude to trigger
instability below the theoretical threshold. Even with this deviation the
numerical curve seems to capture the experimental results well over the whole
range. Based on their analysis of the lee-wave wake at large Reynolds number
Machicoane  _et al._ (2018) postulate that the division between steady and
unsteady flow is of the form $Ro\cdot Re_{\text{crit}}=C$ with the
experimental results giving $C=275$. This curve is included in Figure 14 and
captures the division in the flow regime there well. The table in Figure 13
shows that the relation may not capture the division so well in other flow
regimes.

Ro | 0.3 | 0.5 | 0.75 | 1 | 2 | 3 | 5 | 10
---|---|---|---|---|---|---|---|---
$\text{Re}_{\text{crit}}$ | 1640 | 1050 | 500 | 295 | 260 | 220 | 180 | 170
$\text{Ro}\cdot\text{Re}_{\text{crit}}$ | 492 | 525 | 375 | 295 | 520 | 660 | 900 | 1700

Figure 13:  Top: As figures 7, 9, but for deep flow over a cylinder. Red solid
curve: The smallest value of $Re$, at fixed $Ro$, for which the flow
separates. Blue dotted curve: The largest value of $Re$, at fixed $Ro$, for
which a stationary solution was obtained by the fixed point algorithm. Bottom:
The largest Reynolds number $Re_{\text{crit}}$, for which a stationary
solution was obtained by the fixed point algorithm for fixed $Ro$ and product
$Ro\cdot Re_{\text{crit}}$.  Figure 14:  Comparison of the critical Reynolds
numbers found by numerical simulations (blue curve) against the experimental
results in Machicoane  _et al._ (2018). Green circles indicate combinations of
$Ro$ and $Re$, for which steady solutions were observed experimentally; black
squares indicate unsteady solutions.

##  7 Conclusion

Motivated by the geophysical relevance of rotating flow over topography in the
atmosphere Mason and Sykes (1979) and ocean and by recent theoretical results
and laboratory experiments Machicoane  _et al._ (2018) we have discussed
steady attached flow, boundary layer separation, transition to unsteady flow
and flow patterns in the three numerical examples of shallow and deep flow
over a Gaussian ridge and deep flow over a horizontal cylinder using a
stabilised equal-order finite element method with Continuous Interior Penalty
Stabilisation. The advantage of this technique is that it can obtained steady
solutions at high Reynolds numbers. For shallow flow (with aspect ratio
$\delta=D/L\to 0$) the solutions were compared with integration of the
boundary layer equations valid in the limit $Re\to\infty$. The finite-element
results for $\delta=0.1$ approached the boundary layer results with increasing
Reynolds number and in the region of attached flow the finite-element results
converged for Reynolds numbers of order $10^{6}$. In deep flow over the ridge
the inertial wave wake discussed in Johnson (1982); Machicoane  _et al._
(2018) appeared. The computations for deep flow over a cylinder included the
parameter regime of the experiments of Machicoane  _et al._ Machicoane  _et
al._ (2018) and reproduced many of the features observed there.

In all three cases the flow was steady and attached for sufficiently small
Reynolds number. As the Reynolds number was increased at fixed Rossby number
the flow separated in the lee of the ridge and a steady separation bubble
appeared. At higher Reynolds numbers the fixed point iteration failed to
converge. In flow over a cylinder this coincided closely with the appearance
of unsteady flow in the experiments of Machicoane  _et al._ (2018). The
Reynolds numbers above which separation and, subsequently, unsteadiness appear
decrease by an order of magnitude in moving from shallow to deep flow, as the
ridge slope increases by a factor of 10, and by another order of magnitude in
moving from the smooth ridge to the cylinder, with its infinite slope at the
rear stagnation point. In all cases it appears that sufficiently strong
rotation allows the flow to remain attached at arbitrarily large Reynolds
numbers.

##  Appendix A Boundary-layer integration

The method here follows Prandtl (1938), modified to allow expansion in
Laguerre polynomials. In (31) write $(U,V)=(U_{e}+\tilde{u},V_{e}+\tilde{v})$
so $(\tilde{u},\tilde{v})\to 0$ as $\zeta\to\infty$. Then the $\xi$-momentum
equation (31a) becomes

$U\tilde{u}_{\xi}+W\tilde{u}_{\zeta}=\tilde{u}_{\zeta\zeta}+a\tilde{v}-U_{e\xi}\tilde{u},$
(41)

where the far-field pressure satisfies

$P_{e\xi}=-U_{e}U_{e\xi}+aV_{e}.$ (42)

Substituting (41) into the continuity equation and rearranging gives

$(W/U)_{\zeta}=-(U_{e}U_{e\xi}+a\tilde{v}+\tilde{u}_{\zeta\zeta})/U^{2}.$ (43)

On the solid boundary, $\zeta=0$, $U=0$ for all $\xi$, so $U_{\xi}=0$ and
hence $W_{\zeta}=0$. Thus

$\lim_{\zeta\to 0}(W/U)=\lim_{\zeta\to 0}(W_{\zeta}/U_{\zeta})=0,$ (44)

and (43) can be integrated from zero on the wall to give $W/U$ at fixed $\xi$.
The momentum equations (31a) and (31b) then become

$\displaystyle\tilde{u}_{\xi}$
$\displaystyle=-(W/U)\tilde{u}_{\zeta}+(\tilde{u}_{\zeta\zeta}+a\tilde{v}-U_{e\xi}\tilde{u})/U,$
(45a) $\displaystyle\tilde{v}_{\xi}$
$\displaystyle=-(W/U)\tilde{v}_{\zeta}+(\tilde{v}_{\zeta\zeta}-a\tilde{u}-V_{e\xi}\tilde{u})/U.$
(45b)

It is convenient to obtain the outer flow and to discuss the solutions in the
unstretched coordinate $x$. Multiplying (45) by
$\xi_{x}=\sqrt{1+\alpha^{2}b^{{}^{\prime}2}(x^{\prime})}$ gives

$\displaystyle\tilde{u}_{x}$
$\displaystyle=-\xi_{x}[(W/U)\tilde{u}_{\zeta}+(\tilde{u}_{\zeta\zeta}+a\tilde{v})/U]-U_{ex}\tilde{u}/U,$
(46a) $\displaystyle\tilde{v}_{x}$
$\displaystyle=-\xi_{x}[(W/U)\tilde{v}_{\zeta}+(\tilde{v}_{\zeta\zeta}-a\tilde{u})/U]-V_{ex}\tilde{u}/U.$
(46b)

Equations (46) and (43) are solved spectrally in $\zeta$ by expressing
$\tilde{u}$ and $\tilde{v}$ as the product of a decaying exponential and a
series of Laguerre polynomials through their values at the Laguerre
interpolation points, expressing (43) as a first order differential equation
using a differentiation matrix, inverting to obtain $W/U$ at the interpolation
points, and advancing in $x$ using (46) and a stiff ordinary differential
equation solver.

## References

  * Mason and Sykes [1979] P. J. Mason and R. I. Sykes, Separation effects in Ekman layer flow over ridges, Quart. J. Roy. Met. Soc. 105, 129 (1979).
  * Richards  _et al._ [1992] K. J. Richards, D. A. Smeed, E. J. Hopfinger, and G. C. D’Hiéres, Boundary-layer separation of rotating flows past surface-mounted obstacles, J. Fluid Mech. 237, 343 (1992).
  * Machicoane  _et al._ [2018] N. Machicoane, V. Labarre, B. Voisin, F. Moisy, and P.-P. Cortet, Wake of inertial waves of a horizontal cylinder in horizontal translation, Phys Rev Fluids 3 (2018).
  * Johnson [1982] E. Johnson, The effects of obstacle shape and viscosity in deep rotating flow over finite-height topography, J Fluid Mech 120, 359 (1982).
  * Williamson [1996] C. H. Williamson, Vortex dynamics in the cylinder wake, Ann. Rev. Fluid Mech. 28, 477 (1996).
  * Douglas and Dupont [1976] J. Douglas and T. Dupont, Interior penalty procedures for elliptic and parabolic Galerkin methods, in  _Comput Methods Appl Sci_ , edited by R. Glowinski and J. L. Lions (Springer, 1976) pp. 207–216.
  * Burman and Hansbo [2006] E. Burman and P. Hansbo, Edge stabilization for the generalized stokes problem: a continuous interior penalty method, Comput Methods Appl Mech Engrg 195, 2393 (2006).
  * Burman  _et al._ [2006] E. Burman, M. Fernández, and P. Hansbo, Continuous interior penalty finite element method for Oseen’s equations, SIAM J Numer Anal 44, 1248 (2006).
  * Tong  _et al._ [2022] G. G. Tong, D. Kamensky, and J. A. Evans, Skeleton-stabilized divergence-conforming b-spline discretizations for incompressible flow problems of high Reynolds number, Comput & Fluids 248, 105667 (2022).
  * Moura  _et al._ [2022] R. C. Moura, A. Cassinelli, A. F. da Silva, E. Burman, and S. J. Sherwin, Gradient jump penalty stabilisation of spectral/hp element discretisation for under-resolved turbulence simulations, Comput Methods Appl Mech Engrg 388, 114200 (2022).
  * Fernández [2013] M. A. Fernández, Incremental displacement-correction schemes for incompressible fluid-structure interaction: stability and convergence analysis, Numer Math 123, 21 (2013).
  * Burman  _et al._ [2020] E. Burman, M. A. Fernández, and S. Frei, A Nitsche-based formulation for fluid-structure interactions with contact, ESAIM: M2AN 54, 531 (2020).
  * Jacobs [1964] S. Jacobs, On stratified flows over bottom topography, J Mar Res , 223 (1964).
  * Hecht [2013] F. Hecht, New development in FreeFem++, J Numer Math 20, 251 (2013).
  * Ayachit [2015] U. Ayachit,  _The ParaView Guide: A Parallel Visualization Application_ (Kitware, Inc., USA, 2015).