## Abstract

Recently, the authors presented two numerical studies for capturing the flow structure beneath water waves (Nachbin and Ribeiro-Junior 2014 *Disc. Cont. Dyn. Syst. A* **34**, 3135–3153 (doi:10.3934/dcds.2014.34.3135); Ribeiro-Junior *et al.* 2017 *J. Fluid Mech.* **812**, 792–814 (doi:10.1017/jfm.2016.820)). Closed orbits for irrotational waves with an opposing current and stagnation points for rotational waves were some of the issues addressed. This paper summarizes the numerical strategies adopted for capturing the flow beneath irrotational and rotational water waves. It also presents new preliminary results for particle trajectories, due to irrotational waves, in the presence of a bottom topography.

This article is part of the theme issue ‘Nonlinear water waves’.

## 1. Introduction

The great research interest in surface water waves goes back to past centuries. This problem has attracted attention due to its obvious practical applications in the ocean and open channels, as well as for its many challenging mathematical problems in analysis, asymptotic theory, numerical analysis and scientific computing. Recently, great attention has been paid to the flow in the bulk of the fluid, generated by the propagation of (nonlinear) large-amplitude surface waves. For two-dimensional waves, when the flow is taken to be incompressible, both the theory and the numerical modelling can take advantage of the structure due to harmonic functions, or even of complex-valued analytic functions. Thus, this is all reminiscent of potential theory.

This article reviews the numerical methods presented by Nachbin, Ribeiro-Junior and co-workers [1,2]. In the numerical methods literature, for computing the evolution of the free surface disturbance, variations of potential theory modelling have been around for a while. To name a few, we have the work of Vinje & Brezig [3] where a (complex) Cauchy integral is used for the complex potential, and wave breaking is calculated numerically. Baker *et al.* [4] at the same time were using vortex sheets to compute wave breaking. Again a complex Cauchy-type integral is used with a kernel in the form of a cotangent. Linearization of the free surface leads to a Hilbert transform on the circle. There are many forms to attain a boundary integral formulation of the water wave problem: through Green's third identity, through a Dirichlet-to-Neumann operator, through a non-local formulation that somehow encompasses these ideas, among others. Each approach has its advantages and drawbacks. Based on their computational goals, interested researchers should be aware of the respective features. Recently, Wilkening & Vasan [5] published an article where they compared five different approaches for computing the Dirichlet-to-Neumann operator, which include methods such as a boundary integral formulation (through Green's identity) or the Ablowitz *et al.* [6] non-local formulation.

The flow beneath water waves can be captured by different numerical strategies, as outlined above. The following articles [7–16] are of interest to the problems discussed here. As will be presented herein, we adopt a boundary integral formulation (through Green's third identity) as well as conformal mappings, in order to compute our features of interest with great precision. In one problem presented, the wave elevation is conformally mapped onto a flat line, while in another problem presented, the bottom topography is mapped onto a flat bottom. In both cases, we have a flat strip as our canonical domain, where the harmonic extension of boundary data is more easily done.

Our numerical work has been greatly inspired by results in analysis. Of particular interest we sought closed particle orbits beneath Stokes waves with an opposing current, stagnation points beneath rotational Stokes waves, as well as the related pressure anomalies reported by Teles da Silva & Peregrine [17]. These issues have been numerically investigated by other authors in [11,15] and in [2] we performed a more detailed numerical study.

From the mathematical viewpoint there are still many issues to be investigated. Only recently rigorous theory became available, through several efforts in analysis. Constantin [18] proved in 2006 that particle orbits displayed a drift beneath a propagating irrotational nonlinear surface wave. In 2008, Constantin & Villari [19] used a simpler technique to show that this was also true for linear waves. In 2010, Constantin & Strauss [20] proved that there exists closed particle trajectories when a wave is opposed by a current. Its location and other quantitative information were not provided and this prompted our work in [1]. Concerning the investigation of the flow structure beneath rotational water waves, the literature and results have been growing in recent years. We refer to some recent publications from which the interested reader can explore the literature [21–30]. Laboratory experiments are reported in [31–33].

The paper is organized as follows. In §2, we present the boundary integral formulation following Green's third identity. Using a periodic Green's function [34], we are able to eliminate three boundary segments of the contour integral, thus reducing it to only the free surface. In §3, we present the formulation for rotational waves. We first find our flow domain through a Newton-type continuation method. Namely we deduce equations for finding the wave elevation, according to our choices of parameters. Based on the flow domain obtained, we establish the expressions for the conformal mapping from the canonical domain onto our physical domain. In the simple canonical domain, we are led to Fourier-type expressions for the corresponding harmonic functions and their harmonic conjugates. A preliminary study for particle trajectories beneath a linear wave in the presence of a highly irregular topography is presented in §4. The conformal mapping of the bottom topography, through the Schwarz–Christoffel mapping, plays an important role. Again, through Green's third identity, we reduce the contour integral to only the free surface. Section 5 contains the conclusion.

## 2. Boundary integral formulation of irrotational waves

For the study of flow structure beneath water waves, the Euler equations in two dimensions are the typical model considered. The flow is taken to be incompressible and inviscid. As a starting point for the study of particle trajectories beneath surface waves, the flow is also considered to be irrotational. In a following section, constant vorticity will be considered in the model. The associated problem in partial differential equations (PDEs) has a free boundary at the top of the fluid domain, namely defined by the wave profile. On the other side of the domain, a fixed impermeable bottom leads to a Neumann condition. First, we take the bottom to be flat. Then we will show how a varying topography can be taken into account for the case of irrotational waves.

Let (*u*,*v*) be the velocity vector and *η*(*x*,*t*) be the free surface elevation, namely the wave profile. The Euler system of equations is given by
2.1
and
2.2
where *P*(*x*,*y*,*t*) is the pressure, *ρ* is the constant water density and *g* is the acceleration due to gravity. The flat bottom is taken at depth *d*. We first take the flow to be irrotational, so that
2.3
System (2.1)–(2.3) is completed with the following boundary conditions:
2.4
The atmospheric pressure *P*_{atm} is constant along the free surface [35].

In what follows, we will introduce the potential theory formulation of Euler's equations in order to formulate a boundary integral equation useful for studying particle trajectories beneath a Stokes wave travelling along the free surface. We consider Stokes waves which are smooth solutions to Euler equations. These are nonlinear, periodic travelling wave solutions of spatial period λ and speed *c*. The speed is taken to be positive and, without loss of generality, we study waves propagating to the right. All dependent variables in the above system of equations are λ-periodic and will be written in terms of the travelling variable (*x*−*ct*). The wave profile *η* is taken to have only one crest and one trough per period. The wave profile decreases from crest to trough. Moreover, *η*, *u* and *P* are symmetric with respect to the crest, while *v* is antisymmetric. In the absence of stagnation points in the flow, the symmetry of the free surface (visible in all figures presented in the paper) is guaranteed by the fact that the wave profile is monotonic between crests and troughs. For irrotational waves the reader will find details in the book by Okamoto & Shoji [36], while for rotational waves, details are in the article by Constantin *et al.* [37].

We take the undisturbed water level to be at *y*=0 in such a way that the mean-zero (wave elevation) condition
is satisfied.

We now recast the Euler equations in terms of potential theory. Let *ϕ*(*x*,*y*,*t*) be the velocity potential where (*u*,*v*)=∇*ϕ*. At this stage, the underlying hypothesis is that the flow regime is incompressible and irrotational. As a consequence, *ϕ* is a harmonic function in the domain
and satisfies the following set of equations (see [35]):
2.5
with the free surface (Bernoulli and kinematic) conditions
2.6
together with the impermeability condition
2.7

The trajectory of a particle, initially located at (*x*_{0},*y*_{0}) in the bulk of the fluid, is given by the curve (*x*(*t*),*y*(*t*)) solution of the dynamical system
2.8
In the horizontal component of the vector field, we have incorporated the effect of an underlying current of strength *κ*. For *κ*<0, we have a current in the opposite wave-propagation direction, while for *κ*>0 it is in a favourable direction. We should point out that the average of the horizontal fluid velocity at any depth (beneath the trough) equals *κ*, a fact that explains the interpretation of *κ* as a measure of the strength of the underlying uniform current [21].

### (a) Green's third identity

Now, we take advantage of the potential theory framework in order to evaluate, in a very efficient and accurate fashion, the vector field for the trajectory equations (2.8); namely only using boundary values of the potential. This follows from Green's third identity.

Take a periodic cell in our fluid domain, which will be denoted by . Let be a point beneath the wave in the corresponding periodic cell . A (convenient) Green's function for Laplace's equation will be denoted as *K*_{(x,y)}. Its logarithmic singularity is located at (*x*,*y*). From Green's third identity [38], we have that
2.9
where , with segments of the boundaries labelled with b at the bottom, l at the left and r at the right of the cell. The free surface is denoted by *Γ* and **n** is the outward normal vector along the boundary. The periodic Green's function [34] of our choice is
2.10
Through the method of images we choose as our kernel
This Green's function not only is periodic but also satisfies the Neumann condition at the flat bottom:
2.11
For the contour integral in (2.9), integration along the periodic cell's sides vanishes by cancellation. The Neumann conditions for both *ϕ* and *K* imply that integration along the bottom is zero. Hence the boundary integral, for obtaining the velocity potential beneath the surface wave, has been restricted only to the free surface *Γ*(*t*), which is stationary when considered in a frame moving with the wave speed. Green's third identity reduces to
2.12
We have, for the normal vector, that

### (b) Numerical method for particle trajectories

We now summarize our particle trajectories algorithm [1] based on the boundary integral formulation. We recall that the travelling wave profile does not depend on the value of the current.

#### (i) The algorithm

We have as data the channel depth *d*, the current's intensity *κ*, the initial wave profile *η*_{0}≡*η*(*x*,0) together with its travelling speed *c*, which is defined in the absence of a current. In [1], an approximation for *η*_{0} and *c* is provided according to Fenton's work [39]. The following steps summarize the algorithm:

(A) From the information listed above, a complete set of free surface data is obtained through the stationary Bernoulli's law and kinematic condition. Details are given in [1].

(B) No evolution scheme is needed. Once all initial values

*ϕ*_{x}(*x*,*η*_{0},0),*ϕ*_{y}(*x*,*η*_{0},0) and*ϕ*(*x*,*η*_{0},0) are obtained, their stationary propagation corresponds to translations, namely applying the corresponding rotation in Fourier space [40]. This procedure is spectrally accurate [40]. In the presence of a current, the translation is done for*c*+*κ*.(C) Differentiating the boundary integral representation of the velocity potential provides the vector field for system (2.8). For example, the horizontal speed is given by 2.13 Similarly for

*v*(*x*,*y*,*t*)=*ϕ*_{y}(*x*,*y*,*t*). The kernel's singularity is inside ;*K*is a smooth periodic function and therefore the trapezoidal (integration) rule is spectrally accurate.(D) The fourth-order Runge–Kutta method (RK4) is applied to the dynamical system (2.8). Steps (B) and (C) are then repeated accordingly.

At the top of figure 1, we show three types of particle orbits. Near the free surface (dashed line) we have a trajectory with a forward drift. At the bottom (in a black solid line), we have a particle trajectory with a drift to the left. In the middle, we have several closed orbits as predicted by Constantin & Strauss [20]. The starting point for each closed orbit is different. Therefore, depending on the relative strength between the wave's forward drift and the current's backward motion, we can observe a closed orbit for the particle trajectory. In figure 1, at the bottom, we reproduce all three trajectories, now in the wave's moving frame. We conjecture that all closed orbits lie on the same streamline. This streamline is at a depth where, as mentioned above, there is a balance between the wave action and the current on the particle's trajectory. This balance implies no drift. To the best of our knowledge, there is no theory predicting the precise streamline where this balance takes place. The location of this streamline should clearly exhibit a non-trivial dependence on wave height, on current strength and on other relevant parameters of this problem. For the reader interested in a review which addresses progress in laboratory experiments and field observations of Stokes drift, we suggest [41].

## 3. Rotational waves

In this section, we summarize the work presented in [2] where stagnation points were studied numerically. Let **U**=(*u*(*x*,*y*,*t*),*v*(*x*,*y*,*t*)) denote the velocity field and recall that the Euler equations are
3.1
where **j**=(0,1). The vorticity is defined by
3.2
Take the background *linear shear flow* **U**_{0} (of constant vorticity *ω*(*x*,*y*,*t*)=*γ*) such that
3.3
The velocity field can be written as
3.4
where *ϕ* is the velocity potential produced by an irrotational (wave-induced) perturbation to the background shear flow. The incompressibility condition to be satisfied in the fluid is ∇⋅**U**=*Δϕ*=0 [2].

The equations are non-dimensionalized through *d* and (as the characteristic length and time). The (dimensionless) reference channel has unit depth together with a unit linear long-wave speed.

All dependent variables are taken to be periodic in *x* with wavelength λ and travelling speed *c*. The natural choice for independent variables in the travelling frame is *X*=(*x*−*ct*) and *Y* =*y*. As mentioned before, the *y*-axis is centred so that *η*(*x*,*t*) is mean-zero. The function *η* is considered to be even and *ϕ* to be odd with respect to *X*=0, where the wave crest is centred. In the moving frame, the (stationary) equations for the potential component of the velocity field read as [2]
3.5
and
3.6
where *ψ* is associated with *ϕ* (as its harmonic conjugate) and *B* is the Bernoulli constant.

As presented in [2], a conformal mapping is established from a (rectangular) canonical domain onto the (stationary) fluid domain. Namely we formulate the mapping from a strip in the canonical *w*-plane (*w*=*ξ*+*iζ*) onto the fluid body in the *Z*-plane (*Z*=*X*+i*Y*). In other words, we derive expressions for the harmonic functions and . In particular, represents the wave elevation in terms of the canonical coordinates. It is the Dirichlet condition for the harmonic problem in . Similarly, we obtain Laplace equations (in the canonical coordinates) for and , where and . In the canonical domain ‘linear-like’ Fourier expressions are easily obtained for the respective harmonic functions and their harmonic conjugates. Let **F**_{k}[ *f*] indicate the amplitude of the *k*th Fourier component of a function *f*(*ξ*). The above-mentioned ‘linear-like’ function representations are given as [2]
3.7
3.8
3.9
3.10
where *k*_{j}=2*πj*/λ (for ), , and
Here, *L* is the wavelength in the canonical domain. We have used the Cauchy–Riemann relations and the fact that
3.11
The integral (Hilbert-like) operator acts on the function's respective harmonic conjugate [2]. The boundary value of the velocity potential's harmonic conjugate is
3.12
Through the relation it follows that along the (canonical) free surface and therefore . This relation allows us to compute from the free surface .

Clearly, once the boundary data are obtained, the potential theory problem is solved and one can proceed to solving the trajectory equations in the bulk. Particle trajectories and streamlines are mapped to the physical domain using (3.8) and (3.7). Solving for is equivalent to finding the travelling wave profile compatible with the parameters chosen; namely finding compatible with our choice of vorticity strength, wave height and wavelength as well as channel depth. Having all this in mind, finding is reduced to a Newton-type problem, as presented next.

The free surface conditions from (3.5) can be combined into one expression, in terms of only the free surface elevation :
3.13
As unknowns we have , the wave speed *c*, the Bernoulli constant *B* and the channel depth *D* in the canonical domain. We have chosen to keep the wavelength (hence the spatial period) invariant under the conformal mapping [2]. Therefore, depending on the wave profile, compatible with our choice of parameters, a corresponding channel depth *D* will be found, namely under the restriction that the mapping is conformal. To close our 4×4 system, the above equation is supplemented by three conditions: we impose our target wave height *H* through
3.14
we impose a mean-zero wave profile in the physical domain using
3.15
and finally we have the depth condition [2]
3.16
Equations (3.13)–(3.16) are solved using a Newton continuation method [2].

Using the conformal mapping defined by (3.7), (3.8) and the harmonic potential (3.9), we have completed the information needed for the particle dynamics. The pressure is also readily available through expression (3.6).

In figure 2, we present the streamlines for two different values of the vorticity *γ* and wave steepness *ε*=*H*/λ. In figure 2*a*, we have in grey the recirculation region, displaying three stagnation points, one at the centre of the recirculation zone, as well as two saddle points at the bottom. It was shown in [2] how these stagnation points move as the vorticity value changes. As the vorticity decreases, the three points will collide at the bottom (*x*=0) and leave the fluid domain. Increasing the constant vorticity elongates the recirculation region and the saddle points reach the sides of our periodic cell, becoming one saddle point at *x*=±*π*, for example. Further increasing the vorticity causes the recirculation region to rise, leading to the cat's eye pattern depicted in figure 2*b*. We have a critical layer where, above the recirculation region, the flow is to the left and below that the flow is to the right.

In [2], we also showed how an anomalous pressure arises, manifested as a minimum pressure inside the fluid domain just above the recirculation region. One also observes a maximum pressure along the bottom located below the troughs rather than below the crest. This is due to the saddle points. Both types of anomalies (minimum and maximum at unexpected locations) are due to a ‘hydrofoil-like’ effect of the recirculation region. The pressure is reduced where a fast flow passes above the ‘hydrofoil’. On the other hand, the flow decelerates fast near the saddle point, increasing the pressure. The minimum pressure anomaly had been conjectured by Teles da Silva & Peregrine [17], and was indicated by other authors as pointed out in [2]. Nevertheless, through the conformal mapping strategy we were able to find (in the canonical domain) the precise location of these pressure anomalies, as well as the precise position of the stagnation points.

### (a) Particle trajectories

In the canonical domain, the particle trajectories (*ξ*(*t*),*ζ*(*t*)) are integrated with a Runge–Kutta method through the dynamical system
3.17
Setting *U*=*c* yields the paths in a frame moving with the wave. These are then mapped to the physical domain with .

The stagnation points are identified as points of zero velocity in the wave's moving frame. To find these points, we use the symmetry and the periodicity of *ϕ*, as well as expressions (3.7)–(3.9), to compute *ϕ*_{X} and *ϕ*_{Y}, having in mind (3.4). In the wave's moving frame, as expected, these particle trajectories follow exactly the streamlines computed in the canonical domain and mapped to the physical domain.

### (b) Alternative formulation

One can compute the trajectory equations in the physical domain using only boundary data. The vector field *ϕ*_{X}(*X*,*Y*), *ϕ*_{Y}(*X*,*Y*), for the dynamical system in the bulk of the fluid, is readily available through contour integration, such as
3.18

Here, denotes the parametrization of the free surface elevation, **n** is the outward normal vector and is the periodic Green's function presented before. A similar boundary integral is used for *ϕ*_{Y}(*X*,*Y*). The Neumann data can be computed through a weakly nonlinear Dirichlet-to-Neumann map [42]. Having computed the vector field, we can also find the stagnation points. Nevertheless, the first procedure was more accurate.

## 4. Irrotational waves with topography

In this section, we present preliminary results for irrotational waves interacting with a bottom topography. Namely we are interested in investigating how a topography can affect the flow structure beneath a surface water wave. Linear waves are considered in the form of a wave packet. We consider a localized free surface disturbance together with periodic boundary conditions assigned far away from the region of interest.

The fluid domain may contain a highly variable, non-smooth, bottom topography as displayed in figure 3. Hamiltonian models for this configuration are proposed in [43]. For the two-dimensional problem considered here, conformal mapping provides a perfect change of variables. This is achieved by mapping a canonical (uniform) strip onto the irregular (undisturbed) fluid domain containing, for example, a polygonal-shaped topography, such as in figure 4. The Schwarz–Christoffel mapping [44] ensures the existence of such a conformal map, which is analytic everywhere except at the preimages of the topography's corners. Let *w*=*ξ*+i*ζ* represent a point in the canonical *w*-plane, while *z*=*x*+i*y* is a point in the physical domain. The conformal mapping is represented by the complex function . The Jacobian of the mapping is defined by |*J*|=|*J*|(*ξ*,*ζ*)=|d*z*/d*w*|^{2}. Under this mapping, a Cartesian grid in the canonical domain is mapped onto an orthogonal curvilinear coordinate system in the physical domain by which the Laplacian is invariant. Hence we have a very convenient change of variables.

As previously reported in [42,45–47], under the conformal mapping the dimensionless set of nonlinear potential theory equations can be written in the form
4.1
4.2
4.3
4.4
The parameter *α*=*a*/*h*_{0}, when taken to be small, indicates a weakly nonlinear regime, while the parameter *β*=(*h*_{0}/λ)^{2}, when taken to be small, indicates a weakly dispersive (long-wave) regime. Here, *a*, λ and *h*_{0} are, respectively, the amplitude of the wave, the wavelength and the undisturbed depth in the dimensional physical domain. Note that, under the mapping , *N* is the wave-profile representation in the canonical domain. We point out that now the Laplacian in (4.1) is in terms of *ξ* and *ζ*. The velocity potential is and the wave elevation *η*=*η*(*x*(*ξ*,*ζ*_{FS}),*t*), namely evaluated at points determined along the canonical undisturbed free surface *ζ*=*ζ*_{FS}. For example, the Schwarz–Christoffel toolbox [44] sets (by default) *ζ*_{FS}=1. A straightforward description of the use of the Schwarz–Christoffel toolbox can be found in [46,47].

In the canonical domain, we have a flat bottom boundary condition, as seen by the homogeneous Neumann condition (4.4). Moreover, the physical domain whose length is *L* is mapped onto a strip with length and height . Hence we can use the same periodic Green's function adopted earlier in this article, now given in the canonical coordinates:
4.5
where the kernel
‘erases’ the boundary integral along three sides of the periodic cell. We are again left with integration only along the free surface. The effect of the bottom topography only comes in through the Jacobian |*J*|.

We compute the particle trajectories in the canonical domain through the dynamical system 4.6 4.7 4.8 where (. Note that we have used the chain rule through 4.9 and 4.10 Once these have been integrated with the Runge–Kutta method, the trajectories are mapped to the physical domain using the expressions .

### (a) The evolution scheme

At this first stage of the research, we consider linear waves, and therefore the potential theory formulation becomes 4.11 4.12 4.13 4.14 where .

The wave elevation can be eliminated from the free surface conditions, leading to
It is then recovered from equation (4.13). Our discrete variables are defined by
Take a Taylor series expansion in time, for the velocity potential at *ξ*=*ξ*_{j}:
This yields the first finite difference expression of our scheme:
4.15
where for we have used the Dirichlet-to-Neumann operator (*DtN*):
4.16

By we indicate the use of the fast Fourier transform and its inverse, respectively. Note that through the symbol of this Fourier operator, we have taken into account equations (4.11) and (4.14). Finally, the evolution of the wave elevation is given by 4.17 where (4.16) is used to update , after having computed (4.15). The scheme is therefore explicit.

The particle velocity is computed in the canonical domain. In the present numerical illustration, we have considered only points far away from the free surface. Therefore, we can impose Green's third identity along the undisturbed free surface. In [1], even when considering linear waves, particle trajectories near the free surface had to be computed with Green's third identity imposed on the wave profile in order to avoid the particle from hitting the singularity distribution of the boundary integral formulation. Here, we use the simpler boundary contour (imposed at ) and compute and with 4.18 and 4.19

### (b) Numerical simulations

A domain with a submerged well is depicted in figure 4. Our goal is to show that the particle trajectory can be affected by the presence of the well, namely being protected by its tall side walls. Incident from the left is a wave packet. The problem is linear and the wave height is immaterial. Nevertheless, in figure 4 its height scales with , the parameter used in the nonlinear formulation. The problem is not periodic. We use the localized wave packet to avoid wave activity at the extremes of our computational domain of total length *L*.

We point out that, in the present formulation, the potential theory evolution and respective boundary integral are all restricted to the free surface. Information from the bottom topography is encoded in the Jacobian |*J*|. We consider a regime with very weak dispersion, and we set the initial data and parameters accordingly. In the numerical simulations reported here, we use as parameter values *α*=10^{−3}, *β*=5×10^{−4}, *L*=40, and , according to figure 4.

Our approximation to a right-travelling wave is given as 4.20 The two particles whose trajectories we will follow are located at (left) and at (right).

This profile has been chosen so that the inner part of the packet resembles a periodic wave. The initial velocity potential is obtained from using the dispersion relation . Note that we have a long wave packet in a shallow channel.

In figure 5, we observe the wave packet travelling from the left over a flat bottom and generating two particle trajectories at the same depth. The wave packet displays very little dispersive effects and the particle trajectories are very similar.

In the simulation presented in figure 6, we have placed the second particle (to the right) inside a submerged well with tall side walls rising up to 70% of the total depth. The particle is highly confined inside a very narrow well. From top to bottom in figure 6, we observed the effect of the confinement. Figure 7 has a detail of the bottom inset showing that the confined trajectory is roughly three decades smaller than the corresponding free orbit. For the free particle to the left, we see that its orbit contains outer loops due to the passage of the wave packet, as well as inner loops due to the passage of the (smaller, nearly monochromatic) reflected wave.

## 5. Conclusion

We have reviewed our numerical work from two recent articles [1,2]. Irrotational and rotational waves were considered. The goal was to capture certain flow properties beneath surface waves, in an efficient and accurate fashion. Both articles concern strategies manipulating (mostly) with boundary data for computing wave-induced quantities in the bulk of the fluid.

We have also presented a preliminary study in the presence of a bottom topography. Again we show that we can restrict ourselves to computing particle trajectories with only boundary data. In the future, we hope to extend our numerical simulations to regimes which can provide feedback as well as conjectures for further analysis from researchers in this community.

## Data accessibility

This article has no additional data.

## Authors' contributions

Both authors contributed equally to the work, and gave final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

The work of A.N. was supported, in part, by CNPq under (PQ-1B) 301949/2012-2 and FAPERJ Cientistas do Nosso Estado project no. E-26/201.164/2014. The work of R.R.-J. was supported by Petrobras through the program ANP/PETROBRAS PRH 32 projects nos 6000.0061847.10.4 and 6000. 0069459.11.4.

## Acknowledgements

We are grateful to Prof. Adrian Constantin for several fruitful conversations and also to Prof. Paul Milewski for his important collaboration in [2]. We thank the referee for several comments that improved this article.

## Footnotes

One contribution of 19 to a theme issue ‘Nonlinear water waves’.

- Accepted August 23, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.