## Abstract

Two-dimensional gravity waves travelling under an ice sheet are studied. The flow is assumed to be potential. Weakly nonlinear solutions are derived and fully nonlinear solutions are calculated numerically. Periodic waves and generalized solitary waves are studied.

## 1. Introduction

This paper is concerned with nonlinear waves propagating in a fluid bounded above by a sheet of ice. Such waves can be generated by a load moving on the ice sheet. The research is motivated in part by applications to transport systems in cold regions, where bodies of water are transformed into roads and runways [1] and where air-cushioned vehicles are used to break the ice [2]. Most previous theoretical works on the subject use linear models. Such models are appropriate to describe waves of small amplitude. However, they cannot describe high-amplitude deflections of the ice such as those observed by Marko [3] and Hegarty & Squire [4]. Our aim in this paper is to present nonlinear solutions that describe waves of arbitrary amplitude.

As we shall see the mathematical problem is similar to that of gravity–capillary (GC) waves, but it is more complex because the curvature term in the dynamic boundary condition is replaced by a higher order derivative term.

We restrict our attention to two-dimensional waves. There are then periodic waves, solitary waves with decaying tails and generalized solitary waves. The periodic waves were studied by Forbes [5,6] and the solitary waves with decaying tails by Parau & Dias [7]. Further results on solitary waves in water of infinite depth were recently obtained by Milewski *et al*. [8]. Here, we present numerical computations of generalized solitary waves. In addition, we extend the results of Forbes for periodic waves to larger amplitude and provide numerical evidence for the existence of overhanging waves. Three-dimensional waves will be considered in another paper [9] of this issue.

The problems of waves propagating under an ice sheet and of GC waves are formulated in §§2 and 3. Periodic and generalized solitary waves are discussed in §§4 and 5, respectively. Some concluding remarks are presented in §6.

## 2. Formulation of the ice problem

We consider a train of waves travelling at a constant velocity *c* at the surface of an inviscid and incompressible fluid of constant depth *h*. The flow is assumed to be irrotational. The fluid is covered above by a layer of ice.

We choose a frame of reference moving with the wave, so that the flow is steady. We define Cartesian coordinates with the *y*-axis pointing upwards and *x*=0 at a crest of the waves. Gravity is acting in the negative *y*-direction.

We introduce the potential function *ϕ* and the streamfunction *ψ*. Without loss of generality, we choose *ψ*=0 on the free surface and *ϕ*=0 at the crest of the wave where *x*=0. We denote by −*Q* the constant value of *ψ* on the bottom.

We model the ice by assuming that the pressure *P*_{i} exerted by the ice on the surface of the fluid is given by
2.1
where *D* is the flexural rigidity of the ice and *K* denotes the curvature of the free surface. Relation (2.1) is a simplified version of the model proposed by Forbes [5] when the terms involving the ice thickness are neglected. A typical value of *D* in deep water for an ice thickness of 1.6 m is *D*=1.6×10^{9} Nm.

The mathematical problem is then defined by the equations
2.2
2.3
2.4
2.5
Here, *y*=*η*(*x*) denotes the equation of the (unknown) free surface, *g* the acceleration of gravity and *ρ* the density of the fluid. Equations (2.3) and (2.5) are the kinematic boundary conditions on the free surface and on the bottom. Equation (2.4) is the dynamic boundary condition on the free surface. It is derived by combining equation (2.1) with Bernoulli equation. The expression in the square brackets of equation (2.4) is the curvature *K* of the free surface expressed in terms of *η*(*x*).

We shall restrict our attention to waves which are symmetric about *x*=0. Therefore, we assume that
2.6

Equations (2.2)–(2.6) describe waves propagating under an ice sheet, and therefore we refer to the system defined by equations (2.2)–(2.6) as the ICE system.

## 3. Formulation of the gravity–capillary problem

The problem of GC waves is similar to that of §2, except that there is no ice sheet on top of the free surface. Instead, we take into account the effect of surface tension. The relation (2.1) is then replaced by
3.1
where *T* is the surface tension (assumed to be constant) and *P*_{i} denotes now the jump of pressure across the free surface.

Following the derivation in §2, we find that the mathematical problem is defined by equations (2.2)–(2.6) but with equation (2.4) replaced by 3.2 We shall refer to this system of equations as the GC system.

## 4. Periodic waves for the ICE system

The ICE system for periodic waves is defined by equations (2.2)–(2.6), together with the equations
4.1
4.2
4.3
where *λ* denotes the wavelength of the wave. Equation (4.2) defines the phase velocity *c*, and equation (4.3) fixes the origin of *y* as the mean water level.

Similarly, the GC system for periodic waves is defined by equations (2.2), (2.3), (3.2), (2.5) and (2.6) together with equations (4.1)–(4.3).

### (a) Analytical solutions

We first note that an exact solution of both the ICE system and the GC system is
4.4
The solution (4.4) is simply a uniform stream with velocity *c* and an undisturbed flat free surface.

A natural idea is then to seek a solution as a perturbation of equation (4.4). To achieve this, we introduce a parameter *ϵ* which is a measure of the amplitude of the wave (a precise definition will be given later), so that the wave reduces to the uniform stream (4.4) as . We then write the expansions
4.5
4.6
4.7
4.8

Such expansions were first used by Stokes [10] to study pure gravity waves (i.e. the GC system in the absence of surface tension). Relations (4.5)–(4.8) are often referred to as the Stokes expansion. Stokes’ work was later generalized to the complete GC system by Wilton [11], Schwartz & Vanden-Broeck [12], Chen & Saffman [13] and others. As we shall see *ϵϕ*_{1}(*x*,*y*) and *ϵη*_{1}(*x*) are the classical linear approximations and the terms of order *ϵ*^{2},*ϵ*^{3},… are nonlinear corrections.

Substituting equations (4.5)–(4.8) into the ICE system and the GC system and equating powers of *ϵ* lead to a succession of linear systems for *ϕ*_{1}, *η*_{1}, *c*_{0}, *B*_{0}, *ϕ*_{2}, *η*_{2}, etc. These linear systems can then be solved by separation of variables. Since the calculations are classical only the main results will be summarized. The reader interested in details is referred to Vanden-Broeck [14] for the GC problem and to Forbes [5,6] for the ICE problem.

The first linear system for the ICE problem (i.e. the linear system obtained at the order *ϵ*) is
4.9
4.10
4.11
4.12
4.13
4.14
4.15

Similarly, equating the coefficients of order *ϵ*^{2} for the ICE system gives
4.16
4.17
4.18
4.19
4.20
4.21
4.22

The corresponding equations for the GC problem are obtained by replacing (*D*/*ρ*)*η*_{1xxxx} and (*D*/*ρ*)*η*_{2xxxx} by −(*T*/*ρ*)*η*_{1xx} and −(*T*/*ρ*)*η*_{2xx} in equations (4.11) and (4.18).

The solutions of equations (4.9)–(4.15) for the ICE problem are
4.23
4.24
4.25
provided
4.26
However, if the condition (4.26) is not satisfied for some integer value *m* of *n*, i.e. if
4.27
the solution is
4.28
and
4.29

We note that the relation (4.27) implies that linear waves of wavenumbers *k* and *mk* travel at the same speed.

The corresponding solutions for the GC problem are given by
4.30
4.31
4.32
provided
4.33
However, if the condition (4.33) is not satisfied for some integer value *m* of *n*, i.e. if
4.34
the solution is
4.35
and
4.36

The value of the constant *A*_{1} depends on the definition of the amplitude parameter *ϵ*. One possibility is to define the parameter *ϵ* as
4.37
where *a* is the first Fourier coefficient of *η*(*x*), i.e.
It then follows that *A*_{1}=*λ*. Other definitions of *ϵ* can be used.

The constants *A*_{m} that appear in the solutions (4.28), (4.29) and (4.35), (4.36) when the conditions (4.26) or (4.33) are not satisfied are arbitrary when only terms of order *ϵ* are considered. Their values are determined by solvability conditions at higher order. Such calculations for the GC system are classical and are summarized in Vanden-Broeck [14].

We shall present the details of the calculations of the coefficient *A*_{m} for the ICE problem in the case *m*=2. To simplify and to illustrate the main ideas, we assume that the water is of infinite depth.

Taking the limit in equations (4.23) and (4.27)–(4.29) gives 4.38 4.39 4.40 4.41

For *m*=2, equations (4.39)–(4.41) become
4.42
4.43
4.44
Combining equations (4.38) and (4.42) gives
4.45

Solving equation (4.16) by separation of variables and using equations (4.19) (with ) and (2.6) yield
4.46
where *F*_{n} are constant coefficients.

The periodicity and symmetry conditions (4.20) and (2.6) imply that *η*_{2}(*x*) can be represented by the Fourier series
4.47
where *E*_{n} are constant coefficients.

Substituting equations (4.46), (4.47), (4.43) and (4.44) in equations (4.17) and (4.18) and equating coefficients of and yield 4.48 and 4.49 Similarly, equating the coefficients of and yields 4.50 and 4.51 Solving equations (4.48) and (4.49) gives 4.52 Similarly, equations (4.50) and (4.51) yield 4.53 Using equations (4.38) and (4.42), it can easily be checked that the bracket on the left-hand side of equation (4.53) is equal to zero. Therefore, 4.54 Combining equations (4.52) and (4.54) yields 4.55 and 4.56 Substituting equation (4.55) into equations (4.43) and (4.44), we obtain 4.57 and 4.58

This concludes the determination of the constant *A*_{2}. It is interesting to note that the corresponding calculations for the GC problem yield the same value for *A*_{2} [14].

The above calculations demonstrate the non-uniqueness of periodic waves: we have obtained two different solutions when the condition (4.42) is satisfied (one with the plus sign in equation (4.57) and the other with the minus sign). These solutions are illustrated in figures 1 and 2.

They are known for the GC problem as the Wilton ripples. A similar non-uniqueness occurs when equation (4.39) is satisfied with *m*>2. However, the calculations for larger *m* become quickly intractable because it is necessary to calculate the solution to high order in *ϵ* in order to find the values of *A*_{m}. Such calculations were done numerically by Forbes [6]. They show the existence of many different families of solutions with dimples on their profiles. These solutions are qualitatively similar to those obtained for the GC problem [12–15].

### (b) Numerical procedure

The analytical approach of §4*a* provides useful information for waves of small or moderate amplitudes. However, it does not give information for very steep waves because only a few terms in the expansions can be calculated by hand. Therefore, numerical computations are required to calculate fully nonlinear waves. There are two basic approaches. The first is to calculate numerically many terms in the expansions (4.5)–(4.8) and then to construct solutions by accelerating the convergence of the expansions, for example, by using Padé approximants. This approach was pioneered by Schwartz [16] and later applied by others. The second approach is to compute directly fully nonlinear solutions, for example, by boundary integral equation methods or by series truncations. In this paper, we used the series truncation method. We shall describe this method for the ICE problem in water of finite depth. The corresponding method for the GC problem can be found in Vanden-Broeck [14].

Following Stokes [10], we choose the potential function *ϕ* and the streamfunction *ψ* as independent variables, and we define the complex potential *f*=*ϕ*+i*ψ*. We choose *ψ*=0 on the free surface and *ϕ*=0 at a crest of the wave. The flow domain in the complex *f*-plane is then the strip −*Q*<*ψ*<0. Here, −*Q* denotes the constant value of *ψ* on the bottom.

We satisfy the kinematic boundary condition (2.5) on the bottom by reflecting the flow into the bottom. The value of *ψ* on the image of the bottom is *ψ*=−2*Q* and the extended flow domain becomes then the strip −2*Q*<*ψ*<0.

Next, we map the strip −2*Q*<*ψ*<0 into the complex *t*-plane by the transformation
4.59
The flow domain in the *t*-plane is the annulus , where
4.60

Next, we introduce the complex velocity
4.61
where *z*=*x*+i*y*. Here, *u* and *v* are the *x* and *y* components of the velocity. We define the quantity *τ*−i*θ* by
4.62
One advantage of these new variables is that they yield a very simple formula for the curvature *K* of the free surface, namely
4.63

The function *τ*−i*θ* is an analytical function of *t* in the annulus . Therefore, it can be represented by the Laurent series
4.64
Since *ψ*=−2*Q* is the image of the free surface *ψ*=0 in the bottom, we have
4.65
and it follows from equation (4.64) that
4.66

Substituting equations (4.59) and (4.66) into equation (4.64) gives 4.67

We shall determine the coefficients *a*_{n} in equation (4.67) so that the dynamic boundary condition (2.4) is satisfied. To do so, we need first to rewrite equation (2.4) in terms of *τ* and *θ*. We note that
4.68
so that
4.69
and
4.70

Substituting equations (4.63) and (4.70) into equation (2.4), we obtain
4.71
To evaluate the last term on the left-hand side of equation (4.71), we note that
4.72
and
4.73
Here, equation (4.72) follows from the chain rule and equations (4.73) are the Cauchy–Riemann equations. It then follows that the last term on the left-hand side of equation (4.71) can be expressed in terms of derivatives (up to third order) of *θ* and *τ* with respect to *ϕ*. The calculation is simple but we shall not write the final expression because it is rather long.

Next, we introduce dimensionless variables by taking *c* as the unit velocity and *Q*/*c* as the unit length. We can then rewrite equation (4.60) as
4.74
Here, *l* is the dimensionless wavelength. Furthermore, equation (4.71) becomes in dimensionless form
4.75
where
4.76

We truncate the infinite series in equation (4.67) after *N*−2 terms and determine the *N*+2 unknowns (*a*_{n}, *n*=0,…,*N*−2), *B*, *F* and *l* by collocation. Thus, we introduce the *N* collocation points
4.77
and satisfy the dynamic boundary condition (4.75) at these points. This leads to *N*−1 nonlinear algebraic equations. One more equation is obtained by noting that equation (4.2) implies in dimensionless variables
4.78
An extra equation is obtained by imposing the value of *l* or *F*. The last equation fixes the amplitude of the wave by writing
4.79
where *s* is the given steepness (i.e. the difference of heights between a crest and a trough divided by the wavelength). This system of *N*+2 nonlinear algebraic equations with *N*+2 unknowns is solved by Newton’s iterations with the Jacobian evaluated by finite differences. The free surface profile is then given in parametric form by equations (4.69) and (4.70).

### (c) Numerical solutions in water of infinite depth

The numerical scheme of §4*b* was found to converge rapidly as *N* increases. All the results presented in this paper are independent of *N* within graphical accuracy. Most of the results presented were calculated with *N*=200. In this section, we present results in water of infinite depth (i.e. with *r*_{0}=0). Since as , the choice *Q*/*c* as the reference length is not appropriate for waves in water of infinite depth. Instead, we chose (*D*/*ρU*^{2})^{1/3} as the reference length for the calculations presented in this section.

We checked the accuracy of the procedure by computing numerically the solutions (4.57) with the plus sign (similar calculations can be done with the minus sign). As the amplitude of the wave tends to zero, *ϕ*≈*x* and *η*_{x}≈*θ*. Then a comparison between equations (4.57) and (4.67) (when *r*_{0}=0 and equation (4.45) is satisfied) shows that as the amplitude of the wave tends to zero. To check this analytical result, we computed solutions satisfying equation (4.45) for various values of *s*. We found *a*_{2}/*a*_{1}=0.79 for *s*=0.05, *a*_{2}/*a*_{1}=0.92 for *s*=0.01 and *a*_{2}/*a*_{1}=0.95 for *s*=0.0001. This provides a check on both the numerical procedure and the analytical approach.

As explained in §4*a*, there are many different families with dimples on their free surfaces. They can be calculated numerically by using the procedure of §4*b*. We present some typical profiles in figures 3–5 for *m*=3,4,5.

An interesting question is the determination of the limiting configurations of the waves as their amplitude is increased. In the corresponding GC problem, the waves become overhanging and ultimately approach limiting configurations with a trapped bubble at their troughs. This is clearly shown by the exact solutions derived by Crapper [17] and Kinnersley [18] for pure capillary waves and by the numerical results of Schwartz & Vanden-Broeck [12] for GC waves.

Forbes [6] suggested that overhanging waves should also occur with the ICE system. This is confirmed by our numerical solutions. Typical profiles are shown in figure 6 for a dimensionless wavelength equal to 16. Each profile corresponds to a different value of *s*.

## 5. Generalized solitary waves for the ICE system

Periodic waves in water of finite depth have properties similar to those described in §4*a* for periodic waves in water of infinite depth. In particular, there are many branches of solutions with dimples on their free surface profiles. In this section, we examine the behaviour of periodic waves in water of finite depth for *r*_{0} close to 1 (i.e. in very shallow water).

Figure 7 shows values of 1/*F*^{2} versus the dimensionless wavelength
5.1
for *β*=0.07 and
5.2
The condition (5.2) fixes the difference of heights between the crests and the troughs of the wave.

There are many different families of waves and only two are shown in figure 7 (the solid curves). Two typical free surface profiles are presented in figure 8. The one of shorter wavelength corresponds to a point on the solid curve on the left of figure 7 and the one of longer wavelength to a point on the solid curve on the right of figure 7. As expected, there are ripples on the profiles. For values of *r*_{0} close to 1, these ripples tend to concentrate in the troughs of the waves. The results of figure 8 show that as one moves in figure 7 from one family to another one further on the right, an even number of crests and/or troughs per wavelength are added in the train of ripples. For example, there are 19 crests and 18 troughs per wavelength in the profile of figure 8 of shortest wavelength (i.e. the profile that terminates at a value of *x*<30) and 21 crests and 20 troughs per wavelength in the profile of figure 8 of longer wavelength (i.e. the profile that terminates at a value *x*>30).

For *l* sufficiently large, there are an infinite number of curves similar to the two solid curves shown in figure 7. Each curve can be obtained from a previous one of the left by translating it horizontally by the wavelength of the ripples. For *l* sufficiently large, these parallel curves approximate a family of solitary waves (i.e. a non-periodic wave). These solitary waves are characterized by an infinitely long train of ripples in the far field. They are often referred to as ‘generalized solitary waves’ to contrast them with true solitary waves that have a flat free surface in the far field.

We note that for given values of *β* and *sl*, we have a one parameter family of generalized solitary waves (described by the limit of the solid curves of figure 7 as ). The parameter can be, for example, chosen as the Froude number. It is interesting to contrast this property with what is predicted by the Korteweg de Vries (KdV) equation. The KdV equation describes waves of small amplitude in shallow water. Its derivation can be found in Whitham [19] or in Vanden-Broeck [14]. For the ICE problem, it is found to be the classical KdV equation of pure gravity waves in the absence of the ice sheet. The reason being that the last term on the left-hand side of equation (2.4) does not contribute to the asymptotic expansions up to the order considered in the KdV equation. The corresponding solutions (known as cnoidal waves) do not have ripples on their troughs. The quantity 1/*F*^{2} is then a function of *l*. It is shown in figure 7 by the broken line. In the limit as , the broken line approaches asymptotically the value 1/*F*^{2}≈0.88 and the corresponding profile is a true solitary wave. Therefore, the KdV equation predicts a unique true solitary wave, whereas the complete set of equations predicts a family of generalized solitary wave as . This shows that the KdV equation does not describe accurately the fully nonlinear problem. A similar situation was found before for the GC problem (see [14] and references therein).

## 6. Conclusions

We have presented weakly nonlinear and fully nonlinear solutions for waves under an ice sheet. Periodic and generalized solitary waves were discussed. Our findings complement those in Forbes [5,6] and Parau & Dias [7]. In particular, we have stressed the similarities between the present problem and that of GC waves.

## Acknowledgements

The research presented in this paper was supported by EPSRC.

## Footnotes

One contribution of 13 to a Theme Issue ‘The mathematical challenges and modelling of hydroelasticity’.

- This journal is © 2011 The Royal Society