Solitary flexural–gravity waves in three dimensions

The focus of this work is on three-dimensional nonlinear flexural–gravity waves, propagating at the interface between a fluid and an ice sheet. The ice sheet is modelled using the special Cosserat theory of hyperelastic shells satisfying Kirchhoff's hypothesis, presented in (Plotnikov & Toland. 2011 Phil. Trans. R. Soc. A 369, 2942–2956 (doi:10.1098/rsta.2011.0104)). The fluid is assumed inviscid and incompressible, and the flow irrotational. A numerical method based on boundary integral equation techniques is used to compute solitary waves and forced waves to Euler's equations. This article is part of the theme issue ‘Modelling of sea-ice phenomena’.


Introduction
In this paper, we consider an incompressible and inviscid fluid covered by an ice sheet. Under certain conditions, the floating ice sheets can be modelled as an elastic medium, and their interaction with the fluid makes the resulting hydroelastic problem mathematically challenging [1]. A renewed interest in waves generated by moving loads on top of floating ice sheets has been sparked in the last 40 years by a series of experiments in cold regions (e.g. [2][3][4][5][6]), where ice roads and runways are used during the winter.

Formulation
The Euler equations describing an incompressible, inviscid, irrotational fluid with velocity potential Φ(x, y, z, t) and variable surface ζ (x, y, t) under a sheet of ice in three dimensions are given by (2.1) where ρ is the density of the water, D is the flexural rigidity of ice, g is gravitational acceleration, h is the depth of the water (which could be taken to be infinite), p(x, y, z, t) is an external pressure exerted on the ice and P flex is a term describing the effect of the ice on the surface of water. The schematic for the domain of interest is shown in figure 1, where the Cartesian coordinates Oxyz are defined such that the z-axis is oriented vertically upward and the waves propagate along the x-axis.
The term describing the effects due to the presence of ice is modelled using the Cosserat theory of hyperelastic shells [11]. It assumes that the ice is a thin elastic plate with constant thickness, the ice bends with the water waves and it cannot break. The model neglects friction between the ice and the water. The coefficient for flexural rigidity for ice D is given by where E is Young's modulus (i.e. the modulus of elasticity), ν is the Poisson ratio describing the effects of transverse strain relative to axial strain and d is thickness of the ice. After some algebra, we can express P flex in Cartesian coordinates as where H is the mean curvature and K is the Gauss curvature of the ice-water interface, given by It is worth noting that in previous works on three-dimensional flexural-gravity waves [44,45] a linear model was used for the ice sheet, where P flex has a much simpler form as We analyse the linearization of the nonlinear system (2.1) and look for plane wave solutions of the form e i(k 1 x+k 2 y−ωt) and k = k 2 1 + k 2 2 . The linear dispersion relationship between the frequency ω and the wavenumbers k 1 and k 2 is then  If we only consider two-dimensional waves moving in the x-direction with wavenumber k, then we can further define a phase speed It can be shown that c always has a minimum c min at a finite wavenumber k = k min > 0 for all the values of the physical parameters. At this critical speed, the group and phase speeds are equal.  waves are generated in the far field. However, for U > c min , there are two wavenumbers for which U = c(k) and a more complicated pattern emerges. This analysis generalizes to three dimensions and, in this work, we focus on the case of U < c min and we will look for symmetric solutions in x and y.
Since we restrict our focus to solitary and forced waves, travelling with speed U, it is convenient to use a reference frame moving with the wave by setting x → x − Ut and considering steady-state solutions. We non-dimensionalize the problem using this speed U as the unit of velocity and introduce a unit of length L. This allows us to rewrite the Bernoulli condition at the surface (the third equation in (2.1)) in terms of the non-dimensional parameters, using We introduce β = D/ρU 2 L 3 and for convenience L = 3 D/ρU 2 , which will set β = 1. The nondimensional depth is defined asĥ = h/L, dropping the hat for ease from now on. As mentioned above, there is always a minimum of the dispersion relation, whatever the values of physical parameters. In figure 3, we plot the curve in the F-h plane, which corresponds to this minimum. To solve equations (2.1), we use a boundary integral equations method (e.g. [40,[42][43][44][45]), using Green's functions. We now review the main aspects of the numerical scheme. We note that, for Laplace's equation (the first equation in (2.1)) in three-dimensional free space, Green's function in infinite depth for the points P = (x, y, z) and P * = (x * , y * , z * ) is We use Green's second identity, which states where we can set g = G(P, P * ) and f = Φ − x, which satisfies Laplace's equation. After the proper substitutions and some manipulation, we obtain where n is the normal to the ice-water interface S pointing into the fluid, and P * is a point from S. We define which allows us to rewrite the Laplace equation as well as some of the boundary conditions in terms of a surface integral. The final form of equations to solve for flexural-gravity waves in infinite depth is and We use the notation ζ = ζ (x, y), ζ * = ζ (x * , y * ) and φ * = φ(x * , y * ). Symmetry in the y-direction with ζ (x, y) = ζ (x, −y) and φ(x, y) = φ(x, −y) implies that we can solve the set of equations (2.11) and (2.12) on half of the domain. We use the method of images to account for the symmetry, which introduces two extra terms in each kernel Part of the integral in (2.12) is singular [10] and we remove the singularity by noting that The last term in (2.15) can be computed analytically since it is of the form (1/z) dz = ln z. Moreover, if we consider a fully symmetric solution in x and y, then and this additional symmetry in x will also introduce two more terms to each kernel. This implies that, using Green's function approach, we will have four terms in each kernel due to the method of images. We can thus reduce the computational cost by only integrating over a quarter of the domain. where (2.20) and where, as before, I 2 is singular.

Numerical method
The equations given by (2.11) We evaluate the equations at points (x i+1/2 , y j ) and use finite differences for the derivatives. This gives 2(N − 2)M equations. We obtain 2M equations from symmetry about the y-axis and 2M more equations from decay at the boundaries. To obtain ζ and φ, we use a trapezoid rule. Derivatives are computed using central difference, second-order accurate schemes, except at the border where one-sided schemes are needed. This gives 2NM equations, which can be written as To solve the system, we use Newton's method [46,47], which is summarized as (i) set up an initial guess u 0 (ii) until convergence (a) solve J(u n )δ n = −G(u n ) (b) set u n+1 = u n + δ n (c) test for convergence and repeat steps if not converged.
This method relies on an initial guess u 0 and the Jacobian matrix J. To compute non-trivial solitary waves (ones that are non-zero), we need a good initial guess in order for the numerical method to converge. However, since the set of equations is complicated and very nonlinear, there is no analytical expression for it. In order to obtain solitary waves, the numerical procedure is to compute a branch of forced waves, reduce the forcing, and use that as an initial guess. The particular form of the forcing pressure is not important. The main requirement is that it has a bounded support. Here we choose p(x, y) = p 0 e 1/(x 2 −L 2 0 ) e 1/(y 2 −L 2 0 ) for |x| < L 0 and |y| < L 0 0 otherwise, and, for convenience, we choose L 0 = 1 and p 0 to be small, positive or negative. To obtain fully nonlinear waves, we use a continuation method. We start with a small amplitude, using the flat water as an initial guess with F > F min , where F min corresponds to the minimum of (2.4). We use Newton's method to compute the true solution for the chosen F, and then rescale to use as an initial guess for a larger amplitude solution. This process is continued until we go around a turning point in the bifurcation branch for the forced solutions. Then, decreasing p 0 , we can compute non-trivial solutions to the unforced (free surface with p 0 = 0) problem. We expect bifurcation of solitary waves to start from a point on the curve shown in figure 3 corresponding to the depth h considered, and solitary waves to exist for F > F min . Similar approaches were used before for gravity-capillary waves for the two-dimensional problem [48] and for the three-dimensional problem [42,43].
Another key component to Newton's method is the Jacobian matrix J. Owing to the form of the equations we are solving, (2.11) and (2.17), this Jacobian is hard to compute analytically. Therefore, it is computed numerically via finite differences. For a three-dimensional problem, this matrix can be quite large and it requires a lot of memory to store and also to compute its inverse. In particular, due to the high degree of the nonlinearity from the flexural term (2.2) and the presence of the double integral over the whole space, this Jacobian matrix contains at least two dense quadrants. In general, there are two ways to reduce the computational requirements of Newton's method. One way is to compute an approximate inverse of the Jacobian and another way is to compute an approximate Jacobian [47]. It is known that iterative solvers such as GMRES can compute approximate inverses, but these require preconditioners to converge well. For this problem, preconditioners can be found from previous steps along the bifurcation branch; however, they require a lot of memory to store so this process is inefficient. It is faster to not update the Jacobian at every step along the bifurcation branch, hence we use only an approximate Jacobian. However, for this flexural-gravity wave problem with the nonlinear term due to the Cosserat model used, having an inexact Jacobian was found to result in poor convergence.

Numerical results
In this section, we present the numerical results for forced and free surface waves. For what follows, we use the doubly symmetric equations on a quarter of the domain with n = 80, m = 50, x = 0.5, y = 0.8. The accuracy of the method was checked by varying x and y and the number of points in each direction. The tolerance for convergence of Newton's method was set to 10 −10 . The forcing term (3.3) had p 0 = 1 for depression waves and p 0 = −1 for elevation waves (unforced p 0 = 0). Figures 4 and 5 summarize the computations using a continuation method with a Newton iteration at every step. These figures show the bifurcation branches for both depression and elevation waves with and without a forcing. Dashed lines show the forced waves, which were used as a starting point for the computations of free solitary waves and computed for U < c min , and solid lines are the solitary waves without a forcing. Figure 6 shows the comparison of the branches of solutions of solitary waves, normalized by F min in each case. The dashed line is for infinite depth and the solid line is for h = 1. For infinite depth, it can be shown that in non-dimensional units F min ≈ 0.4725 and for h = 1, F min ≈ 1.0272. These are the values used to normalize the horizontal axis. If these solitary branches bifurcated from F min , then they would be shown to originate fromF = 1 in figure 6. However, this is not the case and the authors believe this is due to numerical errors introduced by truncation or a grid that is too coarse. These numerical difficulties are related to the fact that the solitary waves become more and more oscillatory and less localized as the bifurcation points are approached. Similar difficulties have been encountered previously for gravity-capillary waves (e.g. [42]). Figure 6 (taken together with figures 7-10) shows that solitary waves in finite depth have a smaller free surface slope than those for infinite depth for a similar distance (in F) from the   infinite depth, these branches do not bifurcate from zero amplitude, whereas in shallower water they do. It becomes difficult to follow the branches numerically for large-amplitude solutions. In particular, the elevation branch may have turning points as in two dimensions [49], but more grid points are needed for a good accuracy, and this become prohibitive computationally. Solitary waves are shown in figures 7 and 8 for infinite depth and figures 9 and 10 for h = 1. The full wave profiles are shown on the left and the cross-sections on the right. As anticipated from figure 6 with the infinite depth branches having higher amplitude, we see that, in the wave profiles as well, infinite depth waves show more oscillations in the x-direction, but finite depth waves are less compact in the y-direction. On the elevation branches, the elevation of the wave at the centre point (x, y) = (0, 0) becomes smaller in magnitude than the depressions which are adjacent to it, similar to the gravity-capillary wave problem [42].

Conclusion and discussion
In this work, we employ the boundary integral equations method [40,45] to compute threedimensional solitary flexural-gravity waves for water of finite and infinite depth, covered by an ice sheet, using the model presented by Plotnikov & Toland [11]. This model is very nonlinear with a high number of derivatives, as shown in (2.2), and therefore is numerically difficult to implement. In this case, the usual numerical aids such as iterative methods for matrix inversions and use of inexact Jacobians have not proved to be successful. However, working on a quarter of the domain and using the full symmetry, we are still able to compute the branches of forced and unforced solutions shown in figures 4-6. Sample figures of the wave profiles are also provided, allowing for comparison of solutions in infinite and finite depth as well as depression and elevation waves. These waves are similar to solutions for gravity-capillary waves [43].
One issue is still left to resolve, and that is the asymptotic predictions for the start of the bifurcation branch. It was shown in [36] that flexural-gravity solitary waves in infinite depth bifurcate from a finite-amplitude solution. However, numerically this is difficult to prove. Owing to computational cost, there is a limit to how well these waves can be resolved, resulting in numerical errors due to truncation close to the start of the bifurcation branches, as shown in figure 6, where the branches do not begin at the predicted values of F. In shallow water, such as here with h = 1, the numerical results suggest that the branch starts from zero amplitude, as predicted by the weakly nonlinear model [36]. There has already been work on using highperformance computing techniques to improve the accuracy of similar computations [46], but more methods need to be developed for nonlinear regimes with higher order derivatives such as the model used for flexural-gravity waves.
It is worth noting that some large-amplitude solutions calculated here may become unphysical, as the strain of the ice plate may be higher than the yield strain of ice. In this case, the elastic model for the sheet will become unrealistic and different models should be used (see [32] for a more detailed discussion).
Data accessibility. This article has no additional data. Authors' contributions. All the authors had an equal contribution. All authors gave final approval for publication. Competing interests. We declare we have no competing interests. Funding. This work was partially supported by EP/J019305/1 for E.I.P., EP/J019321/1 for P.M. EP/J019569/1 for J.-M.V.-B and O.T. This work was also supported by: EPSRC grant nos. EP/K032208/1. Parts of the research presented in this paper were carried out on the High Performance Computing Cluster supported by the Research and Specialist Computing Support service at the University of East Anglia.