Three-dimensional time-domain scattering of waves in the marginal ice zone

Three-dimensional scattering of ocean surface waves in the marginal ice zone (MIZ) is determined in the time domain. The solution is found using spectral analysis of the linear operator for the Boltzmann equation. The method to calculate the scattering kernel that arises in the Boltzmann model from the single-floe solution is also presented along with new identities for the far-field scattering, which can be used to validate the single-floe solution. The spectrum of the operator is computed, and it is shown to have a universal structure under a special non-dimensionalization. This universal structure implies that under a scaling wave scattering in the MIZ has similar properties for a large range of ice types and wave periods. A scattering theory solution using fast Fourier transforms is given to find the solution for directional incident wave packets. A numerical solution method is also given using the split-step method and this is used to validate the spectral solution. Numerical calculations of the evolution of a typical wave field are presented. This article is part of the theme issue ‘Modelling of sea-ice phenomena’.

Three-dimensional scattering of ocean surface waves in the marginal ice zone (MIZ) is determined in the time domain. The solution is found using spectral analysis of the linear operator for the Boltzmann equation. The method to calculate the scattering kernel that arises in the Boltzmann model from the single-floe solution is also presented along with new identities for the far-field scattering, which can be used to validate the single-floe solution. The spectrum of the operator is computed, and it is shown to have a universal structure under a special nondimensionalization. This universal structure implies that under a scaling wave scattering in the MIZ has similar properties for a large range of ice types and wave periods. A scattering theory solution using fast Fourier transforms is given to find the solution for directional incident wave packets. A numerical solution method is also given using the split-step method and this is used to validate the spectral solution. Numerical calculations of the evolution of a typical wave field are presented.
This article is part of the theme issue 'Modelling of sea-ice phenomena'.

Introduction
Vast areas of the polar ocean surface are frozen into a thin layer (centimetres to metres thick) known as sea ice. Sea ice plays an important role in the global climate, affecting heat fluxes between the atmosphere and ocean,

Problem formulation
The wave energy balance equation is ∂ t N(x, t, k, θ ) + c g θ · ∇N(x, t, k, θ ) = −c g α s (x, t, k, θ)N(x, t, k, θ) + c g 2π 0 S k (x, t, k, θ, θ )N(x, t, k, θ ) dθ , (2.1a) which is also known as the linear Boltzmann equation. Here, the Cartesian coordinate x = (x, y) denotes location on the ocean surface, ∇ = (∂ x , ∂ y ) is the gradient operator, t is time, N is the so-called wave action density, k is the wavenumber (modulus of the wave vector), θ is the wave direction vector (argument of the wavenumbers in the x-and y-directions), c g is the group velocity, α s is the energy loss and S k is the scattering kernel that redistributes wave energy in different directions. As scattering is an energy-conserving process, the redistribution of energy is exactly the loss of energy in the given wave direction, so that Non-conservative effects could easily be included by modifying the expression for α s . We assume that the scattering is isotropic, in the sense that S(x, t, k, θ, θ ) ≡ S(x, t, k, θ − θ ), anticipating this condition will hold approximately and on average in a typical MIZ over sufficient spatial and/or temporal scales, as ice floes are orientated randomly. Equations (2.1) become respectively. In this study, we achieve isotropic scattering by using modelling the ice floes as being circular. Non-circular models (e.g. [25,26,31]) could be used, but only with the additional burden of averaging the scattering characteristics.

Calculation of the scattering kernel (a) Solution for a circular floe
We briefly summarize the eigenfunction matching solution method proposed by Peter et al. [27] for a circular elastic plate model of an ice floe, floating on a three-dimensional water domain, which has finite depth and stretches to infinity in both horizontal directions. The plate has radius a and uniform thickness h, although varying thickness could be incorporated using the method of Bennetts et al. [32]. Assuming that all motions are time harmonic with radian frequency ω, the velocity potential of the water,φ, can be expressed as is the Laplacian operator in the horizontal plane, k is the open water wavenumber (defined below), φ I is the incident potential, and α = ω 2 /g is a frequency parameter, in which g ≈ 9.81 m s −2 is the constant of gravitational acceleration. The parameters β and γ are where ρ ≈ 1025 kg m −3 is the water density, Y ≈ 6 GPa is Young's modulus of the plate, ν ≈ 0.3 is its Poisson's ratio and ρ i ≈ 992.5 kg m −3 is its density. The incident potential is set to be a plane wave with displacement amplitude A, travelling in the positive x-direction, i.e.
The potential is expanded in its eigenfunctions as where I n and K n are modified Bessel functions of order n. The wavenumbers k m and κ m are the rootsk andκ of the dispersion equationŝ respectively [33]. The first dispersion relation has a unique purely imaginary solution in the upper half of the complex plane, denoted k 0 = ik, and an infinite number of positive real solutions denoted k 1 < k 2 < · · · . The second dispersion relation also has a unique purely imaginary solution in the upper half of the complex plane, κ 0 , and an infinite number of positive real solutions, κ 1 < κ 2 < · · · . It also has two complex roots in the upper-half complex plane, denoted κ −2 and κ −1 , where κ −1 =κ −2 . We define associated vertical eigenfunctions as for the open water region, and As the problem is axisymmetric, we can solve for the amplitudes a mn and b mn of each angular mode n separately, and we also note that the amplitudes for positive and negative n are complex conjugates so that solving for n ≥ 0 is sufficient. For a given value of n, the amplitudes are obtained in the standard way, by matching the potential and its radial derivative at r = a for −H < z < 0, taking inner products with respect to φ m (m = 0, 1, . . .), and truncating the resulting system of equations and all infinite sums to m ≤ M. The system is closed by applying the plate free edge conditions, which are the antepenultimate and penultimate conditions of equation (3.2).

(b) Scattering kernel from the single-floe solution
Meylan et al. [18] gave the following justification for calculating the scattering kernel from the scattering produced by a single floe (based on the work of [29]). Each floe scatters energy, and the energy radiated per unit angle per unit time in the θ direction for a wave incident in the θ direction, E, is given by where D(θ − θ ) is the scattered amplitude. At a large distance from the scatterer, the asymptotic amplitude of the outgoing wave in the θ direction, for an incident wave travelling in the θ direction, is given by The scattering kernel, S k , is found by dividing E by the rate of energy passing under the ice floe. The rate of energy passing under the floe is given by the product of the wave energy density ( 1 2 ρgA 2 ), the effective area occupied by a floe (A f /f i , where A f is the area of the floe and f i is the fraction of the surface area of the ocean which is covered in ice), and the wave group speed (ω/2k). Thus, Using the asymptotic results for Bessel functions with large arguments, the far-field wave amplitude is e n cos(nθ), (3.14) where we have used the property that the solutions for positive and negative n are complex conjugates. For the purposes of calculations in wave prediction code, values of D can be efficiently stored using only a small number of e n values.

(c) Identities
Here we derive identities for the solution of the circular ice floe problem, which are (to our knowledge) not known in the literature, and help validate our numerical results. The identities apply to the far-field calculation for anybody with circular symmetry. They follow from noticing that if the incident wave was a Hankel function of the second kind, then the outgoing wave at infinity would have to be a Hankel function of the first kind with the same amplitude (possibly with a different phase). As the incident wave is written as i n J n (Kr), it follows that   s. This figure shows that the scattering is highly dependent on wave period, noting that from equation (3.13) the square of |D(θ)| controls the scattering of wave energy. As the floe radius increases, the scattering becomes more complicated with respect to angle, but given the variation in floe size in an MIZ this effect will be removed by averaging. The scattering is strongest in the θ = 0 direction, which is related to the so-called shadow zone behind the floe, where the scattered wave almost cancels the incident wave. thickness affect the scattering for small thicknesses, but, above a critical thickness, scattering becomes insensitive to thickness changes. This insensitivity is because the floe bending is negligible above the critical thickness. The critical thickness is period dependent and increases as the period increases. The behaviour is also illustrated in figure 3, where we vary the thickness rather than the period in each subplot.
It is informative to plot α s , as it determines the strength of the scattering. Note that, if all the scattered wave energy was removed from the system, then the wave field would decay as e −α s x . Figure 4 shows α s as a function of wave period for various thicknesses and floe radii. The plot is on a log scale in α s , and shows that the scattering depends strongly on period but only weakly on the other parameters. The sharp dips indicate resonances and would be removed by averaging with respect to the floe properties. The figure shows that the wave period has the strongest control over scattering.

Temporal solution of the half-space scattering problem
Consider the half-space scattering problem for the wave energy balance equation (2.2), in which S k = 0 for x < 0, and S k = S k (t, k, θ − θ ) = 0 for x < 0, i.e. open water exists for x < 0 and a spatially homogeneous MIZ exists for x > 0. Further, we assume that the solution is independent of the y-coordinate.

(a) Transformation to a discrete problem
We approximate the integral operator appearing in the wave energy balance equation (2.2) using 2n angular values where n is odd to ensure that we do not have θ i = ±π/2, which would cause singularities. (A modified solution method for this case is given by Meylan [30].) Equation (2.2a) becomes   is the discrete intensity, χ Γ (x) is the characteristic function, which is zero for x / ∈ Γ and unity for x ∈ Γ , and D is a diagonal matrix, which is self-adjoint but not positive, and S is a positive self-adjoint matrix. The entries of D and S are, respectively, where w i is the weight associated with the chosen quadrature rule (we use the trapezium rule so Symmetry meansα s is independent of i andα s is a numerical approximation of α s consistent with our discretization, ensuring energy is conserved. Further, the matrix S is a circulant matrix.

(b) Non-dimensionalization
Non-dimensionalizing the problem makes the numerical scheme more stable, and, as will be shown, the time and length scales we use control the gross behaviour of the solution. We scale the temporal coordinate with respect to 1/s 11 . For a large number of angles (n 1), we note that s 11 ≈ c gαs ≈ c g α s , which is the exponential decay rate in time in the case that all scattered wave For a large number of discrete angles, x ≈α −1 s ≈ α −1 s , which is the length scale of the exponential decay in space if all the scattered waves are removed. Under this non-dimensionalization, equation (4.1) becomes where the entries of D and S are Thus, the entries of the matrices vary between −1 and 1, and we have s ii = 1. This makes the numerical schemes stable and the same discretization can be used in space and in time. We will drop the star from now on.

(c) Solution using scattering theory
The solution of equation (4.5) can be written in abstract form as is an operator that defines a continuous semi-group, for which Meylan [ how important this is for controlling wave scattering. We note that if μ ∈ spec {L}, then e iμt ∈ spec {exp(iLt)}. We will show how to compute this semi-group using a modified scattering theory, for the case of an incident wave group, which is initially zero in the MIZ (x > 0). The spectrum of L is defined by Separating the spectrum in the two regions x > 0 and x < 0, we obtain where w − j are the eigenvectors of the matrix −iD −1 λ and ζ − j are the corresponding eigenvalues, and w + j are the eigenvectors of the matrix −iD −1 (−iS + λ) and ζ + j are the corresponding eigenvalues. The constants c ± j are at this stage arbitrary. At first sight, it appears that there will be a solution of (4.8) for every value of λ. However, we need to impose the condition that the solution is bounded at x = ±∞. From symmetry, in both cases the eigenvalues come in plus and minus pairs, meaning we can eliminate half the eigenvectors and this leaves us only 2n unknowns. Matching at x = 0, we obtain 2n homogeneous equations, for which there will only exist a trivial solution. However, if either ζ − j or ζ + j is purely imaginary, we will have one extra unknown (as we can keep both the plus and minus solutions). As the number of unknowns is one greater than the number of homogeneous equations, we will have a non-trivial solution. This non-trivial solution is determined up to a constant and this will lead to the spectrum of the operator L. The corresponding functions for the eigenvalues of L are known as generalized eigenfunctions.
For λ ∈ R, the solution for x < 0 gives ζ − j ∈ i R. We denote the corresponding generalized eigenfunctions (with multiplicity n) by where the summation for x > 0 is over eigenvalues with negative real part only, i.e. Re(ζ + j ) < 0. The generalized eigenfunctions v − i are chosen to be non-zero only in a single incoming direction, which makes subsequent calculations easier. The n equations obtained by matching the solution for the first n rows at x = 0 uniquely determine the coefficients c + j . Subsequently, by matching the second n rows at x = 0 the elements of the matrix T ij are found. where I is the identity matrix, for a given μ ∈ R. Equation (4.11) can be solved by finding the eigenvalues of the matrix −μD + iS for μ ∈ R. Therefore, the complex spectrum consists of n lines in the complex plane, which lie in the upper-half plane so that the semi-group is a contraction. The generalized eigenfunctions for the complex spectrum are zero in the incident directions, so that energy is allowed to leave the MIZ only. The eigenfunctions are where, as before, the summation for x > 0 is over eigenvalues with negative real part only, i.e. Re(ζ + j ) < 0. As the ordering of the eigenvectors for x > 0 is arbitrary, we choose, without loss of generality, that eigenvectors corresponding to iμ and −iμ have numbering n and n + 1, respectively. Further, we choose the values of c + j so that the solution is zero in the incoming directions. The elements of the matrix T ij are then found by matching the second n rows at x = 0. Figures 5 and 6 show the branches of the spectrum of the operator L for different values of the angle truncation parameter n. Figure 5 shows the spectrum for a floe of thickness h = 1 m and radius a = 50 m, and wave period of T = 8 s. Figure 6 shows the spectrum for a floe of thickness h = 1 m and radius a = 200 m, and wave period of T = 12 s. These figures are the only plots we present in non-dimensional variables and the reason for this is the following. The imaginary part of this spectrum represents the decay which is proportional to the diagonal element of S since decay in this context is spreading. As the non-dimensional diagonal is always unity, we see that the complex spectrum consists of a vertical branch and lines which all have imaginary part ≈ 1. This pattern of the spectrum implies that the parameters t and x we calculate control the properties of the solution in time and space, respectively. It is interesting to ask what happens in the limit as the number of angles tends to infinity. This limit should allow us to compute the spectrum of the continuous operator which appears in equation (2.2a). Perhaps the spectrum collapses to a single line, or perhaps it fills a region of the complex plane.
We construct the solution to the time-dependent half-space problem using a version of the method developed by Meylan [30], which is based on the generalized eigenfunctions v − i . For an incident wave which is initially zero in the MIZ, it is possible to use only the real part of the spectrum of L and only v − i . Physically, these are the only cases which could be excited, but mathematically we could start with the initial condition non-zero in the MIZ. We assume the initial wave action density in the MIZ is zero, i.e. N 0 (x) = 0 for x > 0, and that it is of the form . . .   (4.13) and, using the properties of the Fourier transform, we deduce that (4.14) Note that the functionsf − j are not exactly the Fourier transform of f j because of the D jj term. Therefore, the solution is The integrals in equations (4.14) and (4.15) are calculated using the fast Fourier transform.

(d) Validation via numerical method
To validate the above spectral method based on scattering theory, we solve equation (4.5) using a so-called split-step method, described in appendix A. Figure 7 shows the evolution of the wave action density calculated using scattering theory and the split-step method. The incident wave field is uni-directional, travelling in the direction of the  are 2n = 6 angles, which is unlikely to be sufficient for practical calculations but is enough to validate the scattering theory. The initial wave profile is the Gaussian The agreement is generally excellent, with only small deviations at large distances, as the fast Fourier transform component of the spectral method is periodic, resulting in a residual solution, evident in figure 7b,c.

Numerical results
Numerical results are given for floe radius a = 50 m and thickness h = 1 m, an incident wave field with period T = 8 s, and shape    . This plot shows that the solution is well converged for 2n = 42, noting that this is more angles than are normally used in global wave prediction codes. On this basis, figure 9 shows the evolution of polar plots for 2n = 42. It is apparent that the amplitude of the scattered field rapidly decreases with distance into the MIZ and that the wave field becomes isotropic. Figure 10 shows corresponding three-dimensional mesh plots of the wave field as functions of wave direction.

Conclusion
We have investigated a Boltzmann equation model for wave scattering in the MIZ. We have shown how to derive the scattering kernel appearing in the Boltzmann equation from the solution for a single circular ice floe and presented new identities to validate it. We developed a mathematical theory to analyse the structure for the half-space problem under discretization in angle. We calculated the spectrum of the linear operator and showed that under a particular nondimensionalization it has universal properties. We showed that a solution can be found using a scattering theory for incident waves, and validated the spectral solution against a split-step method. The method allows computation of the spatial and temporal evolution of wave packets propagating through the MIZ.