Solutions of the nonlinear water wave equations under an ice sheet are computed using a boundary integral equation method. The ice sheet is modelled as a thin elastic plate and the fluid equations are nonlinear. Depending on the velocity of the moving disturbance generating the flow, different types of responses of the floating ice sheet are discussed.
In this paper, three-dimensional flexural–gravity waves generated by a disturbance moving at a constant velocity U on top of a floating ice sheet are considered. The research is motivated in part by applications to transport systems in cold regions, where the frozen water is used for roads or aircraft runaways and where air-cushioned vehicles are used to break the ice [1–3]. There have also been a number of experiments performed with loads moving on ice plates: the McMurdo Sound experiments in deep water in Antarctica  and the Lake Saroma experiments in shallow water in Hokkaido, Japan .
The ice sheet is modelled as a thin elastic plate, and the linear elastic plate equations are used. This model is also useful in studying waves beneath a very large floating structure, such as the one which was constructed in Japan as a prototype for a floating airport (Mega-Float, see ). This approximation is often appropriate , but there are other possible models such as a viscoelastic model  or a finite ice thickness model . In this study, the fluid underneath the plate is of infinite depth and it is assumed to be incompressible and inviscid, while the flow is irrotational. Fully nonlinear equations are used for the flow.
As we shall see, waves under an ice sheet have properties similar to those of gravity–capillary waves. In recent decades, many analytical and numerical results have been obtained for gravity–capillary waves (see  for a review). Some experiments were also recently performed (see  for two dimensions and  for three dimensions) for gravity–capillary waves. However, it should be easier to perform experiments for flexural–gravity waves, as the wavelength of these waves is in metres rather than in millimetres or centimetres, as in the gravity–capillary problem.
There have been some theoretical studies in recent decades on the linear three-dimensional flexural–gravity problem, where the fluid flow is assumed to be linear, as well as on the elastic plate. Davys et al.  considered waves generated by a concentrated point source and obtained wave crest patterns using asymptotic Fourier analysis. More recently, Milinazzo et al.  have computed the linear deflection generated by the uniform motion of a rectangular load using Fourier integrals and the method of residues. The present study complements the findings in these linear studies, by studying the effect of the nonlinearity of the fluid flow.
The type of solutions for this problem depends on the velocity of the moving pressure. When the velocity of the moving pressure on top of the floating ice plate is less than the minimum phase speed cmin predicted by the linear dispersion relation, the disturbance is localized and the waves do not propagate away from the pressure. If the velocity of the pressure approaches this minimum, the amplitude of the waves increases and some oscillations start to appear before and after the moving pressure. For velocities of moving pressure greater than the minimum phase speed cmin, a more complicated pattern emerges, with waves of different wavelengths appearing in front of and behind the moving pressure. Since waves appear both at the front of and at the back of the disturbance in this case, the radiation condition cannot easily be imposed. We adapt, for this nonlinear problem, a technique introduced by Rayleigh  to calculate analytically linear solutions for the gravity–capillary waves problem. This includes a dissipative term in the dynamic boundary condition, known as the artificial Rayleigh viscosity. This technique was successfully employed by Părău et al.  to compute forced nonlinear three-dimensional gravity–capillary waves. It is worth noting that this artificial viscosity does not change the irrotational character of motion [15, art. 242].
Two-dimensional nonlinear hydroelastic waves are considered by Vanden-Broeck & Părău . They complement the weakly nonlinear results obtained by Părău & Dias , who derived a forced nonlinear Schrödinger equation for the ice sheet deflection (see also ), and more recent numerical results by Bonnefoy et al. , which are based on higher order spectral methods. They use a nonlinear plate/nonlinear flow model, in which the ice sheet is modelled by a Kirchhoff–Love plate, and a nonlinear term involving the plate curvature is used in the dynamic boundary condition. In particular, two-dimensional solitary waves with decaying oscillations beneath an ice sheet travelling with velocity U, where 0<cmin−U≪1, have been studied previously by Il’ichev , Părău & Dias  and Dias & Iooss . Very recently, Plotnikov & Toland  have rigorously derived another nonlinear model for hydroelastic waves in two dimensions.
In §2, we describe the governing equations for the problem. The dispersion relation for linear plane periodic waves is studied in §3. The numerical method based on the boundary integral equation method is briefly described in §4. Some of the solutions obtained with this numerical method describing forced steady waves are investigated in §5, with a discussion of the different possible cases.
2. Formulation of the problem
We consider three-dimensional free-surface flows due to a distribution of pressure moving at a constant velocity U on a surface of an ice sheet of thickness h, which lies on top of a fluid of infinite depth. The sketch of the flow is shown in figure 1.
The fluid is assumed to be incompressible and inviscid, and the flow to be irrotational. We introduce Cartesian coordinates x, y, z with the z-axis directed vertically upwards and the x-axis in the direction of the velocity of a moving load U. Gravity is acting in the negative z-direction. We denote by z=ζ(x,y,t) the vertical ice sheet deflection, where z=0 is the undisturbed interface between the ice sheet and the fluid.
In terms of the fluid velocity potential Φ(x,y,t), the equations of the flow are 2.1with the boundary conditions at the interface between the ice sheet and the water 2.2 2.3 and 2.4where ∇4=∂xxxx+∂yyyy+2∂xxyy is the bi-Laplacian operator. Equations (2.2) and (2.4) are the kinematic boundary condition on the free surface and the boundary condition at infinity. Equation (2.3) is the dynamic boundary condition obtained after use of the Bernoulli equation. The dynamic boundary condition depends on the model used to describe the ice sheet. There are many different models for the ice sheet , and here we assume that it satisfies the equations of a thin elastic plate. We neglect the inertia of the thin plate, so the plate acceleration term is not considered here (see  for details), and we assume that there is no initial tension in the plate. The constant D=Eh3/12(1−ν2) is the flexural rigidity of the plate, ρ is the density of the fluid, g is the acceleration due to gravity, h is the ice thickness, E denotes Young’s modulus, ν is the Poisson ratio for ice and p(x,y,t) is the external pressure distribution exerted on the ice sheet. We can choose different formulae for the external pressure, and we will specify it later (see §4).
3. Dispersion relation
By linearizing the problems (2.1)–(2.4) and looking for free waves of the form ei(lx+my−ωt), we obtain the following dispersion relation: 3.1where . By considering uniform plane waves with wavenumber k in the direction of the x-axis, we find that the phase speed c(k)=ω/k is given by 3.2It can be easily shown that the phase speed c(k) tends to for k→0 or and that it always has a minimum given by At this minimum, the group and phase speed are equal. If we consider a moving pressure with U<cmin, the linear solution approaches a uniform flow at infinity since there are no waves with U<cmin (figure 2). For U>cmin, the linear solution is characterized by trains of waves in the far field, with gravity waves of wavenumber kg downstream of the moving pressure and flexural waves of wavenumber kf upstream (figure 2).
In the three-dimensional case, an analysis of the dispersion relation (3.1) was performed by Davys et al.  and steady wave patterns have been studied. The system of steady waves generated by a moving pressure is more complicated than in two dimensions, with different types of waves behind and in front of the pressure. The waves behind the pressure are mainly longer gravity waves and are confined to an angle, while the waves in front of the pressure are shorter flexural waves, curved around the support of the pressure. This angle varies, depending on the ratio U/cmin (see  for more details).
4. Numerical scheme
We are interested in steady solutions of the problem (2.1)–(2.4) and we choose a frame of reference moving with the disturbance. We substitute x−Ut by x, Φ(x,y,z,t) by Φ(x−Ut,y,z)−Ux and ζ(x,y,t) by ζ(x−Ut,y). The pressure distribution that models the disturbance is chosen as 4.1Here, P0 is a constant and L defines the size of the support of the pressure. We non-dimensionalize equations (2.1)–(2.4) by using U as the unit of velocity and L as the unit of length. In the moving frame, the equations of motion and the boundary conditions to be solved become 4.2 4.3 4.4 and 4.5where and ε is the dimensionless magnitude of the pressure.
We can introduce a non-dimensional parameter independent of L, It can easily be shown that the critical speed U=cmin corresponds to .
The numerical scheme used to solve the problem is based on the scheme described by Părău & Vanden-Broeck . It uses a boundary-desingularized integral-equation method introduced initially by Landweber & Macagno  for the problem of uniform flow past an ellipsoid and generalized by Forbes  for three-dimensional gravity free-surface flows past a source.
Only the main points of the formulation and of the numerical procedure are presented here. The reader is referred to Părău & Vanden-Broeck  for details. The formulation involves applying Green’s second identity to the functions Φ−x and the three-dimensional free-space Green function where P=(x,y,z), P*=(x*,y*,z*) and P,P*∈SF on a region V consisting of a large-radius hemisphere bounded above by the interface SF between the ice and the water, except for a small hemisphere around the point P*(x*,y*,z*) on SF (figure 3).
We obtain, after some manipulations, where n is the unit normal vector pointing into the fluid. The integral is made non-singular by adding and subtracting a quantity that can be evaluated in a closed form [24–26] and that is then projected onto the xy-plane. It is convenient to introduce the potential ϕ(x,y)=Φ(x,y,ζ(x,y)) on the ice sheet/water interface.
The discretization involves a regular grid with n points in the x-direction, approximating , and m points in the y-direction, approximating . The uniform mesh sizes on the x- and y-axes are denoted by Δx and Δy. The 2nm unknowns are ζx and ϕx, evaluated at the grid points, and the equations are obtained from the evaluation of the integral and the dynamic boundary condition (4.4) at midpoints. The values of ζ and ϕ are obtained by integrating ζx and ϕx with respect to x by the trapezoidal rule. The values of the dependent variables at midpoints are found by interpolation, and the other derivatives are calculated by using finite differences. In particular, the bi-Laplacian ∇4ζ is discretized by centred finite differences in the x- and y-directions with 13 points . Some truncation and radiation conditions are imposed for the further points downstream (ζx=0,ϕx=1) and for the further points upstream (ζ=0,ϕ=x). On the lateral boundaries, no conditions were imposed, the y-derivatives of the variables being computed by one-sided difference formulae. The 2nm nonlinear equations are solved by a modified Newton’s method. As the initial guess, a flat surface (ζ=0,ϕ=x) everywhere was considered in most cases for forced waves.
In the case U>cmin, the radiation condition cannot easily be imposed, since there are many waves, before and after the moving pressure. We have followed for this case Părău et al.  and introduced a dissipative term in the dynamic boundary condition (4.4). This term is characterized by an artificial viscosity μ>0, known as the Rayleigh viscosity in the case of gravity–capillary waves. It was first introduced by Rayleigh  to obtain a unique solution for the linear problem of gravity–capillary waves by taking the limit μ→0. The dynamic boundary condition (4.4) becomes, in this case, 4.6 In the limit μ→0, we should recover the correct solutions to our problem. However, we cannot take μ arbitrarily small in our numerical simulations, but we can estimate the influence of μ≠0 on the solutions.
The accuracy of the results was checked by varying the number of points n and m and the grid intervals Δx and Δy. An example of the accuracy check is presented in figure 4, where solutions are computed for two values of the grid interval. The solution on the left is calculated with Δx=Δy=0.3 and the solution on the right with Δx=Δy=0.2. A very good agreement was found between the two solutions.
We will present next some computations for forced waves with U<cmin; then, we will discuss some solutions for forced waves with U>cmin. Unless otherwise stated, the results presented here correspond to n=80, m=40 and Δx=Δy=0.6.
We compute here solutions using the dynamic boundary condition (4.4), with parameters fL and β. The minimum phase speed is reached for
For velocities U much smaller than cmin, the deformation of the ice sheet is highly localized near the support of the pressure (figure 5). For symmetric pressures p with respect to x=0, the solutions are also symmetric. The steady displacement of the ice ζ(x,y) is similar to the static deformation of the ice sheet under a load.
When the velocity of the pressure is increased, the depth of the depression under the load starts to increase (figure 5). As the velocity of the pressure U approaches cmin, the area where the displacement of the ice sheet ζ(x,y) is visible increases considerably and symmetric decaying oscillations start to appear in front of and at the back of the perturbation. The solutions still decrease monotonically in the y-direction. This can be observed if we present the centrelines of the solutions in the plane y=0 in figure 6.
Solutions with negative pressure (ε<0) can also be computed. An example is presented in figure 7. The pattern is similar to the positive pressure case, but the amplitude of the solution at (x,y)=(0,0) is now positive.
By keeping the velocity U constant, the influence of the flexural rigidity of the elastic plate on solutions is analysed by varying β. In figure 8, we observe that the amplitude of the solutions increases as β decreases. This is to be expected, as by decreasing β, which corresponds to decreasing the flexural rigidity D, the minimum phase speed cmin decreases as well, so the pressure moves at a velocity closer to this minimum.
In this section, equation (4.6) is used as the dynamic boundary condition, with the artificial viscosity μ included. The influence of μ is analysed first for some solutions with U<cmin. A solution is presented in figure 9, for fL=0.8, β=0.5 and μ=0.1. It can be observed that the solution becomes asymmetric, as a consequence of including the artificial viscosity in the dynamic boundary condition. The elevation wave in front of the pressure is higher than the elevation wave after the pressure. A similar effect was observed when the Rayleigh viscosity was introduced in the dynamic condition for the gravity–capillary waves . This effect is clearer in figure 10, where only the centrelines in the x-direction are plotted for different values of μ. It can also be observed that the artificial viscosity damps slightly the amplitude of the depression under the support of the pressure.
A typical solution for U>cmin is presented in figure 11. Two different wave system are found: one in front of the pressure with a shorter wavelength, and one V-shaped after the pressure with a longer wavelength. The waves behind the pressure are confined in an angle that varies with the velocity of the moving pressure. The waves in front (mainly flexural waves) are curved around the support of the pressure. The highest positive elevation is observed to occur just in front of the pressure.
The influence of varying fL is studied in figure 12. More waves become visible in front of the pressure when fL is decreased. The wavelength of the waves in front of the pressure decreases when we decrease fL. The depression under the pressure becomes wider and shallower. The angle where the waves are confined behind the pressure decreases when fL decreases (figure 13).
As above, we can study the effect of varying the flexural rigidity of the plate, i.e. varying β and keeping the other parameters constant. We observe that, when β increases, the amplitude of the solution under the support of the pressure decreases and the wavelength of the waves in front of the pressure increases (figure 14). The amplitude of the solution will start to decrease again when β reaches the value corresponding to the critical velocity cmin, which for the case presented is β≈3.9.
In the previous graphs, only one or two waves are observed behind the pressure, as the difference in wavelength between the two systems of waves is quite significant. To obtain waves with comparable wavelengths in front of and behind the pressure, solutions with U close to cmin are computed (figures 15 and 16). The waves are nearly two dimensional close to this critical velocity.
It should also be noted that, when the artificial viscosity is present, there is a smooth transition between the region U<cmin and the region U>cmin, as in the gravity–capillary waves problem . This can be observed in figure 17, where the maximum value of |ζ| is plotted against fL. The solution with the maximum amplitude corresponds to U≈cmin. If no artificial viscosity is present (μ=0), the solutions computed for U<cmin cannot be extended in the region U>cmin.
6. Conclusions and future work
We have calculated three-dimensional forced flexural–gravity steady waves for the full nonlinear Euler equations in water of infinite depth. Fully nonlinear boundary conditions have been imposed at the interface between the ice sheet and the fluid. The ice sheet is modelled as a thin elastic plate. The form of the solutions generated by a moving disturbance depends on the value of the velocity U of the moving pressure relative to the minimum phase velocity cmin. In the case U<cmin, a localized symmetric solution is found, with damped oscillations appearing in front of and after the support of the pressure as the velocity U approaches cmin. In the case U>cmin, an artificial viscosity is used in the dynamic boundary condition in order to impose the radiation condition, similar to the Rayleigh viscosity used in the gravity–capillary waves problem. The properties of these solutions are analysed. The solutions are no longer localized near the support of the pressure, with shorter flexural waves travelling in front of the pressure and longer gravity waves behind the pressure.
This study can be extended in the future in a number of directions. One of these is to consider waves in water of finite depth H. The existence of solitary waves with velocities in the region U<cmin needs to be studied. Some preliminary results have been reported in Părău & Vanden-Broeck . Similar three-dimensional solitary waves for a weakly nonlinear model for flexural–gravity waves in finite depth were found by Hărăguş-Courcelle & Il’ichev . Using the nonlinear flow/linear plate equations they derived a three-dimensional generalization of the fifth-order Korteweg–de Vries equation for ice sheet deflection. In the case of gravity–capillary waves, fully localized three-dimensional solitary waves have been computed by Părău et al.  in infinite depth and Părău et al.  in finite depth. They were also studied theoretically for the fully nonlinear gravity–capillary problem by Groves & Sun , and for some weakly nonlinear models by Kim & Akylas  and Milewski . In the related problem of three-dimensional interfacial gravity–capillary waves, they have also been computed by Părău et al. .
In water of finite depth another critical velocity may become important: the long-wave limit , where H is the depth of the water . In two dimensions, its effect was discussed by Schulkes & Sneyd . We expect that it will play an important role in the pattern of solutions, so a boundary integral method based on a Green function that incorporates the boundary condition at the bottom of the fluid needs to be developed .
We also plan to include a more complex model for ice, and in particular to use a nonlinear model for the plate. In two dimensions, Peake  has discussed the different possible combinations: linear/nonlinear plate and linear/nonlinear flow. In three dimensions, there are fewer nonlinear models for the plate that can be used, so an investigation of the different options will be performed.
Another possible extension is to consider the water under the ice sheet to be stratified. The stratification can be the effect of freezing or melting of the ice  and a pynocline or thermocline often exists. The simplest model is to consider two homogeneous layers of fluid of different densities resting on top of each other, with the upper fluid bounded above by an elastic plate. Părău et al.  have considered waves in a similar configuration, but in the absence of the elastic plate.
A time-dependent algorithm should be developed to study the evolution with time of flexural–gravity steady waves. A starting point could be an algorithm developed recently by Părău et al.  for time-dependent gravity–capillary flows.
The research presented in this paper was carried out on the High Performance Computing Cluster supported by the Research Computing Service at the University of East Anglia. The authors also wish to thank the International Centre for Mathematical Sciences, Edinburgh, UK, for their help in the organization of the workshop Mathematical challenges and modelling of hydroelasticity in June 2010, where ideas exchanged during discussion sessions were useful for this study. The work was also supported by EPSRC.
One contribution of 13 to a Theme Issue ‘The mathematical challenges and modelling of hydroelasticity’.
- This journal is © 2011 The Royal Society