# Spontaneous critical layer formation and robustness beneath rotational waves

M. V. Flamarion1, A. Nachbin2 and R. Ribeiro-Junior3

###### Abstract

Non-stationary rotational surface waves are considered, where the underlying
current has constant vorticity. A study is presented on the robustness of a
critical layer in the presence of a bottom topography, as well as on its
spontaneous formation for waves generated from rest. The restriction, from
previous studies, to a traveling-wave formulation is removed leading to a non-
stationary set of equations. In this setting streamlines are not necessarily
pathlines. Particle-trajectories are found evolving the respective submarine
dynamical system with a cloud of tracers. Pathlines are then visualized and
the respective submarine structures identified. Robustness is illustrated
through surface waves interacting with topographic undulations. The respective
Kelvin cat eye structure dynamically adjusts itself over the bottom topography
without loosing its integrity. On the spontaneous formation of a Kelvin cat
eye structure, the surface is initially undisturbed and waves are generated
from either the current-topography interaction or by a surface pressure
distribution suddenly imposed. Under the pressure forcing, an isolated Kelvin
cat eye spontaneously forms. The two extreme critical points of the cat eye
structure are connected with a stagnation segment, thus exhibiting a new form
of a critical layer.

Keywords: Surface gravity waves. Waves/Free-surface Flows. Critical layers.

1 UFRPE/Rural Federal University of Pernambuco, UACSA/Unidade Acadêmica do
Cabo de Santo Agostinho, BR 101 Sul, 5225, 54503-900, Ponte dos Carvalhos,
Cabo de Santo Agostinho, Pernambuco, Brazil.

2 IMPA/National Institute of Pure and Applied Mathematics, Est. D. Castorina,
110, Rio de Janeiro, RJ, 22460-320, Brazil.

3 UFPR/Federal University of Paraná, Departamento de Matemática, Centro
Politécnico, Jardim das Américas, Caixa Postal 19081, Curitiba, PR, 81531-980,
Brazil.

##  1 Introduction and background

Water waves propagating in the presence of a background flow is a problem of
great current interest both regarding applications as well as the mathematical
questions that arise. This topic is broad and therefore difficult to make a
comprehensive overview of contributions. For the interested reader we mention
just a few recent references from which the bibliography can be useful.
Regarding geophysical applications, Soontiens et al. [2010] studied trapped
internal waves over isolated topography, in the presence of a background
shear. The book by Buhler [2009] presents several techniques relevant to GFD
(Geophysical Fluid Dynamics) where an example includes “the interplay between
large-scale Rossby waves and two-dimensional turbulence”. On the other extreme
of the scientific spectrum, we have many contributions in Analysis with
rigorous theoretical results on the existence of surface waves in the presence
of vorticity (namely rotational waves), on the existence of stagnation points
in the wave’s moving frame, the respective Kelvin cat eye recirculation region
below nonlinear (periodic, traveling) Stokes waves, among others. Many recent
references may be found through the work of Constantin and collaborators
(Constantin [2011, 2017]; Constantin et al. [2016]). Numerical studies with
traveling waves and the respective stationary submarine structures, such as
stagnation points, start with the work of Teles da Silva & Peregrine [1988].
More detailed numerical studies on the flow structure beneath traveling waves
appeared recently, as for example that of Vasan & Oliveras [2014] and also
Ribeiro Jr et al. [2017]. A numerical stability study for finite-amplitude
steady rotational surface waves is presented by Francius & Kharif [2017].

Regarding wave-current interaction, special attention has been given to
sheared currents of constant vorticity, a regime which permits the existence
of stagnation points in the wave’s moving frame. These stagnation points are
part of a recirculation region, displaying a Kelvin cat eye formation. Some
theoretical and computational studies, such as those mentioned above, have
considered different issues of interest related to particle trajectories
beneath a rotational traveling wave. Having traveling waves in mind, most
studies are formulated in the wave’s moving frame thus addressing stationary
differential equations. The stationary submarine structure for particle
dynamics, for example, can be visualized through level-curves of the
streamfunction.

In the present study, in order to study the robustness and the spontaneous
formation of the Kelvin cat eye structure and its associated critical layer,
we need to consider a non-stationary set of equations. Even if the initial
disturbance has data related to a traveling wave, the evolution equations are
not “aware” of this class of data and should (hopefully) evolve the initial
wave profile as a pure translation. Unfortunately in the nonlinear case there
is strong evidence that traveling-wave data evolves as a stable wave of
translation only for small amplitudes. Francius & Kharif [2017] studied the
stability of rotational waves and report an interesting result that
“increasing the shear reduces the growth rate of the most unstable sideband
instabilities but enhances the growth rate of these quartet instabilities,
which eventually dominate the Benjamin-Feir modulational instabilities.” To
the best of our knowledge there are no articles with time dependent models and
non-stationary waves studying the submarine Kelvin cat eye structure and its
associated critical layer. The critical layer defines the boundary where the
submarine flow (beneath the wave) is divided into a rightgoing and a leftgoing
stream.

As we remove the restriction of traveling waves, our propagating waves might
change their profiles as time evolves implying that the streamlines are no
longer pathlines. In this case, not only we have to compute the free surface
conditions in time, but also we have to solve the particle-trajectories’
dynamical system for a cloud of tracers in order to visualize the pathlines
and the respective submarine structures.

The robustness of the submarine structure is tested first by allowing an
initial surface disturbance to propagate over the bottom topography. Due to
the problem of nonlinear wave instability we perform our numerical study with
linear waves, which has not yet been explored. We consider a modulated initial
surface disturbance, where the Kelvin cat eye structure is already present at
time $t=0$. The initial wavetrain is chosen so that dispersive effects are
very small and we observe effectively (namely to a good approximation) a wave
of translation. We will show that as the wave interacts with the topographic
undulations, commensurable with its wavelength, the Kelvin cat eye dynamically
adjusts and “slides” over the topography without loosing its integrity. When
the finite-length topography ends, the recirculation region is back to its
original cat eye formation. In the supplemental material a video shows the
entire process.

A rapidly varying topography is also considered. In this case we cannot
identify the Kelvin cat eye over the topography. The material curves, defined
through groups of tracers, end up crossing each other. Nevertheless, it is
remarkable to observe the robust reconfiguration of the submarine structure,
identified by the respective group of tracers, as soon as the bottom becomes
flat again.

The novel results on the spontaneous formation of the Kelvin cat eye structure
are considered through surface waves generated from rest. The free surface is
initially undisturbed when surface waves start being generated by the current-
topography interaction, or by a localized (steady) surface pressure
distribution suddenly applied at time $t=0^{+}$. This case is motivated by the
work of Johnson [2012]. Under the pressure forcing two localized surface
pulses are generated and propagate in opposite directions. An isolated Kelvin
cat eye spontaneously forms under the downstream-propagating pulse. A center
(critical point) appears in the middle of the recirculation region. The two
extreme critical points of the cat eye structure are connected with a
stagnation segment, thus exhibiting a new form of a critical layer. All these
dynamical critical-layer scenarios described above have not been contemplated
in the literature.

The paper is organized as follows. In section 2 we present the mathematical
formulation of the free surface Euler equations in the canonical domain, a
uniform strip where computations are more easily performed. The canonical
domain is defined through a conformal mapping. In section 3 the numerical
method is presented. We introduce the dynamical system for particle
trajectories in canonical coordinates. By not using a traveling-wave
formulation this dynamical system is no longer autonomous and its vector field
must be constantly updated. This update depends on solutions of the Euler
equations and is done through the potential component of the velocity field.
Properties of harmonic functions are used in order to write all “Euler-
information” needed in terms of one-dimensional Fourier expressions. These are
essentially Fourier-type operators acting on the Dirichlet (boundary) data.
This framework leads to the numerical method described in section 3. The
results are presented in section 4 and the conclusions in section 5.

##  2 Mathematical formulation

We have a two-dimensional incompressible flow of an inviscid fluid of constant
density $\rho$. The free-surface Euler equations (Whitham [1974]) is our
reference model:

$\displaystyle\begin{split}&u_{t}+uu_{x}+vu_{y}=-\frac{p_{x}}{\rho},\;\
\mbox{for}\;\ h(x+U_{0}t),<y<\tilde{\zeta}(x,t),\\\
&v_{t}+uv_{x}+vv_{y}=-\frac{p_{y}}{\rho}-g,\;\ \mbox{for}\;\
h(x+U_{0}t),<y<\tilde{\zeta}(x,t),\\\ &u_{x}+v_{y}=0,\;\ \mbox{for}\;\
h(x+U_{0}t)<y<\tilde{\zeta}(x,t),\\\ &p=P(x+U_{0}t),\;\ \mbox{at}\;\
y=\tilde{\zeta}(x,t),\\\ &v=\tilde{\zeta}_{t}+u\tilde{\zeta}_{x},\;\
\mbox{at}\;\ y=\tilde{\zeta}(x,t),\\\ &v=U_{0}h_{x}+uh_{x},\;\ \mbox{at}\;\
y=h(x+U_{0}t).\\\ \end{split}$ (1)

The free surface is denoted by $\tilde{\zeta}(x,t)$, the bottom topography by
$h$ and the forcing pressure distribution over the free surface by $P$. We
consider either a submarine obstacle traveling with speed $U_{0}$ along the
bottom or a pressure distribution traveling along the free surface, also with
speed $U_{0}$. The fluid velocity is denoted by $(u,v)$, the gravity
acceleration by $g$, while the pressure within the fluid body by $p$. All
functions $(u,v,\tilde{\zeta}$, h and $p)$ in equations (1) are considered to
be smooth and periodic in $x$ with period $L$.

We are in a flow regime having a given background current of constant shear.
With this in mind, write the velocity field in the form

$\displaystyle\begin{split}(u,v)=\nabla\tilde{\phi}+(ay,0),\end{split}$ (2)

where $\tilde{\phi}$ is the velocity potential of the irrotational component
of the flow while $-a$ is the given constant vorticity value. Substitute (2)
in (1) to obtain

$\displaystyle\begin{split}&\Delta\tilde{\phi}=0,\;\ \mbox{for}\;\
-h_{0}+h(x+U_{0}t)<y<\tilde{\zeta}(x,t),\\\
&(U_{0}-ah_{0})h_{x}+ahh_{x}+\tilde{\phi}_{x}h_{x}=\tilde{\phi}_{y},\;\
\mbox{at}\;\ y=-h_{0}+h(x+U_{0}t),\\\
&\tilde{\zeta}_{t}+(a\tilde{\zeta}+\tilde{\phi}_{x})\tilde{\zeta}_{x}-\tilde{\phi}_{y}=0,\;\
\mbox{at}\;\ y=\tilde{\zeta}(x,t),\\\
&\tilde{\phi}_{t}+\frac{1}{2}(\tilde{\phi}_{x}^{2}+\tilde{\phi}_{y}^{2})+a\tilde{\zeta}\tilde{\phi}_{x}+\tilde{\zeta}-a\tilde{\psi}=-\frac{P(x+U_{0}t)}{\rho},\;\
\mbox{at}\;\ y=\tilde{\zeta}(x,t),\end{split}$

where $\tilde{\psi}$ is the harmonic conjugate of $\tilde{\phi}$. In a moving
frame, given by $x\rightarrow x+U_{0}t$, denoting
$\overline{\zeta}(x,t)\equiv\tilde{\zeta}(x-U_{0}t,t),\;\
\overline{\phi}(x,y,t)\equiv\tilde{\phi}(x-U_{0}t,y,t),$ these equations read
as

$\displaystyle\begin{split}&\Delta\overline{\phi}=0\;\ \mbox{for}\;\
-h_{0}+h(x)<y<\overline{\zeta}(x,t),\\\
&(U_{0}-ah_{0})h_{x}+ahh_{x}+\overline{\phi}_{x}h_{x}=\overline{\phi}_{y}\;\
\mbox{at}\;\ y=-h_{0}+h(x),\\\
&\overline{\zeta}_{t}+(U_{0}+a\overline{\zeta}+\overline{\phi}_{x})\overline{\zeta}_{x}-\overline{\phi}_{y}=0\;\
\mbox{at}\;\ y=\overline{\zeta}(x,t),\\\
&\overline{\phi}_{t}+\frac{1}{2}(\overline{\phi}_{x}^{2}+\overline{\phi}_{y}^{2})+(U_{0}+a\overline{\zeta})\overline{\phi}_{x}+\overline{\zeta}-a\overline{\psi}=-\frac{P(x)}{\rho}\;\
\mbox{at}\;\ y=\overline{\zeta}(x,t).\end{split}$ (3)

We have thus prescribed the velocity field to have a constant vorticity,
through the given background flow, and to satisfy the Neumann condition (now)
around a fixed submarine obstacle. The surface pressure distribution and the
bottom topography are stationary in this reference frame.

We have the following characteristic scales: $h_{0}$ for distance,
$(gh_{0})^{1/2}$ as a reference long-wave speed, $(h_{0}/g)^{1/2}$ for time
and $\rho gh_{0}$ as a characteristic pressure. Using these characteristic
scales, and abusing of notation, we put (3) in a dimensionless form:

$\displaystyle\begin{split}&\Delta\overline{\phi}=0,\;\ \mbox{for}\;\
-1+h(x)<y<\overline{\zeta}(x,t),\\\ &(F-\Omega)h_{x}+\Omega
hh_{x}+\overline{\phi}_{x}h_{x}=\overline{\phi}_{y},\;\ \mbox{at}\;\
y=-1+h(x),\\\
&\overline{\zeta}_{t}+(F+\Omega\overline{\zeta}+\overline{\phi}_{x})\overline{\zeta}_{x}-\overline{\phi}_{y}=0,\;\
\mbox{at}\;\ y=\overline{\zeta}(x,t),\\\
&\overline{\phi}_{t}+\frac{1}{2}(\overline{\phi}_{x}^{2}+\overline{\phi}_{y}^{2})+(F+\Omega\overline{\zeta})\overline{\phi}_{x}+\overline{\zeta}-\Omega\overline{\psi}=-P(x),\;\
\mbox{at}\;\ y=\overline{\zeta}(x,t),\end{split}$ (4)

where the Froude number is $F={U_{0}}/{(gh_{0})^{1/2}}$ and
$\Omega={ah_{0}}/{(gh_{0})^{1/2}}$. All the variables in (4) are
dimensionless. The dimensionless vorticity is given by $-\Omega$. In the
linear regime equations (4) are of the form

$\displaystyle\begin{split}&\Delta\overline{\phi}=0,\;\ \mbox{for}\;\
-1+h(x)<y<0,\\\ &(F-\Omega)h_{x}+\Omega
hh_{x}+\overline{\phi}_{x}h_{x}=\overline{\phi}_{y},\;\ \mbox{at}\;\
y=-1+h(x),\\\
&\overline{\zeta}_{t}+F\overline{\phi}_{x}=\overline{\phi}_{y},\;\
\mbox{at}\;\ y=0,\\\
&\overline{\phi}_{t}+F\overline{\phi}_{x}+\overline{\zeta}-\Omega\overline{\psi}=-P(x),\;\
\mbox{at}\;\ y=0.\end{split}$ (5)

Since we will study linear waves, from now on each time we mention the Euler
equations we mean system (5).

Either for linear or nonlinear waves, the particle trajectory beneath a
surface wave is governed by the dynamical system

$\displaystyle\begin{split}&\frac{dx}{dt}=\overline{\phi}_{x}+\Omega y+F,\\\
&\frac{dy}{dt}=\overline{\phi}_{y},\\\ &x(0)=x_{0},\;\ y(0)=y_{0}.\end{split}$
(6)

To compute its vector field one needs $\overline{\phi}_{x}$ and
$\overline{\phi}_{y}$ in the bulk of the fluid. These are obtained from the
Euler equations, which in the present study is given by system (5). From the
dispersion relation of system (5) we have that the linear wave speed $c$ is
given by

$c=F-\frac{\Omega\tanh(k)}{2k}\pm\frac{\sqrt{\Omega^{2}\tanh^{2}(k)+4k\tanh(k)}}{2k},$
(7)

where we will consider the mode with the positive sign of the square root.
Suppose we have a traveling wave solution in the form
$\overline{\phi}(x,y,t)=\overline{\phi}(x-ct,y,t)$. If $(x(t),y(t))$
represents a particle trajectory given by the dynamical system (6) then, in
the wave’s moving frame $X=x-ct\mbox{ and }Y=y,$ the trajectories
$(X(t),Y(t))$ satisfy the dynamical system in the form

$\begin{array}[]{l }\dfrac{dX}{dt}=\overline{\phi}_{X}(X,Y,t)+\Omega
Y+F-c=\overline{\psi}_{Y}(X,Y,t)+\Omega Y+F-c,\\\ \vskip 12.0pt plus 4.0pt
minus
4.0pt\cr\dfrac{dY}{dt}=\overline{\phi}_{Y}(X,Y,t)=-\overline{\psi}_{X}(X,Y,t).\\\
\vskip 12.0pt plus 4.0pt minus 4.0pt\cr\end{array}$ (8)

A simple calculation shows that we cannot have waves excessively long while
seeking for a Kelvin cat eye structure. For the cat eye structure to exist we
need the presence of stagnation points, in the wave’s moving frame. In the
linear regime we have small amplitude waves and therefore a weak velocity
potential. Thus to leading order the flow is dominated by the background
current, with $\overline{\phi}_{X},\overline{\phi}_{Y}\approx 0$, and a linear
dynamical system such as

$\begin{array}[]{l }\dfrac{dX}{dt}\approx\Omega Y+F-c,\\\ \vskip 12.0pt plus
4.0pt minus 4.0pt\cr\dfrac{dY}{dt}\approx 0.\\\ \vskip 12.0pt plus 4.0pt minus
4.0pt\cr\end{array}$ (9)

Denote by $(X^{\star},Y^{\star})$ a candidate position for a stagnation point.
This yields the critical point $\Omega Y^{\star}=c-F,$ which by substituting
the wave speed (7) gives

$Y^{\star}=-\frac{\tanh(k)}{2k}\pm\sqrt{\frac{\tanh^{2}(k)}{4k^{2}}+\frac{\tanh(k)}{\Omega^{2}k}}.$
(10)

In the long wave limit with $(k\rightarrow 0)$ we have that

$Y^{\star}=-\frac{1}{2}\pm\sqrt{\frac{1}{4}+\frac{1}{\Omega^{2}}}.$ (11)

No matter how strong is the vorticity of the background flow, we cannot find
$-1\leq Y^{\star}\leq 0$. In other words there are no stagnation points within
our fluid body. For linear waves the stagnation point will appear as a balance
between wavelength, total depth and the vorticity applied.

In the next section we detail the conformal mapping technique for solving the
Euler equations (5) and the dynamical systems (6) or (8) in the canonical
domain, where we have a flat strip.

###  2.1 Conformal mapping

As depicted in figure 1 the conformal mapping $z=f(w)$, from the canonical
domain (a flat strip) onto the physical domain, is defined by

$z(\xi,\eta)=x(\xi,\eta)+iy(\xi,\eta),$

where in the canonical $w$-plane $w=\xi+i\eta$. We have the following boundary
conditions:

$y(\xi,0)=0\;\ \mbox{ and }\;\ y(\xi,-D)=-1+H(\xi).$

We are imposing that the canonical upper boundary $\eta=0$ is mapped onto the
undisturbed free surface $y=0$. We denote by $H(\xi)=h(x(\xi,-D))$ the
topography representation in the $\xi$-variable, running along the bottom of
the flat strip. The flat bottom boundary $\eta=-D$ in the canonical domain is
mapped onto the topography profile, which in the physical domain reads as
$h(x)$. The constant $D$ will be determined by setting that the period in the
canonical domain is the same of the fluid domain.

Figure 1: The inverse conformal mapping. The bottom topography is flattened
out in the canonical domain.

We denote by $\mathbf{X}(\xi)$ and $\mathbf{Y}(\xi)$ the traces of the
respective harmonic functions along $\eta=0$, and by $\mathbf{X}_{b}(\xi)$ the
respective trace along the bottom $\eta=-D$. The $L$-periodic harmonic
function $y$ satisfies

$\displaystyle\begin{split}&y_{\xi\xi}+y_{\eta\eta}=0,\;\ \mbox{in}\;\
-D<\eta<0,\\\ &y=0,\;\ \;\ \mbox{at}\;\ \eta=0,\\\ &y=-1+H(\xi),\;\ \;\
\mbox{at}\;\ \eta=-D.\\\ \end{split}$

Using a Fourier transform $\mathbf{F}$ in the $\xi$-variable, we have that

$\displaystyle\begin{split}&y(\xi,\eta)=\mathbf{F}_{k\neq
0}^{-1}\left[\frac{-\coth(kD)\sinh(k\eta)\hat{H}}{\cosh(kD)}\right]+\frac{\eta(1-\hat{H}(0))}{D},\\\
\end{split}$ (12)

where $k=(2\pi/L)j$, $j\in\mathbb{Z}$, and

$\mathbf{F}_{k}[g(\xi)]=\widehat{g}(k)=\dfrac{1}{L}\int_{-L/2}^{L/2}g(\xi)e^{-ik\xi}d\xi,$
$\mathbf{F}^{-1}[\widehat{g}(k)](\xi)=g(\xi)=\sum_{j=-\infty}^{\infty}\widehat{g}(k)e^{{ik\xi}}.$

The Cauchy-Riemann equation $x_{\xi}=y_{\eta}$ yields

$x(\xi,\eta)=\mathbf{F}_{k\neq
0}^{-1}\left[\frac{i\coth(kD)\cosh(k\eta)\hat{H}}{\cosh(kD)}\right]+\frac{1-\hat{H}(0)}{D}\xi.$
(13)

From (13) we obtain

$\mathbf{X}_{b}(\xi)=x(\xi,-D)=\mathbf{F}_{k\neq
0}^{-1}\left[i\coth(kD)\hat{H}\right]+\frac{1-\hat{H}(0)}{D}\xi.$

Again by the Cauchy-Riemann equation $x_{\xi}=y_{\eta}$ and (12), we obtain an
alternative (and equivalent) expression:

$\mathbf{X}_{b}(\xi)=\frac{1-\hat{H}(0)}{D}\xi+\mathbf{F}^{-1}\left[\frac{i\coth(kD)}{\cosh^{2}(kD)}\widehat{H}\right]+\mathbf{F}^{-1}\left[i\tanh(kD)\widehat{H}\right].$
(14)

We have made the choice that the periodic solutions have the same period in
the canonical and physical domains. Let $L$ and $\lambda$ be the respective
periods, so that

$\mathbf{X}(\xi=L/2)-\mathbf{X}(\xi=-L/2)=\lambda$

with

$\big{<}\mathbf{X}_{\xi}\big{>}\equiv\frac{1}{L}\int_{-L/2}^{L/2}\mathbf{X}_{\xi}(\xi)d\xi=1.$
(15)

From (13) and $\big{<}\mathbf{X}_{\xi}\big{>}=\widehat{\mathbf{X}}_{\xi}(0)$,
it follows that

$1=\frac{\lambda}{L}=\frac{1-\widehat{H}(0,t)}{D},~{}~{}\mbox{where}~{}~{}D=1-\big{<}H\big{>}.$
(16)

This is the depth of the canonical channel.

The velocity potential $\phi$ is a harmonic function in both domains. Changing
variables in the bottom boundary condition $(\ref{F5})_{2}$ yields the
elliptic problem

$\displaystyle\begin{split}&\phi_{\xi\xi}+\phi_{\eta\eta}=0,\;\ \mbox{in}\;\
-D<\eta<0,\\\ &\phi=\mathbf{\Phi}(\xi,t),\;\ \;\ \mbox{at}\;\ \eta=0,\\\
&\phi_{\eta}=(F-\Omega)H_{\xi}(\xi)+\Omega HH_{\xi}(\xi),\;\
\mbox{at}~{}\eta=-D.\\\ \end{split}$

The notation is such that
$\phi(\xi,\eta,t)\equiv\overline{\phi}(x(\xi,\eta),y(\xi,\eta),t)$ and
$\psi(\xi,\eta,t)\equiv\overline{\psi}(x(\xi,\eta),y(\xi,\eta),t)$ is the
harmonic conjugate expressed in the canonical domain. Their traces along
$\eta=0$ are denoted by $\mathbf{\Phi}(\xi,t)$ and $\mathbf{\Psi}(\xi,t)$,
respectively. These time-dependent traces will be updated through the free
surfaces conditions. The harmonic conjugate $\psi$ satisfies

$\displaystyle\begin{split}&\psi_{\xi\xi}+\psi_{\eta\eta}=0,\;\ \mbox{in}\;\
-D<\eta<0,\\\ &\psi=\mathbf{\Psi}(\xi,t),\;\ \mbox{at}\;\ \eta=0,\\\
&\psi=-(F-\Omega)H(\xi)-\frac{\Omega}{2}H^{2}(\xi)+Q,\;\ \mbox{at}\;\
\eta=-D,\\\ \end{split}$

where $Q=Q(t)$. Solving these elliptic problems in Fourier space leads to

$\displaystyle\begin{split}&\phi(\xi,\eta,t)=\mathbf{F}^{-1}\left[\frac{\cosh(k(\eta+D))}{\cosh(kD)}\widehat{\mathbf{\Phi}}(k,t)+\left(\frac{i(F-\Omega)\widehat{H}+i\frac{\Omega}{2}\widehat{H^{2}}}{\cosh(kD)}\right)\sinh(k\eta)\right],\\\
&\psi(\xi,\eta,t)=\mathbf{F}^{-1}\left[\left(\widehat{\mathbf{\Psi}}(k,t)+\frac{(F-\Omega)\widehat{H}+\frac{\Omega}{2}\widehat{H^{2}}}{\cosh(kD)}\right)\frac{\sinh(k(D+\eta))}{\sinh(kD)}-\right.\\\
&-\left.\left(\frac{(F-\Omega)\widehat{H}+\frac{\Omega}{2}\widehat{H^{2}}}{\cosh(kD)}\right)\cosh(k\eta)\right]-\frac{Q(t)}{D}\eta.\end{split}$
(17)

In order to find $Q(t)$, we start with the Cauchy-Riemann equation
$\psi_{\eta}=\phi_{\xi}$ and use the periodicity in $\xi$ to obtain

$\widehat{\psi_{\eta}}(k=0,\eta,t)=\widehat{\phi_{\xi}}(k=0,\eta,t)=\frac{1}{L}\int_{-L/2}^{L/2}\phi_{\xi}(\xi,\eta,t)\,d\xi=0.$

From $(\ref{eq_Pot})_{2}$ we have

$\begin{split}\psi_{\eta}(\xi,\eta,t)=&\mathbf{F}^{-1}\left[\left(\widehat{\mathbf{\Psi}}(k,t)+\frac{(F-\Omega)\widehat{H}+\frac{\Omega}{2}\widehat{H^{2}}}{\cosh(kD)}\right)\frac{k\cosh(k(D+\eta))}{\sinh(kD)}-\right.\\\
&-\left.\left(\frac{(F-\Omega)\widehat{H}+\frac{\Omega}{2}\widehat{H^{2}}}{\cosh(kD)}\right)k\sinh(k\eta)\right]-\frac{Q}{D}.\end{split}$

Imposing $\langle\psi_{\eta}\rangle=0$ in the equation above yields

$\left(\widehat{\mathbf{\Psi}}(0,t)+(F-\Omega)\widehat{H}(0)+\frac{\Omega}{2}\widehat{H^{2}}(0)\right)\frac{1}{D}-\frac{Q}{D}=0,$

which is rewritten as

$Q(t)=\widehat{\mathbf{\Psi}}(0,t)+(F-\Omega)\widehat{H}(0)+\frac{\Omega}{2}\widehat{H^{2}}(0).$

Using $-\phi_{\eta}=\psi_{\xi}$ in (17), and evaluating over $\eta=0$, gives

$\begin{array}[]{l}\mathbf{\Phi}_{\xi}(\xi,t)=\mathbf{F}^{-1}\left[-i\coth(kD)\left(\widehat{\mathbf{\Psi}_{\xi}}(k,t)+\dfrac{(F-\Omega)\widehat{H_{\xi}}+\frac{\Omega}{2}\widehat{\partial_{\xi}H^{2}}}{\cosh(kD)}\right)\right].\end{array}$
(18)

The linearized kinematic condition is obtained from (5) and written as

$\frac{D}{Dt}\big{(}\overline{\zeta}(x,t)-y\big{)}=0,\;\ \mbox{where}\;\
\frac{D}{Dt}\equiv\partial_{t}+F\partial_{x}+\overline{\phi}_{y}\partial_{y}.$

In the canonical domain the wave elevation is given as $N(\xi,t)$, where
$\overline{\zeta}(\mathbf{X}(\xi),t)=y(\xi,N(\xi,t))$. The conformal mapping’s
Jacobian, evaluated along the frees surface, is denoted as
$J(\xi,0)=\mathbf{X}_{\xi}^{2}(\xi)\equiv M(\xi)^{2}$ (Nachbin [2003]). The
kinematic condition infers a material curve where, using the change of
variables $\partial_{\xi}=M(\xi)\partial_{x}$ and
$\partial_{\eta}=M(\xi)\partial_{y}$, we get

$\frac{\mathcal{D}}{\mathcal{D}t}\big{(}N(\xi,t)-\eta\big{)}=0,\;\
\mbox{where}\;\
\frac{\mathcal{D}}{\mathcal{D}t}\equiv\partial_{t}+\frac{F}{M(\xi)}\partial_{\xi}+\frac{\phi_{\eta}}{M(\xi)^{2}}\partial_{\eta}.$

The kinematic condition and Bernoulli law are rewritten from (5) to read as

$\displaystyle\begin{split}&N_{t}+\frac{F}{M(\xi)}N_{\xi}=\frac{{\phi}_{\eta}}{M(\xi)^{2}},\;\
\mbox{at}\;\ \eta=0,\\\
&\mathbf{\Phi}_{t}+\frac{F}{M(\xi)}\mathbf{\Phi}_{\xi}+\overline{\zeta}-\Omega\mathbf{\Psi}=-P(\mathbf{X}),\;\
\mbox{at}\;\ \eta=0,\end{split}$

in the canonical variables.

In a similar fashion as presented in Nachbin [2003], we can reflect our domain
about $\eta=0$ and obtain an odd extension of $y(\xi,\eta)$ on the strip
$-D<\eta<D$. In the physical domain the image is a reflected channel with the
undisturbed free surface ($y=0$) as the axis of symmetry. This does not change
any of the presentation above. Actually the boundary condition $y(\xi,0)=0$,
at the beginning of this section, is automatically satisfied due to the
reflection. Under this reflection, $y(\xi,\eta)$ is an odd function in $\eta$
that can be Taylor expanded as

$y(\xi,\eta)=y_{\eta}(\xi,0)\eta+\frac{y_{\eta\eta\eta}(\xi,0)}{6}\eta^{3}+\cdots,$

which for small wave elevations becomes

$y(\xi,N(\xi,t))\approx M(\xi)N(\xi,t).$

In the canonical domain the kinematic condition and Bernoulli law are

$\begin{array}[]{l}N_{t}=-\dfrac{F}{M(\xi)}N_{\xi}-\dfrac{\mathbf{\Psi}_{\xi}}{M(\xi)^{2}},\\\
\vskip 12.0pt plus 4.0pt minus
4.0pt\cr\mathbf{\Phi}_{t}=-M(\xi)N-\dfrac{F}{M(\xi)}\mathbf{\Phi}_{\xi}+\Omega\mathbf{\Psi}-P(\mathbf{X}(\xi)).\\\
\end{array}$ (19)

We are in position to rewrite the 2D Euler equations only in the canonical
variable $\xi$. The dependence on the $\eta$-variable is implicitly built-in
through the harmonic extension performed by the (Hilbert-type) Fourier
operator containing a $\coth(kD)$ as a multiplier. These operators are defined
below. The 2D Euler system is recast in the form

$\displaystyle\begin{split}&\mathbf{X}_{\xi}=1+\mathcal{C}_{k\neq
0}\left[\mathbf{F}^{-1}\left[\dfrac{\hat{H}_{\xi}}{\cosh(kD)}\right]\right]\\\
&\mathbf{\Phi}_{\xi}=-\mathcal{C}\left[{\mathbf{\Psi_{\xi}}}(k,t)+\mathbf{F}^{-1}\left[\dfrac{(F-\Omega)\widehat{H_{\xi}}+\frac{\Omega}{2}\widehat{\partial_{\xi}H^{2}}}{\cosh(kD)}\right]\right]\\\
&N_{t}=-\frac{F}{M(\xi)}N_{\xi}-\frac{\mathbf{\Psi}_{\xi}}{M(\xi)^{2}},\\\
&\mathbf{\Phi}_{t}=-M(\xi)N-\frac{F}{M(\xi)}\mathbf{\Phi}_{\xi}+\Omega\mathbf{\Psi}-P(\mathbf{X}(\xi)),\end{split}$
(20)

where

$\displaystyle\begin{split}&\mathrm{H}(\xi)=h(\mathbf{X}_{b}(\xi)),\\\
&\mathbf{X}_{b}(\xi)=\xi+\mathcal{C}_{k\neq
0}\left[\mathbf{F}^{-1}\left[\frac{\widehat{H}}{\cosh^{2}(kD)}\right]\right]+\mathbf{F}^{-1}\left[i\tanh(kD)\widehat{H}\right].\end{split}$
(21)

The Fourier operators are:

$\mathcal{C}[\cdot]\equiv\mathbf{F}^{-1}\left[i~{}\coth(kD)\mathbf{F}_{k}[\cdot]\right]\mbox{~{}~{}and~{}~{}}\mathcal{C}_{k\neq
0}[\cdot]\equiv\mathbf{F}^{-1}[i~{}\coth(kD)\mathbf{F}_{k\neq 0}[\cdot]].$

##  3 Numerical Method

Our main goal in this study is to exhibit (numerically) the spontaneous
formation of critical layers and its robustness in overcoming a bottom
topography. In the stationary wave regime, particle pathlines coincide with
streamlines (Ribeiro Jr et al. [2017]). When the surface disturbance is not a
traveling wave we cannot use the streamfunction to depict the pathlines.
Hence, in order to see the formation and respective dynamics of the Kelvin’s
cat eye structure we need to compute several particle orbits. These tracers
will allow for the visualization of the critical layer. Due to the presence of
a bottom topography we adopt a conformal mapping from the canonical (flat
strip) domain onto the physical domain. All the dynamics is computed in the
(simpler) canonical domain and then mapped onto the physical domain for a
proper visualization of the pathlines.

In order to describe our numerical procedure the Euler formulation, presented
earlier, is now connected with the tracer dynamics. Recall that particle
orbits $({x}(t),{y}(t))$ are governed by system (6). In the canonical domain
the pre-image of a trajectory is given by $(\xi(t),\eta(t))$, where the
mapping onto the physical domain yields
$(\tilde{x}(t),\tilde{y}(t))=(x(\xi(t),\eta(t)),y(\xi(t),\eta(t)))$. The
dynamical system for computing the particle pathlines in the canonical domain,
is

$\left\\{\begin{array}[]{l}\dfrac{d\xi}{dt}(t)=\dfrac{\phi_{\xi}(\xi,\eta,t)}{J(\xi,\eta)}+\dfrac{(\Omega~{}y(\xi,\eta)+F)~{}x_{\xi}(\xi,\eta)}{J(\xi,\eta)},\\\
\vskip 12.0pt plus 4.0pt minus
4.0pt\cr\dfrac{d\eta}{dt}(t)=\dfrac{\phi_{\eta}(\xi,\eta,t)}{J(\xi,\eta)}-\dfrac{(\Omega~{}y(\xi,\eta)+F)~{}y_{\xi}(\xi,\eta)}{J(\xi,\eta)},\\\
\vskip 12.0pt plus 4.0pt minus 4.0pt\cr\xi(0)=\xi_{0},\;\ \eta(0)=\eta_{0}.\\\
\end{array}\right.$ (22)

The Jacobian is expressed as
$J(\xi,\eta)=x_{\xi}^{2}(\xi,\eta)+y_{\xi}^{2}(\xi,\eta)$. The particle is
initially located at $(\xi_{0},\eta_{0})$, the pre-image of $(x_{0},y_{0})$.
The physical pathline is obtained from

$\tilde{x}(t)=x(\xi(t),\eta(t))=\mathbf{F}_{k\neq
0}^{-1}\left[\frac{i\coth(kD)\cosh(k\eta)\hat{H}}{\cosh(kD)}\right]+\frac{1-\hat{H}(0)}{D}\xi,$
(23)
$\displaystyle\begin{split}&\tilde{y}(t)=y(\xi(t),\eta(t))=\mathbf{F}_{k\neq
0}^{-1}\left[\frac{-\coth(kD)\sinh(k\eta)\hat{H}}{\cosh(kD)}\right]+\frac{\eta(1-\hat{H}(0))}{D}.\\\
\end{split}$ (24)

Recall that in the canonical variables the topography is denoted by $H(\xi)$.
All Fourier transforms $\mathbf{F}$ are in the $\xi$-variable, computed
numerically through an FFT (Fast Fourier Transform). In computing the vector
field of the dynamical system (22) we use the potential $\phi$ given by its
Fourier representation

$\phi(\xi,\eta,t)=\mathbf{F}^{-1}\left[\frac{\cosh(k(\eta+D))}{\cosh(kD)}\widehat{\mathbf{\Phi}}(k,t)+\left(\frac{i(F-\Omega)\widehat{H}+i\frac{\Omega}{2}\widehat{H^{2}}}{\cosh(kD)}\right)\sinh(k\eta)\right].$
(25)

The dynamical system is not autonomous. Therefore its vector field is updated
through the free surface conditions:

$\left\\{\begin{array}[]{l}N_{t}=-\dfrac{F}{M(\xi)}N_{\xi}-\dfrac{\mathbf{\Psi}_{\xi}}{M(\xi)^{2}}\\\
\vskip 12.0pt plus 4.0pt minus
4.0pt\cr\mathbf{\Phi}_{t}=-M(\xi)N-\dfrac{F}{M(\xi)}\mathbf{\Phi}_{\xi}+\Omega\mathbf{\Psi}-P(\mathbf{X}).\\\
\vskip 12.0pt plus 4.0pt minus 4.0pt\cr\end{array}\right.$ (26)

In other words, we evolve the boundary data for the velocity potential as well
as for the wave profile. Within the fluid body, updated values of $\phi$ and
its derivatives are readily available from expression (25). The boundary data
of its harmonic conjugate is obtained from

$\begin{array}[]{l}\widehat{\mathbf{\Psi}}(k,t)=i\tanh(kD)\widehat{\mathbf{\Phi}}(k,t)-\left(\dfrac{(F-\Omega)\widehat{H}+\frac{\Omega}{2}\widehat{H^{2}}}{\cosh(kD)}\right).\\\
\vskip 12.0pt plus 4.0pt minus 4.0pt\cr\end{array}$ (27)

The term $P(\mathbf{X})$ in the dynamic condition above is used only when a
traveling pressure distribution is considered along the free surface. We use
the fourth order Runge-Kutta method for numerically evolving systems (22) and
(26).

As mentioned all Fourier transforms are approximated by the FFT on a uniform
grid, with all derivatives performed in Fourier space (Trefethen [2001];
Milewski & Tabak [1999]). The computational grid in the canonical domain is
given by $\xi\in[-L/2,L/2)$, with $N$ uniformly spaced points, with grid size
$\Delta\xi=L/N$. This corresponds to a nonuniform grid in physical space. A
typical resolution has $N=2^{13}$ with a time step $\Delta t=0.05$. A Kelvin
cat eye structure is typically captured with a cloud of $90$ tracers.

We have shown that $\mathbf{X}_{b}(\xi)$ and $H(\xi)$ are coupled in a
nontrivial fashion. We need to beforehand compute the topography profile
$H(\xi)$ in the canonical domain. This depends on a nontrivial composition of
the form $h(x(\xi,-D))$. This topographic composition is pre-processed using
an iterative method as presented by Flamarion et al. [2019]. The iterates are
labelled by a superscript $l$. The iterative scheme has the following
structure:

$\displaystyle\begin{split}&\mathbf{X}^{l}_{b}(\xi)=\xi+\mathcal{C}_{k\neq
0}\left[\mathbf{F}^{-1}\left[\frac{\widehat{H^{l}}}{\cosh^{2}(kD)}\right]\right]+\mathbf{F}^{-1}\left[i\tanh(kD)\widehat{H^{l}}\right].\\\
&H^{l+1}(\xi)=h(\mathbf{X}_{b}^{l}(\xi)),\end{split}$ (28)

where the initial step is based on the identities
$\mathbf{X}_{b}^{0}(\xi)=\xi$ and $H^{1}(\xi)=h(\xi)$. The stopping criteria
used is

$\displaystyle\max_{\xi\in[-L/2,L/2)}\big{|}H^{l+1}(\xi)-H^{l}(\xi)\big{|}<\varepsilon,$

where $\varepsilon$ is a given tolerance. In our simulations, we used
$\varepsilon=10^{-12}.$

When we have a traveling wave along the free surface the pathlines are
identical to the streamlines. When this is not the case one can still compute
the streamlines. The potential’s harmonic conjugate $\psi$ is readily
available from the Fourier representation

$\begin{array}[]{r}\psi(\xi,\eta,t)=\mathbf{F}^{-1}\left[\left(\widehat{\mathbf{\Psi}}(k,t)+\dfrac{(F-\Omega)\widehat{H}+\frac{\Omega}{2}\widehat{H^{2}}}{\cosh(kD)}\right)\dfrac{\sinh(k(D+\eta))}{\sinh(kD)}-\right.\\\
-\left.\left(\dfrac{(F-\Omega)\widehat{H}+\dfrac{\Omega}{2}\widehat{H^{2}}}{\cosh(kD)}\right)\cosh(k\eta)\right]-\dfrac{Q(t)}{D}\eta.\end{array}$
(29)

In the moving frame we have that

$X=x-ct\mbox{ and }Y=y,$ (30)

where $c$ is the wave speed given by (7). The streamfunction $\psi_{T}$ is
obtained by putting $\overline{\psi}$ together with the shear flow in the form

$\psi_{T}(X,Y,t):=\overline{\psi}(X,Y,t)+\frac{\Omega Y^{2}}{2}+(F-c)Y.$ (31)

At any instant of time the streamfunction can be evaluated on a uniform grid
in the canonical domain and mapped onto points in the physical domain in order
to trace its level curves there. Yet, at each given time, we denote
$\psi_{T}(t)\equiv\psi_{T}(X,Y,t)$ and consider $\|\psi_{T}(t)\|_{1}$ as its
1-norm evaluated on the computational grid.

##  4 Results

It is well known in the literature that stagnation points and critical layers
can be formed beneath a periodic traveling wave in the presence of vorticity
(Teles da Silva & Peregrine [1988]; Vasan & Oliveras [2014]; Ribeiro Jr et al.
[2017]; Nachbin & Ribeiro Jr [2017]). There are different forms of finding the
initial wave profile, depending whether we have in mind an irrotational
(Nachbin & Ribeiro Jr [2014]) or a rotational surface traveling wave (Ribeiro
Jr et al. [2017]). See the references within these two papers.

In the present study the waves are linear and dispersive, so we do not have a
traveling wave. Nevertheless for weak dispersion the propagating waves change
shape very slowly. A right propagating disturbance is obtained by pre-
processing our desired wave profile over a flat bottom and keeping the right
going mode. The initial wave elevation is usually positioned at a reasonable
distance from the bottom topography so that there the Jacobian is effectively
equal to 1, and therefore $N(\xi,0)=N_{0}(\xi)\equiv\zeta(x,0)=\zeta_{0}(x)$.
Hence the wave elevation representation in the canonical and physical domains
are identical.

Figure 2: A nearly stationary stream function. The small relative error is
given by Error$(t)={\|\psi_{T}(t)-\psi_{T}(0)\|_{1}}/{\|\psi_{T}(0)\|_{1}}$.

In the present study the periodic boundary condition is only for numerical
purposes, namely for the Fourier spectral method. We consider wave
disturbances localized in space which undergo reflection and transmission over
a bottom variation of compact support. Hence, over the time interval of our
simulations, no disturbances should be observed at the endpoints of our
computational domain. We consider a tabletop modulated wavetrain localized in
space which resembles a periodic wave in its central region. The goal is to
observe a segment of pathlines with a Kelvin’s cat eye (critical layer)
structure which resembles that of periodic waves.

###  4.1 Wave-current interaction for a slowly varying wavetrain.

We start with a tabletop wave profile of the form

$N_{0}(\xi)=\alpha e^{-2(\sigma(\xi-\xi_{0}))^{s}}\cos(k(\xi-\xi_{0})),$ (32)

where $\alpha=10^{-3}$, $\sigma=6\cdot 10^{-3}$, $s=8$ and $k=2\pi/50$. The
very first simulation is for an unforced case (with $h(x)=0$ and $P(x)=0$)
where the underlying current is such that $\Omega=-18$ and $F=-9$. The
dispersion is very weak in this case. The wave profile barely changes during
the simulation. As shown in figure 2 the difference between the evolving
profile and the initial disturbance is very small. The error was computed over
a large time interval ($t\in[0,340]$) corresponding to a propagation distance
of about 60 wavelengths. Over this time interval the level curves of the
streamfunction are a good approximation for the pathlines. The vector filed of
equation (8) is effectively autonomous. The initial phase portrait obtained
from the streamfunction $\psi_{T}(0)$ is depicted in figure 3. The central
region is typical of a phase portrait for a periodic traveling wave in the
presence of constant vorticity (Ribeiro Jr et al. [2017]). The novelty here is
that the wavetrain is effectively of compact support, rather than periodic as
in the previous studies. The super-Gaussian envelope, that provides a tabletop
pattern over the central region of the wavetrain, was designed so that we
could observe a Kelvin cat eye structure similar to periodic waves, as
presented by Teles da Silva & Peregrine [1988], Vasan & Oliveras [2014] and
Ribeiro Jr et al. [2017] . This is the case in the region located
approximately within the interval [400,600]. This is a first numerical display
of the cat eye structure dynamically adjusting to the endpoints of the
wavetrain. The wave is propagating to the right and it is remarkable that
beneath the very small (weakly dispersive) oscillatory tail, the method can
detect a number of diminishing cat eyes (recirculation regions) fading to the
left. The small amplitude wave is linear and the Kelvin cat eye structure is
narrow and located near the bottom. Nevertheless, due to a nontrivial
vorticity (Ribeiro Jr et al. [2017]) the cat eye structure is detached from
the bottom exhibiting both types of critical points: saddles and a center. The
cat eye structure and associated critical layer is robust in the present of
the weak dispersion. As the oscillatory tail of wave changes, the
recirculation regions adjust accordingly. In particular, in analogy with the
development of an Airy-type solution, as the lagging oscillatory tail
develops, a small new recirculation (cat eye) region is generated.

Figure 3: Streamlines obtained from $\psi_{T}(0)$. The dynamics is nearly
stationary and the streamlines resemble pathlines. The vertical-axis to the
left has the scales corresponding to the wave elevation, namely the thicker
line along the top. The vertical-axis to the right has values on the depth
scale and corresponds to the scales of the streamlines, indicating that the
critical layer is very close to the bottom.

###  4.2 Critical layer robustness due to wave-current-topography interaction.

Now we introduce a topography located in the middle of our computational
domain. There is no pressure distribution along the free surface. The
topography has compact support and a tabletop configuration, in order to
locally resemble a periodic depth variation. The bottom profile is given by

$h(x)=\delta e^{-2(\sigma(x-\tilde{x}_{0}))^{s}}\cos(k_{b}(x-\tilde{x}_{0})),$
(33)

where $\delta\ll 1$. The right going initial wave profile is the same as
before. We have kept the super-Gaussian (tabletop) modulation the same, for
the wave and for the topography, so that they have comparable lengths. The
physical domain is schematically depicted in figure 4. The initial disturbance
has a Kelvin cat eye structure formed beneath it and will eventually interact
with the topography. At the same time a wave is generated by the current
topography interaction, as studied in Flamarion et al. [2019]. We will examine
the submarine pathline structure in the presence of all these features.

Figure 4: The physical domain is schematically depicted, for the wave-
topography interaction, in the presence of a current of constant vorticity.
Two regimes are considered: $\lambda\approx l$ (comparable wavelengths) and
$\lambda\gg l$ (rapidly varying topography).

Figure 5: Comparable wavelengths regime ($\lambda=2l$). The black dots
indicate the initial position for each tracer, while the white dot indicates
their final position, at the end of the time interval. These are depicted in
the laboratory (fixed) frame. The $y$-axis to the left indicates the scale for
the free surface disturbance. The $y$-axis to the right indicates the depth
values. The tracer and respective structures are very near the bottom.
Nevertheless our numerical method accurately captures the detailed structure
of the tracers’ orbits. An animation (video01) can be found in the
supplementary material of this article.

We consider two wave-topography regimes. First we consider the regime where
the surface wavelength is comparable to that of the topography and we let
$l/\lambda=0.5$. Then we consider the case where the topography is rapidly
varying: $l/\lambda=0.05$. The parameters used for the underlying flow and the
initial disturbance (32) are: $\Omega=-18$, $F=-9$, $\alpha=10^{-3}$,
$\sigma=3\cdot 10^{-3}$, $s=8$, $k=2\pi/50$. For the topography we used
$\delta=2\cdot 10^{-3}$, $k_{b}=2k$ and $k_{b}=20k$ in expression (33). The
amplitude of the topography needs to be small so that the height of the waves
generated by the current-topography interaction is comparable to that of the
initial disturbance. This is clearly seen in the figures that follow.

Recall that in the present (non-stationary) regime the streamlines are not the
same as the pathlines. In order to visualize the pathlines we need to track
the orbits of many tracers, each tracer solving the dynamical system (22).

In figure 5 we consider the $l/\lambda=0.5$ case, first in the (fixed)
laboratory frame. In the top snapshot (at time $t=34$) we see the initial
(right going) disturbance at the left of the topography and the current-
generated wave above the topography. We focus on the initial disturbance and
release a cloud tracers below this wave. The initial cloud (depicted by black
particles near the bottom) is positioned to the left of $x=500$. At time
$t=34$ the cloud has traveled to the right, is positioned below the wave
train, and is depicted by a bunch of white particles. The inset provides a
zoom of the white cloud and we clearly see the tracers in a Kelvin cat eye
formation. The interaction of this structure with the topography is displayed
in the middle picture of figure 5. At time $t=118$ we observe the pathlines in
a terrain-following pattern, moving over the topography. The cat eye structure
is distorted by the bottom undulations but retains its integrity. At this
stage the rightgoing initial wave is interacting with the stationary wave
generated by the current-topography mechanism. There is a nontrivial wave
pattern above the topography, while a wave generated by the current-topography
interaction is moving downstream and out of the picture. The computational
domain is wider then that shown by these pictures. The wave pattern due to the
current-topography interaction is in accordance with that observed in
Flamarion et al. [2019], where there were standing waves above the topography
together with waves moving downstream. In the bottom picture at time $t=170$
the initial rightgoing disturbance is leaving the variable bottom region. The
inset shows the tracers back to (nearly) the same formation as at time $t=34$.

Figure 6: Comparable wavelengths regime ($\lambda=2l$). The black dots
indicate the initial position for each tracer, while the white dot indicates
their final position, at the end of the time interval. These are depicted in
the wave’s moving frame. The $y$-axis to the left indicates the scale for the
free surface disturbance. The $y$-axis to the right indicates the depth values
near the bottom. The respective cat eye structures are very small and very
near the bottom. Nevertheless the accuracy of the numerical method captures in
great detail the Kelvin cat eye structure. It is remarkable to observe the
robustness of the recirculation region, adjusting to depth variations. An
animation (video01) can be found in the supplementary material of this
article.

Figure 7: A rapidly varying topography ($\lambda\gg l$). We track tracers
positioned beneath the initial rightgoing disturbance. The black dots indicate
the initial position for each tracer, while the white dot indicates their
final position, at the end of the time interval. These are depicted in the
wave’s moving frame. For the middle and bottom picture we erased the initial
position for a better visualization of the cat eye. The $y$-axis to the left
indicates the scale for the free surface disturbance. The $y$-axis to the
right indicates the depth values near the bottom. The respective cat eye
structures are very small and very near the bottom. Nevertheless the accuracy
of the numerical method captures in great detail the Kelvin cat eye formation.

The wave profile, due to the initial disturbance (32), changes very slowly.
Hence we proceed as with traveling waves and display particle orbits in the
moving frame associated with the wave speed $c$, given in (7). These will
reveal the Kelvin cat eye flow structure and the respective critical layer. In
figure 6 the Kelvin cat eye structure is evident, with pathlines forming a
closed recirculation region about a critical point (centre), namely the
stagnation point seen in Ribeiro Jr et al. [2017]. At the top picture, the
inset displays tracers that were initially located at the back points, but
have moved along the cat eye contours reaching the white points at time
$t=34$. In the middle picture, at time $t=118$, the inset displays material
curves, where a collection of tracers is following a given contour of the cat
eye structure. In the middle part (in red) we have a closed material curve
where the white points indicate the current position of each tracer. We have
the same particles (cloud of tracers) in all snapshots of figure 6. The top
material curve (in yellow) and the bottom material curve (in blue) eventually
connect to the saddle points, as identified by Ribeiro Jr et al. [2017] and
which were here displayed in figure 3. The cat eye structure is very robust,
as seen deforming to adjust to the topography while practically returning to
its undisturbed shape at time $t=170$, in the bottom picture of figure 6.

We now consider a rapidly varying topography. The initial disturbance and
collection of markers are the same from the previous example. The current-
topography interaction generates rapidly varying wavetrains. In the present
regime there is a stationary wavetrain above the topography and another
wavetrain moving upstream. The initial rightgoing disturbance collides with
the upstream wavetrain not affecting at all the Kelvin cat eye, as shown in
figure 7, at three different times. As the wavetrain propagates over the
rapidly varying topography the robustness feature is modified, in some sense.
In the top part of figure 8 we observe the disintegration of the material
curves, which formed the contours of the cat eye. As soon as the wavetrain
reaches the other side of the depth variations, it is remarkable to observe
that the same group of markers form again the same material curves of the cat
eye (recirculation) region.

Figure 8:  Continuation of figure 7. As the wavetrain crosses over the rapidly
varying topography, the bottom’s fine features destroys the material curves’
well defined formation, as seen at the top picture. Nevertheless as soon as
the wavetrain leaves the variable depth region, the material curves
constituted with the same particles as before, reconfigures itself in (namely)
the same cat eye as observed in figure 7. The robustness of the recirculation
region, in this example, is remarkable.

###  4.3 The spontaneous critical layer formation

The spontaneous critical layer formation is associated with the spontaneous
formation of stagnation points (in the moving frame), and the respective
appearance of the Kelvin cat eye structure. Up to this point we have tracers
released beneath the initial wave disturbance. This initial wave disturbance
immediately establishes the submarine phase-plane picture as seen in figure 3.
In the present case, the free surface starts at rest. We then compute the
submarine orbits for tracers placed beneath the wave generated by the current-
topography interaction. At time $t=0^{+}$ there are no waves on the surface
nor critical points in the bulk. The flow and topographic parameters are the
same as before, but $k_{b}=k$.

Figure 9: The spontaneous formation of the Kelvin cat eye structure. The top
picture shows tracers positioned over a valley of the topography before the
simulation starts. They are positioned in order to capture a cat eye
formation. The middle picture shows oscillating tracers due to their relative
motion with the topography. The black dot indicates the tracer’s initial
position while the white dot its final position at time $t=10$. The bottom
picture, at a later time ($t=82$), shows a wavetrain generated by the
topography, propagating downstream to the right. The tracers’ orbits are
depicted in the wave’s moving frame, as we observe the topography moving to
the left. Again the markers first oscillate vertically, after which we see the
onset of the cat eye formation through the short orbit “filament” connected to
the final position (white dot). The followup for this cat eye formation is
presented in figure 10. The left $y$-axis has the wave scale while the right
$y$-axis the depth scale. An animation (video02) can be found in the
supplementary material of this article.

At time $t=0$ we distribute tracers at positions that would correspond to a
cat eye. This is shown at the top picture of figure 9. Then a standing wave
appears, and at time $t=10$ we observe the tracers moving up and down (middle
picture of figure 9). The color code (see version online) is the same as
before, according to the position along the cat eye. Eventually a downstream
propagating wavetrain is generated by the topography. Note that the cloud of
tracers is always centered about the same position because we are in the
moving frame with respect to this wave. Therefore the topography seems to be
moving to the left. We see at the bottom-picture of figure 9 that, after a
preliminary oscillation, the markers start tracing orbits (at $t=82$) which
eventually will look like a cat eye. Erasing the initial oscillation, this is
clearly seen a bit later in time ($t=100$), at the top of figure 10. In the
other two pictures of figure 10 the Kelvin cat eye is effectively a stationary
structure.

Figure 10: Continuation of figure 9. We have erased the initial positions of
all tracers for a better visualization. As time evolves, from $t=100$ up to
$t=348$, we clearly see the Kelvin cat eye structure and the associated
critical layer. The overall picture for the pathlines is similar to that of
figure 3, but a large amount of markers would be needed to clearly trace them.
An animation (video02) can be found in the supplementary material of this
article.

The most common parameters used in numerical studies for establishing the
presence of critical layers under rotational waves are the vorticity, the
energy and the wavelength (Ko & Strauss [2008a, b]; Nachbin & Ribeiro Jr
[2017]; Ribeiro Jr et al. [2017]; Teles da Silva & Peregrine [1988]; Vasan &
Oliveras [2014]). In order to study the natural emergence of a critical layer,
Johnson [2012] suggests a pressure distribution along the surface.

Figure 11: Wave generated by a localized pressure distribution. Top two
pictures show the markers initially positioned in order to capture the
spontaneous cat eye formation, both in the (fixed) laboratory frame as well as
in the wave’s moving frame. The bottom two pictures are at time $t=4$ at the
onset of the pressure wave formation. A steady depression is observed where
the pressure is applied while two disturbances are about to propagate
outwards. The markers are set to follow the rightgoing pulse-shaped wave. In
the laboratory frame we see a structure moving with the wave. In the moving
frame the markers move very slowly. In figure 12 a large time interval is
considered so that our cloud of markers has fully traced the cat eye
structure. The left $y$-axis has the wave scale while the right $y$-axis the
depth scale. An animation (video03) can be found in the supplementary material
of this article.

With this in mind, at time $t>0$, we apply a localized (bell-shaped) high
pressure distribution

$P(x)=0.2e^{-\sigma x^{2}},$

along the free surface. The bottom is flat ($h(x)=0$) and the free surface
initially undisturbed ($N_{0}(\xi)=0,~{}\mathbf{\Phi}(\xi,0)=0)$. The
underlying shear flow is such that $\Omega=-18$ and $F=-9$.

In figures 11 and 12 we display different stages of the submarine dynamics,
both in the laboratory frame as well as in the wave’s moving frame. The
corresponding wave speed was found by inspection. The top two pictures are at
$t=0$, when the bell-shaped pressure is about to be applied and a cloud of
markers has been prepared in the fluid body. At time $t=4$ the steady pressure
distribution has generated two pulse-shaped waves which are about to propagate
to the left and right accordingly. The tracers’ dynamical systems are set to
follow the rightgoing wave, as depicted by the two bottom pictures. The black
dot is the starting point while the white dot the final position for the
respective time interval. In the bottom picture of figure 11 these are almost
at the same position because the dynamics is relative to the wave position.
Hence in the moving frame each tracer position is slowly varying. The steady
high pressure distribution creates a steady free surface depression.

Note that the waves generated are weakly dispersive. They are similar to a
traveling pulse-shaped wave for a certain time interval as seen at the top of
figure 12. At time $t=99$ we observe a small dispersive tail appearing behind
the leading pulse. Nevertheless the cat eye structure, spontaneously formed,
remained unchanged as the wave propagated for more than 10 pulse widths in
distance. As shown figure 13 the streamfunction is effectively steady during
the time interval $t\in[20,40]$ where, up to a relative error of $O(10^{-8})$,
the initial streamfuntion remains (approximately) unchanged. During this
period of time the streamlines effectively coincide with the pathlines.

Based on the discussion above, in figure 14 we plot the streamlines for
$\psi_{T}(30)$. We have a single cat eye (recirculation) region, traveling
with the wave speed. A critical layer divides a region into a rightgoing and a
leftgoing flow. As in earlier works, the critical layer goes through the
recirculation region. But in this case the critical layer then becomes a sharp
interface emerging from what would have been a saddle point in the usual cat
eye formation. The velocity field is continuous so there is no velocity jump
crossing the stagnation segment. Above the stagnation segment the flow is to
the left while immediately below the interface the flow is to the right. In
the present example we have an isolated stagnation point in the form of a
center, in the middle of the recirculation region, and stagnation segments
connected to the isolated Kelvin cat eye. This is a scenario not reported
before, in particular for a non-stationary wave regime.

Figure 12: Continuation of figure 11. The top two pictures show the tracers
following the expected structure. The bottom picture indicates that after a
large propagation distance our cloud of markers has fully traced the Kelvin
cat eye. An animation (video03) can be found in the supplementary material of
this article.  Figure 13: Effectively-steady streamfunction $\psi_{T}$. The
small relative error, during the time interval $t\in[20,40]$, is computed as
Error$(t)={\|\psi_{T}(t)-\psi_{T}(20)\|_{1}}/{\|\psi_{T}(20)\|_{1}}$. During
an initial small time interval the stream function is quite different from its
reference state $\psi_{T}(20)$.  Figure 14:  Effective phase portrait using
$\psi_{T}(30)$. The left $y$-axis has values on the scale of the free surface,
depicted at the top in a thicker dark line. The right $y$-axis has values on
the depth scale, so that we can see the closed streamlines near the bottom.

##  5 Conclusions

We have studied non-stationary rotational surface waves, in the presence of
constant vorticity. We consider the regime where the Kelvin cat eye structure
can be observed. In the linear regime one has a very narrow structure near the
bottom. The accurate numerical method captured all details of the cat eye
structure and the respective critical layer separating left from rightgoing
streams. Particle-trajectories were computed and visualized by evolving the
respective submarine dynamical system with a cloud of tracers. Waves were
generated from either the current-topography interaction or by a surface
pressure distribution suddenly imposed.

##  6 Acknowledgements

The work of M.F. was supported in part by CNPq-Science without Borders under
(PDE) 200920/2014-6 and by CNPq-Cotas do Programa de Pós-Graduação (GM/GD)
140773/2014-2. The work of A.N. was supported in part by CNPq under (PQ-1B)
304671/2017-7 and FAPERJ Cientistas do Nosso Estado project no.
E-26/203.023/2017. R.R.-Jr is grateful to IMPA for the research support
provided during the Summer Program of 2018 and to University of Bath for
extended visit to the Department of Mathematical Sciences. A.N. and R.R.-Jr
thank Prof. Robin Johnson for his comments during theirs visits to Vienna.
A.N. is grateful to the Erwin Schrödinger International Institute for
Mathematics and Physics, for the support and hospitality during the ESI
Programme on “Mathematical Aspects of Physical Oceanography”. R.R.-Jr thanks
Prof. Adrian Constantin for his support and hospitality during his visit to
the Department of Mathematics of the University of Vienna.

## References

  * Buhler [2009] Buhler, O. 2009 Waves and Mean Flows. Cambridge University Press.
  * Constantin [2011] Constantin, A. 2011 Nonlinear Water Waves with Applications to Water-Current Interactions and Tsunamis. CBMS-NSF Regional Conference Series in Applied Mathematics 81, SIAM.
  * Constantin [2017] Constantin, A. 2017 Nonlinear water waves: introduction and overview, Phil. Trans. R. Soc. A 376.
  * Constantin et al. [2016] Constantin, A., Strauss, W. & Varvaruca, E. 2016 Global bifurcation of steady gravity water waves with critical layers. Acta Math. 217, 195-262.
  * Flamarion et al. [2019] Flamarion, M. V. & Milewski, P. A. & Nachbin A. 2019 Rotational waves generated by current-topography interaction. Stud Appl Math, 142, 433-464.
  * Johnson [2012] Johnson, R. S. 2012 Models for the formation of a critical layer in water wave propagation. Phil. Trans. R. Soc. A 370, 1638-1660.
  * Francius & Kharif [2017] Francius, M. & Kharif, C. 2017 Two-dimensional stability of finite-amplitude gravity waves on water of finite depth with constant vorticity. J. Fluid Mech. 830, 631-659.
  * Ko & Strauss [2008a] Ko, J. & Strauss, W. 2008a Effect of vorticity on steady water waves. J. Fluid Mech. 608, 197-215.
  * Ko & Strauss [2008b] Ko, J. & Strauss, W. 2008b Large-amplitude steady rotational water waves. Eur. J. Mech./B Fluids 27, 96-109.
  * Milewski & Tabak [1999] Milewski, P. A. & Tabak, E. G. 1999 A pseudo-spectral algorithm for the solution of nonlinear wave equations. SIAM J. Sci. Comp. 21 (3), 1102-1114.
  * Nachbin [2003] Nachbin, A.  2003 A terrain-following Boussinesq system. SIAM J. Appl Math. 63, 905-922.
  * Nachbin & Ribeiro Jr [2014] Nachbin A. & Ribeiro Jr, R. 2014 A boundary integral formulation for particle trajectories in Stokes waves 2014 DCDS-A 34 (8), 3135-3153.
  * Nachbin & Ribeiro Jr [2017] Nachbin A. & Ribeiro Jr, R. 2017 Capturing the beneath water waves Phil. Trans. R. Soc. A, 376.
  * Ribeiro Jr et al. [2017] Ribeiro Jr R., Milewski P.A., & Nachbin A. 2017 Flow structure beneath rotational water waves with stagnation points. J. Fluid Mech. 812, 792-814.
  * Teles da Silva & Peregrine [1988] Teles da Silva, A.F. & Peregrine, D.H. 1988 Steep, steady surface waves on water of finite depth with constant vorticity. J. Fluid Mech. 195, 281-302.
  * Trefethen [2001] Trefethen, L. N. 2001 Spectral Methods in MATLAB. Philadelphia: SIAM..
  * Soontiens et al. [2010] Soontiens, N., Subich C. & Stastna, M. 2010 Numerical simulation of supercritical trapped internal waves over topography. Phys. Fluids 22,16605.
  * Vasan & Oliveras [2014] Vasan, V. & Oliveras, K. 2014 Pressure beneath a traveling wave with vorticity constant. DCDS-A 34 (8), 3219-3239.
  * Whitham [1974] Whitham, G. B. 1974 Linear and Nonlinear Waves. New York: John Wiley and Sons.


Critical load arrangement investigation sec:critical_load_arrangements
Polarity sequences sec:polarity_diagrams Influence lines can be used to
identify the critical load arrangements for a given continuous beam system.
The design problem is restricted to positive UDL values only (no uplift) which
can be activated on or off (1 or 0 as activation factors). Therefore, by
integrating the influence line for each individual beam $i$, one can evaluate
the net contribution (positive or negative) a given beam causes in terms of
bending moments and shear forces at the influence line (IL) location when
subjected to a positive, unit UDL. The net-contribution of each beam can be
either positive or negative at the IL location, that is “hogging or sagging”
for bending moments and “clockwise or anti-clockwise” for shear forces,
respectively, which is termed as the polarity of that particular beam. This
procedure is shown in Figure fig:influence_lines_to_polarity_zones.

[width=height=keepaspectratio]figure_3_influence_line_to_polarity.pdf

Figure : An exemplary process of arriving from influence line plots (top row)
to polarity sequences (bottom row) via integrated influence lines (middle row)
for a) major axis bending moment $M_{y}$ and b) major axis shear force $V_{z}$
about the specified influence line (IL) location.

The last row of Figure  therefore reflects a particular polarity sequence for
a given IL location, which can be directly used to identify the critical load
arrangement for that IL location. When all beams of positive polarity are
loaded, then the maximum positive internal forces are generated at the IL
location, and vice-versa, loading the negative polarity members leads to the
maximum negative internal forces.

##  Polarity zones

A rigorous qualitative study of the polarity sequences for different IL
locations and design scenarios revealed 5 unique polarity sequences that occur
along specific segments of a given beam span termed polarity zones, which are
illustrated in Figure  for the central beam highlighted in red.

[width=0.99height=keepaspectratio]figure_4_polarity_zones_and_sequences.pdf

Figure : Polarity zones that occur along various span segments of a $m=5$
homogeneous beam system of equal span and cross-sectional properties. The same
zones and sequences, although at different boundaries, occur in heterogeneous
(varying UDL and span) systems.

These 5 polarity zones are common to all members of both homogeneous (equal
spans and cross-sections) as well as heterogeneous continuous beam systems,
although the exact boundaries between one zone varied depending on the
relative magnitude of spans and cross-section properties. The sequences
identified in Figure  also apply to larger beam systems with the polarity
direction alternating at each successive beam. For example, if the 5-member
system was extended by an additional member on either side of the system (to
give a 7 member system), the left-most member of the Type I polarity sequence
would have a positive polarity, and similarly, the right-most member would
have a negative polarity. The same logic extends to the other four sequences.
Each polarity sequence is indicative of two critical load arrangements that
maximise the positive or negative internal member forces respectively. The
maximum positive load arrangement for Type I is also equal to the maximum
negative load arrangement for Type IV, since these sequences are polar
opposites of each other, which is also true for the Type II and Type IV
sequences. Consequently, these 5 polarity zones correspond to 6 unique
critical load arrangements for a given beam, namely positive Type I, II and
III along with their (negative) polar opposites. The only exceptions occur for
the beams at either end of the spans, named end-span beams, in which the Type
I and Type IV sequences collapse into the Type III sequence (or its polar
opposite) at the left end, and similarly for the Type II and Type V sequences
at the right end, resulting in four unique load arrangements for end-span
beams. Whilst each non-end-span beam has 6 unique critical load arrangements,
it does not mean that the beam system has $6m$ unique load arrangements ($m$
is the number of members in the beam system). This is because, as shown in
Figure , the maximum positive Type V load arrangement for one beam is
identical to the maximum positive Type I load arrangement of the beam
immediately adjacent to (the right of) it. A similar overlap exists between
Type II and Type IV sequences, and the two Type III load arrangements (for
maximum and negative internal forces) are identical for all beams. Through a
process of elimination, it is possible to simplify the actual total number of
potential critical load arrangements to $p_{\mathrm{flex}}=2m$. This set will
be termed the flexural load arrangements set $\mathbf{J_{\mathrm{flex}}}$, and
can be evaluated using Algorithm LABEL:algo:flexural_set provided in
LABEL:app:algorithm_1, with an example output for a $m=5$ system shown in
Figure , grouped in alternating and adjacently loaded arrangements.

[width=height=keepaspectratio]figure_5_polarity_sequences_simplification.pdf

Figure : Polarity sequences are identical for adjacently lying beams
(highlighted in red) for Type I and Type V sequences as shown by Figure a) and
Figure c), as well as Type II and Type IV sequences, as shown by Figure b) and
Figure d).

[width=height=keepaspectratio]figure_6_load_arrangements_example.pdf

Figure : The critical load arrangements set $\mathbf{J_{\mathrm{flex}}}$ of
size $p=2m$ for a 5-member continuous beam system ($p=10$) grouped in
alternating and adjacently loaded arrangements.

The load arrangement set $\mathbf{J_{\mathrm{flex}}}$ of size
$p_{\mathrm{flex}}=2m$ identified here is a literal exponential improvement to
the brute-force approach of analysing and designing against $p=2^{m}$ load
arrangements and for evaluating the influence zone with Equation
LABEL:eq:global_formulation. It needs to be shown though that all critical
load arrangements $\mathbf{J_{\mathrm{crit}}}$ fall within
$\mathbf{J_{\mathrm{flex}}}$ (i.e.
$\mathbf{J_{\mathrm{crit}}}\in\mathbf{J_{\mathrm{flex}}}$).

##  Shear beams and the impact on critical load arrangements

To check if the load arrangement set $\mathbf{J_{\mathrm{flex}}}$ contains all
critical pattern loads, various continuous beam systems of size up to $m=10$
(to facilitate computational feasibility) were numerically generated, using
randomly distributed UDLs and span lengths based on the design constraints
identified in Table LABEL:table:design_sets. By calculating all $2^{m}$ load
arrangements, evaluating the utilisation ratio (based on the moment and shear
force combinations) of each, the load arrangement that caused the worst-case
utilisation ratio could be identified. This set of critical load arrangements
was then compared against the $\mathbf{J_{\mathrm{flex}}}$ set identified by
Algorithm LABEL:algo:flexural_set provided in LABEL:app:algorithm_1. Although
$\mathbf{J_{\mathrm{flex}}}$ tended to cause the critical utilisation ratio in
the majority of cases, there were instances were other load arrangements not
within $\mathbf{J_{\mathrm{flex}}}$ controlled the design. This unexpected
behaviour occurred in cases where short spanning, deep beams were included.
Analysing these special cases in detail indicated that the exact conditions
under which the previously unidentified load arrangements occur were generally
related to the following $L_{\mathrm{shear}}$ span limit quantified by

$\sqrt{\frac{6EI_{yy}}{GA_{z}}}<L_{\mathrm{shear}}$ (1)

where $E$ and $G$ are the Young’s and shear modulus of the material
respectively, and $I_{yy}$ and $A_{z}$ are the major second moment of area and
shear area of the prismatic beam, respectively. Although the
$L_{\mathrm{shear}}$ span limit appears to be related to shear beams, this is
the first time that shear beams have been reported in literature to cause
novel critical load arrangements. As shown in Figure , shear beams appear to
flip the polarity of the immediately adjacent member when measured outwardly
from a given IL location , with all subsequent members alternating the
polarity direction as before.

[width=height=keepaspectratio]figure_7_flex_vs_shear_sequence.pdf

Figure : A schematic demonstrating the impact of a shear beam (highlighted in
yellow) on a standard polarity sequence of a continuous beam system when spans
shorter than the shear span limit $L_{\mathrm{shear}}$ (as identified by
Equation 1) occur. Note the flipped polarity directions of the members on the
right-hand side of the system.

When shear beams (as defined by the $L_{\mathrm{shear}}$ limit) occur, they
introduce new critical load arrangements not found within
$\mathbf{J_{\mathrm{flex}}}$. The increase in terms of the final utilisation
factor of the beams was typically in the range of 4-5%, although larger
increases were also observed. Whilst a thorough analysis of the increase in
utilisation ratio caused by these newly identified load arrangements would be
of interest, it falls outside the scope of this study. Instead, an algorithm
will be presented capable of identifying these new load arrangements a priori,
which is the main objective of this investigation as explained in section
LABEL:sec:key_challenge.

##  Determining shear beam induced critical load arrangements a priori

The principal issue when evaluating the shear beam induced critical load
arrangements a priori, hereafter referred to as the shear load arrangements
$\mathbf{J_{\mathrm{shear}}}$, is the fact that the final material and cross-
sectional properties to evaluate the $L_{\mathrm{shear}}$ limit in Equation 1
are not known until the beam is designed. This creates a causality dilemma and
hence needs to be addressed. In clear opposition to the
$\mathbf{J}_{\mathrm{flex}}$ set, which does not depend on the continuous beam
system properties, the shear load arrangements cannot be established in
universum without some system knowledge. However, by taking advantage of the
design constraints set by Equation LABEL:eq:design_constraints, one can
identify a priori what members are potentially susceptible to cause shear load
arrangements by re-writing Equation 1 as:

$\sqrt{6\left(\frac{E}{G}\right)_{\mathrm{max}}\left(\frac{I_{yy}}{A_{z}}\right)_{\mathrm{max}}}<L_{\mathrm{shear,max}}$
(2)

The above equation groups the maximum material and cross-sectional property
ratios together. By limited the design space to S355 steel and UKB section
sizes as specified in Section LABEL:sec:design_constraints, the maximum
material ratio ($(E/G)_{\mathrm{max}}=2.600$) and cross-sectional property
ratio ($(I_{yy}/A_{z})_{\mathrm{max}}=0.397m^{2}$) can be evaluated.
Consequently, beams shorter than the shear span limit are susceptible to cause
shear load arrangements (in this case $L_{\mathrm{shear,max}}=2.49$m). In
identifying these susceptible members a priori, it is possible to evaluate the
shear load arrangements using Algorithm 2 provided in LABEL:app:algorithm_2.
Algorithm 2 transforms the flexural load arrangement from set
$\mathbf{J_{\mathrm{flex}}}$ based on a list of susceptible shear beams
identified by Equation 2. This is achieved by flipping the on/off activation
factor of the load arrangement if a shear beam is encountered whilst
travelling outwardly in both the left (-1) and right (1) direction from a
start beam index. This operation transforms the flexural load arrangement
based on the behaviour identified visually in Figure , and needs to check four
individual case conditions to account for continuous beam systems that have
multiple, potentially adjacently lying, shear beams. Since every beam system
is of size m, the time complexity of a single pass of Algorithm
LABEL:algo:shear_set is $O(m)$. However, since every flexural load arrangement
($2m$), and every combination of $n$ potential shear beams ($2^{n}-1$
combinations, as the zero set is already considered in
$\mathbf{J_{\mathrm{flex}}}$ by default), and every possible start-index ($m$)
needs to be computed, the time complexity to evaluate the shear set
$\mathbf{J_{\mathrm{shear}}}$ would be $O(m^{3}\,2^{n})$. It should be noted
that this process is computationally expensive. It was observed that passing
every possible start index generated either duplicate shear load arrangements,
or occasionally existing flexural load arrangements. For example, for a given
singular potential shear beam location, the algorithm would result in the same
transformed shear load arrangement for all start-indices starting on the left
and right hand-side of that susceptible shear beam location. Similarly, the
two alternating arrangements from $\mathbf{J_{\mathrm{flex}}}$ would result in
an already existing adjacent arrangement from $\mathbf{J_{\mathrm{flex}}}$ if
only a singular susceptible shear beam exists. Using such logic, it is
sufficient to pass only adjacent arrangements from
$\mathbf{J_{\mathrm{flex}}}$ along with the left-hand (or right-hand) index of
the adjacently loaded spans as the start index for Algorithm
LABEL:algo:shear_set to yield an effective set of potential shear load
arrangements. By not having to evaluate Algorithm LABEL:algo:shear_set for
every possible start index of each load arrangement, the computational
complexity reduces to $O(m^{2}\,2^{n})$. From this, it also follows that since
the alternating load arrangement is never transformed (which leaves only
$2(m-1)$ load arrangements to be passed to the algorithm) and since $2^{n}-1$
possible shear beam combinations can exist, the maximum number of unique
critical shear load arrangements should be of size
$p_{\mathrm{shear}}=2(m-1)(2^{n}-1)$.

##  Validating flexural and shear load arrangement algorithms

A design data set consisting of $32$ UDL and $32$ span values sampled from a
random uniform distribution for a $m=10$ beam system was generated based on
the high-variation design scenario identified in Section
LABEL:sec:design_constraints. Significantly higher variable UDLs
($Q_{k,i}\in[200\text{\,kN/m},400\text{\,kN/m}]$) were applied to increase the
likelihood of deep beams and thereby critical shear load arrangements,
allowing the performance of the algorithm to be stress-tested. This resulted
in $10\times 32\times 32=10240$ individual beam design examples, for which the
critical load arrangement $J_{\mathrm{crit}}$ could be identified. The results
of this validation exercise are illustrated in Figure , which plots the
critical load arrangement index for each design beam example. Every load
arrangement index corresponds to a unique load arrangement out of the naive
set $\mathbf{J}_{\mathrm{naive}}$ of size $p_{\mathrm{naive}}=2^{m}=1024$. The
set $\mathbf{J}_{\mathrm{naive}}$ was ordered so that the load arrangements
for set $\mathbf{J}_{\mathrm{flex}}$ are first, followed by those of set
$\mathbf{J}_{\mathrm{shear}}$, and subsequently all others. The design
examples themselves were sorted twice: first in ascending number of shear beam
occurrences, and subsequently in ascending load arrangement indices. This
results in the gradual increase of the $J_{\mathrm{crit}}$ indices as seen in
Figure .

[width=height=keepaspectratio]figure_8_algorithm_validation.pdf

Figure : Load arrangement index for each design beam example ordered in
increasing number of shear beam occurrences and critical load arrangement
indices. This confirms visually that the critical load arrangement
$J_{\mathrm{crit}}$ for each design beam example from the generated data set
falls within either $\mathbf{J_{\mathrm{flex}}}$ or
$\mathbf{J_{\mathrm{shear}}}$ and are significantly smaller than
$\mathbf{J_{\mathrm{naive}}}$. Figure b) is an enlarged view of Figure a).

Figure  sheds insight on a number of important points. The first is that the
critical load arrangement $J_{\mathrm{crit}}$ for every single beam example
from the 10240 data set occurred within the $\mathbf{J_{\mathrm{flex}}}$ or
$\mathbf{J_{\mathrm{shear}}}$ sets, validating the qualitative analysis based
on the polarity zones and sequences identified previously. This also
emphasises the validity of Algorithm LABEL:algo:flexural_set and
LABEL:algo:shear_set. Furthermore, the set size predictions $p_{flex}=2m$ and
$p_{shear}=2(m-1)(2^{n}-1)$ are also confirmed. For the $m=10$ member system
designed here, $p_{flex}=20$, and depending on the number of shear beam
occurrences of each system, which varied from $n=\\{0,1,2,3,4\\}$, the number
of shear load arrangements varied from $p_{shear}=\\{0,18,54,126,270\\}$. This
corresponded to $p_{total}=\\{20,38,74,146,290\\}$ respectively, as indicated
by the $y$-axis of Figure  b). Figure  a) also emphasises how much smaller
sets $\mathbf{J_{\mathrm{flex}}}$ and $\mathbf{J_{\mathrm{shear}}}$ are in
comparison to $\mathbf{J_{\mathrm{naive}}}$. This will greatly reduce the
number of load-arrangements that need to be analysed, reducing the
computational cost of both optimally designing the continuous beam system and
evaluating the influence zone for system lengths of $m>10$. Further insights
generated by Figure  are discussed in Section LABEL:sec:discussion.

##  Summary of critical load arrangements of continuous beams

By adding the set of critical flexural and shear load arrangements together,
it is possible to explicitly define a priori the set of critical load
arrangements for any continuous beam system under defined design constraints.
For the purpose of the influence zone concept and Equation
LABEL:eq:global_formulation, it is assumed that:
$\,\mathbf{J}\rightarrow\mathbf{J_{\mathrm{crit}}}\in\mathbf{J_{\mathrm{flex}}}\cup\mathbf{J_{\mathrm{shear}}}$.
The results from this systematic critical load arrangement investigation are
summarised in Table .

0.8c cc Set Set Size Algorithm Complexity
Critical load arrangements per internal beam 6 $O(1)$
Critical load arrangements per end-span beam 4 $O(1)$
$\mathbf{J_{\mathrm{flex}}}$ \- Critical flexural arrangements per beam system
$2m$ $O(m)$
$\mathbf{J_{\mathrm{shear}}}$ \- Critical shear arrangements per beam system
$2(m-1)(2^{n}-1)$ $O(m^{2}\,2^{n})$
$\mathbf{J_{\mathrm{naive}}}$ \- Naive load arrangements $2^{m}$ $O(2^{m})$

Table : Load arrangements set summary for $m$ dimensional beam systems
containing $n$ shear beams with associated algorithm complexities.