Stationary water waves on rotational flows of two vortical layers

Stationary waves of constant shape and constant propagation speed on rotational flows of two layers are computed numerically. Two layers are assumed to be of distinct constant vorticity distributions. Three different kinds of waves of finite depth are considered: pure capillary, capillary-gravity, and gravity waves. The problem is formulated as a bifurcation problem, which involves many parameters and produces a complicated structure of solutions. We adopted a numerical method by which waves with stagnation points can be computed, and obtained variety of new solutions. It is also reported that the locations of the stagnation points vary curiously with the prescribed parameters and that they offer an interesting problem.


Introduction
We consider two-dimensional stationary water waves on two layers of different vortices. 'Stationary' implies that the wave profile is of constant shape and the propagation speed is a constant, if it is viewed in a frame of coordinates moving in the same where p 1 is another prescribed constant such that 0 < p 1 < p 0 ; 1 and 2 are prescribed constants which may take any real numbers. Note that = p 0 represents the bottom y = −d , while = 0 represents the free surface. The curve defined by = p 1 is called interface: see Fig. 1. We assume that the wave profile is symmetric in x about the y-axis. This assumption implies that the interface also possesses the same symmetry. This assumption is not trivial since there are non-symmetric waves on irrotational flows, see [17,19,22]. However, we do not study these nonsymmetric waves in the present paper, since = 0 (y = h(x)), = p 0 (y = −d). symmetric waves alone are rich enough. We also assume that the unit of length is so chosen that L∕2 = . Under this convention it is sufficient to consider the flow in 0 < x < , −d < y < h(x) . (Note that d and h are already normalized here.) d is assumed to be finite for the numerical computation's sake. But d = +∞ causes no problem theoretically.
Since the wave profile y = h(x) is a part of the flow domain, and since it is unknown a priori, this is a free boundary problem. Free boundary problems usually accompany mathematical difficulties, and many techniques have been devised to overcome them. In the framework of two-dimensional rotational wave, the following technique is often used (see for instance [17] or [6]): (x, y) is transformed to (q, p) ∶= (x, ) . Then, under certain assumptions, everything can be regarded as a function of (q, p). We then define h by The differential equations are rewritten in terms of h . However, for the notational reason, we use the following symbols instead of h . We define u and U as h in 0 < q(= x) < , p 1 < p < p 0 and 0 < q < , 0 < p < p 1 . They are nothing but restrictions of h . However, in order to keep compatibility with [6,7,14] we use (u, U) instead of h . We now see that the original free boundary problem is converted into the following nonlinear partial differential equations in a fixed rectangular domain (see Fig. 1(right)): The following boundary conditions are imposed: This system of equations appears in [6,7] in the case where the surface tension is neglected. The present general case is due to [14]. This is the mathematical formulation by which we can compute the waves. Here the domain, p 0 and p 1 are given, and u, U, Q are to be sought. The reader should note that the formulation above does not yield solutions which contain stagnation points. They are a priori excluded. In order to compute such solutions as having stagnation points, we must use more complicated set of equations. This will be introduced later in Sect. 4. h(q, p) = d + y. (1) Traditionally, as is described in [17], water waves on irrotational flows have been targets of research. The irrotationality assumption is accepted as a rough approximation. However, in flows of real world, vorticity is inescapable. Water waves on rotational flows are important in real-world applications. For instance, flows near the water surface have vorticity because of wind. Flows near the bottom becomes rotational because of bottom friction. In these cases, it seems to us that a rotational flow of two different vorticity layers is one step closer to the real-world water waves.
We now review some known results. We admit that we are ignorant of recent rapid progress in theoretical studies of rotational waves. Accordingly the following review is far from complete. Rotational water waves with general vorticity distribution F = F( ) were mathematically formulated by Dubreil-Jacotin and Moiseev (see [17]) but their studies did not seem to be pursued by later people, although there have sporadically been some works such as [17,20,21] under the assumption of constant vorticity. Note however that those works were limited to the case where the vorticity is constant throughout the flow domain. Then Constantin and Strauss [6], which triggered a lot of mathematical researches on rotational water waves, considered water waves on general rotational flow from the viewpoint of the bifurcation theory. They proved in [6] the existence of the global branch of solutions for general vorticity distribution. They proved that a branch of nontrivial solutions exists and that it exists just until the solution approaches a flow having a stagnation point in it. In this sense they proved the existence of the 'global' branch. In doing so, they assumed some regularity on the distribution function of vorticity F = F( ) . Later the same authors published [7] where they allowed the discontinuity of the vortex distribution and proved a global branch of solutions. In the papers, however, the authors neglected the capillary force, which was taken into account by Martin and Matioc [14]. We use the formulation of [14] in the present paper.
Simmen and Saffman [20], Teles da Silva and Peregrine [21], and later [17] numerically computed water waves on the flow of constant vorticity of single layer. Recently Ribeiro-Junior and others [18] considered the same problem in more details and obtained many waves having more than one stagnation point. On the other hand, Ko and Strauss [10,11] numerically computed rotational water waves not only of constant vorticity but also of variable vorticity, including the case of two layers of different vorticity. They however neglected capillary force. Recent works [4,5] by Constantin, Kalimeris, and Scherzer are communicated to the authors by one of the reviewer. Those include interesting results for vortical flows of one-layer.
In what follows, we use the mathematical formulation of [14,15] to numerically compute the capillary-gravity water waves on rotational flows. New contributions of the present paper are: (1) capillary-gravity rotational waves are computed; (2) waves having stagnation points beneath the water surface are computed. Our works are closely related to [10,11]. But, since we included capillarity, we have studied solutions with more varieties of parameters and some of our water waves with stagnation points may be new. Existence/non-existence of stagnation points is important, but it is also important to see where they may appear-whether they appear near the bottom, or near the free surface. We will compute solutions with these questions in mind.
The present paper consists of 8 sections. In Sect. 2, we explain trivial solutions and bifurcation points around them. In Sect. 3, we explain our method of numerical computations of solutions without a stagnation point. In Sect. 4, we introduce an idea, due to [14,15], to compute solutions having a stagnation point. In Sect. 5 we present our numerical results for the pure capillary waves. The capillary-gravity waves are explained in Sect. 6 and the gravity waves in Sect. 7. The final section is devoted to discussions.

Trivial solution and bifurcation points
We begin our analysis by computing the 'trivial' solution. The trivial solution is, by definition, the one where all the streamlines including the free surface are straight horizontal lines. This can be phrased as: the functions u and U are independent of q. The trivial solution (u 0 , U 0 ) is explicitly written as follows: where we have used the following notation Hereafter, as in [14,15], we use instead of Q as the bifurcation parameter whenever is more convenient than Q. In the case of 1 = 0 or 2 = 0 , appropriate limit must be taken in the formulas above. For instance, if 2 = 0 , U 0 (p) becomes U 0 (p) = u 0 (p 1 ) + (p 1 − p)∕

√
. There is no a priori reason that > 0 . However, the formula (3.3) of [14] shows that > 0 is necessary if a non-trivial solution exists. Therefore we do not lose generality by assuming the positivity of .
Nontrivial solutions bifurcate from this trivial solution. The bifurcation point can be found by linearizing the equation at the trivial solution. We then see that the linearized equation has nontrivial solution if and only if the following equation in is satisfied: where we have set Equation (6) is nothing but the dispersion relation in the terminology of fluid mechanics under disguise of the bifurcation theory. We computed this complicated relation by Mathematica and compared the result with those in [14] to see the complete agreement.
Once the roots of (6) are obtained, we resort to a path-continuation method to compute bifurcating solutions.

Discretization
We discretize u and U by a combination of the finite difference method and the spectral method. The functions are expanded in q variable into a cosine series. Recall that the wave profiles are assumed to be symmetric about the y axis. This amounts to the evenness of u and U. To the variable p, we apply the finite difference method. Namely we divide the interval [p 1 , p 0 ] into N 1 subintervals of equal length. We then set The unknowns are now and are taken as bifurcation parameters. The differential operators in p variable are discretized by central differences and give us nonlinear equations, which we will solve by the path-continuation method. The path-continuation method is a standard method for numerically computing bifurcating solutions: see [17] or [1], for instance.
As we shall see soon, the computation of gravity waves is considerably more difficult than that of capillary-gravity waves. However, even capillary-gravity waves, where we do not have a sharp crest, require many grid points. We found that in some cases stagnation points appear in solutions of relatively small amplitude. Thus we have realized that a new method is needed to more fully understand the global picture of rotational waves.
In the following computations, we tested several values of N 1 , N 2 and M. If the solutions are close to the trivial solution, N 1 = N 2 = 16, M = 100 were enough. But when we traced the branch to compute solutions of larger amplitude, we needed larger numbers. For instance we used N 1 = N 2 = 16 and M = 200 .
Because of limitations of our computers, we were forced to use smaller N 1 = N 2 if larger M were needed. For instance we used N 1 = N 2 = 16 and M = 400 in some cases.

Formulation for waves with stagnation points
As we have noted, stagnation points appear rather soon and the formulation in Sect. 3 is no longer applicable. In this case we employ the following tactics: we compute, by the method above, the solutions near the trivial solution and trace the branch of the solutions. If the solution approaches a solution with stagnation point, we can recognize it by the occurrence of a flow region of slow speed. We then switch the numerical code to the following new one, which is adopted from [15]. The flow region is mapped from a rectangular region by the following mappings: , respectively. See Fig. 2. Note that h(x) denotes the free surface, and f(x) denotes the interface where the vorticity becomes discontinuous. Recall that d 1 and d 2 are related to p 0 and p 1 by (7).
Before using this mapping, the unknowns are u and U. After being mapped in this way, the unknowns are f , h, . The Poisson equations − △ = 1 and − △ = 2 are written in the following way. The boundary conditions arise from the requirement that the boundaries are streamlines, and are written as follows: In addition, is required to be continuously differentiable at the interface. Namely we require The Neumann boundary conditions are imposed on the side boundaries: Finally we must impose the Bernoulli condition, which is written in the following way: Applying a discretization just analogous to the one in Sect. 3, we can compute the water waves with stagnation points.

Results: pure capillary waves
Before we go into the detailed study of the capillary-gravity waves, we warm up ourselves by computing the pure capillary waves on the flow of a constant vorticity. 'Pure capillary' implies that the gravity is neglected and the capillary force acting on the free surface is the only force. Pure capillary waves on irrotational flows can be written explicitly by elementary function (the case of infinite depth-Crapper's waves) or by the elliptic functions (the case of finite depth-Kinnersley's waves), and accordingly can be used for the check of computer program. See, e.g., [9,17]. The case of pure capillary waves on rotational flows do not seem to have been considered in the references. Since the waves of extremely short wavelength (e.g., a few centimeters) may be approximated by this pure capillary waves, pure capillary waves may not be so odd things. We begin our analysis with this case.

Single-layer pure capillary waves
In this subsection we denote the depth and the vorticity by d and , respectively. Formulation and discretization conform to the way described in the previous sections. We used N 1 = 32, N 2 = 0 and M = 100, 200, 400 . (Note that N 2 = 0 since we consider flows of a single layer.) The Poisson equation and the Bernoulli condition are written as follows: The solutions of (10), (11) are drawn as the diagrams in Fig. 3. Here the trivial solutions lie on the -axis. Solutions with a stagnation point are drawn in red, and those without a stagnation point are drawn in green. Since there is a symmetry x ↦ x + , with which the governing equations are equivariant, we drew only one side of the branch. The solutions on the other side of the branch are obtained by the symmetry.
The solutions of the maximal a 1,0 in Fig. 3 are shown in Fig. 4. In both cases, we set = 4 and d = 1 . The figures on the top are perspective views (3D views) of the stream functions. The middle figures show streamlines and the lower ones do velocity vectors. When = −1 , a stagnation point appears in the upper half and vortex flow is clockwise as shown in (a). When = 1 , it appears in the lower half and vortex flow is counterclockwise as in (b). Stagnation points of two types are observed: elliptic and hyperbolic. The elliptic stagnation point appears on the y-axis inside the flow region, while the hyperbolic stagnation points appear on the free surface (Fig. 4a) or on the bottom (Fig. 4b). Figure 5 shows those cases of = 0.7, ±0.5 and −0.1 , while d and are the same as above. As | | becomes smaller, the stagnation region (those parts of flows consisting of closed streamlines) becomes smaller as is seen in Fig. 5a, c. At = 0.5, −0.1 , a stagnation point does not appear in any solution of either branch: See Fig. 5b, d. Experiments for some other (d, ) showed similar results. Namely no stagnation point appears if | | is small. The diminution of the region of closed streamlines complies with the fact that such a region does not exist in the case of irrotational flow, i.e., = 0 . See [17].
The solutions in Figs. 4 and 5 are those on the extreme points of the bifurcation diagrams-those from which we cannot trace the branch any further: The iteration process does not converge if we try to compute a solution of larger amplitude. The reason that solutions cannot be traced seems to be explained by the existence of a steep slope on the free boundary. In the case of the capillary water waves on irrotational flows, we do know those water waves whose free surfaces are partly overhanging (see [9,17]), namely y = h(x) is not single-valued. We suspect that such solutions exist in the present case. However, unfortunately, our formulation excludes such a solution. (See also Fig. 11 and its explanation for the same phenomenon.)

Two-layer pure capillary waves.
We now consider pure capillary waves on flows of two layers. There are five parameters , d 1 , d 2 , 1 and 2 . Therefore we gave up exhaustive experiments and computed only several cases. Since the bifurcation diagrams are qualitatively similar to those in the preceding subsection, we may safely omit them.  Fig. 7b. We now show another pattern of closed streamlines. Figure 8 shows the solution of d 1 = 0.5, d 2 = 0.1, = 4, 1 = −1, 2 = 1 and A 1,0 = 0.8227 : Here, unlike the previous cases, there exist two stagnation points appearing simultaneously just above and below the interface. Such phenomenon that two stagnation points exist simultaneously is interesting, but it occurs in a very narrow range of parameters in our computations. See also [18], where similar phenomena are presented.
In the same way, we carried out experiments in other parameter regions. In some parameter region, a solution with a stagnation point exists. In another region, a In the present section we consider the case where both capillarity and gravity are taken into account. We used N 1 = N 2 = 16 and M = 100, 200, 400 . A larger M was necessary to compute solutions of large amplitude. If both gravity and capillary forces are taken into account, there exist double bifurcation of mode m and mode n with different integers (wave numbers) m ≠ n , whence bifurcation structures become more complicated than those of pure capillary waves or the gravity waves. We must therefore pay attention not only to individual solutions but also to bifurcation structures among them. Because of the richness of the solutions, we concentrate ourselves to those solutions of mode 1 and mode 2, i.e., those where the coefficients of cos q and cos 2q are large and others are relatively small. In the upper-left diagram of Fig. 9, the horizontal axis represents , and the vertical one does 0.1 × A 1,0 + A 2,0 , where A 1,0 and A 2,0 are the Fourier coefficients of h defined in (9). The trivial solutions lie on the -axis. The upper-right diagram is a view from the -direction. There exist three primary bifurcation points of mode 1 and mode 2 on the horizontal -axis. The mode 1 branch, which emanates from a trivial solution at the rightmost point of the green line, extends to the left, forming a closed loop and finally joins the mode 2 branch, which is drawn in red. The mode 2 branch also forms a closed loop which emanates from a trivial solution and joins another trivial solution. The bottom figures are plots of streamlines of the solution at the points A, B, and C on the mode 1 branch. It shows a transition from mode 1 solutions to mode 2 solutions. These solutions have no stagnation point. Figure 10 shows the case where a stagnation point appears. Bifurcation diagrams are drawn in a way similar to Fig. 9. The bottom figures show streamlines of solutions at the points from A to E on the mode 1 branch. A stagnation point emerges somewhere between A and B (transition from green to blue). The stagnation region, i.e., the region enclosed by the closed streamlines, becomes larger as the solution approaches point D. If we go on further from D, then the stagnation region gradually becomes smaller and moves to the right as shown in E, which is a solution of the mode 2 branch. In this case, the mode 1 branch joins the mode 2 branch, too. The mode 2 branch forms a closed loop connecting a trivial solution and another trivial solution. The dotted part on the mode 2 branch shows where we could not The reason that we could not compute can be seen if we look at the solution in Fig. 11, where we plotted the wave profiles of the solution from which the numerical continuation of the solutions is no longer possible. We see that some parts of the free boundary are quite steep. If we trace the solution further we are sure that the overhanging wave profile or even self-intersecting waves would appear just as in the case of capillary-gravity waves on irrotational flows [17]. It is unfortunate that our formulation does not allow an overhanging wave profile: In order to obtain a complete bifurcation diagram, we need to invent a numerical method that permits both stagnation points and overhanging wave profiles. This is left for the future work. Note that all are incorporated if the problem involves one and only one value of vorticity. See [4,5,10,11,20,21]. It seems to us that the methods in these papers do not apply to the waves of overhanging profiles with a stagnations point on two layers.
In Fig. 12, the diagrams of mode 1 branches are drawn altogether for various values of 1 .

Gravity waves
This section is devoted to gravity waves, which are computed by the same method as in the previous sections. Setting = 0 in (6), the bifurcation point of mode n at the trivial solution is given by: Rotational gravity waves of finite depth on one layer were studied in [10,11,18,21] and some numerical examples of stagnation points are obtained. Our results are consistent with them too. Notwithstanding, in order to understand our results on twolayers, it would be better to examine our results on single-layer gravity waves,

Single-layer gravity waves
Numerical treatments are the same as in the previous sections. We used N 1 = 32, N 2 = 0 and M = 100, 200, 400 . We know that, in the case of gravity waves, stagnation points never appear when the vorticity is negative. Figure 14 shows the diagram of mode 1 branch of d = 0.5 , = 0.5, 1, and 2. As before, red/green parts represent solutions with/without a stagnation point, respectively. When = 0.5 , the mode 1 branch is rather short and no stagnation point appears. For smaller , the branch is shorter and no stagnation point exists. As increases, the branch becomes longer and a stagnation point appears in more solutions. Similarly to the case of capillary-gravity waves, the stagnation point appears in the lower side of the flow region. With d = 0.5 , streamlines for = −0.5, 1, 2, and 6.4 are shown in Fig. 15. As increases, the bulge of the closed streamlines around a stagnation point becomes more conspicuous.
In Fig. 16 the size of the closed streamlines increases then decreases. Here d = 0.8 and = 0.8 , and we plotted streamlines of three solutions, which are marked as A, B, and C in the bifurcation diagram. A stagnation point appears a little after the point A. The region of closed streamlines is the largest at B. If we trace the branch further, it gradually becomes smaller and almost on the verge of disappearing in C. Though we could not compute the solutions until it actually disappears, the disappearance seems to be very likely.

Two-layer gravity waves
Since we did not find a solution with a stagnation point if vorticity is negative, we experiment mainly in the cases where 1 > 0 & 2 = 0 , and, 1 = 0 & 2 > 0 .  Figure 17 shows the case of d 1 = d 2 = 0.5 and 2 is fixed to 0, the diagram of mode 1 branches are shown on the top. Plots of the highest wave solution of each branch are shown in the bottom, for 1 = 1, 1.5 and 2. When 1 = 1 , a stagnation point does not appear. The larger 1 is, the larger the closed streamlines are. Figure 18 shows the cases of different ratio of d 1 and d 2 , (A) d 1 ∶ d 2 = 1 ∶ 5 and (B) d 1 ∶ d 2 = 5 ∶ 1 . From these results, we do not see a correlation between the size of closed streamline and the ratio of d 1 to d 2 .
These and other experiments show that a stagnation point appears if 1 > 0 and 2 = 0 . On the other hand, it does not appear in any case where 1 = 0, 2 > 0 . Finally, Fig. 19 shows the results when 2 ≠ 0 for the example of Fig. 17c.

Discussions
The problem of water waves on two layers of piecewise constant vorticity involves many parameters, and it is not easy to classify the structure of bifurcation systematically. We accordingly chose some parameters, , (d 1 , d 2 ), ( 1 , 2 ) , which we hope to represent typical cases. It turned out that a stagnation point, if it exists, appears on the upper side of flow region with negative vorticity and on the lower side of the flow region with positive vorticity. On the other hand, whether it appears in the upper layer (the 2 -layer) or lower layer (the 1 -layer) depends on the combination of depth and vorticity. We could not find any simple law concerning their locations. It does not necessarily appear in the thicker layer. Nor does it always appear on the side with larger vorticity. When 2 = 0 and 1 < 0 , no stagnation point was found either in the gravity or the capillary-gravity waves. In other cases, their locations depend on the combination of (d 1 , d 2 ) and .
Although we succeeded in obtaining solutions with stagnation points, our method has a drawback: it does not allow us to compute wave of overhanging profiles. We need to come up with a new method for computing such solutions. If the vorticity is constant throughout the flow domain, then there are mathematical formulations which incorporate both overhanging profiles and stagnation points (see [3,8,13,14,17,18]). If vorticity is not a constant, we may have formulations which include stagnation points, like ours in this paper. Others may include overhanging waves but stagnation points are excluded. We do not know a mathematical formulation which includes all of (i) a stagnation point, (ii) overhanging part of the free surface, (iii) nonconstant vorticity. Note that our problem in this paper has two values of vorticity and is considered to be the case of nonconstant vorticity.
Ours seems to be the first attempt to numerically compute capillary-gravity water waves on two vortical layers with a stagnation point. Although we obtained new solutions, some problems such as the location of the stagnation point are left for future study.  Fig. 17c