Flow dynamics for coupled surface and internal deep-water waves

This article determines the fluid motion underlying coupled linear internal and surface waves in a deep-water two-fluid-layer model (with the lower layer being of infinite depth). A detailed Eulerian description of the wave-field kinematics for coupled linear travelling waves is achieved using phase-plane analysis. The qualitative motion of individual fluid particles is elucidated through analysis of the relevant nonlinear dynamical systems from the Lagrangian viewpoint.


Introduction
Determining the underlying fluid motion generated by a wave propagating on an interface is an intriguing area of mathematical research which has important practical implications in the broad field of fluid mechanics. For surface water waves, it has been established in the past decade or so that fluid particle paths are uniformly nonclosed throughout the fluid domain, both in the (approximate) linearised setting [6,11,16,20,21,26,32], and for exact solutions of the fully nonlinear governing equations [2,4,7,19,22,23,31,35]. These results relate to periodic travelling surface waves which, from a theoretical perspective, already present substantial theoretical challenges (cf. the discussions and references in [3,5,30]).
The scale of these theoretical challenges is greatly increased when one also considers an additional (unknown) free interface. Internal water waves, which arise where there is 2 Linear wave solutions

Governing equations for coupled wave motion
We consider the two-dimensional motion of a stratified fluid, denoting horizontal and vertical coordinates by x and y, respectively. The fluid is assumed inviscid and incompressible, with an external restoration force due to gravity.
The physical regime consists of two vertically stratified fluid layers of differing (but constant) densities separated by a sharp internal interface denoted by y = (x, t) , which fluctuates about the mean water level y = 0 ; hence, ∫ ℝ (x, t)dx = 0 . The lower-fluid layer (denoted Ω l ) is of infinite extent, comprising the region −∞ < y < (x, t) , and the velocity field in the lower-fluid layer is expressed as (u, v). The upper-fluid layer (denoted by Ω u ) lies in the region (x, t) ≤ y ≤ h 1 + 1 (x, t) , where the a priori unknown free-surface boundary 1 represents fluctuations around the undisturbed surface water level y = h 1 , that is ∫ ℝ 1 (x, t)dx = 0 . The velocity field in the upper layer Ω u is denoted (u 1 , v 1 ) . Here, h 1 > 0 is a physical constant, which determines the mean depth of the upper-fluid layer. We assume that the wave profiles and 1 are such that and max | (x, t)| + max | 1 (x, t)| < h 1 , which precludes any intersection of the surface and internal wave interfaces. We assume stable stratification, with the upper layer being less dense than the lower layer, in which case we denote the density of the upper layer by and the lower layer by (1 + r) , where r > 0 is constant. In practice, r ≪ 1 ; for instance, in an oceanographical context, the value r = O(10 −3 ) may be taken as reasonable [9,27].
The equations of motion for an inviscid and incompressible fluid are the Euler equation, which is expressed in the lower layer Ω l by together with the equation of continuity The scalar function P(x, y, t) represents the internal fluid pressure, and g denotes the gravitational acceleration constant. In the upper-fluid layer Ω u , The kinematic boundary condition in the lower layer Ω l expresses the fact that the fluid is motionless at large depths: The dynamic and kinematic boundary conditions at the internal interface take the form The dynamic boundary condition (3b) ensures that the pressure is always continuous throughout a fluid. The kinematic boundary conditions (3c) and (3d) ensure that the normal components of the respective velocity fields match, and are continuous, at the interface. For inviscid fluid motion, this need not be true for the tangential velocity components. Finally, at the free surface, the governing equations (2) in the upper layer Ω u have the associated dynamic and kinematic boundary conditions The fluid is assumed to be irrotational in each fluid layer separately which, in two dimensions, corresponds to (1a) u t + uu x + vu y = − P x (1 + r) , v t + uv x + vv y = − P y (1 + r) − g, (3a) (u, v) → 0 as y → −∞.
This implies the existence of velocity potentials (x, y, t) and 1 (x, y, t) for which The governing equations (1) and (2) and boundary conditions (3) can be reformulated in terms of the velocity potentials as follows. It follows from definition (5), coupled with the incompressibility equations (1b) and (2b), that the velocity potentials are harmonic functions: Furthermore, formulating the kinematic boundary conditions (3a), (3c), (3d) and (3f) in terms of the velocity potentials leads to where n = (− � (x), 1)) is the normal vector to the interface exterior to the lower domain, and n 1 = (− � 1 (x), 1)) is the exterior normal vector to the free surface. The pressure-matching dynamic boundary condition (3b) formulated in terms of velocity potentials becomes while the dynamic surface condition (3e) becomes the Bernoulli condition

Linearised equations
The governing equations and boundary conditions can be linearised by invoking nondimensionalisation and scaling procedures (cf. [25] for full details in the setting where the lower-fluid layer has finite depth). Let denote a characteristic wavelength for the water waves being considered, and let be a characteristic depth scale (the choice = h 1 is natural in the present context). Let a be a characteristic amplitude of the internal wave, with a 1 a characteristic amplitude for the surface wave. Then, the nondimensional parameters = a∕ , 1 = a 1 ∕ (which measure the magnitude of the wave amplitudes relative to the characteristic vertical depth scale) and = ∕ (a 'shallowness' parameter which measures the magnitude of the characteristic vertical depth scale relative to the wavelength) are naturally introduced into the governing equations and boundary conditions. Furthermore, linearisation requires that the wave-steepness parameters = ak = 2 (a∕ ) (for the internal wave) and 1 = a 1 k = 2 (a 1 ∕ ) (for the surface wave) are small. The linear wave approximation is valid under the assumption that these nondimensional parameters satisfy The linearisation process eliminates all product terms in the governing equations (1) and (2) and boundary conditions (3). Additionally, an important consequence of linearisation is that the boundary conditions (3) are now evaluated at the constant mean levels y = 0, h 1 , as opposed to on the unknown interfaces y = , h 1 + 1 , respectively. The linearisation of (6) results in the following Neumann boundary value problem for the velocity potential in Ω l , while in the upper-fluid layer Ω u the velocity potential 1 must solve the Neumann boundary value problem The linearised dynamic boundary conditions (6f) and (6e) take the form and

Linear wave solutions
Seeking periodic travelling wave solutions of (8) which have a functional x, t dependence of the form kx − t , where is the wave frequency, k = 2 ∕ is the wavenumber, and c = ∕k is the wave phase speed, suggests the Ansatz: Here, a 1 and a are the amplitudes of the free surface and interface, respectively. Unless otherwise stated, in the following we assume that a, a 1 ≠ 0 , thereby implying a nontrivial coupling of wave motions at the free surface and the interface. The physical set-up illustrated in Fig. 1 represents waves where the crests (and the troughs) of the surface and internal interfaces coincide, which corresponds to the ratio a∕a 1 > 0 in (9). Ansatz (9) permits (7) , 1 , , 1 ≪ 1.
an alternative configuration, whereby the crests (respectively, the troughs) of the surface wave coincide with the troughs (respectively, the crests) of the internal wave. This scenario arises when a∕a 1 < 0 and is represented schematically in Fig. 2.
The coupled waves represented in Fig. 1, with a∕a 1 > 0 , will be referred to as being 'in-phase', whereas the coupled waves represented in Fig. 2, with a∕a 1 < 0 , will be referred to as being 'out-of-phase'. We note that |a| + |a 1 | < h 1 must hold for out-ofphase waves. It will be established in relation (22) that there exists precisely one set of coupled wave solutions which are in-phase, and one set which are out-of-phase, for a given fixed wavelength (or, alternatively, fixed frequency). Implementing Ansatz (9b) in (8b), system (8a)-(8c) has the solution  Invoking Ansatz (9b) in (8e) and (9a) in (8f), system (8d)-(8f) can be solved for where the nondimensional 'wave amplitude ratio' parameter A is given by It is convenient to re-express (12) in terms of the amplitude ratio The value of A is determined by two nondimensional parameters, namely the ratio of wave amplitudes a∕a 1 and the shallowness parameter kh 1 = 2 ⋅ h 1 ∕ .
Remark Note that the limiting case a 1 → 0 (no surface wave) in the upper-fluid layer corresponds to the classical 'rigid-lid' model of internal wave motion, with the associated velocity potential The boundary conditions (8b) and (8f) ensure that the velocity potentials and 1 prescribe matching linearised normal velocities at the internal interface. Evaluating (8g) for the wave surface profile (9b) and velocity potential 1 (11) leads to a dispersion relation at the free surface Note that this is not a dispersion relation for the surface wave alone: there is an intrinsic coupling between the surface and internal waves induced by the wave amplitude ratio a∕a 1 . An immediate consequence of (15) is that the parameter A must be positive, which confers a bound on the internal wave amplitude for in-phase waves, namely: Evaluating (8h) for the velocity potentials (10) and (11) leads to the following dispersion relation at the interface: As remarked above, this is not a dispersion relation for the interface alone: this relation prescribes coupling between surface and internal waves through terms involving the ratio of wave amplitudes.

Dispersion relations for the coupled waves
Dispersion relations are formulae which specify the linear wavespeed c in terms of various physical parameters. They are so-called since, if the wavespeed c varies with respect to some parameter, then the waves are dispersive: waves corresponding to different parameter values will travel at different speeds. Assuming we know (either by measurement, or prescription) the mean depth h 1 of the upper-fluid layer, Eqs. (15) and (16) feature three unknown parameters, namely the wavenumber k = 2 ∕ , the wavespeed c = ∕k (or, alternatively, the frequency ) and the ratio of wave amplitudes a∕a 1 . Using (13), we can re-express (15) as and (16) as another quadratic for c of the form: The evident coupling between relations (17) and (18) illustrates that the wavespeed c depends on the ratio of wave amplitudes, a∕a 1 , and vice versa. Together, Eqs. (17) and (18) constitute a system of dispersion relations prescribing the wavespeed c of coupled internal and surface waves, and the ratio a∕a 1 of associated wave amplitudes, in terms of the wavenumber k. To obtain a single dispersion relation involving the wavespeed c, we substitute (17) into (18) to obtain the classical dispersion relation [28,36]: The nondimensional wave amplitude ratio parameter A plays a key role in subsequent analysis of the underlying fluid motion induced by coupled waves. Hence, it is beneficial for our purposes to reformulate relation (19) as a dispersion relation in terms of A by substituting kc 2 ∕g = 1∕A: This constitutes a dispersion relation for the ratio of wave amplitudes since, given solutions of (20), the ratio of the corresponding wave amplitudes a∕a 1 can be established directly from (13) (in terms of k). Denoting the polynomial (20) as P(A) = 0 , one root is A 1 = 1 (since P(1) = 0 ) and the other is Typically r ≪ 1 , with r = O(10 −3 ) constituting a reasonable value in the ocean [9,27], in which case A 2 ≫ 1 by (21). Relation (13) implies from which we conclude that Coupled wave motions at the interface and free surface are in-phase for solutions corresponding to the root A = A 1 = 1 , and out-of-phase for A = A 2 .
Remark Tangential fluid velocities at the internal interface are continuous for the coupled wave solutions corresponding to A = A 1 = 1 , and discontinuous for A = A 2 . This can be seen by comparing x and 1,x (corresponding to linearised tangential velocities in the respective layers) at y = 0 . By (10) and (11), these match only if a = a 1 − sinh kh 1 + A cosh kh 1 which, by (16), holds only if A = 1 . For coupled wave solutions corresponding to A = A 2 , although the fluid is irrotational in each fluid layer separately, there must be a vortex sheet located at the internal interface. In practice, viscosity (neglected in this model) acts to blur the sheet into a vortex film.
The roots of (19), the dispersion relation for the wavespeeds, can be determined from c 2 = g∕kA to get Each root c 2 i in (23) represents a given wavespeed, with the choice of sign +c i (or −c i ) corresponding to right-moving (or left-moving) coupled waves, respectively. It follows from (23) that c 1 > c 2 > 0 ; indeed, c 2 ≪ c 1 for r ≪ 1 . Solutions corresponding to the wavespeeds c 1 and c 2 are referred to as barotropic and baroclinic, respectively, cf. [36]. In the ocean, it is observed that internal waves are typically much slower than surface waves, with significantly greater amplitudes. This is due to the restoring force at the internal interface being substantially less than at the free surface. Internal waves in the ocean can have periods ranging from tens of minutes to several hours, with wavelengths ranging from hundreds of metres to tens of kilometres, and their height can often exceed 50 m (cf. [17,18,29,36]). In contrast, for ocean surface gravity waves the period ranges from 1 to 25 s, with ocean swell having a typical wavelength that is greater than 260 m (up to a maximum of approximately 900 m) with a period larger than 13 s (up to a maximum of 24 s) (cf. [37]). We expect these characteristic properties of surface and internal waves to be reflected in the solutions of the linear dispersion relations above. Since A 1 = 1 for the larger wavespeed c 1 , from (13) we have and so a 1 > a : the surface wave is more prominent than the internal wave for the coupled solution corresponding to the faster wavespeed c 1 . This distinction is more pronounced for shorter wavelengths (larger k) and greater upper-layer depths h 1 . Taking indicative values h 1 = 120 m for the equatorial Pacific thermocline depth (cf. [9]), and k = 2 ∕300m −1 , relation (24) gives a ≈ 0.08 a 1 . For coupled waves propagating with the slower wavespeed c 2 , Hence, the internal wave is (much) more prominent than the surface wave for the coupled wave solution corresponding to the slower wavespeed c 2 . Taking the indicative values h 1 = 120 m, k = 2 ∕5km −1 and r = 2 × 10 −3 , relation (25) gives |a| ≈ 508 |a 1 | . Of course, relations (24) and (25) come with the caveat that they relate only to linear coupled waves satisfying conditions (7).

Remark
Letting r → 0 simplifies the physical model, which now consists of one homogeneous fluid layer of infinite depth. Taking r → 0 gives c 2 = 0 and c 2 1 = g∕k : this is the classical dispersion relation for a surface wave over a fluid of infinite depth.

Remark
The 'rigid-lid' model describes wave motion at the interface separating two fluid layers which are bounded above a by rigid horizontal wall. If the mean thickness of the upper layer is h 1 , and the lower layer is infinitely deep, the dispersion relation for the rigidlid model takes the form Relation (26) can be formally derived from (18) by letting a 1 → 0 (with a ≠ 0 ). Alternatively, it arises from ensuring continuity at the internal interface of normal velocities prescribed by velocity potentials (10) and (14).
Remark When the depth of the upper-fluid layers becomes very large ( h 1 → ∞ ), the influence of the upper boundary diminishes, and the dispersion relation (26) for the interface becomes: This matches the wavespeed c 2 obtained in (23) in the limit kh 1 → ∞.

Dynamical systems formulation
If (x(t), y(t)) is the path of a particle in the lower-fluid layer Ω l , then the motion of the particle is described by the nonlinear dynamical system for −h < y < 0 , with initial data (x 0 , y 0 ) . In the upper-fluid layer Ω u , particle trajectories (x(t), y(t)) are determined by the nonlinear dynamical system for 0 < y < h 1 , with initial data (x 0 , y 0 ) . The mean level of the oscillating internal wave interface y = is located at y = 0 , whereas the free surface y = h 1 + 1 oscillates about the mean level located at y = h 1 . The right-hand sides of the differential systems (27) and (28) are smooth; therefore, the existence of unique local smooth solutions for both (27) and (26) .
is ensured by the Picard-Lindelöf theorem [33]. Furthermore, since the right-hand sides of (27) and (28) are bounded in the respective fluid layers, these unique solutions are defined globally [33]. The right-hand sides of both (27) and (28) are nonlinear, and thus, such systems cannot be solved explicitly. The fluid layers are separated by an impermeable interface y = (x, t) ; therefore, we perform a phase-plane analysis of system (27) in the lower-fluid layer Ω l and system (28) in the upper-fluid layer Ω u , separately, and then piece together the information to obtain the motion of the entire two-layer body.
As we are dealing with travelling waves it is possible to transform to moving frames where the motion is steady by the changes of variables The mapping (29) transforms system (27) to the autonomous system with initial data (X(0), Y(0)) = (kx 0 , ky 0 ) . We denote M ∶= ak = ≪ , where = ak ≪ 1 is the wave-steepness parameter for the internal wave (cf. Sect. 2.2). Note that (31) can be expressed as a Hamiltonian system for the Hamiltonian function If (X(t), Y(t)) is a solution of (31), then dt dt H(X(t), Y(t)) = H XẊ + H YẎ = 0 , and so, H is constant along the phase curves. The mapping (30) transforms (28) to the autonomous system (12) and the parameter M 1 ∶= a 1 k = 1 ≪ (since the wave-steepness parameter for the surface wave satisfies 1 ≪ 1 : cf. Sect. 2.2). System (33) also has a Hamiltonian structure for the Hamiltonian function Since (31) and (33) (30) reflects vertical coordinates through the line y = h 1 with the effect that, when constructing phase-plane diagrams in terms of the transformed (X, Y 1 ) variables, wave crests in the physical system (28) correspond to troughs in the transformed system (33), with a similar correspondence between troughs for (28) and crests for (33). Additionally, the streamline denoting the surface wave lies beneath that representing the internal wave in the resulting phase portraits.
In subsequent phase portrait analysis of the dynamical systems (27) and (28), the parameters k and will be regarded as fixed constants, with the nondimensional wave amplitude ratio parameter A taking the values A 1 = 1 and A 2 > 1 for the upper-fluid layer system (28). While these parameters are fixed for a given coupled wave motion, in general they are not free-parameters; rather, they are determined by the dispersion relations (15) and (16).

Phase portrait analysis: lower-fluid layer
The velocity field for the lower-fluid layer is qualitatively identical to one which describes fluid motion in a single homogenous (uniform density) fluid layer whose upper interface separates the fluid from a source of constant pressure (such as the atmosphere), cf. [6]. Indeed, we note phase-plane approaches have previously proven successful in revealing the underlying flow-structure of a variety of surface water waves (cf. [6,11,15,16,20,21]. The physical influence of the upper-fluid layer is conveyed implicitly by way of the dispersion relations (17) and (18). The autonomous system (31) meets standard regularity assumptions for the uniqueness of the Cauchy problem [33]; therefore, its trajectories do not intersect. We note that the right-hand side of (31a) is an even function in X, while the right-hand side of (31b) is an odd function in X; therefore, the trajectories of (31) have a mirror symmetry with respect to the X−axis. Without loss of generality, we choose a > 0 (hence M > 0 ) throughout this subsection.
The 0-isocline is defined to be the set where dY∕dt = 0 , and the ∞-isocline is the set where dX∕dt = 0 . Therefore, the 0-isocline is given by the lines X = 0, ± . The ∞-isocline is given in the region X ∈ − 2 , 2 by the curve (X, (X)) , where is defined by (X) = ln ( ∕M cos X) : (X) ∈ [Y * , ∞) for Y * = ln ( ∕M) = ln (1∕ ) , where = ak is the wave-steepness parameter. The even function is smooth, it takes on its infimum Y * at X = 0 , and it satisfies the limiting condition lim X→± ∕2 (X) = ∞ . The only singular point of system (31) ln (1∕ )) . Note a unique positive solution > 0 exists since ≪ 1 . This implies that the streamline that represents the internal wave lies beneath the separatrix, as required, since the wave trough is located at X = and Y = − < . For X ∈ ( 2 , ) , we have dX∕dt < 0, dY∕dt > 0 . If X ∈ (0, 2 ) , then dX∕dt < 0 below the curve of (X) and is positive above it, while dY∕dt remains positive in this region. The corresponding signs for X ∈ (− , 0) are obtained using symmetry with respect to the Y-axis. A phase portrait for the lower-fluid layer is given in Fig. 3.

Phase portrait analysis: upper-fluid layer
Phase-plane analysis for the upper-fluid layer proceeds along the lines of that for coupled waves propagating on a finite fluid domain over a flat bed [25], with simplifications resulting in the present setting due to the explicit form of the wave solution parameters A = A 1 , A 2 (given by (21)). For completeness, we present the salient aspects of the phase-plane analysis which are germane to the present context, referring to [25] for full details.
The autonomous system (33) meets the standard regularity assumptions for the uniqueness of the Cauchy problem [33]; therefore, its trajectories do not intersect. Moreover, since F(X, Y 1 ) is an even function and G(X, Y 1 ) an odd function, with respect to X, any trajectory of system (33) is symmetric with respect to the Y 1 -axis when viewed as a curve in the (X, Y 1 )− phase plane. Re-express the right-hand sides of (33) as for the functions Without loss of generality, we fix M 1 > 0 ( a 1 > 0 ) throughout subsequent considerations, with the sign of a matching that of the ratio a∕a 1 as prescribed by A in (22), namely: In order to investigate the phase portrait of system (33), we must begin by investigating the 0-isocline (defined as the set of points where the vector field is horizontal, Ẏ 1 = 0 , and so G(X, Y 1 ) = 0 in (35a)) and the ∞-isocline (defined as the set of points where the vector field is vertical, Ẋ = 0 , and so F(X, Y 1 ) = 0 in (35a)).

System (33) with A = A 1 = 1
For A = 1 , relations (35b) reduce to f (Y 1 ) = e −Y 1 = −g(Y 1 ) . Accordingly, the phase portrait for positive values of Y 1 rapidly converges to a series of flat, horizontal lines as Y 1 increases, since the velocity field converges exponentially fast to the uniform system Ẋ ≡ − and Ẏ 1 ≡ 0 . We have G(X, Y 1 ) < 0 for 0 < X < and vanishes at X = 0 and X = ; hence, these vertical lines comprise the 0-isocline for system (33) in [0, ] . The ∞ -isocline consists of points Y 1 where cos(X)e −Y 1 = 1∕ 1 , and for X ∈ 2 , , there is no such solution. For each X ∈ 0, 2 , there is a solution, and for X = 0 , we have e −Y * is a singular point, and since the Hessian of the Hamiltonian (34) is it follows immediately that Q 1 is a saddle point which lies at the confluence of four separatrices. The maximum value of the ∞-isocline is located at Q 1 , and it decreases monotonically to −∞ as X → 2 . The phase portrait for system (33) when A = 1 is given in Fig. 4.
Note the magnitude of the amplitude of the surface wave is greater than that of the internal wave, in accordance with (24). The surface wave streamline must be located above the separatrix for physically relevant solutions of (33), and this is the case if − 1 > Y * 1 , that is, 1 e 1 < 1 : this holds since 1 ≪ 1.
For 2 ≤ X ≤ , F(X, Y 1 ) is strictly positive since M 1 > 0 ; hence, there are no ∞-isocline or singular points in this region. In the region 0 ≤ X < 2 , the ∞-isocline is given by points (X, Y 1 ) where cos(X)f (Y 1 ) = ∕M 1 = 1∕ 1 ; therefore, since the left-hand side is maximised at X = 0 , an ∞-isocline exists in this region if f (Y * 1 ) = 1∕ 1 for some value of Y * 1 . It follows from the schematic in Fig. 5 We denote these singular points of system (33) by Q 1 = ( ,Ỹ * 1 ) , Q 1 = ( , Y * 1 ) , and let us first examine the point Q 1 . By virtue of the properties of f (Y 1 ) , the ∞-isocline is given in the interval X ∈ 0, 2 by 1 (X) = +∞ and, in the interval 0, 2 , 1 is an increasing function. Given the Hamiltonian (34) for system (33), the Hessian at Q 1 is and arguments similar to those above involving Morse theory imply that Q 1 is a saddle point which lies at the intersection of four separatrices: two that reach Q 1 in infinite time in the future and two that need infinite time backwards to reach the saddle point. Analogous reasoning shows that the singular point Q 1 is also a saddle point, and the ∞-isocline emanating from Q 1 is qualitatively similar to that described for Q 1 , except 'flipped' vertically. Having characterised the 0-and ∞-isoclines, we can easily determine the signs of the two components F(X, Y 1 ) and G(X, Y 1 ) of the vector field given by system (33), and employing symmetry properties of (33), we infer the phase portrait behaviour in the region − ≤ X ≤ 0 . The complete phase portrait of system (33) is given in Fig. 6.
The surface wave profile has a mean water level Y 1 = 0 , which corresponds to y = h 1 , while the internal wave profile has a mean water level Y = kh 1 , which corresponds to y = 0 . From the phase portrait in Fig. 6, we see that there are qualitatively different fluid motions possible in the upper-fluid layer. Streamlines located beneath Y 1 =Ȳ 1 in Fig. 6 are in-phase with the surface wave; the horizontal line Y 1 =Ȳ 1 is itself a (flat) streamline; streamlines located above Ȳ 1 in Fig. 6 are out-of-phase with the surface wave: the internal wave with mean water level Y = kh 1 lies in this region. Note that the magnitude of the internal wave amplitude is significantly greater than that of the surface wave, as prescribed by (25). In order for system (33) to describe physically relevant solutions, the singular points Q 1 and Q 1 must lie outside the upper-fluid layer; hence, we must have Ỹ * 1 < − 1 and kh 1 + < Y * 1 or, equivalently (cf. Fig. 5), max f (− 1 ), f (kh 1 + ) < 1∕ 1 . Since , 1 ≪ 1 in the linear regime, and f (− 1 ) ≈ A + 1 , for this condition to breakdown requires A to be sufficiently large that 1 f (− 1 ) ≈ 2 1 + A 1 ≈ A 1 > 1 . Otherwise, the singular points lie outside the upper-fluid layer.

Particle trajectories
Phase portraits have been constructed for the nonlinear dynamical systems (31) and (33), which describe fluid motion in the lower-and upper-fluid layers, respectively. This analysis is performed in terms of the (X, Y)− and (X, Y 1 )-variables, in (moving) reference frames for which the flow is steady. For steady fluid motion, the dynamical systems (31) and (33) are autonomous and, hence, particle trajectories coincide with the streamlines. Accordingly, a complete qualitative picture of fluid motion may be ascertained in terms of the transformed coordinates. From the phase portraits in Figs. 3, 4, 6, it is also possible to ascertain which streamlines correspond to physically admissible fluid trajectories in each fluid layer. This global picture of the qualitative behaviour of the dynamical systems (31) and (33) provides an Eulerian description of the fluid motion.
In this section, we pursue a Lagrangian description of fluid motion by ascertaining qualitative features of the motion of specific fluid particles described in terms of the physical variables (x(t), y(t)) . We prove that there does not exist any closed particle trajectory in either the lower, or upper, fluid layers, and prove that all particles experience a forward drift. Regarding particle drift, it bears remarking that the Stokes' drift phenomenon (whereby fluid particles experience a mean net drift velocity in the direction of wave motion) is intrinsically nonlinear (cf. [24,38]): the existence of a forward drift induced by linear internal waves for all fluid particles, first established in [25], is apparently new. Finally, we establish monotonicity properties for this forward drift.

Lower-fluid layer
The significant difference between systems (27) and (31) arises due to the latter being autonomous, while the former is nonautonomous and hence considerably more difficult to analyse. Nevertheless, it is possible to use qualitative properties of the fluid motion induced by (31), established in the phase portrait analysis of Sect. 3.1, to construct a qualitative description of fluid motion in terms of physical variables as prescribed by (27). System (27) is qualitatively similar to one which describes the motion of a surface wave propagating on a single homogeneous fluid layer of finite depth, cf. [6]. Physical variables are obtained by reversing the coordinate transformations (29) by defining in the lower-fluid layer. Suppose (X(t), Y(t)) describes a streamline in the lower-fluid layer (in the moving reference frame) such that (X(0), Y(0)) = ( , Y 0 ) : in the lower-fluid layer focus is restricted to streamlines for which Y 0 ∈ (−∞, k(h − a)] . Let t Y 0 (− ) denote the time it takes for the particle to intersect the line X = − . The next lemma is stated for motion in the lower-fluid layer, but applies similarly to fluid motion in the upper layer.

Lemma 4.1 If the particle trajectory prescribed by (x(t), y(t)) is a closed path with period
, then, necessarily, we have = 2 . Conversely, suppose t Y 0 (− ) = 2 , then the particle path prescribed by (x(t), y(t)) is closed.
Proof The proof follows from the periodicity of system (31) with respect to X, together with the definition (37). ◻ The main result concerning motion in the lower-fluid layer is stated in the following. Proof The proof follows, bearing in mind Lemma 4.1, by proving that For such Y 0 , along the streamlines we have dY∕dt > 0 for X ∈ (0, ) , and dY∕dt < 0 when X ∈ (− , 0) . If this streamline intersects the line X = 2 at the value Y = Y , then (X(t), Y(t)) lies below the line Y = Y for X(t) ∈ [− , − 2 ) ∪ ( 2 , ] and lies above the line for X(t) ∈ (− 2 , 2 ) . Thus, Introducing the differential equation with (0) = , it follows immediately from (38) and the fact that X(0) = (0) = that

Proposition 1
The horizontal drift experienced by a fluid particle in the lower fluid layer over one wave period decreases strictly with depth, and furthermore, the drift goes to zero at great depths.
Proof The horizontal drift experienced by a fluid particle over one wave period is given by

which can be expressed from (31a) as
If Y 1 < Y 0 , with Ỹ =Ỹ(X) denoting the streamline with Ỹ ( ) = Y 1 , then which, in the limit Y 1 → Y 0 , has the same sign as Thus, D(Y 0 ) > D(Y 1 ) , and the particle drift is decreasing with depth. It is clear from (40) that particle drift goes to zero at great depths (as Combining these results with the phase-plane analysis of system (31) undertaken in Sect. 3.1, facilitates a qualitative description of physical particle motion in the lower-fluid layer as prescribed by (27). Assume a given fluid particle is initially at its lowest depth y(0) = y 0 : we label this position A. This corresponds to X(0) = , and since Ẋ < 0 along streamlines it follows that, in the moving frame, the variable X(t) decreases continuously from to − , and we have: ̇x < 0 , ̇y > 0 for X(t) ∈ ( ∕2, ) ; ̇x > 0 , ̇y > 0 for X(t) ∈ (0, ∕2) ; ̇x > 0 , ̇y < 0 for X(t) ∈ (− ∕2, 0) ; ̇x < 0 , ̇y < 0 for X(t) ∈ (− , − ∕2) . The particle returns to its lowest position in the fluid layer (with depth y = y 0 ) at time t = t Y 0 (− ) > 2 ∕ , having experienced a positive horizontal drift: we label this position B. A representation of this qualitative particle motion is illustrated in Fig. 7.
Remark The primary qualitative difference between fluid particle motion in a lower-fluid layer which has infinite depth, rather than finite depth over a flat impermeable bottom, is that particle trajectories are (almost) circular in the present setting, as opposed to being (almost) elliptical in the finite depth case, when fluid motion is small. To see this, observe that making the approximation (x, y) ≈ (x 0 , y 0 ) on the right-hand side of (27) enables system (27) to be integrated, leading to the solution which implies that Hence, assuming fluid motion is small, the fluid particles follow approximately circular trajectories, with centres located at x * 0 = x 0 + ae ky 0 sin kx 0 and y * 0 = y 0 − ae ky 0 cos kx 0 , and radii ae ky 0 that decrease with increasing depth.

Upper-fluid layer
Physical variables in the upper-fluid layer are obtained by reversing the coordinate transformations (30) by defining Note the reflection involved in the vertical coordinate transformation in (42) reverses the vertical orientation of fluid motion when expressed in terms of the physical coordinate y, as opposed to the Y 1 coordinate. If (X(t), Y 1 (t)) is a solution to (33) in the upper-fluid layer with (42)

Fig. 7
Schematic of a typical particle trajectory in the lower-fluid layer, representing its location at time: for which a∕a 1 < 0 ). The analysis of particle motion in the upper-fluid layer follows along the lines of [25], although with more specificity due to the nature of the solutions A = A 1 = 1 and A = A 2 given by (21). For the sake of completeness, we present the main results and pertinent qualitative features, referring the reader to [25] for full details.

Upper-fluid layer with
There are no closed particle paths in the upper-fluid layer whose motion is governed by system (28). That is, system (28) has no solutions (x(t), y(t)), which are periodic.
Proof With reference to Fig. 6, streamlines of the upper-fluid layer are restricted to the phase portrait region for which Y 0 1 ∈ [ 1 , kh 1 + ] . Although there is a change in the qualitative behaviour of streamlines depending on whether Y 0 1 <Ȳ 1 or Y 0 1 >Ȳ 1 , where Ȳ 1 is prescribed by (36), it can be shown (cf. [25]) that in either case the inequality holds: Noting similarities between inequalities (38) and (43), the arguments used in the proof of Theorem 4.2 apply here to the differential equation For the case Y 1 =Ȳ 1 , in which the streamline is the flat 0-isocline, we note that f (Y 1 ) ≡ f (Y 1 ) ≡ f (Ȳ 1 ) and relation (43) becomes an equality, giving ◻ Combining these results with the phase-plane analysis for system (33) presented in Sect. 3.2 facilitates a qualitative description of physical particle motion in the upper-fluid layer for A = A 2 . Note that, along streamlines in the moving frame, Ẋ (t) < 0 for solutions of system (33): X(t) decreases continuously from X = to reach X = − in time t Y 0 1 (− ) . For streamlines with Y 0 1 ∈ [ 1 ,Ȳ 1 ) in this region (which corresponds physically to the top section of the upper-fluid layer) we have: ̇x < 0 , ̇y > 0 for X(t) ∈ ( ∕2, ) ; ̇x > 0 , ̇y > 0 for X(t) ∈ (0, ∕2) ; ̇x > 0 , ̇y < 0 for X(t) ∈ (− ∕2, 0) ; ̇x < 0 , ̇y < 0 for X(t) ∈ (− , − ∕2) . Fluid particles in this region return to their lowest position in the fluid layer, with depth y = y 0 , say, after the time t Y 0 1 (− ) , having experienced a forward horizontal drift (− ) − 2 ∕k > 0 . This particle motion is captured in schematic (a) of Fig. 8.
For streamlines with Y 0 1 ∈ (Ȳ 1 , kh 1 + ] (which corresponds physically to the bottom region of the upper-fluid layer), we get: ̇x < 0 , ̇y < 0 for X(t) ∈ ( ∕2, ) ; ̇x > 0 , ̇y < 0 for X(t) ∈ (0, ∕2) ; ̇x > 0 , ̇y > 0 for X(t) ∈ (− ∕2, 0) ; ̇x < 0 , ̇y > 0 for X(t) ∈ (− , − ∕2) . In the upper-fluid layer, the horizontal drift experienced by a fluid particle over one wave period is given by ∕k > 0 , which can be expressed from (33) as has the same sign as The following result applies for particle drift in the upper-fluid layer (cf. [25]). lim . The horizontal drift experienced by a fluid particle in the upper-fluid layer over one wave period decreases with depth for fluid motion of the type illustrated in schematic (a) of Fig. 8, while it increases with depth for fluid motion of the type illustrated in schematic (c) of Fig. 8. Accordingly, the minimum horizontal drift experienced by fluid particles occurs at Y =Ȳ 1 , as illustrated in schematic (b) of Fig. 8.

Remark
Fluid particle trajectories in the upper-fluid layer are (approximately) elliptical when fluid motion is small since making the approximation (x, y) ≈ (x 0 , y 0 ) on the righthand side of (28) enables this system to be integrated. This leads to the solution for Hence, for small fluid motion in the upper-fluid layer, particle trajectories approximately follow the elliptical orbits prescribed by where x * 0 = x 0 + a 1 F(y 0 ) sin kx 0 and y * 0 = y 0 − a 1 G(y 0 ) cos kx 0 . The depth-dependent parameters F(y 0 ) and G(y 0 ) determine the size of the major and minor axes, respectively. Note that F(y 0 ) does not vanish in the upper-fluid layer for A = A 2 , while G(y 0 ) vanishes for ȳ 0 = h 1 −Ȳ 1 ∕k : this value corresponds to particle trajectories prescribed by the flat line in schematic (b) of Fig. 8. The change in the sign of G(y 0 ) that occurs at ȳ 0 corresponds to the change in the orientation of particles trajectories, which move clockwise in schematic (a) and anticlockwise in schematic (c) of Fig. 8.

Upper-fluid layer with A = A 1 (= 1)
Due to the straightforward form that the functions f (Y 1 ), g(Y 1 ) assume in the setting A = A 1 = 1 , a comprehensive qualitative description of fluid motion can be achieved.

Theorem 4.4
Let A = 1 . There are no closed particle paths in the upper-fluid layer whose motion is governed by system (28).
Proof For A = 1 , we have f (Y 1 ) = e −Y 1 = −g(Y 1 ) , and it can easily be shown that the analogue of inequality (43) holds for streamlines in the region Y 0 1 ∈ [ 1 , kh 1 + ] .

Proposition 3 ([25])
If A = A 1 , the horizontal drift experienced by fluid particles in the upper-fluid layer over one wave period decreases with depth.
Fluid particle motion matches that illustrated in schematic (a) of Fig. 8, with all fluid particles experiencing a forward drift, except now fluid trajectories are approximately circular for A = 1 (rather than elliptical) since F(y 0 ) = G(y 0 ) = 2e k(y 0 −h 1 ) in (45). Hence, for A = 1 equation (46) prescribes approximately circular particle trajectories for small fluid motion in the upper-fluid layer.

Conclusions
In this article, the fluid motion underlying coupled linear internal and surface waves is studied, from both the Eulerian and Lagrangian viewpoints, for a deep-water two-fluidlayer model (with the lower layer being of infinite depth). Although the linearised governing equations are studied, the resulting dynamical systems which prescribe the fluid motion are themselves nonlinear. In contrast to the setting of two fluid layers of finite depth (recently studied in [25]), when the lower layer is of infinite extent the dispersion relations simplify, and the nondimensional wave amplitude ratio parameter A for the upper-fluid layer system (28) is restricted to either of the values A 1 = 1 and A 2 > 1.
A detailed Eulerian description of the wave-field kinematics for coupled linear travelling waves is achieved using phase-plane analysis, a by-product of which is an elegant graphical verification that coupled surface and internal waves are in-phase for the parameter value A = A 1 = 1 , while they are out-of-phase when A = A 2 > 1 . Furthermore, a complete qualitative Lagrangian description of the fluid particle trajectories in both fluid layers is obtained. In the deep-water regime, the wave-field kinematics in the lower-fluid layer are necessarily quite different to that of the finite depth case featured in [25]. Particle trajectories in the lower-fluid layer all exhibit a forward drift which decreases with increasing depth, and they are approximately circular. In the case of a lower-fluid layer of finite depth (as analysed in [25]), the particle motion is approximately elliptical (with similar forward drift properties).
A significant benefit of the deep-water setting under consideration here is that we can particularise results from [25] concerning the dynamics of the upper-fluid layer, giving a complete qualitative picture of fluid particle motion in the upper-fluid layer. All particles exhibit a forward drift: for A = A 1 = 1 the particle motion is approximately circular, moving clockwise with a drift which decreases with increasing depth. For the parameter value A = A 2 > 1 , the particle motions are approximately elliptical, with upper particles moving clockwise and lower particles moving anticlockwise. A reverse in the orientation of particle trajectory motion occurs at the flat streamline located at Y 1 =Ȳ 1 , which is also the location of minimum forward particle drift: particle drift decreases with depth for fluid particles located between the surface and this streamline, while it increases for fluid particles located between this streamline and the internal wave.