## Abstract

The evolution of the free surface of a three-dimensional conducting fluid in the presence of gravity, surface tension and vertical electric field due to parallel electrodes, is considered. Based on the analysis of the Dirichlet–Neumann operators, a series of fully nonlinear models is derived systematically from the Euler equations in the Hamiltonian framework without assumptions on competing length scales can therefore be applied to systems of arbitrary fluid depth and to disturbances with arbitrary wavelength. For special cases, well-known weakly nonlinear models in shallow and deep fluids can be generalized via introducing extra electric terms. It is shown that the electric field has a great impact on the physical system and can change the qualitative nature of the free surface: (i) when the separation distance between two electrodes is small compared with typical wavelength, the Boussinesq, Benney–Luke (BL) and Kadomtsev–Petviashvili (KP) equations with modified coefficients are obtained, and electric forces can turn KP-I to KP-II and vice versa; (ii) as the parallel electrodes are of large separation distance but the thickness of the liquid is much smaller than typical wavelength, we generalize the BL and KP equations by adding pseudo-differential operators resulting from the electric field; (iii) for a quasi-monochromatic plane wave in deep fluid, we derive the cubic nonlinear Schrödinger (NLS) equation, but its type (focusing or defocusing) is strongly influenced by the value of the electric parameter. For sufficient surface tension, numerical studies reveal that lump-type solutions exist in the aforementioned three regimes. Particularly, even when the associated NLS equation is defocusing for a wave train, lumps can exist in fully nonlinear models.

## 1. Introduction

Electrohydrodynamics (EHD) deals with the interaction between electric field and fluid flow field via the Maxwell stress tensor which plays a crucial role in the coupling. EHD enjoys a wide usage in chemistry, biology and engineering. Recent applications of EHD include, but are not limited to, the bilayer patterning process which is induced by EHD instability and the coupling of kinetics and thermodynamics [1]; EHD conducting pumping, which is used to enhance the heat transfer capability of cooling systems [2]; the industrial coating process that is employed to manufacture a vast number of different products [3]; electrospray ionization, an essentially useful technique in converting solution ions into highly charged gas-phase icons of macromolecules [4]. Owing to the practical significance, in-depth knowledge of the characteristics of EHD, especially the electrohydrodynamic interfacial waves, is therefore of great importance.

The research of electrohydrodynamic interfacial waves was initiated by Taylor & McEwan [5], whose theoretical and experimental results showed that the net force induced by the normal electric field can destabilize the interface between conducting and non-conducting fluids, and Melcher & Schwarz [6], who conducted the linear stability analysis of interfacial waves under a tangential electric field which produces a dispersive regularization for short waves. Recent theoretical studies on electrohydrodynamic interfacial waves for inviscid liquid sheets and layers were able to capture the nonlinear features either based on multi-scale models or via direct numerical simulations for the primitive equations. Examples include the touch-down singularity of liquid sheets under electric fields [7,8], the control and suppression of the Rayleigh–Taylor instability using horizontal electric fields [9], nonlinear electrohydrodynamic Kelvin–Helmholtz and Rayleigh–Taylor instabilities [10,11] and arbitrary amplitude electrocapillary travelling waves in the full Euler equations [12]. All the above-mentioned works are confined to two-dimensional problems, but three-dimensional electrohydrodynamic waves have barely begun to be studied, since the primitive equations are not easily amenable to analysis, and the free boundary nature adds complexity and precludes many of the tools used in other wave problems.

As a first step towards a comprehensive understanding of electrohydrodynamic waves in three dimensions, we provide in the present paper a number of reduced models for a relatively simple physical setting. The physical set-up is sketched in figure 1: two horizontal parallel electrodes separated by a distance with a layer of conducting fluid (such as impure water or mercury) attached to the lower electrode. We investigate the disturbances of the air–fluid interface, and particular attention is paid to the incompressible, inviscid and irrotational flows. The competing forces resulting from gravity, surface tension and normal electric field are all taken into account. The fully nonlinear models derived in this study afford the remarkable simplification over the primitive equations in computations, while the weakly nonlinear models are expected to be more suitable for theoretical explorations.

Two-dimensional free-surface waves propagating on perfect conducting inviscid fluids under normal electric fields have been studied by different authors, and the case that the thickness of the conducting layer is much smaller than the typical wavelength is most commonly examined in the literature. When the thickness of the gas layer is of the same order as that of the conducting layer, Easwaran [13] and Perel’man *et al.* [14] derived the Korteweg–de Vries (KdV) equation and modified KdV equation, respectively, with coefficients depending on the electric parameter. As the top plane electrode is placed far from the free surface, the KdV–Benjamin–Ono equation was derived by Gleeson *et al.* [15] which combines the KdV equation with an additional linear term involving the Hilbert transform. These works were further extended to the full range of scalings for the thickness of the gas layer by Hammerton and his collaborator [16,17].

There has been considerable research devoted to the modelling of non-electrical gravity-capillary waves propagating at the liquid–gas or liquid–liquid interface in three dimensions. In the long-wave approximation and in the weakly nonlinear regime, of note is the work of Kadomtsev & Petviashvili [18], who proposed the so-called KP equation to study the transverse instability of solitons of the KdV equation, Berger & Milewski [19], who studied the generation of fully localized three-dimensional waves by a pressure forcing in the modified Benney–Luke (BL) equation, and Kim & Akylas [20], who investigated the existence, generation and properties of interfacial lumps between two immiscible fluids based on the two-dimensional Benjamin equation with weak transverse variations. In the short-wave limit, Ablowitz & Segur [21] derived the Benney–Roskes–Davey–Stewartson system governing the envelope and induced mean flow, and the system reduces to the cubic nonlinear Schrödinger (NLS) equation in deep water as the mean flow vanishes in the limiting situation. In this paper, we generalize the existing models to include an applied normal electric field, discuss their various limits, explore novel physical differences arising from the electric field and investigate the existence of fully localized electrohydrodynamic travelling waves for sufficiently strong surface tension.

The purpose of this paper is to develop fully nonlinear models and weakly nonlinear theory for gravity-capillary waves propagating on a three-dimensional conducting fluid under a normal electric field. Since the Dirichlet–Neumann operator (DNO) appearing in the Hamiltonian is crucial to free-surface water-wave problems in modelling and computations [22–24], we prove in §2 that the system with electric field is also a Hamiltonian but involves two distinct DNOs, and fully nonlinear models can be obtained via expanding and truncating the DNOs either in the Hamiltonian formulation or in the free-surface boundary conditions. In the weakly nonlinear regime, we derive model equations to describe the evolution of long waves in §3 and the dynamics of the envelope of short waves in §4, respectively. The set of equations derived here is based on systematic asymptotic expansions of the DNOs for different wavelength and amplitude scalings. When surface tension is large enough, locally confined travelling waves propagating on three-dimensional fluids, which are known as lumps, are numerically computed with different models using Fourier pseudo-spectral methods. The last section provides our conclusions.

## 2. Formulation

### (a) Governing equations

Figure 1 shows the definition of the problem. We consider an incompressible, inviscid and perfect conducting fluid, flowing irrotationally on a solid electrode. The Cartesian coordinate system (*x*,*y*,*z*) is introduced, such that *x* is the wave propagating direction, *y* is the transverse direction and *z*-axis rises from the unperturbed fluid level *z*=0 and is positive upwards, the bottom being at *z*=−*h*^{−}. A vertical electric field is imposed and bounded above by a second horizontal electrode at *z*=*h*^{+} demanding *V* =*V* _{0}. We remark that the upper layer can be extended to an infinite height by replacing the Dirichlet boundary condition with *E*_{0} is a constant. We assume that the fluid is a perfect conductor so that the electric strength vanishes within the fluid. If we denote the displacement of the free surface of the fluid by *z*=*η*(*x*,*y*,*t*), then the voltage potential above the fluid is governed by
_{xx}+∂_{yy} is the Laplace operator acting on the horizontal coordinates. The motion of the liquid sheet is also governed by Laplace's equation. If the velocity potential is denoted by *ϕ*(*x*,*y*,*z*,*t*), then
_{x},∂_{y})^{⊤} is the gradient operator acting on horizontal variables. In this paper, the restoring forces due to gravity, surface tension and electric field are taken into account; therefore, equality of the normal stress on the free surface gives the dynamic equation (e.g. [25])
*ρ* is the density of the fluid, *g* accounts for the gravitational acceleration, and *σ* represents the surface tension coefficient of the liquid. Let **n**⋅*Σ*⋅**n** arises from the interfacial electric stress given by the Maxwell stress tensor (*Σ*_{ij})_{i,j=1,2,3}:
*ϵ*_{p} is called the dielectric constant. Noticing *V* _{x}+*η*_{x}*V* _{z}=*V* _{y}+*η*_{y}*V* _{z}=0 at *z*=*η*(*x*,*y*,*t*), a straightforward calculation shows that

We introduce a modified voltage potential *W* by defining *V* =*V* _{0}(*W*+*z*)/*h*^{+}. Therefore, *W* satisfies Laplace's equation, together with the boundary condition *W*=0 at *z*=*h*^{+}. After redefining the velocity potential to absorb constant, the dynamic boundary condition reads
*W* satisfies
*V* _{0}/*h*^{+} should be replaced with *E*_{0}.

### (b) Dirichlet–Neumann operators

The DNO which maps the Dirichlet boundary condition to normal derivatives on the boundary via solving Laplace's equation is essential to address the present problem. Recalling that *W*=−*η* at the free surface, the DNO for the electric potential, denoted by *G*^{+}, can be defined by
*W* and *η* satisfy the periodic boundary conditions in horizontal variables, and *W* solves the system (2.6). Similarly, if we define the surface velocity potential as *G*^{−}, can be expressed as
*h*^{±} in subsequent analyses for simplicity of notations. Following Craig & Sulem [23], by use of the DNOs *G*^{±}(*η*), the kinematic and dynamic boundary conditions can be recast in terms of *ξ* and *η* as follows:
*η*(*x*,*y*,*t*) and the surface velocity potential *ξ*(*x*,*y*,*t*). In this formulation of EHD surface-wave problems, the DNOs play a central role. It was proved by Craig *et al.* [26] that *G*^{−}(*η*) is an analytical operator if the *C*^{1}-norm of the free-surface displacement *η* is considerably small for three-dimensional water-wave problems; therefore, they can be expanded as convergent Taylor series
*G*^{−}_{j}(*η*) is homogeneous of order *j* in *η*, and the first three terms of the series are given as follows:
*D*=−*i*∇, |*D*|^{2}=−Δ and *h*=*h*^{−} for the current problem. As *G*_{0} reduces to |*D*|.

In order to determine the expression of *G*^{+}(*η*), we assume that *G*^{+}(*η*) has a Taylor series representation *G*^{+}_{j}(*η*) (*j*=0,1,2,…), and apply the *operator expansion method*. To proceed, we denote **x**=(*x*,*y*) and **k**=(*k*,*l*), which are the wavenumber with its components along the propagation direction (e.g. *x*) and in the transverse direction (e.g. *y*), respectively. Upon noting that the function *h*=*h*^{+} for the present problem), one obtains the following formula by using the definition of *G*^{+}(*η*):
*η*=0 yields
*η*, we can write down the explicit expressions of *G*^{+}_{n}(*η*) (*n*=0,1,2,…) which are recursive. However, upon noting that the DNO is self-adjoint, Nicholls & Reitich [24] suggested to use the adjoint formulae in computations since the numerics is substantially faster and memory efficient than the original formulae given by Craig & Sulem [23]. Following Nicholls & Reitich's work, we obtain, for *n*>0 even,
*n* odd,
*G*^{+}_{0} simplifies to |*D*|. It is obvious that if *G*^{+}_{0}=*G*^{−}_{0}, *G*^{+}_{1}(*η*)=*G*^{−}_{1}(−*η*) and *G*^{+}_{2}(*η*)=*G*^{−}_{2}(*η*). In fact, it is not difficult to obtain that under this assumption *G*^{+}_{j}(*η*)=*G*^{−}_{j}(−*η*) for arbitrary *j*.

Substituting the expansions of *G*^{±} into the kinematic and dynamic boundary conditions (2.9)–(2.10) and truncating at certain order gives two evolution equations for two unknowns *η* and *ξ*. This approximation has reduced the three-dimensional problem to a two-dimensional one involving only the variables on the surface, which is computationally reasonably simple. In a doubly periodic setting, each term can be efficiently computed using a pseudo-spectral method and the fast Fourier transform (FFT). We consider this approximation as a computational model since it provides an efficient way to compute EHD surface waves propagating on three-dimensional fluids. When the electric field is absent, similar computational methods have been used in a variety of water wave problems (e.g. [23,24]).

### (c) Hamiltonian formulation and fully nonlinear models

In the absence of the electric field (i.e. *E*=0), many attempts were made to prove that water wave problem (2.9)–(2.10) is a Hamiltonian system, but it was Zakharov who finally published the Hamiltonian structure in [27]. Here we show that as *E*≠0, (2.9)–(2.10) is also a Hamiltonian system, and the Hamiltonian functional is the total energy defined by
*W*=−*η* at *z*=*η*, which implies the following relation at the free surface
*η* and *ξ* as

Based on the Hamiltonian formulation, we can construct fully nonlinear models by expanding the DNOs in the Hamiltonian and truncating at certain order. For example, we expand the energy in power of *η* and *ξ* in the Hamiltonian as

Linearizing the system (2.26)–(2.27) and seeking solutions in the form *e*^{i(kx+ly−ωt)} leads to the dispersion relation
*ω* is the wave frequency and *c*_{p} is called the phase speed. The flow is linearly stable if *c*^{2}_{p} is non-negative, so it is obvious that short waves are always stable, whereas long waves are stable only when

## 3. Long-wave approximation

### (a) Shallow (fluid layer)–deep (gas layer) limit

#### (i) Boussinesq scaling

In this section, we consider the case that the typical wavelength is large in comparison with the mean depth of the fluid *h*^{−}, but small in comparison with the thickness of the gas layer *h*^{+}. First of all, we choose *h*^{−}, *ε*. If we further assume that waves are isotropic in *x*- and *y*-directions, then |*D*|=(−Δ)^{1/2}=*O*(*ε*) and
*ε*:
*O*(*ε*^{8}) and apply the formula (2.24). After some algebra, one obtains

The BL-type equation, which is a single equation describing bidirectional waves, can be derived from (3.6)–(3.7). First of all, equation (3.7) implies that *η*=−*ξ*_{t}+*O*(*ε*^{4}). Let us differentiate equation (3.7) with respect to time,
*O*(*ε*^{5}) yields
*η* with −*ξ*_{t}, we arrive at
*θ* and *ξ*, where *θ* is defined by

### Proof.

Calculating the variational derivatives with respect to the canonical variables gives

#### (ii) KP-type models

For the inhomogeneous case, where the variation in the *y*-direction is much slower than that in the *x*-direction, we assume ∂_{y}=*O*(*ε*^{2}) while all other variables are of the same scales as listed in (3.2). It follows that
*O*(*ε*^{6}):
*Q* is a pseudo-differential operator defined as *η*:
*O*(*ε*^{7}) and returning back to the original variables, one can then neglect the small parameter *ε*,
*x* and letting *p*=*ξ*_{x}, one obtains a KP-like equation

It is clear that the *p*_{xxx} term in equation (3.21) vanishes at *EQp*_{x} can be used to balance the nonlinear effects. However, if *E* is much smaller (say, of order *ε*^{3}), we need to proceed to the next step to derive a higher order equation. To obtain a balance between dispersion and nonlinearity near *E*=0, we choose the following scales

### Remark 3.1

For two-dimensional fluids (namely, the governing equations are *y*-independent), the model equations (3.21) and (3.23) were derived by Hammerton & Bassom [17].

### Remark 3.2

For the limiting case *i* sgn(*k*), hence *et al.* [25] using standard asymptotic techniques. Equation (3.24) is also called the two-dimensional Benjamin equation.

### Remark 3.3

If we confine ourselves to one-dimensional surface waves, equation (3.24) is the celebrated Benjamin equation, a reduced model written for long interfacial waves between two immiscible fluids when the upper layer is bounded above and the lower layer is of infinite depth, with additional conditions that the densities of two fluids are nearly equal and the interfacial tension is large (see Benjamin [29]). And it was first derived by Gleeson *et al.* [15] in the context of EHD surface waves.

#### (iii) Lumps

In the absence of electric fields, it is well known that fully localized travelling-wave solutions exist in free-surface/iterfacial water waves when the surface tension is strong. This fact has been confirmed in the long-wave limit by the exact lump solutions of the KP-I equation, and the numerical solutions of the BL equation [19] and two-dimensional Benjamin equation [20]. In this section, based on the model equation (3.10), we investigate the existence and generation mechanism of lumps when the vertical electric field is also taken into account. The numerical scheme for constructing fully localized travelling-wave solutions is an extension of Petviashvili's method [30]. The basic idea is to perform the iteration in the Fourier space supplemented by a normalization factor upon the degree of nonlinearity. Assuming a lump propagates with speed *c* in the *x*-direction, and applying the Fourier transform to equation (3.10) leads to
**k**=(*k*,*l*) is the vector of wavenumber, and *R* is taken to be infinity throughout the computation for simplicity. In order to prevent the unlimited growth or reduction, a multiplier needs to be introduced in every iteration step. Following [31], we propose the numerical scheme as
*α*_{n} is given by

In figure 2, we present the bifurcation curves and typical profiles of solitary waves in equation (3.26) for *E*=1/2 and *a*) shows the relation between *c* and *cξ*_{x} which indicates that both plane solitary waves (dashed curve) and lumps (solid curve) bifurcate from free stream at *c*≈0.866. For positive *c*, only depression solitary waves which feature a negative displacement at their centre were found. Typical examples of two- and three-dimensional solitary-wave solutions are shown in (*b*) and (*d*), respectively, while *x*- and *y*-cross-sections of the lump (*d*) are plotted in (*c*). It is beyond the scope of this paper to study the full bifurcation of solitary waves in equation (3.26).

Kim & Akylas [20] showed that plane solitary waves in the two-dimensional Benjamin equation (3.24) are unstable subject to transverse perturbations of sufficiently long wavelength, and the long-time dynamics of this instability results in the formation of lumps. The new model (3.10), which can be viewed as the bidirectional and isotropic counterpart of (3.24), is expected to have the similar property. We first sketch the numerical scheme for time-dependent simulations. Taking the Fourier transform of the system (3.13)–(3.14) in *x*- and *y*-yields
*Ω* is defined by
*θ* and *ξ* are real. Thus, the problem is reduced to solving
*θ* and *ξ* can be recovered by
*Ξ*(*x*) which is shown in figure 2*b* and perturb it in the transverse direction using a long cosine function, namely
*c*≈0.593 (figure 3*a*). Figure 3*c*,*d* demonstrates the comparisons of *x*- and *y*- cross-sections between the resultant lump (solid lines) and the exact solution (dashed line), which show a remarkable agreement. This also provides us with the validation of our numerical codes. Finally, we remark that in some sense the isotropic equation (3.10) is better than the two-dimensional Benjamin equation (3.24), since the choice of the initial condition for the latter equation requires a careful check. If we assume the wave is locally confined in the *x*-direction, the quantity *y*, and the numerical scheme adopted should persevere this property.

### (b) Shallow–shallow limit

Throughout this subsection, we make the assumption that both *h*^{−} and *h*^{+} are small compared with the typical wavelength λ. The non-dimensionalization is done in the same way as in the shallow-deep case, where *h*^{−} is taken as the length scale. The isotropic Boussinesq regime includes the additional expectation of small-amplitude motions. One chooses the following scaled variables
*ε* as *E* and *R*. Finally, we remark that the KdV equation, which is a celebrated reduced model in classical water-wave problems, was derived by Easwaran [13] and Hammerton [16] in the context of EHD surface waves using conventional perturbation method. Here the KdV equation can be derived by dropping the *y*-dependence in equation (3.45).

## 4. Short-wave approximation

The cubic NLS equation provides a canonical description for the envelope dynamics of a quasi-monochromatic plane short wave propagating in a weakly nonlinear dispersive medium. Let us now consider the derivation of the NLS equation from the primitive system when the electric field and conducting fluid are both of infinite depth. Here we are following a packet of nearly one-dimensional waves, travelling in the *x*-direction, with an identifiable wavenumber *G*^{−}_{j}(*η*) takes the form
*G*^{−}_{0}=|*D*|. We truncate the DNO expansions in the Hamiltonian (2.22), and take the variational derivatives to obtain the Euler–Lagrange equations. We omit the tedious calculations, but list below are the quadratic, cubic, quartic and quintic models (

### (a) Normal form analysis

The NLS equation is a conventional tool to investigate wavepacket solitary waves. In order to derive the cubic NLS equation, we expand the kinematic and dynamic boundary conditions about *z*=0 and retain quantities valid up to the third order terms in *ξ* and *η* which are assumed to be proportional to a small parameter *ε*:
*x*- and *y*-directions, but whose fast oscillation is only in the propagation direction *x*. To derive the governing equation for the wave envelope, we introduce variables of multiple scales *X*=*εx*, *Y* =*εy*, *T*=*εt* and *τ*=*ε*^{2}*t*, and choose e^{i(kx−ωt)} as the carrier wave. One then seeks a solution of the form
*Θ*=*kx*−*ωt*. It is noted that *A*_{j} and *B*_{j} include all the harmonic modes up to *j*, namely *A*_{j}=(⋅)*e*^{0}+(⋅) *e*^{iΘ}+⋯+(⋅) *e*^{jiΘ}. The wave envelope *A* can then be found to satisfy the NLS equation. We omit the details of the derivation, and just state the results (the interested reader is referred to Wang & Milewski [32])
*c*_{g} is equal to the phase velocity *c*_{p}. Then wavepacket solitary waves may bifurcate from infinitesimal periodic waves at this minimum as long as the associated NLS equation is focusing. When *k*=1, *c*_{p} attains its minimum _{1} and λ_{2} are positive when *E*<2, while *γ* is also positive as

### (b) Lumps

However, it is more interesting to search for lump solutions in the primitive equations or in the fully nonlinear models when the associated NLS equation is defocusing. That is because the weakly nonlinear theory only works for small-amplitude waves, but does not rule out lump solutions of finite amplitudes. Hence, we focus our attention on the case *E*=0.5 in the subsequent computations. The fully nonlinear model (4.2)–(4.3) with *m*=5 is numerically solved to seek for travelling lumps with symmetry in both *x*- and *y*-directions, based on an accurate pseudo-spectral method. We write the travelling solution as
*a*_{m,n}=*a*_{−m,n}=*a*_{m,−n} and *b*_{m,n}=*b*_{−m,n}=−*b*_{m,−n}, where *ζ*=*x*−*ct* with the translating speed *c*. *η* and *ξ* are periodic functions on the rectangle of size 2*π*/*k* and 2*π*/*l*. *a*_{m,n} and *b*_{m,n} are unknown Fourier coefficients which need to be solved using Newton's method.

It is presented in figure 4 the bifurcation diagram and typical profiles of electrohydrodynamic lumps which feature oscillatory decaying tails in the propagation direction but monotonic decaying tails in the transverse direction akin to non-electric gravity-capillary lumps. As implied by the defocusing NLS, no small-amplitude lumps can be found below the minimum of the phase speed *E* as the bifurcation parameter. Typical computations shown in figure 4 use *M*=64 an *N*=128.

## 5. Conclusion

In this paper, we have presented the modelling of nonlinear electrohydrodynamic surface waves propagating on three-dimensional perfect conducting fluids under vertical electric fields. The problem was proved to be a Hamiltonian system involving two distinct DNOs. Fully nonlinear models were derived via formally expanding the DNOs as series of pseudo-differential operators, either in the free boundary conditions or in the Hamiltonian functional.

When the Bond number is of order one, both shallow- and deep-fluid configurations were considered, and different weakly nonlinear models were proposed for long- and short-wave approximations based on the systematic expansions of the DNOs in the Hamiltonian. It turns out that the electric field has a great effect on the physical system: (i) when the separation distance between two electrodes is small compared with typical wavelength, the Boussinesq, BL and KP equations with modified coefficients are obtained, and electric forces can turn KP-I to KP-II and vice versa; (ii) as the parallel electrodes are of large separation distance but the thickness of the liquid is much smaller than typical wavelength, we generalize the BL, KP and fifth-order KP equations by introducing linear pseudo-differential operators arising from the electric field; (iii) for a quasi-monochromatic plane wave in deep fluid, we derive the cubic NLS equation whose type (focusing or defocusing) is strongly affected by the value of the electric parameter.

As surface tension is strong enough, locally confined travelling waves in three dimensions were numerically computed in equation (3.10), and the transverse instability of plane solitary waves was further examined, and found to be a mechanism of generation of lumps. Lump solutions also exist in the fully nonlinear model in deep fluids, even though the associated NLS equation is of defocusing type. The stability and dynamics of electrohydrodynamic lumps, which are beyond the scope of this paper, are of great interest, and we leave for future studies. Finally, it must be pointed out that when gravity and surface tension are equally important, the viscous effect, which we have neglected in the modelling, plays an important role due to small typical wavelength. Therefore, in practice, we can derive models from the Navier–Stokes equations when the viscous effect is strong (e.g. [33,17]), or use the visco-potential theory for weakly viscous fluids [34].

## Competing interests

I declare I have no competing interests.

## Funding

This work was supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (grant no. XDB22040203) and the Key Research Program of Frontier Sciences, CAS (grant no. QYZDB-SSW-SYS015).

## Acknowledgements

The author is grateful to all referees for their valuable suggestions.

- Received November 2, 2016.
- Accepted March 1, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.