Model of attenuation of long waves under continuous ice layer

In this work new mathematical model of long wave propagation on water surface with ice cover is proposed. The model of thin elastic plate is used to describe ice layer movement. Equation for ice cover contain additional term to takes into account dissipation effects in the ice cover to explain wave attenuation. Proposed model was reduced to one nonlinear evolution equation for water level perturbation. The expression for wave energy was obtained under assumption of long waves. Proposed model is numerically studied, energy of system is computed. Obtained results are compared with results of suggested before model that takes into account the flow law of Glen.


Introduction
The process of attenuation of surface wave under ice layer is well studied at the present time. There are some models for description energy loss in the system ice-water under wave propagation. Almost all proposed models consider system under linear theory approximation.
The one model was proposed by Wadham [1]. He takes into account the Glen law to describe energy lose in ice layer due to creep. He consider the propagation of a one-dimensional monochromatic wave propagating at right angle to the ice edge. As a result, he obtain that energy of wave is directly proportional to expression 1 (n − 1)Sx + 1/A n−1 i 2 (1) where S is some constant depending on system and wave properties, A i is initial amplitude of wave and n is power in flow law of Glen. Wadham find that attenuation rate of wave is best fitted by exponent n = 3.
Another approach to problem was proposed by Weber [2]. He assume that ice floes very small and treat ice field as a viscous Newtonian fluid. Moreover he studies case of infinite rotating ocean. His model fit well field data from the marginal ice zone.
Considering only linear effects gives good approximation for wave propagation over few tens of wavelengths [3]. But waves in ocean can propagate over longer distance and nonlinear effects will accumulate and affect the wave form. So a nonlinear model for describing considered problem is needed. The model that takes into account Glen flow law proposed by Wadham can not be directly applied in this case.
We consider a problem of long wave propagation under ice layer. In this work the new phenomenological model is proposed. Thin plate model is used for ice layer. Proposed model In nature the waves in open ocean usually are long, so perturbation method can be applied to reduce origin system of equations to simpler problem. In section II the nonlinear evolution equation that describe water level perturbation is derived. The expression for the total energy of wave in terms of water level perturbation written in non-dimensional variables is obtained. Section III contain description of numerical algorithm for solving boundary value problem for derived equation. In section IV boundary value problem is solved numerically. The behaviour of solitary wave is studied, time dependence of the wave energy is calculated.

Equation derivation
Let us consider the potential incompressible flow of ideal fluid under thin elastic plate over flat bottom. This motion is described by the following equations [4,5] Here ϕ is the velocity potential of the liquid, h 0 is depth of liquid, η is the deviation of liquid from its equilibrium position (h(x, t) = h 0 + η(x, t)), ρ i and h i are the density and thickness of the ice, g is the gravitation acceleration, P 0 is the atmospheric pressure, ρ w is the density of water, P is the external load, E is the elastic modulus and σ xx = const is the component of stress tensor. We introduce the dissipative term −κη xxx to the right part of the equation of plate motion, so κ is the phenomenological coefficient describing dissipation of wave energy under ice cover due to bend. After changing of variables in (2) where l is a wavelength, a is a characteristic amplitude and c 0 = √ gh 0 is a phase velocity in the linear approximation, we obtain the following problem Here ε, µ, δ, β, τ, and γ are small parameters of the same orders defined by the following expressions Here we suppose that all parameters have the same order. Expanding ϕ(x, y, t) in a Taylor series in powers of y we obtain, from the first equation and first boundary condition of problem (4) Substitution of (6) into the second and third boundary conditions of (4) yields system of nonlinear partial differential equation, which with an accuracy up to terms of order O(ε 2 ) looks where w( In the zero order approximation system (7) takes a form and for compatibility of this system we have to assume Applying perturbation method in the form given in [5] to (7) we obtain differential equation for describing nonlinear waves of liquid under continuous ice layer Values of constant coefficients c k are defined through small parameters (5) by following expressions

Small parameters of equation
Values of coefficients (10) were analyzed for typical parameters of model: the length of wave is from 100 m to 1 km, the wave amplitude is of the order of 1 m, water depth is from 10m to 50m, the thickness of ice is of the order of 1 m. Direct calculation of all coefficients in this range shows that the biggest coefficient is c 01 . The another coefficients in equation is usually two to three orders less than c 01 .The coefficients c 3 , c 03 , c 12 , c 001 can be of the order 10 −3 -10 −4 . The remaining coefficients has the lower order. The order of coefficient c 4 depend on value of parameter κ. We want to take into account dissipative effects in the model, so we keep it in equation.
Waves in ocean can penetrate several hundreds of kilometers. It corresponds to time 10 3 -10 4 in nondimensional variables for wave length from 100 m to 1 km. So we can omit all coefficients in equation having the order 10 −5 and less: c 5 , c 6 , c 7 , c 8 , c 9 , c 04 , c 22 , c 13 , c 05 , c 14 , c 23 .

Energy of long wave
The total energy of long wave consist of kinetic and potential energy. Elastic energy of ice layer we do not take into account. We suppose that potential energy of nonperturbed system is zero. The energy of unit-width wave can be calculated by the formula (potential energy can be negative for the wave with negative amplitude)

In nondimentional variables (3) energy has the form
Here E 0 is potential energy of unit-width unperturbed water reservoir measured with respect to bottom y = 0. Note, that ϕ 2 x = O(1), while ϕ 2 y = O(µ 2 ). So the kinetic energy related to vertical motion has the higher order of smallness than the energy related to horizontal motion.
Taking into account expansion (6) for potential and the expansion (8) expression for the energy (11) can be written in terms of η(x, t) Therefore the first correction to the total energy of nonperturbed system has the order O(ε 2 ).

Numerical modeling
As was said above, we can omit following coefficients in Eq. Here L is big enough for wave was localized in domain [0; L] for all t < T . Finite difference method is used to solve problem (13), (14). Eq.(13) can be written in the form For numerical modeling we use implicit finite difference scheme, which leads to system HereL is finite difference operator. All differential operators in L are approximated with central differences to a second order of accuracy. The truncation error of this scheme is O(τ 2 ) + O(h 2 ). Numerical experiments shows that this scheme is unconditionally stable. A couple of numerical algorithms for solving KdV equations is presented in work [6]. The suggested methods can be generalized for solving derived is Section II equation. We have a system of nonlinear algebraic equations on each time layer, which is solved by Newton's method. On the each iteration of Newton's method solution is computed by formula Exact form of Jacobin matrix J can be found in Appendix. For the first iteration of algorithm the solution from the previous time layer is used L ∞ -norm is used to control convergence of Newton's method So, a five-tridiagonal system of equations is required to be solved at each step of Newton's method. We use DGBSV solver from LAPACK library to solve this system of equations.

Results
To verify the proposed mathematical model several cases of problem was numerically solved. In table 4 the values of model parameters and corresponding values of coeffcieints in Eq.(13) are presented. The following initial condition was used in numerical simulation .
Here x 0 is some constant. In the first two cases value x 0 = 10 was used. In the third case  (14) is presented for all cases from table 4. As it can be seen, the wave shape is various for all cases and heavily  depend on the parameters of the model. The attenuation of the wave energy was numerically calculated with taking into account expression (12). Results presented on the figure 4. The shape of plot for the energy is similar to results, obtained by Wadham [1]. The value of parameter τ was fixed in all three cases. The additional tests for other values of τ shows that the shape of energy curve does not crucially depend on it exact value. Additional investigation is required to explore the dependence of energy attenuation on it value. However obtained curve of energy against time can not be properly fitted by the expression (1), which has good agreement with field data [1]. There some reasons for this. The proposed model is phenomenological and require additional argumentation. Furthermore proposed model does not take into account some processes. As it can be seen on figure 4, the curve of the solution in some points is very bended. In practice, this leads to the ice break up and subsequent energy release. The model must be changed to take into account such processes. The bottom surface of the ice layer is not plain for real ice floe and turbulent mixing must be considered.