{'fraction_non_alphanumeric': 0.030902321978092923, 'fraction_numerical': 0.0144182142933999, 'mean_word_length': 4.830614805520702, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Progressive flexural–gravity waves with constant vorticity

Abstract This paper is concerned with the interaction of vertically sheared currents with two-dimensional flexural–gravity waves in finite depth. A third-order Stokes expansion is carried out and fully nonlinear computations are performed for symmetric, steadily travelling waves on a linear shear current. For upstream periodic waves, two global bifurcation mechanisms are discovered. Both branches bifurcate from infinitesimal periodic waves, with one stopping at another infinitesimal wave of different phase speed, and the other terminating at a stationary configuration. Generalised solitary waves are found for downstream waves. More surprisingly, the central pulse of the generalised solitary wave can become wide and flat as the computational domain is enlarged. This provides strong evidence for the existence of wave fronts in single-layer free-surface waves. Particle trajectories and streamline structures are studied numerically for the fully nonlinear equations. Two patterns, closed orbits and pure horizontal transport, are observed for both periodic and solitary waves in moving frames. The most striking phenomenon is the existence of net vertical transport of particles beneath some solitary waves due to wave–current interactions. The streamline patterns alternate between net vertical transport and a closed orbit, resulting in the formation of a series of nested cat's-eye structures.


Introduction
There has been a long-standing scientific interest in flexural-gravity waves (also called hydroelastic waves in the literature) due to their importance for marine structures and sea transport. Flexural-gravity waves resulting from the interaction between moving fluids and deformable sheets have a wide range of applications in the polar regions where large floating ice sheets are used as roadways and landing strips. More recently, very large floating structures, such as floating airports (e.g. the Maga-Float project in Tokyo) and ultra-large merchant container vessels were thought to be environmentally friendly and self-sustained for converting ocean waves into energy. A better understanding of large-scale fluid-structure interactions, including flexural-gravity waves, is of fundamental importance in the process of building and utilising these engineering structures.
The study of flexural-gravity waves was initiated by Greenhill (1886Greenhill ( , 1916, who obtained the dispersion relation of this new type of wave and proposed the first practical application. Thereafter, flexural-gravity waves and their applications in marine engineering received growing attention from the scientific community. Research findings before the 1990s, most of which are based on linear theories, are summarised in the monograph by Squire et al. (1996). Observations of intense waves-in-ice events reported by Marko (2003) indicate that nonlinearity may play an important role. The first mathematical result on nonlinear flexural-gravity waves was the computation of large-amplitude periodic waves carried out by Forbes (1986) via a high-order series truncation method.
However, a key motivation to study nonlinear flexural-gravity waves is the generation of waves as a load moves on an ice cover. The linear theory shows that there exists a critical speed c min such that, if the speed of the moving load is close to c min , energy can hardly radiate away from the load. Although the linear theory identifies the critical speed, it fails to describe accurately the wave phenomenon near c min , since it predicts unlimited growth of wave amplitudes. Pȃrȃu & Dias (2002) first performed the weakly nonlinear normal-form analysis for both free and forced problems near c min , which showed the existence of envelope solitons in shallow fluids that are qualitatively similar to the experimental measurements carried out at Lake Saroma in Hokkaido (Takizawa 1988). However, their analysis and numerics cannot be generalised to the deep-water case, though a similar critical phenomenon was observed at McMurdo Sound in Antarctica by Squire et al. (1988). Milewski, Vanden-Broeck & Wang (2011) revisited the problem and showed numerically from the full Euler equations that, even though small-amplitude localised travelling wave solutions are not predicted to exist in deep water by standard perturbation analyses, they do occur along a new type of bifurcation branch.
All the aforementioned nonlinear analyses are based on the nonlinear Kirchhoff-Love plate theory, which has been widely used but does not have a clear conservation form for the elastic potential energy. Most recently, Toland (2007) proposed a model for plates based on the Cosserat theory of hyperelastic shells satisfying Kirchhoff's hypotheses, with the elastic energy being the total squared curvature. Since then, nonlinear flexural-gravity waves with the Toland elastic model have attracted intensive attention. Of interest are the following works: Guyenne & Pȃrȃu (2012, 2014 searched for hydroelastic solitary waves for the full Euler equations using the boundary integral method and performed unsteady simulations by truncating the Dirichlet-Neumann operator in arbitrary depth; Gao & Vanden-Broeck (2014) investigated the elevation generalised solitary waves in finite depth; Gao, Wang & Vanden-Broeck (2016) studied the stability and dynamics of solitary waves for the fully nonlinear equations via a time-dependent conformal mapping technique; and Trichtchenko et al. (2019) carried out the linear spectral analysis for periodic waves using the Fourier-Floquet-Hill method and compared the results with those obtained by a modulational instability analysis.
Satellite measurements of ice cover displacements induced by moving vehicles reported by van der Sanden & Short (2017) stress the need for continued efforts in research on three-dimensional fully localised flexural-gravity waves, known as lumps. The nonlinear elastic model in three dimensions was proposed by Plotnikov & Toland (2011) using the Willmore functional (namely the total squared mean curvature). Milewski & Wang (2013) derived the Benney-Roskes-Davey-Stewartson system in the vicinity of the minimum of the phase speed to predict the existence and elucidate the bifurcation mechanism of hydroelastic lumps. Recently, hydroelastic lumps were found numerically for the full Euler equations by Trichtchenko et al. (2018) using a boundary integral equation method.
The results mentioned above were obtained for irrotational flows. However, sea surface waves are commonly accompanied by underlying currents, and sometimes the current speed varies with depth (e.g. tidal currents and wind-driven currents). Early numerical works in this direction were carried out by Simmen & Saffman (1985), Teles Da Silva & Peregrine (1988), Milinazzo & Saffman (1990) and Vanden-Broeck (1994), who computed surface gravity waves for the fully nonlinear equations with constant vorticity. Recently, the cubic nonlinear Schrödinger equation was derived by Thomas, Kharif & Manna (2012) for pure gravity waves and by Hsu et al. (2018) for capillary-gravity waves to investigate the modulational instability of wave trains propagating on a linear shear current. Curtis, Carter & Kalisch (2018) expanded the primitive equation to the next asymptotic order to obtain the Dysthe equation with constant vorticity and investigated the motion and mean properties of particle paths. In addition, Hsu et al. (2016) extended the Stokes expansion to capillary-gravity waves and paid particular attention to the effect of vorticity on the phase velocity, wave profile and Wilton-type waves. On the theoretical side, wave-current interactions have received great attention since the pioneering work by Constantin & Strauss (2004) on local and global bifurcations of periodic gravity waves propagating on an arbitrary vorticity distribution. Subsequent research has focused on particle paths and flow structures beneath free-surface waves in the presence of vorticity. The interested reader is referred to Ehrnström & Villari (2008), Wahlén (2009) and Matioc (2014) and references therein for more details.
There are relatively fewer studies on the interaction between an underlying current and flexural-gravity waves. This is a special kind of wave-current-structure interaction. Peake (2001Peake ( , 2004 considered the nonlinear stability and the dynamics of a fluid-loaded elastic plate interacting with a mean flow using the method of multiple scales. He showed that the interaction gives interesting phenomena, including negative-energy waves and convective instability. Xia & Shen (2002) studied flexural-gravity waves in river channels in the presence of a mean flow, and derived the fifth-order Korteweg-de Vries (KdV) equation in the weakly nonlinear shallow-water regime. Bhattacharjee & Sahoo (2009) examined the effect of underlying shear currents on flexural-gravity waves in the linear shallow-water approximation. Wave scattering and trapping by jet-like shear currents were both analysed.
For hydroelastic waves with constant vorticity, a cubic nonlinear Schrödinger equation was derived by Gao, Wang & Milewski (2019), together with equations in the resonant cases. Fully nonlinear computations of solitary waves and the study of Benjamin-Feir instabilities were carried out to validate the weakly nonlinear models and to extend their results. Based on the same physical setting (a schematic is shown in figure 1), we focus in the present paper on asymptotics, global bifurcations, new steady solutions and particle trajectories in both periodic and solitary waves, as well as flow structures beneath solitary waves. The numerical computations used in this paper to solve the fully nonlinear equations rely on a conformal mapping method. This is similar to the approaches of Choi (2009), who studied the influence of a linear shear current on the Benjamin-Feir instability of a gravity wave train, and of Ribeiro, Milewski & Nachbin (2017), who investigated the flow structure as multiple stagnation points appear due to wave-current interactions.
The outline of the paper is as follows. The mathematical formulation of the problem is given in § 2. The results for symmetric, periodic and steadily travelling waves, including the third-order Stokes expansion, global bifurcations, generalised solitary waves and wave fronts, and particle trajectories, are presented in § 3. The numerical results for solitary waves are shown in § 4, with particular attention paid to particle paths and streamline patterns beneath the free surface. Finally, § 5 contains our conclusions and further remarks.

Governing equations
We consider a two-dimensional incompressible and inviscid fluid of finite depth h covered by an elastic sheet that provides a restoring force through its bending deformation. We introduce Cartesian coordinates with the x-axis along undisturbed elastic sheet and the z-axis directed vertically opposite to gravity. In non-perturbed states, the flow is assumed to be a shear current varying linearly in z, namely U(z) = U 0 + Ω 0 z, where U 0 and Ω 0 are both constants and Ω 0 is called the vorticity strength. There exists a frame of reference where the velocity vanishes at the undisturbed free surface; therefore, without loss of generality, we can let U 0 = 0 throughout the paper. Flow perturbations superimposed on the shear are assumed to be irrotational, with a potential function φ. Therefore, the governing equations of the problem are where η(x, t) is the elevation of the elastic sheet, ρ the density of the fluid, p the pressure at the top surface and g the acceleration due to gravity. We denote by ψ the streamfunction satisfying the Cauchy-Riemann equations ψ z = φ x and ψ x = −φ z . Following Toland (2007), the pressure across the elastic sheet is assumed to be where P a is atmospheric pressure, D the flexural rigidity, κ = η xx /(1 + η 2 x ) 3/2 the curvature of the sheet and s the arclength parameter with ∂ s = ∂ x / 1 + η 2 x . We study longitudinal progressive waves translating at a constant wave speed c. Waves become steady in a reference frame moving with the speed c, and the kinematic boundary condition can be written as The pressure equation at z = η can be recast as where B is the Bernoulli constant. It is noted that the bulk equation (2.1) and the bottom condition (2.4) are unchanged by changing the reference frame.

Periodic waves
3.1. Stokes expansion An exact solution of (2.1), (2.4), (2.6) and (2.7) is which is simply a uniform stream with velocity c and an undisturbed flat free surface. Non-trivial travelling waves can be obtained by perturbing the solution (3.1a-d). To achieve this, we introduce a small parameter , which is a measure of the amplitude of the wave, and write the expansions where B 0 = c 2 0 /2 and a precise definition of will be given later. This expansion was pioneered by Stokes (1847) for pure gravity waves, and now bears the name of the Stokes expansion. It was generalised to capillary-gravity waves by Wilton (1915) and to flexural-gravity waves by Vanden-Broeck & Pȃrȃu (2011). In the presence of a linear shear current, the Stokes expansion was carried out by Kishida & Sobey (1988) for gravity waves and by Hsu et al. (2016) for capillary-gravity waves.
The difficulty due to the unknown free surface can be overcome by writing the potential function on the free surface as a Taylor series, and expanding the kinematic and dynamic boundary conditions around z = 0. Substituting the various expansions into (2.1), (2.4), (2.6) and (2.7) and equating the powers of leads to a succession of linear systems. The linear system obtained at order reads (3.10) If we assume a periodic surface elevation with a fundamental wavenumber k, the solution to this system takes the form and c 1 will be determined at the next order in . At the second order in , we have the following sequence of equations: Eliminating the secular term yields c 1 = 0 and hence B 1 = c 0 c 1 = 0. Solving for other modes yields where (3.19) with It is noted that the Stokes expansion is valid provided that the denominators of A 22 and C 22 are non-zero. However, if then the expansion needs to be modified to include two modes: k and jk (the interested reader is referred to Vanden-Broeck & Pȃrȃu (2011) for more details). This was first achieved by Wilton (1915) for capillary-gravity waves, who provided evidence for the non-uniqueness of periodic water waves.
In the same vein, by collecting the terms of O( 3 ) we obtain Eliminating the secular term yields (3.28) Solving for non-resonant modes gives (3.31) (3.32) In addition, B 3 = c 0 c 3 and the next order of gives B 3 = c 3 = 0. Upon noticing that we denote by a the first Fourier coefficient of η(x), i.e. (3.37) Following Vanden-Broeck (2010), if we define the parameter as = a/λ, it then follows that A 11 = λ. In the discussion above, we carried out the Stokes expansion to the third order, which results in a correction to the linear speed of wave propagation since c 1 = 0 and c 2 / = 0. This procedure was previously applied by Kishida & Sobey (1988) and Hsu et al. (2016) in different contexts. As a check, we compare c 2 with the value from Hsu et al. (2016) when only gravity is considered (i.e. taking p = const. in (2.3)). It turns out that, under this circumstance, (3.28) can be reduced to where σ = tanh(kh) and we have used 1/k, Ω 0 √ gk and c 0 √ g/k to replace A 11 , Ω 0 and c 0 , respectively, for comparison purpose. Equation (3.38) is exactly the same as (3.32) in Hsu et al. (2016) when surface tension is neglected, providing a partial verification for our calculation.

Validation
We start with comparing the asymptotic results of the Stokes expansion with numerical solutions of the full Euler equations. To seek travelling waves for (2.1)-(2.4), we first non-dimensionalise the system by choosing as the units of length, time and potential, respectively. Therefore, the coefficients g and D/ρ are equal to 1 in (2.7), and Ω 0 in (2.6) and (2.7) can be replaced by a non-dimensional vorticity strength Ω defined as Following Gao et al. (2019), the problem can be handled by using a conformal transformation that maps the physical fluid domain to a strip in a complex plane. For travelling waves, after the transformation the unknown free surface can be parametrised by η(ξ ), which satisfies an integro-differential equation ξ is the Jacobian of the map, and the curvature κ in the new plane is of the form It is noted that the translating speed c and the Bernoulli constant B are also unknowns and need to be determined together with η(ξ ). The detailed derivation of (3.41) can be found, for example, in Gao et al. (2019). The pseudo-differential operator T is defined as where λ is the wavelength andh is the thickness of the fluid in the transformed plane defined ash We approximate η(ξ ) by the truncated Fourier series where a n are real and a n = a −n due to symmetry. We introduce collocation points uniformly distributed along the ξ -axis. This provides discrete algebraic equations by projecting (3.41) onto each Fourier mode. The equations are solved by Newton's method with the classical pseudo-spectral algorithm. To obtain solutions with high accuracy, The parameter space consists of four dimensionless parameters: the vorticity strength Ω, the fluid depth h, the wave speed c, and the wavelength λ (or, equivalently, the wavenumber k). To obtain more solutions or complete bifurcation curves, we use continuation methods where one previously computed solution is used as an initial guess to compute a new solution for slightly perturbed values of the parameters. The validity and accuracy of such schemes have been checked by several groups, and the interested reader is referred to Vanden-Broeck (1994), Choi (2009) and Gao et al. (2019) and references therein for more details.
We now use both numerical results for the full Euler equations and asymptotic predictions from the Stokes expansion to compare periodic travelling wave solutions. The asymptotic solutions can be obtained by substituting the expressions of φ i , η i , c i and B i into (3.2)-(3.5), choosing A 11 = λ, and varying the value of . Figure 2 illustrates the comparison between the asymptotic and numerical solutions for different values of . As expected, the difference between asymptotic predictions and numerical results increases as the wave steepness increases. As opposed to downstream waves (Ωc > 0, shown in figure 2a,c), steeper waves exist in the upstream case (Ωc < 0, shown in figure 2b,d), and the third-order approximation still works well for moderate-amplitude waves. It is observed that the downstream waves have a flat and wide crest while the upstream waves feature a narrow crest.
To second order, the Stokes expansion predicts that the wave speed is independent of the wave amplitude since c 1 = 0. At third order, the asymptotic expansion gives a correction c 2 to the linear phase velocity. In figure 3, wave speeds predicted by the third-order Stokes theory (dashed lines) and the fully nonlinear numerical computations (solid lines) are compared. Six solution branches for various values of Ω are presented. It is shown that the asymptotic wave speeds match the numerical results very well for smalland moderate-amplitude waves. We stop continuing the branches when Newton's method diverges or oscillates and fails to reach the desired accuracy. Wave profiles corresponding to the right endpoints of these curves are presented in figure 4, where the numerical results are shown as solid lines and the asymptotic predictions are shown as dashed lines. In general, the Stokes expansion works well for moderate-amplitude waves except for the case of zero vorticity, where the third-order solution is a poor approximation of the numerical solution, as shown in figure 4(c). If we check the denominator of A 22 , it is found that 1 + 16k 4 − 2 coth(2kh)c 2 0 ≈ −0.05 and the situation is very close to resonant harmonics or Wilton ripples, which explains the poor performance of the Stokes expansion in such a case.

Global bifurcation
In this subsection we study the global bifurcation of periodic waves of the fully nonlinear equations. To avoid taking into account the effects of all the parameters, we fix Ω, h and k, where Ω is chosen to be non-zero to include vorticity effects, and explore the wave speed-amplitude relationship. We start from the flat state and compute a branch of solutions until it returns to a trivial or stationary solution, termed a global bifurcation in the current paper. It is obvious that changing the signs of Ω and c simultaneously in (3.41)  results in the same equation, hence we only need to consider positive Ω. Furthermore, numerical evidence shows that the profiles and the bifurcations are much richer for upstream waves (waves propagating against the background shear). Therefore, we choose Ω > 0 and c < 0 in the following numerical calculations.
A small-amplitude monochromatic wave with propagating speed satisfying the linear dispersion relation is used as the initial guess for Newton's method. After iterating to a solution of the nonlinear integro-differential equation (3.41) within a desired tolerance, a continuation method is used to search for more solutions along the same branch by perturbing the previously computed solution by a small amount in some bifurcation parameter. The translating speed c, the value of the centre point of the free-surface displacement η(0), and the wave amplitude, which is defined as are used as bifurcation parameters to complete the bifurcation curves. It is noted that solutions with an overhanging structure are allowed in our computations due to the formulation of water waves in holomorphic coordinates. We stop the computation when the solution reaches the boundary of the speed-amplitude diagram (i.e. the wave becomes either a free stream or a stationary profile). An example of a global bifurcation branch of waves, which terminates in one endpoint at c = 0 and another at H = 0, is shown in figure 5. We start from small-amplitude waves and increase H to trace the branch. Then η(0) is used as the continuation parameter to traverse the very sharp turning point labelled 2 on the curve. Finally, we vary the wave speed to complete the bifurcation diagram. The most striking phenomenon observed in this example is the closing of the overhanging structure and its reopening. Akers et al. (2016) and Akers, Ambrose & Sulon (2017) numerically investigated the global bifurcation of interfacial hydroelastic waves based on the Birkhoff-Rott integral and arclength parametrisation. Their formulation for water-wave problems also allows multivalued wave profiles. However, they had to stop the computation at a limiting configuration where a self-intersecting point appears, enclosing a pendant-shaped bubble. In this situation, in contrast to their results, we can still continue the branch by changing some bifurcation parameter (speed or amplitude), taking advantage of conformal mapping, even though solutions with multiple intersecting points are non-physical (see 3 in figure 5).
A global bifurcation curve can connect two trivial solutions with different propagating speeds. It is shown in figure 6 that a branch of waves with a wavelength of π (dashed line) intersects with the 4π-period branch (solid line) at point 4 on the curve. Following the path 1 -4 and starting from point 1 , figure 7 provides a sequence of profiles that illustrate how a profile with four crests gradually emerges. We compute the solution with the period of π at the same speed (circles in figure 7d), which is exactly on top of the 4π-period solution (solid line in figure 7d), confirming our observations. It is also found in figure 6 that the 4π-period branch forms a closed loop. Waves on one half of the branch . Bifurcation diagrams of hydroelastic periodic waves with h = 5 and Ω = 1, but with different wavelengths. The 4π-period solution branch is shown by solid curves, while the π-period branch is plotted as a dashed curve. These two branches both bifurcate from zero amplitude (labelled with stars) and intersect each other at point 4 , where they have exactly the same wave profiles.
are the same as those on the other half but with a phase shift of 2π. One-and-a-half periods of solutions 3 and 5 are plotted together in figure 7(c) to clearly demonstrate the phenomenon of phase shift.

Generalised solitary waves and fronts
Generalised solitary waves are nonlinear non-periodic travelling waves with a central core similar to a classical solitary pulse and a non-decaying train of ripples extending up to infinity. Generalised solitary waves were previously computed by, among others, Hunter & Vanden-Broeck (1983) and Champneys, Vanden-Broeck & Lord (2002) for capillary-gravity waves and by Gao & Vanden-Broeck (2014) for flexural-gravity waves; all the computations were carried out in periodic domains (in other words, the generalised solitary waves were approximated by long periodic waves in numerics). Proofs of the existence of generalised capillary-gravity solitary waves were provided by Beale (1991) and others.
We now compute periodic waves with non-decaying oscillatory tails akin to generalised solitary waves for downstream flexural-gravity waves. We take a long small-amplitude cosine function (k ≈ 0.1 say) as the initial guess for the Newton-Raphson iteration and use the amplitude or speed as the bifurcation parameter. As the amplitude increases, the solution gradually approaches the configuration of a solitary pulse in the middle with several periodic waves in the tails (see figure 8). Further numerical experiments show that the algorithm appears to converge well if we add more and more oscillations to the obtained profile as the initial guess. It is shown in figure 8(a) that the profile obtained for . Typical wave profiles that correspond to 1 -5 shown in figure 6. The solid curves in (a-d) correspond to 1 -4 , respectively. Point 5 is shown by a dashed curve in (c), which is the same as 3 except for a phase shift of 2π. The π-period solution on the dashed curve at point 4 in figure 6 is shown by circles in (d), which exactly matches the 4π-period solution.
the wavelength λ = 96.96 is almost exactly on top of the profile for λ = 127.19, which provides strong evidence for the existence of true generalised solitary waves. On the other hand, when we use the wavelength as the bifurcation parameter, the broadening phenomenon of the central core is found. It is observed in figure 8(b) that the central core is flat and becomes broader as λ increases. This solution can even serve as a good approximation for a true solitary wave if we split the wave profile down the middle and glue the two endpoints together, considering the periodic nature of the computational domain (see figure 9a,b). As we enlarge the domain, the coexistence of the two phenomena, namely, the increase in the number of periodic waves in the tails and the broadening of the main core, implies the existence of wave fronts, which were previously only found in interfacial waves in multilayer fluid systems.
Wave fronts in hydrodynamics often occur in the flow of contiguous homogeneous fluids of different densities, which are usually called internal fronts. Interfacial gravity solitary waves under the rigid-lid approximation were computed by Turner & Vanden-Broeck (1988), and the most striking feature found in their work is the broadening of the wave, namely, the midsection of the interface develops a plateau, which becomes infinitely long when the wave speed approaches a limiting value. This numerical result provides evidence for fronts, since broad solitary waves can be viewed as the superposition of two fronts. Dias & Vanden-Broeck (2003) later showed that, in the limiting configuration, the flow in the far field and the flow in the middle can be referred to as parallel conjugate flows, and the wave indeed becomes a front. Fochesato, Dias & Grimshaw (2005) proposed a coupled KdV system for multilayer fluids to combine generalised solitary waves and fronts, and they found that ripples can appear on one side of the wave front. According to the above argument, we consider these solutions as solitary waves with many oscillations in the central core. In figure 9, the possibility of an increase in the number of ripples in the midsection of solitary waves is shown for both small-and moderate-amplitude solutions (see figure 9a,b). More obviously, if we consider half of the solutions, figure 9(c,d) provides strong evidence for the existence of wave fronts, since both flat tails and oscillations in the centre can be extended. The wave profiles shown in figures 8 and 9 are qualitatively similar to the travelling dispersive shock waves (TDSWs) found by Hoefer, Smyth & Sprenger (2019) and Sprenger & Hoefer (2017, 2020 in the fifth-order KdV equation and related models. As pointed out in their papers, TDSWs feature a partial non-monotonic solitary wave at the trailing edge connected with a periodic travelling wave train, and can be interpreted as nonlinear resonance between different types of nonlinear waves moving with the same speed. However, to the best of the authors' knowledge, this is the first time that wave fronts in the free-surface Euler equations have been reported.

Particle trajectories
In this subsection, we calculate particle trajectories numerically for the fully nonlinear equations and compare the results with asymptotic approximations. The particle paths under nonlinear and periodic water waves were initially considered by Stokes (1847) in his seminal work for irrotational flows and pure gravity waves. He showed that, in contrast to the linear theory, the particle path in the laboratory frame is not a closed loop but features a slight forward drift in the horizontal direction after a wave period has elapsed. This is known as the Stokes drift. A recent breakthrough on the theoretical side of this topic was made by Constantin (2006), who generalised Stokes' asymptotic work to steep waves to gain a qualitative understanding of the transport properties of arbitrary Stokes waves based on a rigorous mathematical argument. However, considering the boundary layers at the bottom and at the free surface, recirculation orbits (a phenomenon similar to the particle dynamics of a Gerstner wave) were shown to exist by the asymptotic work of Longuet-Higgins (1953) and the experimental research of Grue & Kolaas (2017). The existence of closed orbits was also proved by Constantin & Strauss (2010) when an underlying mean flow exists, which is intuitively understandable since the forward Stokes drift can be eliminated by the counter-propagating uniform current.
When a linear sheared current is added, theoretical studies on particle paths under periodic gravity waves have only been carried out for waves of arbitrarily small amplitude (Ehrnström & Villari 2008;Wahlén 2009). On the numerical side, Ribeiro et al. (2017) investigated particle trajectories under nonlinear periodic travelling waves with multiple stagnation points in a frame moving with the wave speed, and two dynamic behaviours of fluid particles were observed: periodic transport trajectory and closed orbit.
We start by describing the numerical scheme used to trace the trajectory of a fluid particle. In the frame of reference moving with the wave, particle trajectories can be obtained by solving the following ordinary differential equations: For the fully nonlinear equations, the physical space can be conformally mapped to the ξ -ζ plane and the equations in the new plane read Solving for ξ t and ζ t yields Equations (3.50a,b) can be integrated numerically by the fourth-order Runge-Kutta method and particle positions off the grid points can be obtained by interpolating the mesh grid data. Particle trajectories in a laboratory frame can be calculated by a simple Galilean transformation. For the first numerical computation, we compare particle paths resulting from calculations of the full Euler equations with those obtained from the asymptotic expansion. The particle paths based on the Stokes approximation are also computed numerically using a fourth-order Runge-Kutta method but in the physical space (3.47a,b). The velocity potential is given in (3.2) with appropriate φ i . Two examples of particle trajectories for upstream waves with = 0.02 (H = 0.25) and = 0.04 (H = 0.5) are shown in figure 10(a) and (b), respectively. Numerical simulations (solid curves) and theoretical predictions (dashed curves) show a very good agreement for small-and moderate-amplitude waves, partially demonstrating the validity of the numerical algorithm and confirming the asymptotic findings. It is noted that, in the shallow-water regime, long-wave models can also be used to reconstruct the velocity field beneath the free surface so as to compute particle trajectories (see Borluk & Kalisch (2012) for particle dynamics in the KdV approximation for irrotational gravity waves). Figure 11 shows numerical results of particle trajectories of a downstream wave in figure 11(a) and of an upstream wave in figure 11(b) in the moving frame. Two different patterns, periodic transport trajectories and closed loops, are both found. However, the figure illustrates the difference in the location of closed orbits between downstream and upstream waves. If we consider waves in the laboratory frame, this phenomenon actually indicates that, when Ω and c are of the same sign, the wave carries fluid particles beneath its crests, while, on the other hand, the upstream wave moves forward with fluid particles near the bottom. It is remarked that the frequency of the particle motion in the periodic transport regime is not uniform due to the shear current, which provides a non-uniform velocity distribution over water depth.

Solitary waves
The existence and stability of hydroelastic solitary waves propagating on a linear sheared current were investigated by Gao et al. (2019). Solitary waves were computed numerically and envelope equations near to and away from resonance were derived to assist stability analyses. We follow their numerical techniques, and focus on the flow structure beneath solitary waves. The algebraic decay of hydroelastic solitary waves usually requires a long computational domain with a large number of grid points. Long periodic waves with flat tails are usually considered to be a good approximation of solitary waves. Here we take λ = 200 as the domain side, and 4096 Fourier modes are used in most of the computations to achieve a sufficient accuracy.
Two fundamental branches of symmetric solitary waves, including one family of waves with a positive free-surface elevation at the centre (denoted waves of elevation) and the other family of waves with a negative free-surface elevation at the centre (denoted waves of depression) are shown in figure 12 for downstream waves and in figure 13 for upstream waves. All the branches presented bifurcate from the minimum of the phase speed shown by a vertical dashed line where the group velocity is equal to the phase velocity. The interested reader is referred to Gao et al. (2019) for more details on the bifurcation mechanism and its connection to the nonlinear Schrödinger equation. It is noteworthy that a smooth transition between upstream and downstream solitary waves along the same bifurcation curve is possible as the wave speed passes through zero (see figure 12a). Typical wave profiles are found to be similar to other wavepacket-type solitary waves (e.g. gravity-capillary waves in deep water). Hydroelastic solitary waves with constant vorticity are also characterised by oscillatory decaying tails. For hydroelastic solitary waves propagating against a non-uniform flow, particle trajectories show interesting patterns in the frame of reference moving with the wave, the most notable being the net displacement in the vertical direction indicating that the interaction between a solitary wave and a linear shear current can result in vertical mass transport (see figure 14). This is not observed in periodic waves. We take the particle on the left of the wave which experiences a net vertical displacement as an example. First the particle is chased by the solitary wave and swept downwards. Because the shear velocity increases with the water depth, the particle moves faster in deeper water and finally leaves the solitary wave far behind. It is worth noting that particles can also be swept upwards and never catch up with the solitary wave. In addition, figure 14 shows two other possibilities of particle trajectory patterns: closed orbits and pure horizontal transport. We remark that the calculations are carried out for large time periods. The uppermost curve without a circle is the free-surface displacement, while the others represent particle trajectories with different starting points denoted by circles. Three trajectory patterns are observed in the frame of reference moving with the surface wave: closed orbit, pure horizontal transport, and net vertical displacement.
For a better understanding of the structure of the flow field, we plot a set of streamlines which represent different particle paths in a time-independent system. We show the flow structure resulting from an upstream elevation solitary wave in figure 15(a) with parameters c = −2.95, h = 2 and Ω = 3. Only the active part of the horizontal domain near the main pulse is shown for better visibility. It turns out that, in the moving frame, the flow field can be divided into three layers in which streamlines and the dynamic behaviours of fluid particles are significantly different. For convenience, thick solid lines are used to show important streamlines separating different areas, while other streamlines are plotted as dashed curves.
In the layer between curve A and the free surface, particles overall keep moving from left to right and oscillate when they are swept by the solitary wave, and the path profile is in general qualitatively similar to the wave profile. In the region between curve B and the bottom, particles are less affected by the free surface and move from right to left due to the large horizontal speed of the shear current. In the region between curve A and curve B, particles can either move along closed orbits or move vertically when they are swept by the solitary wave but finally move in the opposite horizontal direction. Under each crest of surface oscillations, there is a family of closed trajectories which is bounded by another family of vertical-transport curves and eventually form a series of cat's-eyes nested from large to small. These families of closed orbits are located in a fluid layer where the shear speed is nearly equal to the wave speed, and the cat's-eye structure gradually shrinks when it stays away from the middle pulse of the solitary wave. Similar structures and trajectories can also be found in depression solitary waves, and a typical example is shown in figure 15(b) with parameters c = −3.03, h = 3 and Ω = 3. It is remarked that closed streamlines were also computed in weakly nonlinear models for rotational gravity waves at much lower computational costs (see, for example, the numerical studies on the Benjamin equation by Segal et al. (2017)). However,richer flow structures, such as nested cat's-eyes, can be expected in the fully nonlinear equations and beneath complicated wave profiles.
In the subsequent analyses, we explore the conditions under which the vertical-transport layer exists. Since there is no wave in the far field, the shear velocity should coincide with the wave speed in the horizontal centreline of the vertical-transport layer. It follows that the existence of vertical-transport zones requires a critical depth h c such that −Ωh c = c. This condition is twofold: For positive Ω, c should be negative and no vertical-transport zone exists if c < −Ωh, which provides the two boundaries shown in figure 16, namely c = 0 and c = −Ωh (dashed lines). On the other hand, for fixed Ω and h, there exists c min < 0 such that hydroelastic solitary waves of the wavepacket type can exist only for c > c min , which gives another boundary shown as solid lines in figure 16 for different values of h. In the Ω-c plane for fixed h, these boundaries form a semi-infinite region where vertical-transport zones exist (namely the right-hand side of the solid/dashed curve shown in figure 16). Figures 17-19 show how the vertical-transport layer varies when the parameter set approaches each boundary for h = 5. Since h c decreases along with the absolute value of c, the vertical-transport layer moves upwards. It is observed in figure 17(a) that the layer of right-going trajectories first disappears, and the closer c is to zero, the thinner the vertical-transport layer becomes (compare figure 17a and 17b). It turns out that the vertical-transport layer totally vanishes as c becomes zero. Figure 18 compares the flow structures between two parameter sets sitting on both sides of the dashed line, which correspond to the black dot and asterisk in the inset of figure 16. As h c increases, the vertical-transport layer moves towards the bottom (figure 18a) and completely disappears when h c exceeds the upper bound h ( figure 18b). Finally, figure 19 shows the trend of the vertical-transport layer as the amplitude of the solitary wave decreases (or, equivalently, c approaches c min ). The layer stays in the middle since h c ≈ −Ω/c min , while the thickness of the layer narrows as the surface wave decreases from 0.2 to 0.08 in amplitude. Whether the middle layer will completely disappear depends on the bifurcation mechanism of hydroelastic solitary waves. More precisely, if the associated nonlinear Schrödinger equation is of focusing type at c min , indicating that hydroelastic solitary waves bifurcate from infinitesimal periodic waves (Gao et al. 2019), then the vertical-transport layer will vanish as the wave speed reaches c min , but not vice versa.

Concluding remarks
In this paper the Stokes expansion up to third order has been carried out for flexural-gravity waves with a constant vorticity so that the nonlinearity manifests itself not only in the generation of higher-order harmonics but also in the correction of translating speeds. The full Euler equations were solved numerically using a conformal mapping technique, and travelling wave solutions, including periodic waves, bright solitary waves and generalised solitary waves, were computed. The Stokes expansion was used to validate the numerical algorithm by comparing periodic wave profiles, as well as particle trajectories, and very good agreements were found.
Further numerical calculations for the fully nonlinear equations focused on three topics: global bifurcation mechanisms of periodic waves, the existence of wave fronts, and flow structures beneath solitary waves. For upstream periodic waves, we showed that the global bifurcation includes a curve joining two infinitesimal periodic waves of different phase speeds, and a curve starting from an infinitesimal periodic wave and ending with a stationary state (c = 0). For downstream waves, the key finding of the broadening of the middle table-top structure of generalised solitary waves strongly suggests the existence of wave fronts characterised in the far field by a uniform state on one side and a train of waves on the other. To the best of the authors' knowledge, it is the first example of wave fronts discovered in the full Euler equations in single-layer fluid problems. For particle trajectories beneath solitary waves, in the frame of reference moving with the wave, three patterns, including pure horizontal transport, net vertical displacement and closed orbit, are possible due to wave-current interactions. For upstream waves, a nested cat's-eye structure of streamlines was observed for both elevation and depression solitary waves. Our numerical results raise further questions. A natural question is whether or not there are other global bifurcation mechanisms. On the theoretical side, the global bifurcation of pure gravity waves with arbitrary vorticity was initially investigated by The zone shrinks as the amplitude of the free surface decreases. Constantin & Strauss (2004). They showed three bifurcation mechanisms. Apart from two cases presented in § 3.3, an unbounded bifurcation is also a possibility. Akers et al. (2016) and Akers et al. (2017) investigated the global bifurcation of interfacial capillary-gravity waves and interfacial hydroelastic waves, respectively, using analytical and numerical tools. They provided some numerical evidence for the existence of unbounded bifurcation curves on which the wave amplitude increases without limit. Therefore, one can ask whether or not there are unbounded branches of hydroelastic periodic waves propagating on a linear shear current. The discovery of wave fronts also introduces questions. Since all the wave-front solutions were found in downstream waves in the present paper, a first question is if we can also find wave fronts for upstream waves. With respect to the asymptotic models for this phenomenon, the fifth-order KdV equation with non-convex dispersion, which admits TDSWs (Sprenger & Hoefer 2017), is a reduced model for flexural-gravity waves in the shallow-water regime for potential flows (Xia & Shen 2002). Therefore, it is also expected to be an appropriate model in the presence of a linear shear current so as to explain the wave-front phenomenon found in this paper. A comparative study of wave-front solutions between the asymptotic model and the primitive equations, as well as numerical simulations of the generalised Riemann problem for the full Euler equations, will be reported elsewhere in the near future. Another direction of extension related to TDSWs is to generalise the Whitham modulation equations (see Whitham (1974) for details) to waves with vorticity. In order to use Whitham's argument for slowly varying wave trains, the Lagrangian formulation of the problem is a prerequisite. Therefore, it is of great interest to find an explicit Lagrangian density for unsteady water waves with vorticity and extend Whitham's 'averaged variational principle' to these waves.