Spontaneous critical layer formation and robustness beneath rotational waves

TL;DR

This paper investigates the formation and robustness of critical layers beneath rotational waves with constant vorticity, revealing how Kelvin cat eye structures adapt to bottom topography and spontaneously form under certain conditions.

Contribution

It introduces a non-stationary framework for analyzing rotational waves, demonstrating spontaneous Kelvin cat eye formation and their robustness over complex bottom topography.

Findings

Kelvin cat eye structures adjust dynamically over topography.

Spontaneous formation of Kelvin cat eye occurs under pressure forcing.

Critical layers exhibit a new form with connected stagnation points.

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.

Equations

$$\begin{array}{rl}& u_t+u{u}_x+v{u}_y=-\frac{p_x}{\rho },\text{mbox}forh(x+U_{0}t),<y<\tilde{\zeta }(x,t),\\ & v_t+u{v}_x+v{v}_y=-\frac{p_y}{\rho }-g,\text{mbox}forh(x+U_{0}t),<y<\tilde{\zeta }(x,t),\\ & u_x+v_y=0,\text{mbox}forh(x+U_{0}t)<y<\tilde{\zeta }(x,t),\\ & p=P(x+U_{0}t),\text{mbox}aty=\tilde{\zeta }(x,t),\\ & v={\tilde{\zeta }}_t+u{\tilde{\zeta }}_x,\text{mbox}aty=\tilde{\zeta }(x,t),\\ & v=U_0{h}_x+u{h}_x,\text{mbox}aty=h(x+U_{0}t).\end{array}$$

$$\begin{array}{r}(u,v)=\nabla \tilde{\varphi }+(ay,0),\end{array}$$

$$\begin{array}{rl}& \Delta \tilde{\varphi }=0,\text{mbox}for-h_0+h(x+U_{0}t)<y<\tilde{\zeta }(x,t),\\ & (U_0-a{h}_0)h_x+ah{h}_x+{\tilde{\varphi }}_x{h}_x={\tilde{\varphi }}_y,\text{mbox}aty=-h_0+h(x+U_{0}t),\\ & {\tilde{\zeta }}_t+(a\tilde{\zeta }+{\tilde{\varphi }}_x){\tilde{\zeta }}_x-{\tilde{\varphi }}_y=0,\text{mbox}aty=\tilde{\zeta }(x,t),\\ & {\tilde{\varphi }}_t+\frac{1}{2}({\tilde{\varphi }}_x^2+{\tilde{\varphi }}_y^2)+a\tilde{\zeta }{\tilde{\varphi }}_x+\tilde{\zeta }-a\tilde{\psi }=-\frac{P(x+U_{0}t)}{\rho },\text{mbox}aty=\tilde{\zeta }(x,t),\end{array}$$

$$\begin{array}{rl}& \Delta \overline{\varphi }=0\text{mbox}for-h_0+h(x)<y<\overline{\zeta }(x,t),\\ & (U_0-a{h}_0)h_x+ah{h}_x+{\overline{\varphi }}_x{h}_x={\overline{\varphi }}_y\text{mbox}aty=-h_0+h(x),\\ & {\overline{\zeta }}_t+(U_0+a\overline{\zeta }+{\overline{\varphi }}_x){\overline{\zeta }}_x-{\overline{\varphi }}_y=0\text{mbox}aty=\overline{\zeta }(x,t),\\ & {\overline{\varphi }}_t+\frac{1}{2}({\overline{\varphi }}_x^2+{\overline{\varphi }}_y^2)+(U_0+a\overline{\zeta }){\overline{\varphi }}_x+\overline{\zeta }-a\overline{\psi }=-\frac{P(x)}{\rho }\text{mbox}aty=\overline{\zeta }(x,t).\end{array}$$

$$\begin{array}{rl}& \Delta \overline{\varphi }=0,\text{mbox}for-1+h(x)<y<\overline{\zeta }(x,t),\\ & (F-\Omega )h_x+\Omega h{h}_x+{\overline{\varphi }}_x{h}_x={\overline{\varphi }}_y,\text{mbox}aty=-1+h(x),\\ & {\overline{\zeta }}_t+(F+\Omega \overline{\zeta }+{\overline{\varphi }}_x){\overline{\zeta }}_x-{\overline{\varphi }}_y=0,\text{mbox}aty=\overline{\zeta }(x,t),\\ & {\overline{\varphi }}_t+\frac{1}{2}({\overline{\varphi }}_x^2+{\overline{\varphi }}_y^2)+(F+\Omega \overline{\zeta }){\overline{\varphi }}_x+\overline{\zeta }-\Omega \overline{\psi }=-P(x),\text{mbox}aty=\overline{\zeta }(x,t),\end{array}$$

$$\begin{array}{rl}& \Delta \overline{\varphi }=0,\text{mbox}for-1+h(x)<y<0,\\ & (F-\Omega )h_x+\Omega h{h}_x+{\overline{\varphi }}_x{h}_x={\overline{\varphi }}_y,\text{mbox}aty=-1+h(x),\\ & {\overline{\zeta }}_t+F{\overline{\varphi }}_x={\overline{\varphi }}_y,\text{mbox}aty=0,\\ & {\overline{\varphi }}_t+F{\overline{\varphi }}_x+\overline{\zeta }-\Omega \overline{\psi }=-P(x),\text{mbox}aty=0.\end{array}$$

$$\begin{array}{rl}& \frac{dx}{dt}={\overline{\varphi }}_x+\Omega y+F,\\ & \frac{dy}{dt}={\overline{\varphi }}_y,\\ & x(0)=x_0,y(0)=y_0.\end{array}$$

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

$$\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}$$

$$\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}$$

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

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

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

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

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

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

$${\mathbf{F}}_k[g(\xi )]=\hat{g}(k)=\frac{1}{L}{\int }_{-L/2}^{L/2}g(\xi )e^{-ik\xi }d\xi ,$$

$${\mathbf{F}}^{-1}[\hat{g}(k)](\xi )=g(\xi )=\sum _{j=-\infty }^{\infty }\hat{g}(k)e^{ik\xi }.$$

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

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

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

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

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

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

$$\begin{array}{rl}& {\varphi }_{\xi \xi }+{\varphi }_{\eta \eta }=0,\text{mbox}in-D<\eta <0,\\ & \varphi =\mathbf{\Phi }(\xi ,t),\text{mbox}at\eta =0,\\ & {\varphi }_{\eta }=(F-\Omega )H_{\xi }(\xi )+\Omega H{H}_{\xi }(\xi ),\text{mbox}at\eta =-D.\end{array}$$

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

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

$$\widehat{{\psi }_{\eta }}(k=0,\eta ,t)=\widehat{{\varphi }_{\xi }}(k=0,\eta ,t)=\frac{1}{L}{\int }_{-L/2}^{L/2}{\varphi }_{\xi }(\xi ,\eta ,t)d\xi =0.$$

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

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

Full text

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

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

[TABLE]

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

[TABLE]

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

[TABLE]

[TABLE]

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

[TABLE]

From (13) we obtain

[TABLE]

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

[TABLE]

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

[TABLE]

with

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

which is rewritten as

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

which for small wave elevations becomes

[TABLE]

In the canonical domain the kinematic condition and Bernoulli law are

[TABLE]

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

[TABLE]

where

[TABLE]

The Fourier operators are:

[TABLE]

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

[TABLE]

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

[TABLE]

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

In the moving frame we have that

[TABLE]

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

[TABLE]

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.

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

[TABLE]

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.

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

[TABLE]

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.

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$.

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.

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$.

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.

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.

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

[TABLE]

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.

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.