On rotational flows with discontinuous vorticity beneath steady water waves near stagnation

We numerically investigate the flow structure of periodic steady water waves of fixed relative mass flux propagating on rotational flows with piece-wise constant vorticity. We show that for wave solutions along the global bifurcation diagram the stagnation point can first occur internally or at the bottom, and then again occurs at the crest with further increase in wave height. We observe that the bifurcation diagram has a new branch which is not connected to the trivial solution. Furthermore, we present an in-depth discussion of the results on pressure distributions and particle trajectories beneath large-amplitude steady waves near stagnation. We also expand on previous results concerning the amplitude and mass flux of steady water waves traveling on rotational flows with discontinuous vorticity.


Introduction
The interaction between surface water waves and underlying currents affects the wave properties significantly (Longuet-Higgins & Stewart 1961), as also confirmed by experiments, e.g., Swan et al. (2001). The interaction is hence of practical interest for design of offshore structures, see recent comparisons of structural responses with and without this interaction considered (Chen & Basu 2018a,b) where linear wave-current interaction models were used though; for coastal engineering applications (Kirby & Dalrymple 1983;Dong & Kirby 2012) and in geophysical fluid dynamics (Kirby 1984;Kirby & Chen 1989). Of most interest is the problem of two-dimensional (2D) periodic water waves propagating on steady current. A large number of theoretical studies and numerical investigations have been dedicated to this problem. However, it is still far from being fully understood. A comprehensive introduction to nonlinear wave-current interactions can be found in Constantin (2011).
In fact, it is only recently that the existence of the large-amplitude wave solutions for rotational flows with general vorticity (thus allowing arbitrary current profiles) in finite water depth has been proved by Constantin & Strauss (2004). Prior to this, only waves of small-amplitude were known to exist in the presence of vorticity (Dubreil-Jacotin (1934)). Hur (2006) further extended the proof of Constantin & Strauss (2004) to the case of infinite-depth with the additional restriction that the vorticity distribution is positive and monotone. Hur (2011) later completed the proof of existence of large-amplitude waves with infinite depth removing the mentioned restrictions on vorticity. Moreover, Constantin & Strauss (2011) and Strauss (2012) addressed the existence of periodic traveling waves of small and large amplitude in a flow with an arbitrary bounded but discontinuous vorticity.
Following the work by Constantin & Strauss (2004), rotational flows with constant and non-constant vorticity have received considerable attention. Notably, many numerical investigations, on both irrotational and rotational flows, appeared much earlier. For example, for irrotational flow in finite water depth, Chappelear (1961) expanded the velocity field and the wave profile in Fourier series and the Fourier coefficients were numerically determined from the Bernoulli equation by using the least squares method. The direct numerical simulation results were compared with the fifth-order Stokes theory. Dalrymple (1974a) presented a numerical model for steady water waves propagating over a linear shear current (constant vorticity) where stream function was expressed as the superposition of a linearly varying current profile and an irrotational wave field. The results thereof showed that the vorticity has a significant impact on the wave height. Simmen & Saffman (1985) and Teles da Silva & Peregrine (1988) studied periodic steady water waves with constant vorticity for infinite and finite water depth, respectively. The numerical methods for irrotational flows and rotational flows with constant vorticity are generally based on the boundary integral method or the Fourier expansion. Recently, Clamond & Dutykh (2018) developed a novel conformal mapping method for irrotational flows, where the governing equation is rendered of Babenko kind which is then solved using the Petviashvili method. The model is very efficient and is able to solve for large amplitude waves with quite large steepnesses over arbitrary depth. Dyachenko & Hur (2019b) further modified the Babenko equation to permit constant vorticity and applied Newton-GMRES method following the use of the Fourier expansion. There are also other methods that have been developed. For instance, Ablowitz et al. (2006) proposed a nonlocal formulation which was further applied to consider non-flat bottom and moving non-flat bottom boundaries (Fokas & Nachbin 2012;Francius & Kharif 2017). Wahlén (2009) and Ehrnström et al. (2011) devised a novel transform for small-amplitude waves with constant and linear vorticity, i.e., a vertical scaling transformation.
For waves on rotational flows with non-constant vorticity, early numerical investigations date back to Dalrymple (1973Dalrymple ( , 1977 where the author considered the effects of various current profiles on the finite-amplitude steady waves over finite water depth. The method was based on the Dubreil-Jacotin (DJ) transformation and finite difference discretization of the transformed rectangular domain. For a given boundary condition, the solution to internal points is improved by a number of sweeps; the error on the free surface boundary condition is then evaluated; subsequently, the solution on the free boundary is updated and the procedure is repeated until convergence. Using this method, Dalrymple & Cox (1976) showed the influence of the vorticity on wavelength and crest elevation of the wave, considering that the currents vary as trigonometric and hyperbolic sine and cosine functions of the depth. Thomas (1990) performed experimental studies for a finite-amplitude wave train interacting with a steady current containing an arbitrary distribution of vorticity, and moderate agreement was found between the measurements and the results obtained using the model of Dalrymple (1977). Ko & Strauss (2008b) presented a numerical continuation method to find largeamplitude 2D periodic steady water waves with arbitrary vorticity, which was used for discussing the relationship between amplitude, hydraulic head, depth and mass flux for variable vorticity cases (Ko & Strauss 2008a). This method neither uses truncated Fourier expansion nor makes any assumption on the water depth or wave amplitude. Further numerical analyses of water waves on rotational flows with discontinuous vorticity can be found in Strauss (2010), Strauss (2011), andConstantin (2012a). The ethos of the approaches towards wave computation is that the local bifurcation from a laminar flow to a small-amplitude wave solution is predicted by a dispersion relation, the algorithm beginning with the linearized solution. The wave solution is then continued with increasing wave height until the solution becomes close to a wave in a flow with stagnation point (global bifurcation). The dispersion relation has been derived for small-amplitude periodic water waves traveling on rotational flows for various vorticity cases, e.g., see Constantin & Strauss (2011);and Constantin (2012a) and Martin (2015) for flows with two layers. Karageorgis (2012) discussed dispersion relations for several types of non-constant vorticities. High-order approximations for the constant vorticity case were discussed in Constantin et al. (2015a). Amann & Kalimeris (2018) expanded upon the method of Ko & Strauss (2008b), by allowing nonuniform grid points and enabling the continuation with vorticity as a bifurcation parameter. Using the latter scheme they successfully discovered a new part of the bifurcation curve for some critical constant vorticity values (see also Kalimeris (2018)). Furthermore, Constantin et al. (2015b) proposed a penalization method to solve for large-amplitude waves. Note that using the method of Ko & Strauss (2008b), the mean water depth is recovered from the solution, which leads to a small variation in the value of mean water depth along the bifurcation curve (Ko & Strauss 2008a). The mean-depth can be held constant by following the reformulation proposed by Henry (2013a,b).
Considerable understanding of water waves on irrotational flows and flows with constant vorticity has been established recently. For irrotational flows, Constantin (2006) proved analytically that the particle trajectories are not closed on Stokes waves and Henry (2006) extended the proof to deep water. Constantin & Villari (2008) and Matioc (2010) showed that there are no closed particle orbits in linear periodic water waves. Borluk & Kalisch (2012) reconstructed the velocity fields and the particle paths of long-crested waves on irrotational flows using the KdV approximation. These particle paths were found to be consistent with the Stokes drift but with increasing amplitude. Carter et al.
(2020) studied the particle paths in the nonlinear Schrödinger models. Grue et al. (2014) experimentally studied the velocity field below large-amplitude periodic water waves and the results are generally consistent with nonlinear theories. Grue & Kolaas (2017) further investigated the particle paths in steep waves at finite water depth (without background current) experimentally, focusing on the Stokes drift. It was shown that for steep waves on irrotational flow, two closed particle paths occur at two vertical positions, in contrast to the theoretical results (Constantin 2006). The lower closed particle paths occur at a position quite close to the bed due to the boundary layer effect, whereas the theoretical study by Constantin (2006) assumes inviscid water. It is also clear from Grue & Kolaas (2017) that the thickness of the boundary layer is about 5 percent of the water depth, which is significantly affected by viscosity effect. Paprota & Sulisz (2018) developed a theoretical model for studying the particle paths of nonlinear waves generated in a closed flume and the Stokes drift velocity obtained from experimental data (Paprota et al. 2016) were found to be quite consistent with theories. For a summary of the present analytical knowledge on the properties of irrotational flows, including the velocity field, pressure and particle trajectories, see Constantin (2015). For waves on a linear-shear flow, Ehrnström (2008) studied the pattern of streamlines and particle trajectories and pointed out that for non-positive vorticity the particles display a forward drift. Using an asymptotic expansion for the stream function, Ali & Kalisch (2013) constructed the pressure beneath steady long gravity water waves traveling on flows with constant background vorticity. Curtis et al. (2018) investigated the effect of linear background shear on the properties of wavetrains, particularly on particle paths, in deep water using nonlinear Schrödinger models. Numerically, for understanding the particle paths or trajectories, Nachbin & Ribeiro-Junior (2014) proposed a method for particle trajectories in Stokes waves in the boundary integral formulation for irrotational and constant vorticity cases. A summary of the flow structure beneath rotational flows with constant vorticity and the corresponding numerical method can be found in Nachbin & Ribeiro-Junior (2017). Lately, Dyachenko & Hur (2019a) recovered various characteristics of wave profiles in flows with constant vorticity.
For rotational flow with general vorticity but subject to certain constraints, Constantin & Strauss (2007) studied the points of maximal horizontal velocity of waves near stagnation. Further understanding on the velocity field has been established by Varvaruca (2008) and Basu (2019). Vasan & Oliveras (2014) discussed the pressure beneath rotational flows with constant vorticity pointing out interesting applications in wave profile recovery (Oliveras et al. 2012;Basu 2017b,a).
For surface waves on rotational flows, there are still a number of open problems yet to be solved (Strauss 2012). However, it becomes extremely difficult to capture the properties of the surface wave and the flow structure analytically for complex cases. As manifested in the literature, numerical investigations can provide us with important insights for such cases. To deal with the otherwise intractable problem of rotational flows with arbitrary voriticity distribution, contrary to the widely assumed linear shear flows, we focus on the shear flow with piece-wise constant vorticity which is the simplest case deviating from constant vorticity. This is a valid approximation for waves traveling on depth-varying currents (Dalrymple 1974b;Swan et al. 2001). In particular, we are concerned with shear flows with two layers. It is known that wind typically induces vorticity in the layer near the surface, while near the ocean shore, vorticity can be generated at the bottom, e.g., by tidal action (Ko & Strauss 2008a). In this study, we report on a new branch of the bifurcation curve for waves on rotational flows with discontinuous vorticity and investigate the pressure beneath the waves for flows with nearly approaching stagnation points. Further, we develop a method for computing particle paths based on the exact velocity field solved numerically. The particle paths are of increasing interest since they are related to the Lagrangian transport and the Stokes drift in ocean (Van den Bremer & Taylor 2016; Van den Bremer & Breivik 2017).
The rest of this paper is structured as follows. In Section 2 and Section 3, we briefly introduce the problem and the numerical methods respectively. The numerical results are then presented in Section 4. We close the study by a brief summary in Section 5.

Statement of the problem
We consider two-dimensional steady, periodic, rotational gravity traveling surface waves, of period 2L (L > 0 being a constant) propagating over water flows bounded below by the flat bed y = −d, where d > 0 is a constant. The water is assumed to be inviscid. In the numerical investigation, we shall place our focus on the case with discontinuous vorticity of piece-wise constant type.

Governing equations
The time-dependent water wave problem is presented in the Cartesian (X, Y )coordinate system with the origin on the mean water level, having the X-axis pointing in the horizontal direction and the Y -axis pointing upwards. The steady character of the problem we consider allows us to repeal the time-dependence by using the change of coordinates (X − ct, Y ) → (x, y), where c is the constant speed of the traveling waves and t denotes time. Hereafter, we will be working in the moving frame (x, y) with the origin placed under the crest, except for the computation of particle paths. In this setting, the motion of the water flow is governed by the Euler (or momentum) equations, the equation of mass conservation supplemented by the kinematic boundary conditions on the two boundaries of the flow and the dynamic boundary condition on the free surface. The latter condition at the free surface decouples the motion of the water from the motion of the air above, while the former kinematic conditions convey the fact that the bed and the free surface are impermeable boundaries, cf. Constantin (2011).
Aiming at a simplification of the analysis we introduce the stream function, denoted by ψ, and defined (up to a constant) by (2.1) Here, (u, v) is the velocity field in the (x, y)-coordinates. The mass flux relative to the uniform flow at speed c (normalized with respect to density) is further referred to as relative mass flux. It turns out from mass conservation equation and from the kinematic boundary condition on the surface that p 0 is a constant. From (2.1) and the kinematic condition on the free surface we infer that ψ is constant on the free surface. Therefore, we choose ψ = 0 on the surface and thus, by (2.2), we have ψ = −p 0 at the bottom y = −d. From the first relation in (2.1), we have Since we assume that u < c throughout the fluid (Constantin & Strauss 2004), p 0 is always negative. Following Constantin & Strauss (2004), the fluid vorticity is defined as Moreover, the assumption of no-stagnation (since u < c) implies, cf. Constantin (2011), that there is a function of one variable, which we denote by γ, such that ω(x, y) = γ(ψ(x, y)) for all (x, y) in the water domain. Also from formula (2.1), we have that ω = ∆ψ. Hence, the fluid vorticity satisfies Here, γ is the function already found by the previous argument involving the nostagnation condition, as presented in the book by Constantin (2011). The density normalized hydraulic head E is expressed as which is a constant owing to conservation of energy. In the expression, P is the pressure (normalized with respect to density) in the fluid and at the surface P = P atm , the atmospheric pressure which also normalized with respect to density. Therefore, for a given water depth, we can introduce another constant Q = 2(E − P atm + gd), and the boundary condition on the surface is given by Eventually, we obtain the fixed domain boundary-value problem from the free-surface wave problem by applying the following DJ transformation, i.e. q = x and p = −ψ(x, y) (2.8) to the relations (2.5) and (2.7), and including the boundary condition at p = p 0 . The equations are where h(q, p) = y + d is the height function and R is the rectangle defined by R = {(q, p) : −L < q < L, p 0 < p < 0}.

Recovery of wave properties
We will solve the problem (2.9) for the height function using a numerical method that will be introduced in the next section. For characterizing the waves, we shall recall the recovery of wave properties from the height function. Since in this formulation, the mean water depth d is not fixed, it needs to be recovered first from h by Afterwards, the water profile is determined by The wave height a is obtained as as we are concerned with waves of one crest and one trough per period. Further, the velocity field is recovered by recalling the DJ transformation (see Constantin (2011) for details). The pressure in the fluid domain can be obtained from the relation (2.6) and here we define a relative pressure, P = P − P atm , for the following discussions as (2.14) For computing water particle paths or trajectories, the traveling speed of the wave is recovered (Constantin 2011) by using where κ is the average horizontal current strength on the bed. We then need to solve a set of differential equations (Constantin & Strauss 2010; Nachbin & Ribeiro-Junior 2014) where the prime denotes the derivative with respect to time. Here, we use a numerical integration method to compute the particle paths once the velocity field is solved, as will be detailed in the next section.

Numerical methods
In this section, we briefly introduce the numerical methods to solve (2.9), for performing numerical continuation to arrive at the solution of waves in flows with closely approaching stagnation points, and for computing the particle trajectories from the velocity field in the Eulerian setting. The dispersion relations and the linearized solutions for local bifurcation from the laminar flow solutions to small-amplitude wave solutions are also reviewed for the sake of completeness.

Discretization and implementation
Following the methods in Ko & Strauss (2008b) and Amann & Kalimeris (2018), we use a finite difference scheme to discretize the computational domain. Taking advantage of the symmetry of h function with respect to q = 0 , we set the computational domain as a fixed rectangle [0, L] × [p 0 , 0] to improve the computation efficiency. We first discretize the domain using uniformly distributed grid points and then refine the grid size near the surface and the bottom under the crest (q = 0). In the case of discontinuous vorticity, we additionally refine the grid size near the location where the vorticity changes abruptly. We apply the central difference scheme at interior grid points, and at left and right boundary points by mirroring h along q = 0 and q = L, considering the periodicity and symmetry properties of the waves in absence of stagnation points. At the boundaries p = 0 and p = p 0 , we apply the backward and forward difference schemes, respectively. Note that the refinement is achieved by reducing the grid size smoothly so that the central difference still has a second-order accuracy, approximately.
For implementation, we use the open-source library Eigen (Guennebaud et al. 2010) for solving the nonlinear algebraic equations subsequent to the finite difference discretization. We use the predictor-corrector method detailed in Amann & Kalimeris (2018) for numerical continuation and write our own code on the basis of Eigen. In the following numerical analysis, we discretize the domain first by an equidistant 250 × 750 grid with additional refinement at locations where the stagnation point may occur. Following the refinement, the smallest grid size achieved is 1/8 th of the initial size prior to the mentioned refinement.
Grid points Particle Figure 2. Interpolation strategy for instantaneous particle velocity. In the figure, u = [u v] is the velocity vector; l1, l2, l3 and l4 are respective distances of the particle to the four neighboring grid points.

Dispersion relations and linearization solutions for local bifurcation
As previously mentioned, we focus on the discontinuous vorticity. The vorticity function is expressed as where p 1 corresponds to the value of the streamline at the interface y = y 1 , and γ 1 and γ 2 are the constant vorticity values of the surface layer and the layer adjacent to the bed. The setting of the problem is illustrated in Figure 1. For the case of γ 1 = γ 2 = 0, the problem reduces to irrotational flow, and the dispersion relation and the linearized solutions are available in Constantin (2011). For the situation of a rotational flow over an irrotational and vice versa, the dispersion relations are explicitly derived in Constantin & Strauss (2011) and Constantin (2012a), respectively. For steady water waves on two rotational layers, see Martin (2015). The dispersion relations are nonlinear equations. Here, we solve them using the Newton method. For large vorticity, we apply a simple continuation scheme by using the solution for smaller vorticity for initialization (Seydel 2009).

Computation of the particle trajectories
We numerically compute the particle paths in the physical domain, i.e. in the fixed (X, Y )-coordinate. For a particle located at (X(t), Y (t)) at time t, its velocity is interpolated from the velocities at the neighboring four grid points (following appropriate translation) using a similar interpolation method as in Umeyama (2012). The interpolation scheme is illustrated in figure 2. Note that in our case, the cell is of trapezoidal shape in the physical domain. Subsequently, we use numerical integration to solve the first-order differential equations (2.16). This simple interpolation scheme is found to be sufficiently accurate.

Results
Following Ko & Strauss (2008b,a) and Amann & Kalimeris (2018), we assign g = 9.8, 2L = 2π and p 0 = −2. We focus on the new part of the bifurcation curves, pressure distribution, and particle paths under flows with nearly approaching stagnation points. Note that in the present system all the variables are expressed in SI unit and a physical system of a different wavelength can be scaled to the 2π-periodic system analyzed in  this paper by scaling g and γ using the wave number (Constantin 2011). However, the scaling has no influence on the laminar flow solution. In the following, as will be seen, the current velocity at the crest (under zero-amplitude waves) in most of the cases is less than 2.5 m/s, which is physically realistic in oceans.

Critical vorticity and the new part of the branch of bifurcation diagram
For the case of constant vorticity, it is known that when the vorticity is below a critical value γ crit , the stagnation point first appears at the bottom and at that stage numerical continuation breaks down (Ko & Strauss 2008b). However, Amann & Kalimeris (2018) and Kalimeris (2018) found that when γ < γ crit but is very close to γ crit , the branch of the bifurcation diagram (a − Q plot) has another part which is bounded below and above by two values of the bifurcation parameter Q. These bounds correspond to two different waves in flows with stagnation point at the bottom and at the surface. On the new part of the branch, the wave height is further increased, and it is thus referred to as the upper part. To reach this upper part of the solutions, we need to overcome the gap of the bifurcation curve corresponding to waves in flows with stagnation points. For this purpose, Amann & Kalimeris (2018) resorted to the γ−continuation scheme where a solution on the bifurcation diagram corresponding to γ 0 (γ 0 > γ crit > γ) is used to initialize local bifurcation to arrive at one solution on the upper part corresponding to γ. Once one solution is obtained, the continuation in the solution set corresponding to γ can be performed again with Q or the wave height as the bifurcation parameter. In the following, we apply this scheme when needed.
Proceeding from the observations of Amann & Kalimeris (2018), we investigated the case of discontinuous vorticity and found that the upper part also exists. Furthermore, for the case of discontinuous vorticity, the upper bound on this part always corresponds to waves in flows with the stagnation point at the surface, while its lower bound corresponds to waves in flows with stagnation point either at the bottom or at an internal point where vorticity jump occurs. We first show an example of the former situation. An irrotational flow (γ 1 = 0, p 1 = −0.5) was considered on a rotational flow with vorticity γ 2 . We varied γ 2 to find out the critical value γ 2,crit below which the stagnation point first occurs at the bottom. It was found that γ 2,crit is slightly less than −3.22. Note that with the vorticity values γ 1 = 0 and γ 2 = −3.23 and hence a shearing laminar flow with profile shown in figure 4a over a total water depth of 0.787 m, we have the surface velocity, u = −2.07 m/s (the current flow direction is opposite to the wave propagation direction). Eq. (2.5) with a bifurcation value of Q = 26.85 (see figure 3a) leads to the relative surface velocity with respect to the wave speed as c − u = 3.38 m/s, since u < c. The wave speed c is  1.31 m/s as also verified by the dispersion relation for this case. In the case of positive vorticity in the bottom layer, e.g., γ 1 = 0 and γ 2 = 3, the water depth is 0.74 m, the current velocity at the surface is around 1.53 m/s, and at the local bifurcation the wave speed is c = 3.69 m/s. These values of the physical parameters indicate that the setting considered in this paper are realistic.
We then further decreased γ 2 and found the upper part, as shown in figure 3. For the case of γ 2 = −3.23, the streamlines of wave solutions at the four endpoints (marked with cross in figure 3) of the two parts of the bifurcation curve are shown in figure 4. Figure 4 clearly shows that the stagnation point first occurs at the bottom and then at the crest. velocities c − u at the bottom in these two cases are comparable, as in figure 5a, while the wave in figure 4c has a larger height, a sharper crest and smaller horizontal velocity c − u at the surface (see figure 5b). It is seen from figure 5 that the velocity c − u at the bottom is first decreasing and then increasing along the bifurcation curve, with the presence of a gap for γ 2 < γ 2,crit , while the velocity c − u at the crest is decreasing monotonically. Indeed, at the uppermost end of the bifurcation diagram in all those three cases, the velocities c − u at the surface are still larger than those at the bottom, the stagnation point is expected to occur at the crest because the velocity at the surface is decreasing while at the bottom is increasing with respect to the increasing wave height.
We further consider a layer (p 1 = −0.7) of rotational flow with vorticity γ 1 on an irrotational flow (γ 2 = 0), following the setting in Ko & Strauss (2008a) and Constantin & Strauss (2011). We found that the critical vorticity γ 1,crit is slightly below −10.41 (a value refined from −10 in Ko & Strauss (2008a)). By using γ 1 as the bifurcation parameter,   we also found the upper part of the branch of the bifurcation diagram when γ 1 is less but very close to γ 1,crit . Figure 6 shows the bifurcation diagrams for γ 1 = −10.41, −10.42, and −10.45. For γ 1 = −10.42, the streamlines of the waves corresponding to the four endpoints of the bifurcation curve are shown in figure 7, which clearly illustrates that the stagnation point first occurs internally and later occurs at the crest with further increased wave height. Figure 8 shows the variation of the velocities c − u at p = p 1 and at the surface along the bifurcation curves. The velocity c − u at the position where vorticity jumps, is first decreasing and then increasing with respect to the wave height. Notably, the velocity decreases slightly again when close to the wave in flows with stagnation point at the surface. Again, the velocity c − u on the crest is decreasing monotonically along with the increasing wave height.
To summarize, in this subsection, we have observed new features of the bifurcation curves with fixed mass flux, for large amplitude waves on flows with discontinuous vorticity. Specifically, we show that along the bifurcation curve the stagnation point can first occur at the bottom or internally and then occur again at the surface for some particular vorticity distributions.

Pressure distribution
We now present the pressure distribution for waves with stagnation points at the bottom, at an internal position and at the surface, for the case with discontinuous vorticity.
For an irrotational flow on a rotational flow with p 1 = −0.5 and γ 1 = 0, we consider the   figure 11 and figure 12 shows the pressure along the bottom and the crest line. Two features of the pressure distributions are consistent with the pressure beneath a Stokes wave as proved by Constantin & Strauss (2010) and those can be stated as ''the pressure in the fluid strictly decreases horizontally away from a crest line toward its neighboring trough lines, and the pressure strictly increases with depth''. Even though, when γ 2 is  close to the critical value, the pressure along the bottom displays a "plateau" behavior, as previously pointed out by Amann & Kalimeris (2018) for the case of constant vorticity, the pressure is still decreasing away from the crest line, slowly.
For a rotational flow on an irrotational flow, we analyze the cases for γ 1 = 5, −5, and −10.42. Figure 13 shows the pressure distribution in flows close to having a stagnation point at the crest, and variations of the pressure along the bottom and along the crest line are illustrated in figure 14. The pressure distribution in flows when internal stagnation points are almost developed is illustrated in figure 15. Correspondingly, figure 16 shows the pressure along the streamline ψ = −p 1 and along the crest line. From these figures, it is seen that none of the facts regarding the pressure distribution for an irrotational Stokes wave (Constantin & Strauss 2010) hold, as also observed by Ribeiro et al. (2017) for rotational flows with constant vorticity. When the vorticity γ 1 is below a critical value, the pressure along the streamline ψ = −p 1 achieves a local minimum under the crest, as opposed to the maximum there when γ 1 is large. For smaller γ 1 , the relative pressurẽ P at this position becomes negative (i.e. sub-atmospheric), which is the minimum along the streamline. The pressure attaining a value less than the atmospheric pressure has also been observed by Ali & Kalisch (2013) for the case of constant vorticity using a model based on an asymptotic expansion for the stream functions. When γ 1 is below a critical value, the pressure below but close to the layer where the vorticity jump occurs also exhibits the ''plateau'' behavior, internally, see figure 15. Figures 14-16 show that when γ 1 is much larger than the critical value, the pressure is increasing with depth along the crest line. Interestingly, when γ 1 is quite close to γ 1,crit , the pressureP along the crest line first increases slightly, and then decreases, and increases again when close to the position where vorticity jump occurs; the overall pressure in the surface layer is close to the atmospheric pressure, see figure 14b for γ 1 = −10.42. With a further decrease in the value of γ 1 below the critical value, the relative pressure in the surface layer becomes negative, first decreasing from the crest and then increasing near the position where the vorticity jump takes place. To understand the reason for the development of a negative relative pressureP in the surface layer near the crest line, we recall the expression (2.14). For the cases illustrated in figure 14, the term p 0 γ(−s)ds in (2.14) is positive in the surface layer and increases rapidly with depth since |γ 1 | is large, the term g(y + d) is also large considering the large wave height, and the velocity c − u is also not small, as shown in figure 17; henceP in the surface layer becomes negative. Further decreasing γ 1 (with increasing |γ 1 |), it seems that the minimal pressure along the crest line will keep decreasing.

Particle trajectories
We proceed to compute the particle trajectories under the waves corresponding to flows with nearly developed stagnation points, for the previously mentioned cases. To recover the wave speed using expression (2.15), we assume that the current strength at the bed is zero, i.e., κ = 0. Note that the choice of a different setting (adding an arbitrary uniform current) may lead to variations of the resulting particle paths. For example, when studying particle paths under waves on linear shear flow, Ribeiro et al. (2017) set the current strength at the mid-depth to be zero (to ensure a flow with no net mass flux in the case of zero-amplitude waves). However, we here focus on the relative horizontal displacements of the particles along the depth, which are not affected by the presence of a uniform background current.
In the following computation and illustration of the particle trajectories, we set that: (i) the initial positions of the particles are under a zero-crossing point of the wave profile (Paprota & Sulisz 2018); (ii) the particles trajectories are plotted for the waves traveling rightwards for two wave periods, i.e., for a time of 4L/c = 4π/c, and the final position of the particles under zero-amplitude waves are also indicated in the graphs for comparison. In the graphs, six particle trajectories are plotted. Three of those particles are located at the surface, at the interface, and at a position very close to the bottom but not exactly at the bottom. The other three particles include one located in the top layer and two located in the bottom layer. The initial positions are denoted by dots and the final positions are indicated by square markers.
First of all, we plot the computed particle trajectories in the moving frame, as illustrated in figure 18, which theoretically are the corresponding streamlines. The streamlines passing through the initial particle locations are plotted for comparison, and a good agreement is seen, suggesting the accuracy of the method for computing particle trajectories.
For an irrotational layer on a rotational layer, the particle trajectories in flows close to having a stagnation point at the crest are plotted in figure 19 and those in flows close to having a stagnation point at the bottom are plotted in figure 20. We show that as theoretically predicted by Constantin (2006) when a particle reaches the crest of the wave in flows having a stagnation point at the crest, it does not pause there but moves further (figure 19a). However, for waves in flows close to having an internal stagnation point or  (ii) In the top irrotational layer, the horizontal displacement (with respect to the background current induced displacement) or dirft is larger when the particle is closer to the surface than when it is near the interface of the two layers, which is consistent with experimental and approximate theoretical results (Paprota & Sulisz 2018;Curtis et al. 2018). In the bottom rotational layer, when the vorticity is positive, the wave induced drift is larger when the particle is close to the interface as compared to the drift induced  (iii) For flows close to having a stagnation point at the crest, as shown in figure 19, the negative vorticity in the bottom layer clearly increases the drift of the particles in the top irrotational layer relative to the background current, particularly for particles at the surface.
(iv) For flows close to having a stagnation point at the bottom, as in figure 20, a stronger negative vorticity in the bottom layer tends to decrease the relative drifts of the particles at different vertical positions in the top layer, as observed by comparing figure  19c and figure 20c.
With a shear layer on an irrotational layer, figure 21 shows the particle trajectories for flows close to having a stagnation point at the crest and figure 22 demonstrates the particle trajectories for flows close to having a stagnation point at the location when the vorticity jump occurs. Again, we can see from both figures that the wave-induced drifts are positive in the direction of the propagating waves. Some conspicuous features are observed as described below: (i) In the bottom irrotational layer, the drift is always larger when the particle is close to the interface, consistent with particle paths in irrotional flow (Paprota & Sulisz 2018).
(ii) In the top rotational layer, when the vorticity is positive, the wave-induced drift is larger when the particle is near the surface as comapred to the drift when the particle is close to the interface of the two layers; while with negative vorticity, for flows close to having a stagnation point at the interface, the particle with the largest drift is at the interface.
(iii) A negative vorticity in the top layer tends to increase the wave-induced drifts of particles in both layers.

Conclusions
In this paper, we investigate the behavior of flow along a new bifurcation curve, and the pressure and particle trajectories under waves traveling on a rotational flow with discontinuous vorticity (close to having a stagnation point). For this purpose, we have developed a method for computing the particle paths using the velocity field obtained from a numerical continuation method based on DJ transform of the Euler equations for periodic waves on rotational flow. The continuation method generates a group of waves of fixed relative mass flux until a point where a flow with a stagnation point is reached.
By using the previously mentioned methods, we have recovered new branches of the bifurcation curves not connected to the trivial (laminar) solution, for waves on rotational flow with discontinuous vorticity. Particularly, we show that along the bifurcation curves, the stagnation point can first occur at the position where the vorticity jumps and then again occur at the surface. On the pressure distribution, we show that for the case of shear layer on an irrotational layer, the pressure in the top layer under the crest could be less than the atmosphere pressure and the pressure in the top layer near the interface displays a "plateau" behavior when the vorticity is large with a negative sign (i.e., adverse). This "plateau" behavior of the pressure has only been observed along the bottom for linear shear flow. We further show the particle trajectories in flows close to having a stagnation point for various cases and qualitatively discuss the influence of vorticity of the rotational layer on the drifts of the particles in both layers.

Declaration of Interests
The authors report no conflict of interest.