Nonlinear modeling of wave-topography interactions, shear instabilities and shear induced wave breaking using vortex method

TL;DR

This paper develops a vortex method-based framework to simulate fully nonlinear stratified shear instabilities and wave interactions, capturing key features like wave breaking and cusp formation beyond linear theory.

Contribution

It extends the classical vortex method with multiple interfaces and background shear, enabling nonlinear simulations of complex shear instability and wave-topography interaction problems.

Findings

Successfully modeled wave breaking and cusp formation.

Captured nonlinear features consistent with experiments.

Simulated diverse phenomena like Holmboe and Taylor-Caulfield instabilities.

Abstract

Theoretical studies on linear shear instabilities often use simple velocity and density profiles (e.g. constant, piecewise) for obtaining good qualitative and quantitative predictions of the initial disturbances. Furthermore, such simple profiles provide a minimal model for obtaining a mechanistic understanding of otherwise elusive shear instabilities. However, except a few specific cases, the efficacy of simple profiles has remained limited to the linear stability paradigm. In this work we have proposed a general framework that can simulate the fully nonlinear evolution of a variety of stratified shear instabilities as well as wave-wave and wave-topography interaction problems having simple profiles. To this effect, we have modified the classical vortex method by extending the Birkhoff-Rott equation to multiple interfaces, and furthermore, have incorporated background shear across a density interface. The latter is more subtle, and originates from the understanding that Bernoulli's equation is not just limited to irrotational flows, but can be modified to make it applicable for piecewise velocity profiles. We have solved diverse problems that can be essentially reduced to the multiple interacting interfaces paradigm, e.g. spilling and plunging breakers, stratified shear instabilities like Holmboe and Taylor-Caulfield, jet flows, and even wave-topography interaction problem like Bragg resonance. Free-slip boundary being a vortex sheet, its effect can also be effectively captured using vortex method. We found that the minimal models capture key nonlinear features, e.g. wave breaking features like cusp formation and roll-ups, which are observed in experiments and/or extensive simulations with smooth, realistic profiles.

Equations

$$\left[\frac{\partial {\varphi }_0^{\prime }}{\partial t}+\frac{1}{2}\left(u_0^{{}^{\prime }2}+w_0^{{}^{\prime }2}\right)-{\Omega }_0{\psi }_0^{\prime }+g{\eta }_1\right]+\frac{P_0^{\prime }}{{\rho }_0}=0,$$

$$\left[\frac{\partial {\varphi }_1^{\prime }}{\partial t}+\frac{1}{2}\left(u_1^{{}^{\prime }2}+w_1^{{}^{\prime }2}\right)-{\Omega }_1{\psi }_1^{\prime }+g{\eta }_1\right]+\frac{P_1^{\prime }}{{\rho }_1}=0.$$

$$P_0^{\prime }=P_1^{\prime }+\sigma \kappa ,$$

$$\kappa =-\tilde{\nabla }\cdot \hat{\mathbf{n}},$$

$$\left[\frac{\partial {\varphi }_0^{\prime }}{\partial t}+\frac{1}{2}\left(u_0^{{}^{\prime }2}+w_0^{{}^{\prime }2}\right)-{\Omega }_0{\psi }_0^{\prime }\right]-\left[\frac{\partial {\varphi }_1^{\prime }}{\partial t}+\frac{1}{2}\left(u_1^{{}^{\prime }2}+w_1^{{}^{\prime }2}\right)-{\Omega }_1{\psi }_1^{\prime }\right]+\frac{P_0^{\prime }}{{\rho }_0}-\frac{P_1^{\prime }}{{\rho }_1}=0.$$

$$\begin{gathered} \left[\frac{\partial {\varphi }_0^{\prime }}{\partial t}+\frac{1}{2}\left(u_0^{{}^{\prime }2}+w_0^{{}^{\prime }2}\right)-{\Omega }_0{\psi }_0^{\prime }+g{\eta }_1\right]+ \\ \left[\frac{\partial {\varphi }_1^{\prime }}{\partial t}+\frac{1}{2}\left(u_1^{{}^{\prime }2}+w_1^{{}^{\prime }2}\right)-{\Omega }_1{\psi }_1^{\prime }+g{\eta }_1\right]+\frac{P_0^{\prime }}{{\rho }_0}+\frac{P_1^{\prime }}{{\rho }_1}=0. \end{gathered}$$

$$\gamma \equiv \left({\mathbf{u}}_0^{\prime }-{\mathbf{u}}_1^{\prime }\right)\cdot \hat{\mathbf{s}},$$

$$\left\{\left[\frac{D_0{\text{u}}_0^{\prime }}{Dt}-{\Omega }_0\nabla {\psi }_0^{\prime }\right]-\left[\frac{D_1{\text{u}}_1^{\prime }}{Dt}-{\Omega }_1\nabla {\psi }_1^{\prime }\right]+\frac{\nabla P_0^{\prime }}{{\rho }_0}-\frac{\nabla P_1^{\prime }}{{\rho }_1}\right\}\cdot \hat{\mathbf{s}}=0.$$

$$\left\{\left[\frac{D_0{\text{u}}_0^{\prime }}{Dt}-{\Omega }_0\nabla {\psi }_0^{\prime }+g\hat{\mathbf{k}}\right]+\left[\frac{D_1{\text{u}}_1^{\prime }}{Dt}-{\Omega }_1\nabla {\psi }_1^{\prime }+g\hat{\mathbf{k}}\right]+\frac{\nabla P_0^{\prime }}{{\rho }_0}+\frac{\nabla P_1^{\prime }}{{\rho }_1}\right\}\cdot \hat{\mathbf{s}}=0.$$

$$\begin{gathered} \frac{d\gamma }{dt}=2At\frac{d{\text{U}}^{\prime }}{dt}\cdot \hat{\mathbf{s}}+\frac{\beta +At}{4}\left(\frac{\partial {\gamma }^2}{\partial s}\right)-(1+\beta At)\gamma \frac{\partial {\text{U}}^{\prime }}{\partial s}\cdot \hat{\mathbf{s}}+2Atg\frac{\partial {\eta }_1}{\partial s}-\frac{\sigma }{{\rho }_{av}}\frac{\partial \kappa }{\partial s} \\ -\left(At-1\right){\Omega }_0\nabla {\psi }_0^{\prime }\cdot \hat{\mathbf{s}}-\left(At+1\right){\Omega }_1\nabla {\psi }_1^{\prime }\cdot \hat{\mathbf{s}}. \end{gathered}$$

$$\mathbf{q}\equiv {\mathbf{U}}^{\prime }+\frac{1}{2}\gamma \beta \hat{\mathbf{s}}.$$

$${\text{U}}^{\prime }\left(s,t\right)=-\frac{1}{2\pi }P.V.{\int }_{-\infty }^{+\infty }\frac{\hat{\mathbf{j}}\times \left(\text{X}\left(s,t\right)-\text{X}\left(\tilde{s},t\right)\right)}{|\text{X}(s,t)-\text{X}(\tilde{s},t){|}^2}\gamma (\tilde{s},t)d\tilde{s}.$$

$$u^{\prime}-iw^{\prime}=\frac{i}{2\lambda}\int_{0}^{\lambda}\tilde{\gamma}\cot\bigg{[}\frac{\pi(\chi-\tilde{\chi})}{\lambda}\bigg{]}d\tilde{s}.$$

$$u^{\prime}_{l}-iw^{\prime}_{l}=\sum_{k=1}^{M}\frac{i}{2\lambda}\int_{0}^{\lambda}\tilde{\gamma}_{k}\cot\bigg{[}\frac{\pi(\chi_{l}-\tilde{\chi}_{k})}{\lambda}\bigg{]}\,d\tilde{s}_{k}.$$

$$\frac{d\mathbf{X}}{dt}=\mathbf{q}+U_B\hat{\mathbf{i}},$$

$$u_{il}^{\prime }=\frac{1}{2\lambda }\sum _{k=1}^M\sum _{j=1,j\ne i}^N{\Gamma }_{jk}\frac{\sinh [\alpha (z_{il}-z_{jk})]}{\cosh [\alpha (z_{il}-z_{jk})]-\cos [\alpha (z_{il}-z_{jk})]+{\delta }^2},$$

$$w_{il}^{\prime }=-\frac{1}{2\lambda }\sum _{k=1}^M\sum _{j=1,j\ne i}^N{\Gamma }_{jk}\frac{\sin [\alpha (x_{il}-x_{jk})]}{\cosh [\alpha (z_{il}-z_{jk})]-\cos [\alpha (z_{il}-z_{jk})]+{\delta }^2}.$$

$${\Gamma }_{jk}={\gamma }_{jk}\Delta s_{jk},$$

$$\Delta s_{jk}=\frac{1}{2}\sqrt{(x_{j+1,k}-x_{j-1,k}{)}^2+(z_{j+1,k}-z_{j-1,k}{)}^2}.$$

$${\gamma }_b=\left({\mathbf{u}}_0^{\prime }-{\mathbf{u}}_1^{\prime }\right)\cdot \hat{\text{s}}={\text{u}}_0^{\prime }\cdot \hat{\text{s}},$$

$${\mathbf{U}}^{\prime }=\frac{{\mathbf{u}}_0^{\prime }+{\mathbf{u}}_1^{\prime }}{2}=\frac{{\mathbf{u}}_0^{\prime }}{2}.$$

$${\gamma }_b=2{\mathbf{U}}^{\prime }\cdot \hat{\text{s}}.$$

$$u^{\prime}_{b}-iw^{\prime}_{b}=\frac{i}{2\lambda}\int_{0}^{\lambda}\tilde{\gamma}_{f}\cot\bigg{[}\frac{\pi(\chi_{b}-\tilde{\chi}_{f})}{\lambda}\bigg{]}d\tilde{s}_{f}+\frac{i}{2\lambda}\int_{0}^{\lambda}\tilde{\gamma}_{b}\cot\bigg{[}\frac{\pi(\chi_{b}-\tilde{\chi}_{b})}{\lambda}\bigg{]}d\tilde{s}_{b}.$$

$$\Omega(z)=\left\{\begin{array}[]{cc}\Omega_{0}&\quad h\leq z\\ \Omega_{1}&\quad-h\leq z<h\\ \Omega_{2}&\quad z<-h\end{array}\right.\qquad\qquad\bar{\rho}(z)=\left\{\begin{array}[]{cc}\rho_{0}&\quad h\leq z\\ \rho_{1}&\quad-h\leq z<h\\ \rho_{2}&\quad z<-h.\end{array}\right.$$

$$P_j^{\prime (1)}=-{\rho }_j\left[\frac{\partial {\varphi }_j^{\prime }}{\partial t}+U_1\frac{\partial {\varphi }_j^{\prime }}{\partial x}+g{\eta }_1-{\Omega }_j{\psi }_j^{\prime }\right],j=0,1,$$

$$P_j^{\prime (2)}=-{\rho }_j\left[\frac{\partial {\varphi }_j^{\prime }}{\partial t}+U_2\frac{\partial {\varphi }_j^{\prime }}{\partial x}+g{\eta }_2-{\Omega }_j{\psi }_j^{\prime }\right],j=1,2.$$

$$P_0^{\prime (1)}=P_1^{\prime (1)}+\sigma \kappa ,$$

$$P_1^{\prime (2)}=P_2^{\prime (2)},$$

$${\nabla }^2{\varphi }_j^{\prime }=0.$$

$${\varphi }_j^{\prime }\left(x,z,t\right)={\tilde{\varphi }}_j\left(z\right)e^{i\left(\alpha x-\omega t\right)},$$

$${\eta }_j\left(x,t\right)={\tilde{\eta }}_j{e}^{i\left(\alpha x-\omega t\right)}.$$

$$\frac{d^2{\tilde{\varphi }}_j}{d{z}^2}-{\alpha }^2{\tilde{\varphi }}_j=0.$$

$${\tilde{\varphi }}_0(z)=C_1{e}^{-\alpha z},$$

Full text

Vortex methods, stratified shear instabilities, wave breaking, vortex sheets

††preprint: AIP/123-QED

I Introduction

Velocity shear, with or without density stratification, are ubiquitous in environmental and industrial flows. Examples of shear flows include pycnocline region in lakes and oceans, Antarctic circumpolar current, Stratospheric polar night jet, density currents, fuel-air mixing layer in internal combustion engines. Stratified shear flows usually exhibit region(s) of strong gradients in background velocity and density, and are often unstable, leading to various types of shear instabilities. When gradients are very sharp, it leads to the formation of interfaces. A density interface is marked by a sharp change in background density (e.g. air-sea interface), whereas a vorticity interface is produced due to a sharp change in background vorticity. A crucial feature of interfaces is that it can support progressive waves. A density interface can support two oppositely propagating gravity waves (e.g. the waves we observe on the ocean surface), whereas a vorticity interface supports a vorticity wave. Wave breaking is a hallmark nonlinear feature of shear instabilities, and is crucial from application point of view. While surface gravity wave breaking is pivotal in understanding air-sea coupling Melville (1996), internal gravity wave breaking at the pycnocline plays a major role in understanding ocean/lake mixing and the subsequent biogeochemical processes Smyth and Moum (2012). In the atmosphere, Rossby wave (vorticity waves in a rotating frame) breaking is key to the understanding of ozone mixing, and the fate of the Antarctic ozone hole Guha et al. (2016).

Stratified shear instabilities being reasonably complicated, it is often helpful to theoretically simplify the background flow into piecewise linear and/or constant profiles. In the linear regime (i.e., when the wave-like perturbations are infinitesimally small), a fairly accurate qualitative and quantitative description of shear instabilities can be obtained from simplified piecewise profiles Drazin and Reid (2004). A pertinent question to ask is - do piecewise profiles also provide reasonably accurate description in the nonlinear regime as well? Since nonlinearity brings about additional complexities, it is not unreasonable to expect that the simplistic nature of piecewise background profiles may fail to capture the essential nonlinear features of shear instabilities emanating from continuous profiles. Fortunately the agreement is reasonably good for homogeneous flows (density is constant). A numerical technique known as contour dynamicsPozrikidis and Higdon (1985); Dritschel (1989); Pullin (1992); Guha, Rahmani, and Lawrence (2013), developed around forty years ago Deem and Zabusky (1978), captures the nonlinear evolution of shear instabilities emanating from piecewise linear velocity profiles (vorticity constant in each layer). The results obtained from contour dynamics simulations match well with that of shear instabilities ensuing from continuous shear layers. This breakthrough technique has been used to study a variety of homogeneous shear flow problems, including but not limited to geostrophic turbulence Scott and Dritschel (2012).

The success of the applicability of contour dynamics is however limited to homogeneous flows. The primary objective of this paper is to devise a generalized framework such that density stratified shear instabilities emanating from simple piecewise profiles can be extended to the nonlinear regime. In other words, the goal is to simulate the nonlinear evolution of density interface(s) sandwiched in a shear layer. It is important to note here that the numerical strategy behind the nonlinear evolution of stable or unstable density interfaces (in absence of background velocity shear) is already well established, and is efficiently and accurately implemented using vortex methods Saffman (1992). This technique capitalizes on the fact that density interfaces are basically vortex sheets, hence can be discretized into a series of point vortices. Point vortices interact with each other; each point vortex moves with the resultant velocity imposed on it by all other vortices at any given time. Thus the interface defined by the point vortices also evolves in time. A basic assumption here is that the flow is inviscid, and is therefore suitable in mimicking high Reynolds number flows occurring in geophysical settings. The computation is also numerically efficient since the equations of motion are solved only at the interfaces, thereby reducing the dimensionality of the problem by one. In the past, vortex methods have been successfully implemented to accurately simulate the dynamics of a single unstable interface, e.g. Rayleigh-Taylor instability Baker, Meiron, and Orszag (1980); Tryggvason (1988), Richtymer-Meshkov instability Sohn (2004), Kelvin-Helmholtz instability (KH) Sohn, Yoon, and Hwang (2010), as well as that of a stable interface, e.g. propagation of surface gravity waves Baker, Meiron, and Orszag (1982).

Before proceeding further we briefly summarize the advantages and shortcomings of the two above-mentioned techniques, namely contour dynamics and vortex method. While contour dynamics allows multiple vorticity interfaces, it is only applicable to homogeneous flows. Vortex method, however, can capture the evolution of a single density interface (hence can model stratified flows) provided there is no finite background shear. The method can only capture infinite shear (jump in background velocity) with or without density jump, for example, the classic KH instability set-up Sohn, Yoon, and Hwang (2010). This is a major drawback, especially for simulating stratified shear instabilities, since in such flows, one or more density interfaces are usually embedded in a shear layer. The most general scenario of piecewise density stratified shear layer is shown in figure 1, which even includes a bottom topography. The interface at $z=\eta_{2}$ is the infinite shear interface, and as already mentioned, only the evolution of this interface can be captured using vortex method. The fact that restricts vortex method to shearless (or irrotational) background flow is the unsteady Bernoulli’s equation, which governs the evolution of a vortex sheet. The fundamental assumption behind the derivation of Bernoulli’s equation is that * the flow must be irrotational*. Hence our aim in this paper is threefold :

(i) To find the evolution equation of a vortex sheet in presence of shear. In other words, we need to generalize the well known unsteady Bernoulli’s equation, which is only applicable to irrotational background velocity profile, to profiles that are piecewise linear.

(ii) To find the interaction between multiple vortex sheets. This is because shear instabilities arise from the interaction between multiple interfaces in presence of background shear (rather, waves present at these interfaces)Caulfield (1994); Baines and Mitsudera (1994); Heifetz and Methven (2005); Guha and Lawrence (2014). A point worth mentioning is that free-slip boundary is also a vortex sheetBaker, Meiron, and Orszag (1982), hence an algorithm that captures the interaction between multiple vortex sheets would automatically include the interaction between a free interface and a solid boundary.

 (iii) Once the above two conditions are fulfilled, it remains to be shown whether interacting density interfaces embedded in a shear layer replicate the key nonlinear features of shear instabilities, at least qualitatively.

The organization of the paper is as follows. In section II we discuss the mathematical formulation and the numerical technique. To show the versatility of the technique in capturing wave nonlinearities, including breaking under very different scenarios, we broadly discuss three kinds of settings - flows with finite shear (sections III - IV), strongly sheared flows (section V), and flows without shear (section VI). In section III we capture surface wave breaking due to shear for both short waves (in the presence of surface tension) as well as long waves (where bottom topography becomes crucial). These are respectively known as the spilling and plunging breakers. In section IV we analyze the nonlinear evolution of two well known types of stratified shear instabilities that occur inside lakes and oceans - Holmboe instability and Taylor-Caulfield instability. In section V we discuss liquid sheets or jet flows that is susceptible to KH instability. Unlike previous studies, here KH instability occurs in more than one interfaces. In section VI we highlight the importance of bottom topography by switching off background velocity. Here we show a kind of wave-triad interactions, known as Bragg resonance, occurring between two oppositely propagating surface waves and the rippled bottom topography. The paper is summarized and concluded in section VII.

II Mathematical formulation and numerical procedure

We consider an inviscid, incompressible flow in the $x-z$ plane with piecewise linear background velocity profile $\bar{U}(z)$ (parallel to the streamwise direction $x$) and piecewise constant background density profile $\rho(z)$; see figure 1. Unless otherwise stated, the vertical axis $z$ points upwards, i.e.  opposite to the direction gravity. In this system there are $N$ material interfaces $z=\eta_{1},\,\eta_{2},\ldots\eta_{N}$, which are material curves across which at least one of the flow variables - background density, vorticity (defined as $\Omega=-d\bar{U}/dz$), or velocity, may be discontinuous. In figure 1 we have density discontinuity across all interfaces. Density stratification is always assumed to be stable, i.e. $\rho_{0}\leq\rho_{1}\leq\rho_{2}\leq\ldots\leq\rho_{N}$. The interface at $z=\eta_{2}$ is special; there is a discontinuity in $\bar{U}$, hence that interface is a vortex sheet by definition. Such an interface can become unstable and yields KH instability Drazin and Reid (2004). However, all other interfaces are also vortex sheets since discontinuity in either $\Omega$ or $\rho$ yields discontinuity in the perturbation velocity’s tangential component Caulfield (1994); Baines and Mitsudera (1994). Below we obtain the evolution equation of vortex sheet strength for $\Omega$ and/or $\rho$ discontinuity under a very generalized setting.

II.1 Dynamics of an interface in presence of shear - the “shear-modified” Bernoulli’s equation

For simplicity we focus on a single interface across which both $\Omega$ and $\rho$ are discontinuous (e.g. $z=\eta_{1}$ in figure 1). Since $\bar{U}$ is assumed continuous and piecewise linear, $\Omega$ is piecewise constant. Under such circumstances the perturbed flow is irrotational in each layer111Here we emphasize that for irrotationality of the perturbed flow, it is not mandatory for $\Omega$ to be identically zero in each layer. $\Omega$ being piecewise constant also leads to irrotational perturbed flow in each layer.. Hence we can integrate the Euler’s equation in each layer, just above and below the interface, to obtain

[TABLE]

Here $g$ is the acceleration due to gravity and $\rho_{i}$, $\phi^{\prime}_{i}$, $u^{\prime}_{i}$, $w^{\prime}_{i}$, $\Omega_{i}$, $\psi^{\prime}_{i}$ and $P^{\prime}_{i}$ are respectively the background density, perturbation velocity potential (defined since perturbation velocity is irrotational), perturbation horizontal velocity, perturbation vertical velocity (same as total vertical velocity), background vorticity, perturbation streamfunction and perturbation pressure in the $i$-th layer, where $i=\,0$ or $1$. In absence of background shear, (1) and (2) basically represent unsteady Bernoulli’s equation. The extra term, $-\Omega_{i}\psi^{\prime}_{i}$, appears due to the presence of background shear. To the best of our knowledge, (1)-(2) is the most general description of the “shear-modified” Bernoulli’s equation. A special case of these equations was formulated in Simmen (1984), which considered the effect of constant shear on a free surface gravity wave (air considered as a passive fluid of zero density) in a steady frame. The formulation of Simmen has been used in many surface gravity wave problems that includes linear shear current in the water region and a passive air above it Kishida and Sobey (1988); Constantin and Strauss (2004).

The dynamic boundary condition is usually obtained by subtracting (1) and (2) and equating the pressure, which has to be continuous across an interface in absence of surface tension. However, when surface tension is present, pressures above and below an interface are related by the Young-Laplace equation

[TABLE]

where $\sigma$ is the coefficient of surface tension and $\kappa$ is the curvature of the interface. The curvature is calculated by taking divergence of the vector normal to the interface:

[TABLE]

where $\tilde{\nabla}\equiv\nabla-\left(\mathbf{\hat{n}}\cdot\nabla\right)\nabla$. On subtracting (2) from (1) we get:

[TABLE]

Similarly, on adding (1) and (2) we get:

[TABLE]

Our objective is to find the time evolution of vortex sheet strength $\gamma$, which is defined as follows:

[TABLE]

where $\mathbf{u^{\prime}_{0}}=(u^{\prime}_{0},w^{\prime}_{0})$, $\mathbf{u^{\prime}_{1}}=(u^{\prime}_{1},w^{\prime}_{1})$ and $\mathbf{\hat{s}}$ is the unit tangent vector to the interface.

For converting the formulation to velocities instead of potentials, we operate $\nabla$ on (5), then project it along the tangential direction (by taking an inner product with $\mathbf{\hat{s}}$) to obtain

[TABLE]

Performing the same operation on (6) yields

[TABLE]

Here $D_{i}/Dt\equiv\partial/\partial t+\mathbf{u^{\prime}_{i}}\cdot\nabla$ is the (perturbation) material derivative in the $i^{\mathrm{th}}$ layer; $i=0,1$. Furthermore $\nabla=\mathbf{\hat{s}}\,\partial/\partial s$, where $s$ is the arc-length coordinate. Eliminating the pressure terms in (8) and (9) and using (3) we obtain the final evolution equation of vortex sheet strength:

[TABLE]

Here $\rho_{av}=(\rho_{0}+\rho_{1})/2$ is the average density, $At=(\rho_{1}-\rho_{0})/(\rho_{1}+\rho_{0})$ is the Atwood number and $\mathbf{U^{\prime}}=(\mathbf{u^{\prime}_{0}}+\mathbf{u^{\prime}_{1}})/2$ is the average perturbation velocity at the interface. The term $2Atg\partial\eta_{1}/\partial s$ denotes the baroclinic generation of vorticity. The operator $d/dt\equiv\partial/\partial t+\mathbf{q}\cdot\nabla$ denotes (perturbation) material derivative following the interface, where $\mathbf{q}$ is the interfacial perturbation velocity defined as follows:

[TABLE]

In the above equation, $\beta=[-1,1]$ is a weight parameter such that $\beta=1\,(-1)$ implies that the interface is moving with the upper (lower) fluid. We have chosen $\beta=0$ in this study, implying that the interface travels with the average velocity of the two layers.

The vortex sheet evolution equation, (10), is in its most general form. The last two terms of it denote the effect of piecewise linear background velocity (or piecewise continuous background vorticity). We would emphasize here that these terms were absent in the previous derivations of vortex sheet evolution; see (2.15) of Baker, Meiron, and Orszag (1982), (15) of Tryggvason (1988), or (9) of Sohn, Yoon, and Hwang (2010), since these authors chose the base flow to be irrotational. Furthermore, our derivation strategy has also been different - we have obtained (10) from unsteady Bernoulli’s equation modified by background shear; see (1)-(2).

II.2 Evolution of multiple interfaces - the modified Birkhoff-Rott equation

The interface variables can be parameterized using the arc-length coordinate $s$. The average perturbation velocity $\textbf{U}^{\prime}$ of an infinitely long interface, which is also a vortex sheet, can be evaluated using the Biot-Savart integral:

[TABLE]

Here X is the position vector representing an arbitrary point on the interface, and is given by $\textbf{X}=x(s,t)\mathbf{\hat{i}}+z(s,t)\mathbf{\hat{k}}$. P.V.  in front of the integral implies principal value (hereafter, we don’t put P.V. in front of Biot-Savart integrals, assuming it is expected), and $\mathbf{\hat{j}}$ is the unit vector perpendicular to the plane containing the vortex sheet.

The Biot-Savart law for periodic boundary condition (as is present in the interfaces we consider) reduces to the Birkhoff-Rott equation:

[TABLE]

Here $u^{\prime}$, $w^{\prime}$ are the horizontal and vertical components of $\mathbf{U^{\prime}}$, $\chi$ is the complex position of the interface: $\chi=x(s,t)+iz(s,t)$, $\lambda$ is the wavelength, and $\tilde{\gamma}\equiv\gamma(\tilde{s},t)$. In this paper we have extended the Birkhoff-Rott equation to multiple interfaces, which can be mathematically written as follows:

[TABLE]

Here we are considering a system with $M$ interfaces; the indices $k$ and $l$ represent the $k^{\mathrm{th}}$ and $l^{\mathrm{th}}$ interfaces respectively.

II.3 Numerical technique

An efficient and accurate numerical method for solving the Biot-Savart or the Birkhoff-Rott equation is by discretizing each interface into a series of point vortices. These vortices are Lagrangian markers, given $\gamma$ at a particular time $t$ an interface evolves following

[TABLE]

where $U_{B}$ is the background velocity. In general we have $U_{B}=\bar{U}(z)$, i.e. $U_{B}$ does not evolve in time. However, if contour dynamics is implemented (e.g.  the plunging breaker problem in III.2), then $U_{B}$ would evolve according to (32).

Since we are using Lagrangian formalism, q and $\mathbf{U^{\prime}}$ are both Lagrangian velocities, which are related by (11). The average perturbation velocity $\mathbf{U^{\prime}}$ can be obtained from the discretized version of (14), i.e.  Birkhoff-Rott equation for multiple interfaces:

[TABLE]

Here $u^{\prime}_{il}$ and $w^{\prime}_{il}$ respectively denote the horizontal and vertical components of the average perturbation velocity $\mathbf{U^{\prime}_{il}}$ induced by the $i^{\mathrm{th}}$ point vortex located at the $l^{\mathrm{th}}$ interface. The wavenumber $\alpha=2\pi/\lambda$; $\Gamma_{jk}$ denotes the circulation strength of the $j^{\mathrm{th}}$ point vortex located at the $k^{\mathrm{th}}$ interface. Similarly $x_{il}$ and $x_{jk}$ respectively denote the $x$ coordinates of the $i^{\mathrm{th}}$ and $j^{\mathrm{th}}$ point vortices located at the $l^{\mathrm{th}}$ and $k^{\mathrm{th}}$ interfaces. The same notation applies for the $z$ coordinates. The circulation strength of point vortices is given by

[TABLE]

where $\Delta s_{jk}$ is the arc-length about the $j^{\mathrm{th}}$ point vortex at the $k^{\mathrm{th}}$ interface, and is given by

[TABLE]

Point vortices being Lagrangian parcels can cluster in certain region(s) of the interface, which may lead to singularities. Such issues can be avoided by introducing a desingularization parameter $\delta$, known as the Krasny parameter Krasny (1986), in the denominators of (16) and (17).

Once each interface moves to its new location following (15) in a time interval $\Delta t$, their vortex sheet strengths are evolved in time according to (10). For solving these equations, we require initial conditions for vortex sheet strength and initial shape of the interfaces. Once the initial shape is provided, the initial velocities of the point vortices are found using (16)-(17). Since vortex strength and velocities are both coupled together we will need to solve for them iteratively. For this we have followed the method outlined in Sohn, Yoon, and Hwang (2010).

We have also followed Sohn, Yoon, and Hwang (2010) for handling resolution issues of the vortex sheets. Unless otherwise stated, in all our calculations we initially discretize one wavelength ($2\pi$ distance) of a vortex sheet by $257$ point vortices. As mentioned previously, with the progress of time, an interface can develop poor resolution in certain areas and clustering in others. At late times, point clustering and diverging occur irregularly along the interfaces. To deal with this issue we use an insertion-deletion scheme. We define a threshold arc-length $\Delta s_{\mathrm{thresh}}$, if arc-length at any location $i$ becomes less than the threshold, i.e. $\Delta s_{ik}<\Delta s_{\mathrm{thresh}}$, then we delete that vortex; vice-versa we insert a point vortex. The $x$ coordinate of the inserted vortex $\left(x_{p},y_{p}\right)$ is taken to be the average of the $x$ coordinates of $i^{\mathrm{th}}$ and $(i+1)^{\mathrm{th}}$ vortices: $x_{p}=\left(x_{i}+x_{i+1}\right)/2$; likewise we interpolate for the $z$ coordinate and the vortex strength $\gamma_{p}$. We have validated our code with the published results of Baker, Meiron, and Orszag (1980) and Sohn, Yoon, and Hwang (2010). Our code is found to accurately replicate the late time KH roll-up, as well as surface gravity wave breaking.

II.4 Modeling fixed boundary as a vortex sheet

Up to this point we have discussed the mathematical formulation of interaction between free interfaces, which are basically vortex sheets. By free interface we mean interfaces that can freely evolve. Examples of such interfaces include air-sea interface, the fresh water - salt water interface (pycnocline). However, a bottom topography, e.g. the sea bed, is also an interface. The difference is that it is fixed (cannot evolve freely in time). The key point to note here is that a fixed interface is also a vortex sheet, and therefore can be understood following similar mathematical principles discussed earlier. The reason behind a bottom topography acting like a vortex sheet can be understood as follows. Above the topography (which satisfies free-slip boundary condition owing to the fact that the fluid is inviscid), the fluid has a velocity (that is generated by the point vortices on the free interfaces), while below the topography the velocity has to be zero. Hence there is a velocity discontinuity across the topography, making it a vortex sheet. The prime differences between a free and a fixed interface are that for the latter (i) the normal velocity has to be zero (impenetrability condition), and (ii) there is no separate evolution equation for the vortex sheet strength. The vortex sheet strength of a fixed interface evolves with the vortex sheet strength of the free interfaces, hence the problem must be solved in an iterative fashion. Our formulation in this regard is similar, but simpler than the formulation used by Baker, Meiron, and Orszag (1982) for simulating breaking nonlinear surface gravity waves.

II.4.1 Mathematical formulation

Our system can be composed of multiple free interfaces and multiple fixed interfaces (i.e. boundaries). The evolution of the vortex sheet strength of free interfaces is governed by (10). The vortex strength of the fixed bottom is evaluated using the fact that the velocity below the interface is zero and the velocity above the interface is only present in the direction tangential to the bottom interface:

[TABLE]

where $\bf{u^{\prime}_{0}}$ is the perturbation velocity above the fixed interface, while $\bf{u^{\prime}_{1}}=\bf{0}$ is the perturbation velocity below it. Since the Birkhoff-Rott equation evaluates the average perturbation velocity at the interface we have

[TABLE]

Hence the vortex strength of the fixed interface becomes

[TABLE]

The principal vortex velocity $\textbf{U}^{\prime}$ is itself calculated from the vortex strengths of both free and fixed interfaces, and is given by (in complex representation):

[TABLE]

Here the subscripts $f$ and $b$ denote the free interface and the fixed bottom respectively. Since the equation for the vortex strength of the bottom interface is implicit in nature, we solve for it iteratively.

II.4.2 Numerical Technique

The numerical technique for solving the system of equations governing the interactions between free and fixed interfaces is similar to that discussed in section II.3. The system of governing equations (10), (14), (15), (20) and (21) represent an initial value problem and hence the knowledge of the initial shape of the interfaces and their vortex strengths is imperative. The only difference here is that at every time level we iterate until the acceleration of the free surface and vortex sheet strength of the bottom converge. The vortex sheet strength of the bottom is updated using equation (20).

III Flows with finite shear - surface wave breaking problems

III.1 Spilling breaker

Deep water gravity wave breaking is characterized as either spilling or plunging. In a plunging breaker the wave grows, steepens, forms a horizontal jet which later falls under the effect of gravity onto the front face of the crest of wave. In numerical simulations of a plunging breaker the wave is forced by an asymmetric pressure distribution. Longuet-Higgins and Cokelet (1976) did simulations of a plunging breaker using boundary integral methods. A spilling breaker is different in that the wave grows, becomes nonlinear, then a bulge forms at the crest of the wave, and small capillary waves originate at the toe of this bulge. A theoretical model of spilling breaker was proposed by Longuet-Higgins (1994), we have used this model for our simulations. In this section we present a detailed description of the initial conditions and then present the numerical simulation results. The initial condition derivation for the system is quite general in the sense that the air-water interface (located about $z=h$) has density jump, vorticity jump, as well as surface tension. A schematic of the system is shown in figure III.1. The interface (located about $z=-h$) has only a vorticity jump. The base vorticity and density profiles are given by

[TABLE]

For a spilling breaker $\Omega_{0}=\Omega_{2}=0$ and $\Omega_{1}$ is a constant base shear. Furthermore $\rho_{0}=0$ (density of air assumed to be zero) and $\rho_{2}=\rho_{1}$. However we will proceed without making these assumptions at this point so as to arrive at a general initial condition.

III.1.1 Initial conditions

Following (1)-(2), linearized Bernoulli’s equations (in the presence of shear) above and below the interfaces are given by

[TABLE]

Note here that the subscript $j=0,1,2$ refers to the different fluid layers in our three-layered system. Furthermore, $P_{j}^{\prime\,(k)}$ implies perturbation pressure of the $j^{\mathrm{th}}$ layer in the neighborhood of the $k^{\mathrm{th}}$ interface. Also $U_{1}\equiv\bar{U}(z=h)$ and $U_{2}\equiv\bar{U}(z=-h)$ are basically the background velocities at the respective interfaces. Applying normal stress balance at the interfaces, we obtain

[TABLE]

where $\sigma$ is the surface tension coefficient at the air-water interface. The perturbation velocity potentials satisfy the Laplace equation

[TABLE]

Now we assume normal mode forms for all the perturbation quantities:

[TABLE]

On substituting the normal mode solutions into the Laplace equation we get:

[TABLE]

This implies that the eigenfunctions of the potentials in the $z$-direction are of exponential form. Using the fact that far away from the interfaces these potentials must decay to zero (evanescent condition), we can eliminate a few terms. Finally we are left with four constants appearing in eigenfunctions of velocity potentials and the two interfacial displacements ($\tilde{\eta_{1}}$ and $\tilde{\eta_{2}}$):

[TABLE]

We need six equations to solve for this system. Hence we write four kinematic equations, just above and below each of the two interfaces, and two normal stress balance conditions across each of the interfaces.

[TABLE]

Finally we formulate the problem as a system of linear homogeneous equations $\mathcal{A}\mathcal{X}=0$. The matrix $\mathcal{A}$ in the most general form is given by

[TABLE]

where

[TABLE]

The variable vector $\mathcal{X}$ is given by

[TABLE]

On applying the condition for a non-trivial solution we get the dispersion relation for the system, $\mathcal{X}$ being the null vector of the matrix $\mathcal{A}$.

In case of modal instabilities $\omega$ obtained from the dispersion relation will be complex, and the null vector corresponding to that frequency will give us the phase locked configuration that results in exponentially growing modes. Phase locking is a state in which the relative phase between the waves does not change. This occurs when the intrinsic phase speed of the waves is counter to the base flow at the respective interfaces as well as counter to each other (counter-propagating waves)Heifetz and Methven (2005); Carpenter et al. (2013); Guha and Lawrence (2014). Now that we have solved for all the variables, the vortex sheet strength of the interfaces are given by

[TABLE]

III.1.2 Results

The formation of spilling breaker can be understood in terms of the resonant interaction of counter-propagating waves - it forms as a result of resonance between the leftward propagating surface capillary-gravity wave (at the free surface) and the rightward propagating vorticity wave (at the vorticity interface). The nonlinear structure of a spilling breaker from our simulation is given in figure III.1.2. The values of various parameters are as follows: $\sigma=0.072$ $\mathrm{Nm}^{-1}$, $g=10$ $\mathrm{ms}^{-2}$, $\Omega_{0}=\Omega_{2}=0$ $\mathrm{s}^{-1}$, $\Omega_{1}=10$ $\mathrm{s}^{-1}$, $\rho_{0}=0$ $\mathrm{kgm}^{-3}$, $\rho_{1}=\rho_{2}=1000$ $\mathrm{kgm}^{-3}$, $\alpha=1$ $\mathrm{m}^{-1}$. The initial amplitude of the surface wave is $\tilde{\eta}_{1}=10^{-3}$ m, and $h=0.5$ m.

The initially linear surface and interfacial waves lock in phase, resulting in an exponential growth that eventually gives rise to finite amplitude waves. These finite amplitude waves become steep at some location after which, instead of forming a horizontal jet that falls under the influence of gravity (as in the case of plunging breakers), the crest forms a bulge. At this point of time smaller scale secondary instabilities arise near the toe of the bulge, which occur due to the vorticity shed by the primary instability, as postulated by Longuet-Higgins (1994). We can see in figure III.1.2 that at the toe of the bulge there are some smaller scale curls which represent these secondary instabilities. Furthermore, our simulation results show good qualitative agreement with Duncan et al. (1994); see figure III.1.2.

III.2 Plunging breaker

In the previous section we implemented vortex method to simulate spilling breaker, an essentially deep water phenomenon occurring at capillary scales. Here we aim to show the applicability of vortex method to capture wind shear induced breaking of long surface gravity waves - the plunging breaker phenomenon. A linear surface gravity wave in shallow water travels with a constant phase speed of $\sqrt[]{gh}$, where $h$ is the water depth. Due to the effect of wind blowing over the ocean surface, the wave grows exponentially until it reaches an amplitude of $\mathcal{O}(h)$, at which point it interacts with the ocean bed and breaks. The effect of wind has been taken into account by considering the “boundary layer” over the free surface as a vortex patch of constant vorticity. The schematic of the system is shown in figure III.2.

In our simulation, at every time we first evolve the vortex patch and then we impose the velocity field due to the vortex patch on the surface gravity wave. The vortex patch evolves following contour dynamics techniquePozrikidis and Higdon (1985). The velocity induced by a periodic patch of constant vorticity $\Omega$ is given by

[TABLE]

where $C$ denotes the contour of a single patch that is periodically repeated in the $x$ direction with wavenumber $\alpha$. This interaction makes the surface wave grow exponentially. After the wave’s amplitude becomes comparable with the ocean depth, we remove the effect of the shear layer and let the wave interact with the bottom. The interaction between a surface wave and bottom boundary has been outlined in section II.4.

The base vorticity and density profiles for this case are given by

[TABLE]

For this particular case $\Omega_{0}=\Omega_{2}=0$ and $\Omega_{1}$ is a constant base shear. Furthermore $\rho_{0}=0$ (density of air assumed to be zero) and $\rho_{0}=\rho_{1}$.

III.2.1 Initial conditions

In this problem we initially start with a linear shallow water surface gravity wave which interacts with the constant vortex patch in the wind. The perturbation velocity potentials satisfy Laplace’s equation. On substituting the normal mode form, the eigenfunctions in the $z$ direction become

[TABLE]

Impenetrability at the bottom boundary yields

[TABLE]

We have four variables to solve - the three arbitrary constants from the velocity potentials ($C_{1},C_{2}$ and $C_{3}$), and the interfacial elevation $\tilde{\eta}_{1}$. There are two kinematic conditions, one just above and the other just below the free surface. Also we have Bernoulli’s equation across the free surface, and impenetrability condition at the bottom. Note here that we have assumed the density of air to be negligible in comparison to fluid density. So the matrix $\mathcal{A}$ for this case becomes:

[TABLE]

The variable vector $\mathcal{X}$ is as follows:

[TABLE]

We apply the condition for a non-trivial solution and obtain the dispersion relation and then find the null vector of the system. The vortex sheet strengths are given as

[TABLE]

III.2.2 Results

We start with a small amplitude surface gravity wave of $\tilde{\eta}_{1}=10^{-3}$ m. Above it there is a vortex patch of thickness $2h$ where $h=0.5$ m, whose lower contour matches with the surface wave. The nonlinear structure that the wave forms at some later time instant is known as a plunging breaker, where initially a horizontal fluid jet is formed at the crest of the wave which then falls down under the influence of gravity. Figure III.2.1 shows the structure of a plunging breaker. The values of the various parameters involved are: $\rho_{0}=0$ $\mathrm{kgm^{-3}}$, $\rho_{2}=1$ $\mathrm{kgm^{-3}}$, $g=1$ $\mathrm{ms}^{-2}$, $\Omega_{0}=0$ $\mathrm{s}^{-1}$, $\Omega_{1}=10$ $\mathrm{s}^{-1}$, $\Omega_{2}=0$ $\mathrm{s}^{-1}$, $\alpha=1$ $\mathrm{m}^{-1}$, and the depth of the sea is $d=0.5$ m. These parameter values, except that of $\Omega_{1}$, have been taken from Baker, Meiron, and Orszag (1982). Figure III.2.1 qualitatively and quantitatively matches very well with that of Baker, Meiron, and Orszag (1982), even though their mechanism of nonlinear wave generation is very different from ours.

IV Flows with finite shear - internal wave breaking problems

IV.1 Holmboe instability

Holmboe instability occurs as a result of the interaction between two counter-propagating vorticity and interfacial gravity waves (Holmboe, 1962; Carpenter et al., 2013; Guha and Lawrence, 2014). In this configuration the waves phase lock with each other and show exponential growth, finally forming nonlinear structures. Holmboe instability is not as strong as KH, where complete overturning of the density interface occurs. Instead, cusp like structures that eject from the density interface are the signature of Holmboe instability (Carpenter et al., 2010).

A schematic of the system has been shown in figure IV.1. The top interface shows a vorticity discontinuity and therefore supports a vorticity wave, while the bottom interface has a density discontinuity and supports an interfacial gravity wave. The initial condition for the problem has been derived from the linear theory. The base vorticity and density profiles read

[TABLE]

where $\Omega_{0}=0$ and $\Omega_{1}$ is a constant shear.

IV.1.1 Initial conditions

The derivation is similar to that of the spilling breaker except that here we neglect surface tension effects. So, the matrix $\mathcal{A}$ for this case becomes

[TABLE]

where we define:

[TABLE]

$U_{1}$ is the velocity at the upper interface while $U_{2}$ is the same at the lower interface. The variable vector $\mathcal{X}$ is given by

[TABLE]

On applying the condition for a non-trivial solution we obtain the dispersion relation as well as the null vector $\mathcal{X}$ of the matrix $\mathcal{A}$. For exponentially growing modes, the frequencies $\omega$ will be complex. The vector $\mathcal{X}$ corresponding to that particular frequency will give us the phase locked configuration. The vortex sheet strength of the interfaces would be given by (30)-(31).

IV.1.2 Results

We initialize the problem with an interfacial gravity wave of amplitude $\tilde{\eta}_{1}=10^{-3}$ m. The initial amplitude of the vorticity wave can be calculated from the null vector $\mathcal{X}$. The values of the various parameters involved are: $\rho_{0}=0.99$ $\mathrm{kgm^{-3}}$, $\rho_{2}=1$ $\mathrm{kgm^{-3}}$, $g=1$ $\mathrm{ms}^{-2}$, $\Omega_{1}=0.2$ $\mathrm{s}^{-1}$, $\alpha=1$ $\mathrm{m}^{-1}$, and $h=0.5$ m. In figure IV.1.2 we have shown the nonlinear structure of the density interface (the lower interface in the schematic). We can see that at later times there is ejection of fluid from the interface, just as observed by Carpenter et al. (2010) in their Direct numerical simulations as well as experiments. For comparison we have shown experimental results on wave ejections associated with Holmboe instability in figure 2.

IV.2 Taylor-Caulfield instability

A stable density stratification in the presence of a linear velocity profile gives rise to unstable modes even though both the density stratification and the velocity profile are stable by themselves. Taylor (1931) proposed that an instability can arise due to the resonant interaction between two counter-propagating interfacial gravity waves, respectively located at $z=h$ and $z=-h$ in figure IV.2.1. These waves lock in phase and eventually give rise to nonlinear structures.

The base vorticity and density profiles are as follows:

[TABLE]

Here $\Omega_{0}$ is some constant base shear.

IV.2.1 Initial conditions

The derivation of initial conditions is similar to the case of Holmboe instability but here the background shear is uniform throughout. So, the matrix $\mathcal{A}$ for this case becomes

[TABLE]

where we define:

[TABLE]

The variable vector $\mathcal{X}$ is given by

[TABLE]

We look for complex frequencies $\omega$ and the corresponding null vector $\mathcal{X}$ for the nonlinear structure. Again, the vortex strength of the interfaces are given by (30)-(31).

IV.2.2 Results

The values of the various parameters involved are: $\rho_{0}/\rho_{1}=0.99$, $\rho_{1}/\rho_{2}=0.99$, $g=1$ $\mathrm{ms}^{-2}$, $\Omega_{0}=0.13$ $\mathrm{s}^{-1}$, $\alpha=1$ $\mathrm{m}^{-1}$. The initial amplitudes of the two waves are $\tilde{\eta}_{1}=10^{-3}$ m and $\tilde{\eta}_{2}=10^{-3}$ m, and $h=0.5$ m. In figure IV.2.1 we have shown the late time nonlinear structure of the Taylor-Caulfield instability. For better representation we have chosen to depict three wavelengths. The thick lines represent the envelope of the instability (looks similar to Kelvin’s cat eyes); these are the late time structures of the vortex sheets shown in the schematic figure IV.2.1. The dashed inner contour represents a material line in the intermediate density region which displays a Kelvin-Helmholtz like roll up. Our results are qualitatively similar to Lee and Caulfield (2001); see figure IV.2.1. It is to be noted that our simulation is limited to temporally unstable waves, while DNS has no such restrictions. Therefore figure IV.2.1 shows spatially growing Taylor billows.

V Strongly sheared flows

V.1 Jet stream instabilities

Fluid flowing in the form of a jet stream is encountered in both natural flows like polar jet stream, as well as in industrial flows like fan spray nozzles, fuel injectors in diesel engines, etc. The system consists of a fluid of a given density flowing through a fluid of a different density; see figure V.1. There are density and velocity jumps across the two interfaces, which causes the latter to roll up like KH instability. Surface tension effects would be important for jets emanating from nozzles, but in this study those effects have been neglected. Presence of two interfaces leads to two modes of instabilities, viz. sinuous ([math] phase shift between interfaces) and varicose ($\pi$ shifted interfaces) mode. We first derive the linear theory of the system. Then we show numerical simulations for various density jumps across the interfaces, and discuss its effect on the roll up of the interfaces.

The base velocity and density profiles are as follows:

[TABLE]

We emphasize here that there is no gravity in the system (flows occurring in the horizontal plane).

V.1.1 Initial conditions

Like previous cases, we formulate the problem as a system of linear homogeneous equations $\mathcal{A}\mathcal{X}=0$. The $\mathcal{A}$ matrix is given by

[TABLE]

where

[TABLE]

The variable vector $\mathcal{X}$ is given by

[TABLE]

Similar to previous studies, we look for complex $\omega$ and the corresponding $\mathcal{X}$ for obtaining the eigenstructure. There are two possible initial configurations - sinuous mode and the varicose mode. The vortex strength of the interfaces are given by (30)-(31).

V.1.2 Results

The nonlinear structures for a jet stream flowing in a stationary media or even in a media having a different velocity than the stream look similar to the KH roll ups. From linear theory we get two possible unstable modes - the sinuous mode and the varicose mode. If we start with one of these modes, the system starts to grow exponentially and yields nonlinear structure. Figure V.1.2 shows nonlinear structures corresponding to each of these modes. The various parameters used for these simulations are: $\rho_{0}/\rho_{1}=0.9$, $U_{0}=0$ $\mathrm{ms}^{-1}$, $U_{1}=1$ $\mathrm{ms}^{-1}$, $\alpha=2\pi$ $\mathrm{m}^{-1}$, and $h=0.5$ m. Our results are very similar to the Direct numerical simulation results of Minion and Brown (1997) and Hashimoto et al. (2013).

On increasing the density ratio, the interfaces show reduced tendency for KH like roll-ups; see figure V.1.2. In this case $\rho_{0}/\rho_{1}=0.33$, the rest of the parameters are same as before.

VI Wave breaking in absence of shear: wave-topography interactions

VI.1 Bragg resonance

In resonant triad wave interactions, energy exchange occurs between three nonlinear waves. These resonances occur when some special conditions are satisfied - the wavenumber vectors as well as the corresponding wave frequencies add up to zero (Ball, 1964; Craik, 1988). Bragg resonance is a special case of such resonant triad interactions; it occurs when one of the three waves forming the resonant triad is a stationary undulated bottom topography (implying it is a wave of zero frequency). The undulated bottom acts as a mediator for the exchange of energy amongst the other two waves present in the system. The specific case considered here is a single layer fluid above a rippled bottom topography. Here Bragg resonance can occur between the two surface gravity waves (which are oppositely propagating with a speed of magnitude $c=\sqrt{g\tanh(\alpha_{s}h)/\alpha_{s}}$, where $h$ is the mean depth of the fluid and $\alpha_{s}$ is the wavenumber) present at the air-water interface and the bottom topography, as shown in figure VI.1. If we initialize with the configuration that a surface wave of wavenumber $\alpha_{s}$ is propagating to the right over a rippled bottom of wavenumber $\alpha_{b}$, then another wave of wavenumber $\alpha_{s}$ traveling towards left would be generated. This resonance occurs under the condition

[TABLE]

In figure VI.1, we have a surface wave of elevation $\eta_{1}$ and wavenumber $\alpha_{s}$, and a rippled bottom with elevation $\eta_{b}$ and wavenumber $\alpha_{b}$. The condition for a single layer Bragg resonance has been shown in figure VI.1. Here individual waves are represented by vectors in the $\alpha-\omega$ plane. This representation is known as the Ball’s diagram (Ball, 1964).

VI.1.1 Initial conditions

The simulation of Bragg resonance using vortex method requires the initial interface elevation and vortex sheet strength. The initial conditions are exactly the same as the plunging breaker problem presented in section III.2.1, since we initialize the system with a shallow water surface gravity wave which interacts with the rippled bottom topography to generate an oppositely propagating surface gravity wave. The matrix $\mathcal{A}$ and the variable vector $\mathcal{X}$ are therefore the same as given in section III.2.1. Finally, the vortex strengths of the surface wave and the bottom are again given by (34) and (35).

VI.1.2 Results

The values of various parameters involved in the simulation are: $\alpha_{b}=4$ $\mathrm{m}^{-1}$, $\alpha_{s}=2$ $\mathrm{m}^{-1}$, $\rho_{0}=0$ $\mathrm{kgm}^{-3}$, $\rho_{1}=0.7877$ $\mathrm{kgm}^{-3}$ , $g=1$ $\mathrm{ms}^{-2}$, $h=0.1$ $\mathrm{m}$ and initial amplitude of the surface wave $\tilde{\eta}_{1}^{(0)}=0.0016$ $\mathrm{m}$.

Initially, there is only a rightward traveling wave of wavenumber $\alpha_{s}$ present at the free surface. Resonant interaction with the rippled bottom gives rise to an oppositely traveling surface wave having the same wavenumber. The amplitude of the generated wave is initially zero and slowly increases with time, as shown in figure VI.1.1. Since the frequency of both the original surface wave and the generated wave are same, they show rapid constructive and destructive interference to produce an amplitude variation similar to that observed for a standing wave. For a standing wave, the amplitude fluctuates rapidly between maximum value (constructive interference) and zero (destructive interference), and the envelope of the amplitude fluctuation stays constant. However in this case, the envelope of amplitude fluctuation increases and then decreases because the amplitudes of the two oppositely propagating waves are unequal. It can be seen that around $t=100$ $\mathrm{s}$ the envelope is maximum (varies from [math] to $2\tilde{\eta}_{1}^{(0)}$). Hence at this instant the amplitudes of the two waves are equal. After this, the amplitude of the generated wave becomes larger than the original surface wave, and the difference continues to increase until the amplitude of the original surface wave becomes zero and then the cycle repeats.

VII Summary and conclusions

Various problems involving waves and hydrodynamic instabilities are conventionally studied in the linear regime using piecewise base velocity and density profiles. The advantage of using simplified profiles is that, on many occasions, they provide reasonably accurate qualitative and quantitative predictions of instabilities (as well as stable wave fields) in the initial stages. In fact, one can write analytical expressions for the range of instability, most unstable eigenvalue, corresponding eigenfunction, etc.  in terms of the physical parameters of the problem Drazin and Reid (2004). While this strategy has been very useful, it has largely remained confined to linear waves and instabilities (initial stages of the instability development). Previous works on extension of simplified profiles to the nonlinear regime were limited in the sense that they were applicable to the following cases: (i) homogeneous shear flows (contour dynamics technique)Pozrikidis and Higdon (1985), (ii) infinite background shear (jump in background velocity) with or without density jumpRosenhead (1931); Sohn, Yoon, and Hwang (2010), or (iii) density jump without any background velocityBaker, Meiron, and Orszag (1982); Tryggvason (1988). In density stratified shear layers, it is quite common to find one or more density interfaces embedded in a background shear flow. The aim of this paper has been to find a strategy that can model such flows. We have shown that the classical vortex method under suitable modifications can capture the fully nonlinear structures emanating from stratified shear flows. The main mathematical difficulty in this case was whether one can write unsteady Bernoulli’s equation in the presence of background shear. We have formulated a “shear-modified” Bernoulli’s equation, see (1)-(2), which is the corner stone in deriving the vortex sheet evolution equation (10). To the best of our knowledge, (10) is the most general vortex sheet evolution equation known to us. Once the evolution of a vortex sheet embedded in background shear has been known, the next objective, in a mathematical sense, was to find the interaction equation between multiple vortex sheets (i.e. interfaces). In this regard we have derived Birkhoff-Rott equation for multiple interfaces, see (16)-(17).

The mathematical model has been numerically implemented following the techniques outlined in Sohn, Yoon, and Hwang (2010). In this numerical technique, known as vortex method, interfaces are represented by a discrete array of point vortices that interact with each other. We have explored a wide variety of nonlinear phenomena - spilling (involves capillary-gravity waves) and plunging (involves long surface gravity waves) breakers, two classic stratified shear instabilities (Holmboe and Taylor-Caulfield), jet flows, as well as Bragg resonance (interaction of surface/interfacial gravity waves with rippled bottom topography, the latter being treated as a stationary wave). We observe that simple piecewise velocity and density profiles, when extended to the nonlinear regime, capture the essential nonlinear features, e.g. cusp formation and roll-ups observed in wave breaking, observed in experimental and/or involved numerical simulations with smooth, realistic profiles.

Finally, we concentrate on the usefulness and applications of this work. Since vortex method provides an exact (or near exact) solution of the (stratified) Euler equations, it is able to capture the correct nonlinear wave-mean feedback. This is of prime importance in the understanding of transition to turbulence in shear layers. Direct numerical simulations (DNS) are limited in this sense since they are only applicable to continuous profiles (no discrete interfaces allowed) and modest Reynolds numbers. Geophysical flows occur at very high Reynolds numbers, DNS of which would be extremely computationally expensive, perhaps impossible. Hence the “shear-modified” vortex method outlined in this paper would be useful, not only in capturing the detailed structures of high Reynolds number stratified shear flows, but also in understanding the mechanistic picture of wave-wave and wave-mean interactions. We emphasize here that there is a growing body of literature that tries to pinpoint the mechanisms behind stratified shear instabilities (e.g. works of Heifetz and co-workers, see Heifetz, Bishop, and Alpert (1999); Heifetz and Methven (2005); Heifetz, Harnik, and Tamarin (2009); Heifetz and Mak (2015)). This is because shear instabilities are often elusive and non-intuitive. Using piecewise profiles, these instabilities can be understood in a simplified manner in terms of resonant interaction of interfacial waves. While this wave interaction interpretation has been able to provide deep insights into shear instability mechanism, it has however remained limited to the linear regime. Our technique can very well extend the wave interaction interpretation to fully nonlinear regime, and can therefore help in providing a mechanistic understanding of the transition to turbulence.