On the attenuation of waves through broken ice of randomly-varying thickness on water of finite depth

TL;DR

This paper extends previous models of wave attenuation through broken ice to account for water of finite depth, providing explicit formulas and validating them with numerical simulations, revealing frequency-dependent attenuation behaviors.

Contribution

It introduces a theoretical model for wave attenuation in non-shallow water with randomly-varying ice thickness, validated by numerical simulations and applicable to real-world scenarios.

Findings

Attenuation proportional to frequency to the eighth power at low frequencies

Good agreement between theoretical predictions and numerical simulations

Identification of a roll-over effect at higher frequencies

Abstract

The recent work of Dafydd and Porter [2024] on the attenuation of waves propagating through floating broken ice of random thickness is extended to consider water of non-shallow depth. A theoretical model of broken floating ice is analysed using a multiple scales analysis to provide an explicit expression for the attenuation of waves as they propagate from a region of constant thickness ice into a semi-infinite region of ice whose thickness is a slowly-varying random function of distance. Theoretical predictions are shown to compare well to numerical simulations of scattering over long finite regions of ice of randomly-varying thickness computed from an approximate depth-averaged model derived under a mild-slope assumption. The theory predicts a low-frequency attenuation proportional to the eighth power of frequency and a roll-over effect at higher frequencies. The relationship between the results and field measurements are discussed.

Equations

$$\nabla^{2}\phi=0\qquad\mbox{on $-h<z<-d_{pc}(x)$}$$

$$\phi_{z}=0\qquad\mbox{on $z=-h$}.$$

$$\phi_{z}=-\mathrm{i}\omega\zeta_{i}\qquad\mbox{on $z=-d_{i}$, $x_{i}<x<x_{i+1}$},$$

$$-\omega^{2}\rho\ell_{i}d_{i}\zeta_{i}=-\rho g\ell_{i}\zeta_{i}+\mathrm{i}\omega\rho\int_{x_{i}}^{x_{i+1}}\phi(x,-d_{i})\;\mathrm{d}x\quad\mbox{on $x_{i}<x<x_{i+1}$}$$

$${\varphi }_x=0,x=x_i,-\text{mbox}max(d_{i-1},d_i)<z<-\text{mbox}min(d_{i-1},d_i).$$

$$x^{\prime }=Kx,z^{\prime }=Kz,\varphi (x,z)=-(\mathrm{i}Ag/\omega ){\varphi }^{\prime }(x^{\prime },z^{\prime }),{\zeta }_i(x)=A{\zeta }_i^{\prime }(x^{\prime }),$$

$$\nabla^{2}\phi=0\qquad\mbox{in $-h<z<-d_{pc}(x)$}$$

$$\phi_{z}=0,\qquad\mbox{on $z=-h$}$$

$$\phi_{z}=\zeta_{i},\qquad\mbox{on $z=-d_{i}$, $x_{i}<x<x_{i+1}$}$$

$$(1-d_{i})\zeta_{i}=\frac{1}{l_{i}}\int_{x_{i}}^{x_{i+1}}\phi(x,-d_{i})\,dx,\qquad\qquad\mbox{on $x_{i}<x<x_{i+1}$}.$$

$$l_i{\zeta }_i={\int }_{x_i}^{x_{i+1}}{\varphi }_z(x,-d_i)dx={\int }_{x_i}^{x_{i+1}}{\varphi }_n(x,-d(x))\sqrt{1+[d^{\prime }(x){]}^2}dx$$

$${\zeta }_i=\frac{1}{l_i}{\int }_{x_i}^{x_{i+1}}(d^{\prime }(x){\varphi }_x(x,-d(x))+{\varphi }_z(x,-d(x)))dx\approx d^{\prime }({\bar{x}}_i){\varphi }_x({\bar{x}}_i,-d({\bar{x}}_i))+{\varphi }_z({\bar{x}}_i,-d({\bar{x}}_i))$$

$$(1-d(x)){\varphi }_n(x,-d(x))=\varphi (x,-d(x)).$$

$$\nabla^{2}\phi=0\qquad\mbox{on $-h<z<-d(x)$}$$

$$\phi(x,z)\sim\left\{\begin{array}[]{l}(\mathrm{e}^{\mathrm{i}k_{0}x}+R_{\infty}\mathrm{e}^{-\mathrm{i}k_{0}x})Z_{0}(z),\qquad x\to-\infty,\\ 0,\qquad x\to\infty\end{array}\right\}$$

$$Z_0(z)=N_0^{-1/2}\cosh k_0(z+h),\text{mbox}where{N}_0=\frac{1}{2}\left(1+\frac{\sinh 2k_0(h-d_0)}{2k_0(h-d_0)}\right)$$

$$k_0\tanh k_0(h-d_0)=\kappa \equiv 1/(1-d_0).$$

$$\phi(x,z)\sim\left\{\begin{array}[]{l}(\mathrm{e}^{\mathrm{i}k_{0}x}+R_{L}\mathrm{e}^{-\mathrm{i}k_{0}x})Z_{0}(z),\qquad x\to-\infty,\\ T_{L}\mathrm{e}^{\mathrm{i}k_{0}x},\qquad x\to\infty\end{array}\right\}$$

$$d(x)=d_0(1+\sigma r(x))$$

$$⟨r(x)⟩=0,⟨r^2(x)⟩=1,$$

$$⟨r(x)r(\check{x})⟩={\mathrm{e}}^{-(x-\check{x}{)}^2/{\Lambda }^2}.$$

$$\varphi (x,z)={\varphi }_0(x,X,z)+\sigma {\varphi }_1(x,X,z)+{\sigma }^2{\varphi }_2(x,X,z)+\dots$$

$$\left({\nabla }^2+2{\sigma }^2\frac{{\partial }^2}{\partial x\partial X}+{\sigma }^4\frac{{\partial }^2}{\partial X^2}\right)\left({\varphi }_0+\sigma {\varphi }_1+{\sigma }^2{\varphi }_2+\dots \right)=0.$$

$$\frac{\partial }{\partial z}\left({\varphi }_0+\sigma {\varphi }_1+{\sigma }^2{\varphi }_2+\dots \right)=0,\text{on}z=-h$$

$$(1-d_0-d_0\sigma r(x))\left(\frac{\partial }{\partial z}-\sigma d_{0}r(x)\frac{{\partial }^2}{\partial z^2}+{\sigma }^2\frac{d_0^2{r}^2(x)}{2}\frac{{\partial }^3}{\partial z^3}+\sigma d_0{r}^{\prime }(x)\frac{\partial }{\partial x}-{\sigma }^2{d}_0^{2}r(x)r^{\prime }(x)\frac{{\partial }^2}{\partial x\partial z}+\dots \right)$$

$$\left({\varphi }_0+\sigma {\varphi }_1+{\sigma }^2{\varphi }_2+\dots \right)$$

$$=\left(1+{\sigma }^2\frac{d_0^2[r^{\prime }(x){]}^2}{2}+\dots \right)\left(1-\sigma d_{0}r(x)\frac{\partial }{\partial z}+{\sigma }^2{d}_0^2{r}^2(x)\frac{{\partial }^2}{\partial z^2}+\dots \right)\left({\varphi }_0+\sigma {\varphi }_1+{\sigma }^2{\varphi }_2+\dots \right)$$

$${\varphi }_0(x,0,z)\sim ({\mathrm{e}}^{\mathrm{i}{k}_{0}x}+R_{\infty }{\mathrm{e}}^{-\mathrm{i}{k}_{0}x})Z_0(z),\frac{\partial }{\partial x}{\varphi }_0(x,0,z)\sim \mathrm{i}{k}_0({\mathrm{e}}^{\mathrm{i}{k}_{0}x}-R_{\infty }{\mathrm{e}}^{-\mathrm{i}{k}_{0}x})Z_0(z),x\to -\infty$$

$$\nabla^{2}\phi_{0}=0,\qquad\mbox{in $-h<z<-d_{0}$}$$

$$\frac{\partial\phi_{0}}{\partial z}=0,\qquad\mbox{on $z=-h$}$$

$$\frac{\partial\phi_{0}}{\partial z}=\kappa\phi_{0},\qquad\mbox{on $z=-d_{0}$}.$$

Full text

Key words: sea ice, finite depth flows, wave scattering, ice floes, wave energy attenuation

1 Introduction

In a recent paper, Dafydd and Porter (2024) considered how waves attenuate as they pass through regions of broken ice floating on water. Their study was conducted under a shallow water assumption allowing a basic model to be developed which took the form of a second order ODE whose coefficients encoded the variation of the ice thickness. Subsequently, a multiple scales analysis, inspired by the approach adopted by Grataloup and Mei (2003), Mei and Li (2004) and Mei et al. (2005), but incorporating significant adjustments to the methodology, was used to predict the decay of wave energy through a region of randomly-varying ice thickness as a function of distance. Theoretical predictions were shown to compare well with numerical simulations.

The model of Dafydd and Porter (2024) is energy conserving with no source of physical damping incorporated into the dynamics. The mechanism responsible for the decay of wave energy propagating through a random environment results from multiple scattering and is related to the phenomenon of localisation, an effect first described by Anderson (1958) in quantum systems.

The main purpose of the work of Dafydd and Porter (2024) and of the present paper is to consider whether randomness and multiple scattering is a candidate for describing the data from field measurements of attenuation of wave energy collected in polar regions. This study adds to a number of similar studies which share the same purpose in which randomness is introduced in a number of different ways. In Kohout and Meylan (2008), Squire et al. (2009) and Mokus (2023) two-dimensional scattering is considered through floes whose length is the order of the wavelength and where randomness is introduced through the length of the floes. In Bennetts et al. (2010), Montiel et al. (2016), and Meylan et al. (2021) three-dimensional studies were conducted into random configurations of randomly-sized circular floes. Many researchers have also considered different mechanisms by which attenuation occurs as a result of physical damping such as viscoelastic effects (Wang and Shen (2010), Keller (1998)), roughness, and a range of phenomological effects (as in Vaughan et al. (2009)). None of these studies have been able to provide a definitive account of the observed data.

A number of field studies over many decades including Wadhams et al. (1988), Liu et al. (1992), Kohout et al. (2014), Doble et al. (2015), Kohout et al. (2020), Montiel et al. (2022) have sought to establish an increasingly refined description of the attenuation with a particular emphasis on characterising its dependence on wave frequency, $\omega$. For low to mid frequencies the data is suggestive of a power-law dependence of $\omega^{n}$ where $n$ lies between 2 and 4 (see, for example, Meylan et al. (2014), Meylan et al. (2018), Montiel et al. (2022)). However the evidence for extending the power-law fit to very low frequencies appears much less clear and Meylan et al. (2021) have suggested that models with $n=8-10$ may capture the behaviour. In fact, the work of Herman (2024) suggests that the $2<n<4$ power-law dependence may be a flawed approach to evaluating models that address ice-related attenuation and in fact models with a relatively high $n$ should not necessarily be disregarded. For higher frequencies most data sets available indicate a roll-over effect (see Li et al. (2017) for a comprehensive discussion) although a more recent study by Thomson et al. (2021) provides good evidence to suggest roll over is a statistical effect due to the methods used for sampling noisy data. The analysis of Rogers et al. (2021) and Montiel et al. (2022) uses this work to filter out these effects.

In Dafydd and Porter (2024) the basic shallow water model, theory and numerical simulations resulted in an attenuation coefficient which behaved as $\omega^{2}$ for low frequencies with a roll-over effect – a peak in attenuation followed by exponential decay – at higher frequencies. Dafydd and Porter (2024) assumed floes of small size and introduced randomness in a different way to previous authors via the ice thickness. A number of simplifying modelling assumptions were introduced in order to reduce the complexity of a fluid interacting with broken ice to a 2nd order ordinary differential equation which could be subjected to analysis. Since the shallow water assumption was the most unrealistic modelling assumption, here we have extended that work and report on the case of water of finite depth albeit still working in a two-dimensional setting.

In this paper we assume an ideal fluid undergoing small amplitude motion and hence described by a potential satisfying Laplace’s equation. Nevertheless, the multiple-scales analysis applies in a similar way to Dafydd and Porter (2024), with our methodology following Mei et al. (2005) and Bennetts et al. (2015) who investigated wave localisation over variable bathymetry in the presence of an unloaded surface. Significantly we modify the approach to incorporate the steps made in Dafydd and Porter (2024) which corrects an overprediction of the attenuation. Numerical simulations of Bennetts et al. (2015) highlighted the differences between the decay from averaging individual realisations of wave scattering over variable bathymetry with the theoretical predictions based on the method of Mei et al. (2005) and attributed this to coherent phase cancellation in the averaging process in the multiple scales analysis. In Dafydd and Porter (2024) two modifications were made to address this: (i) a scattering problem was formulated from the outset with incident waves entering a semi-infinite region of randomness; (ii) in the averaging process, terms that could be identified as contributing to decay via coherent phase cancellation rather than by multiple scattering were removed from the computation of the attenuation coefficient. The resulting numerical simulations showed a similar character to those of Bennetts et al. (2015) and were represented accurately by the theory.

In Section 2 we formulate a continuum model for floating broken ice based on the same assumptions employed in Dafydd and Porter (2024). Thus ice completely covers the surface of the fluid and is constrained to move vertically and independently of its neighbours. The floes are sufficiently narrow in width for the governing equations and boundary conditions to be transformed into a continuous description of the ice, equivalent to the mass loading model used by, for example, Keller and Weitz (1953), Mosig et al. (2017) but modified for variable ice thickness, equivalent to the model derived by Porter and Porter (2004) without bending forces. In Section 3, a multiple-scales analysis is employed to the two-dimensional boundary-value problem the result of which is an explicit expression for the attenuation coefficient. In Section 4 we derive a mild-slope equation for broken ice (MSEBI) on the assumption that the ice thickness is a slowly-varying function of space starting from a variational principle for the governing equations. Results of numerical solutions to the MSEBI are compared with the theory in Section 5. Section 5 also includes a discussion regarding the results of the model and how it might relate to to data from field measurements. Further potential work and limitations will also be presented.

2 Modelling floating broken ice

Cartesian coordinates, $(x,z)$, are used where $z$ is directed upwards from an origin located in the horizontal rest position of the water surface in the absence of an ice cover. Fluid of density $\rho$ is bounded below by a flat rigid bed at $z=-h$ and above by a continuous ice cover, broken into rectangular blocks/floes by narrow vertical cracks along $x=x_{i}$. The depth of submergence of the floe lying between $x_{i}<x<x_{i+1}$ is $d_{i}$ and its width, $l_{i}=x_{i+1}-x_{i}$, is assumed small compared to the wavelength of waves that propagate across the ice-covered water surface. On account of the Archimedes principle the thickness of the ice floe in $x_{i}<x<x_{i+1}$ is given by $d_{i}\rho/\rho_{ice}$, where $\rho_{ice}$ is the density of ice.

Assuming inviscid incompressible irrotational flow, the velocity of the fluid is given by the gradient of a potential $\Phi(x,z,t)$. Small amplitude motion is assumed implying that the governing equations and boundary conditions are linearised. We factorise a harmonic time dependence of angular frequency $\omega$ from the motion such that $\Phi(x,z,t)=\mathrm{Re}\left(\phi(x,z)\mathrm{e}^{-\mathrm{i}\omega t}\right)$ where $\phi(x,z)$ is a complex potential satisfying

[TABLE]

in the fluid where $d_{pc}(x)=d_{i}$ for $x_{i}<x<x_{i+1}$ is a piecewise constant function. On a flat bed the impermeability condition is expressed by

[TABLE]

At the fluid-ice interface the kinematic boundary condition is expressed by

[TABLE]

where $\mathrm{Re}\{\zeta_{i}\mathrm{e}^{-\mathrm{i}\omega t}\}$ is the heave amplitude of the ice block with left- and right-hand edges along $x=x_{i}$ and $x=x_{i+1}$. This function is also coupled to $\phi$ through a dynamic condition for the constrained heave motion of each floe by

[TABLE]

where $g$ is gravity.

Additionally there is no flow across the vertical faces of the steps in the ice in contact with the fluid and this is expressed by

[TABLE]

We non-dimensionalise these governing equations by writing

[TABLE]

where $A$ is some characteristic amplitude of motion of the ice and $K=\omega^{2}/g$. This results in (after dropping primes)

[TABLE]

with

[TABLE]

and

[TABLE]

with

[TABLE]

and $h$, $x_{i}$ and $d_{i}$ are all now dimensionless variables and $l_{i}\ll 1$ for each $i$. The condition (2.8) remains unchanged after non-dimensionalisation.

The next step is to replace the discrete description of the ice with a continuum model. In $x_{i}<x<x_{i+1}$ we write

[TABLE]

using the divergence theorem in the small region $-d(x)<z<-d_{pc}(x)$ between the discrete description of the ice cover and a continuous version formed by $z=-d(x)$ where $d(x)$ is defined such that $d(x_{i})=\max\{d_{i},d_{i-1}\}$ and $d^{\prime}(x)$ exists everywhere (practically, we define $d(x)$ and infer $d_{pc}(x)$). In deriving (2.11) we have used (2.5), (2.9) and (2.1) in the fluid region $-d(x)<z<-d_{pc}(x)$; the operator $\partial_{n}=(1+[d^{\prime}(x)]^{2})^{-1/2}(d^{\prime}(x)\partial_{x}+\partial_{z})$ is the derivative normal to the curve $z=-d(x)$ directed out of the fluid. Thus

[TABLE]

for $\bar{x}_{i}=\frac{1}{2}(x_{i}+x_{i+1})$ with error $O(l_{i}^{2})$. Combining with (2.10) and making approximations consistent with those made above results in

[TABLE]

The governing equation in the fluid is

[TABLE]

and holds in the domain bounded above by the continuous function $d(x)$ Thus (2.14), (2.13) and (2.8) represents the continuum model governing equations for floating ice. We note that it coincides with the mass-loading model for floating ice in the case $d(x)$ is a constant, assumed by Weitz and Keller (1950), Keller and Weitz (1953), Wadhams and Holt (1991), Squire et al. (2009), Mosig et al. (2017), etc.

We are concerned with the solution of a scattering problem where waves are incident from $x<0$ through a region of broken ice of constant submergence $d(x)=d_{0}$ into a region in $x>0$ where the ice thickness/submergence is varying randomly. In this case we pose far field boundary conditions

[TABLE]

where it is assumed that randomness gives rise to wave attenuation and where $R_{\infty}$ is the reflection coefficient satisfying $|R_{\infty}|=1$ is required to conserve energy whilst

[TABLE]

is the depth variation associated with propagating waves which is found by separating variables; $k$ is the positive real root of

[TABLE]

In Section 4 the extent of this region is finite and for $x>L$ (where $L\gg 1$ implies the region is large compared to the wavelength) $d(x)=d_{0}$ once again and the boundary conditions at infinity are replaced by

[TABLE]

where $R_{L}$, $T_{L}$ are reflection and transmission coefficients satisfying the energy relation $|R_{L}|^{2}+|T_{L}|^{2}=1$.

In the region of variable ice thickness we define

[TABLE]

where $r$ is a random function with zero mean and unit variance, expressed as

[TABLE]

whilst $r(x)$ is chosen to satisfy a Gaussian correlation relation

[TABLE]

We assume $\sigma\ll 1$ as the RMS of the vertical variations of $d(x)$ with respect to $d_{0}$ and $\Lambda=O(1)$ ensures that significant variations in $x$ are on the scale of the wavelength. We also ensure that the $r(0)=r^{\prime}(0)=0$ so that the ice thickness joins the constant values in $x<0$ smoothly. The same applies in $x>L$ by also letting $r(L)=r^{\prime}(L)=0$.

3 A multiple scales expansion to determine attenuation coefficient

We consider the limit $L\to\infty$. The method used to calculate the attenuation rate is based on the the approach documented in Mei et al. (2005) and altered to address random variations in ice thickness over a constant depth fluid in place of random variations in depth beneath a fluid free surface. Additionally, we adopt the variation in the approach developed in Dafydd and Porter (2024) which corrects for the overprediction in attenuation rate, as identified in Bennetts et al. (2015), by considering waves entering a semi-infinite region of randomness and eliminating any terms contributing to wave attenuation that are not associated with multiple scattering.

We introduce a slow variable $X=\sigma^{2}x$ where $\sigma\ll 1$ and write

[TABLE]

Consequently, $\partial_{x}\to\partial_{x}+\sigma^{2}\partial_{X}$ and so (2.14) is

[TABLE]

The bed condition (2.8) becomes

[TABLE]

and the condition (2.13) is expanded as

[TABLE]

and evaluated on $z=-d_{0}$. We assume that the terms of $O(l_{i}^{2})$ neglected in formulating the continuum description of the floating ice are smaller than $O(\sigma^{2})$. Thus, the theory described below formally applies for sufficiently narrow ice floes.

We are also required to match to (2.15) in $x<0$ which results in

[TABLE]

whilst $\phi_{1}(x,0,z)$, $\phi_{2}(x,0,z),\ldots$ should behave as outgoing waves into $x<0$ (i.e. proportional to $\mathrm{e}^{-\mathrm{i}k_{0}x}$ as $x\to-\infty$). As $x\to\infty$ we require $\phi_{i}\to 0$, $i=0,1,2,\ldots$.

3.1 Order $1$

At leading order we have

[TABLE]

with

[TABLE]

and

[TABLE]

Solutions in $x>0$ are therefore given by

[TABLE]

where $k$, $Z_{0}$ are defined by (2.17), (2.16) and where $A$ and $B$, representing the macroscale modulation to the wave amplitude, are to be found. We also have from (2.15),

[TABLE]

and so $B(0)=R_{\infty}$ and $A(0)=1$.

Following the arguments in Dafydd and Porter (2024), $|R_{\infty}|=1$ requires there to be a local balance of energy $|A(X)|=|B(X)|$ and so $B(X)$ only varies from $A(X)$ by a phase factor. This means the second term in (3.11) can be written as the complex conjugate of the first term, modulo a phase. In order to simplify the presentation henceforth we present the calculations with the first term and it is to be assumed that the second term should be added to this.

3.2 Order $\sigma$

At next order we have

[TABLE]

with

[TABLE]

and

[TABLE]

Using (3.11) this last equation can be rewritten as

[TABLE]

(recalling that we are henceforth suppressing the term proportional to $B(X)$). Additionally, $\phi_{1}(x,0,z)$ is outgoing as $x\to-\infty$ where $r(x)=0$ and $\phi_{1}\to 0$ as $x,X\to\infty$ for $-h<z<d_{0}$. The solution to the problem above is formulated using a Green’s function $G(x,\check{x},z)$ defined in $-\infty<x,\check{x}<\infty$ and $-h<z<-d_{0}$ satisfying

[TABLE]

We also impose a radiation condition on $G$ that waves are outgoing from the source, $x=\check{x}$. Employing Green’s Identity with $\phi_{1}$ and $G(x;\check{x},z)$ we have

[TABLE]

where $\partial D$ is the boundary of the area defined in the left-hand integral. Using the conditions of the problem for $\phi_{1}$ and those placed upon $G$ results in the left-hand side of (3.20) being zero and only the line $z=-d_{0}$, $0<x<\infty$ contributes from the boundary integral so that

[TABLE]

Using Fourier transforms (see e.g. Mei et al. (2005)), it is readily found that

[TABLE]

where $\Gamma$ represents the contour along the real axis of the complex $\alpha$-plane with indents below the pole at $\alpha=k_{0}$ and above the pole at $\alpha=-k_{0}$. By residue calculus, it can be shown that, along $z=-d_{0}$,

[TABLE]

where $\alpha=\pm\mathrm{i}\gamma_{n}$ are the location of poles on the imaginary axis and $\gamma_{n}$ denotes the positive real roots of

[TABLE]

Thus, using the symmetry in $x,\check{x}$ of the Green’s function, the solution to (2.16) on the boundary $z=-d_{0}$ is given by

[TABLE]

3.3 Order $\sigma^{2}$

The governing equations for the second-order potential $\phi_{2}$ are

[TABLE]

[TABLE]

and

[TABLE]

We can use (3.11), (3.13), (3.14) to rewrite (3.28) as

[TABLE]

Next we ensemble average the governing equations for $\phi_{2}$ using (3.11) to give

[TABLE]

and

[TABLE]

and

[TABLE]

using (2.20) and noting that $\langle r(x)r^{\prime}(x)\rangle=0$ and $\langle r^{\prime}(x)^{2}\rangle=2/\Lambda^{2}$. The first two terms on the right-hand side need to be considered carefully before addressing the problem for $\langle\phi_{2}\rangle$ further. Thus, we first note from (3.25) that

[TABLE]

since $\partial_{x}{\mathcal{G}}=-\partial_{\check{x}}{\mathcal{G}}$ is Cauchy singular and the integral is thus momentarily interpreted as a Cauchy-principal value integral. Integrating by parts and using the vanishing of $r$ and $r^{\prime}$ at the limits of the integral gives

[TABLE]

In exactly the same manner we find

[TABLE]

The next step is to multiply (3.34) and (3.35) by $r^{\prime}(x)$ and $r(x)$, respectively, and take ensemble averages according to (2.20) and (2.21). This results in

[TABLE]

and

[TABLE]

As described in Dafydd and Porter (2024), we choose to discard contributions that lead to “fictitious decay” through phase cancellations in the averaging process and these are identified as the integrals over $0<\check{x}<x$. Thus, the remaining integrals from $\check{x}>x$ contain contributions that arise from scattering events upwave of the source point only, meaning that we proceed with the expressions

[TABLE]

and

[TABLE]

We remark that a later consequence of this choice is that there is local conservation of energy, $|A(X)|=|B(X)|$, with both $A$ and $B$ decaying as $X\to\infty$. In (3.32) we note the particular combination of (3.38) and (3.39) featured leads to the following simplification

[TABLE]

Using the substitution $\xi=(\check{x}-x)/\Lambda$ allows us to write

[TABLE]

We can see that the integral is independent of $x$ and the right-hand side is proportional to $\mathrm{e}^{\mathrm{i}k_{0}x}$. In order to proceed, we insert the expression for ${\mathcal{G}}$ into the integral above. Thus, we have from (3.23)

[TABLE]

The integrals above are calculated by integrating by parts and using the results found in Mei et al. (2005) to give

[TABLE]

and

[TABLE]

where $\mbox{erfc}(\cdot)$ and $\mbox{erfi}(\cdot)$ are the complementary error function and the imaginary error function respectively. Returning to the equations (3.32), (3.41) we now note that $\mathrm{e}^{\mathrm{i}k_{0}x}$ is a common factor on the right-hand side of the equations governing $\langle\phi_{2}\rangle$. Thus we employ the ansatz $\langle\phi_{2}\rangle=\mathrm{e}^{\mathrm{i}k_{0}x}F(X,z)$ from which it follows that $F$ satisfies

[TABLE]

and

[TABLE]

and

[TABLE]

The problem for $F(X,z)$ does not need to be solved in order to progress. Instead, following Mei et al. (2005), we apply Green’s identity to $F(X,z)$ and $Z_{0}(z)$ over $-h<z<-d_{0}$

[TABLE]

and it follows that

[TABLE]

where (3.23), (3.24) have been used along with $Z_{0}^{\prime\prime}(z)=k_{0}^{2}Z_{0}(z)$, $Z_{0}^{\prime}(-h)=0$, $Z_{0}^{\prime}(-d_{0})=\kappa Z_{0}(-d_{0})$ and

[TABLE]

according to the definition of $Z_{0}(z)$ in (2.16). Bringing the results above together gives

[TABLE]

where

[TABLE]

We can confirm that the infinite series converges since $\gamma_{n}\sim n\pi/(h-d_{0})$, $Z_{n}^{2}(-d_{0})\sim 2$ and

[TABLE]

as $n\to\infty$. The solution of (3.51) with $A(0)=1$ is given by

[TABLE]

In parallel to the calculations we imagine having also made calculations for terms proportional to $B(X)$ whose functional dependence of $x$ is the conjugate of the terms multiplying $A(X)$. This means that $B(X)=R_{\infty}\mathrm{e}^{-\lambda*X}$. Thus, the leading order solution is, from (3.11)

[TABLE]

where $*$ denotes complex conjugation and since $|R_{\infty}|=1$ we can write $R_{\infty}=\mathrm{e}^{\mathrm{i}\delta}$. Thus, the theoretical prediction for the attenuation coefficient is

[TABLE]

The first term in (3.52) is the dominant term in the attenuation coefficient and originates from multiple scattering from propagating waves; the infinite series in (3.52) has a real part that also contributes to a much lesser extent (as evidenced by numerical results) and originates from evanescent wave interactions. Thus, a reasonable approximation to make is

[TABLE]

In deep water, $(h-d_{0})$ is large, $Z_{0}(-d_{0})\sim\sqrt{2k_{0}(h-d_{0})}$ and so

[TABLE]

In shallow water, $(h-d_{0})$ is small, $Z_{0}(-d_{0})\sim 1$ and so

[TABLE]

and this agrees with Dafydd and Porter (2024).

These formulae are unchanged upon redimensionalisation of variables. We note that since $k_{0}\propto\omega$ for shallow water, $\langle\langle k_{i}\rangle\rangle\propto\omega^{2}$ for small frequencies. In deep water $k_{0}\propto\omega^{2}$ and so $\langle\langle k_{i}\rangle\rangle\propto\omega^{8}$ for small frequencies.

In all depths, there is a peak in the value of $\langle k_{i}\rangle$ at intermediate values of frequency before exponential decay at higher frequencies. The peaks are at $k_{0}\Lambda\approx\sqrt{2}$ and $k_{0}\Lambda\approx 1$ for deep and shallow water limits respectively. The theory has been developed under the assumption that $k_{0}\Lambda\not\gg 1$ and so we may expect the validy of the theoretical results to be in question for wave frequencies far beyond the peak although there is no reason to suppose that the existence of a peak attenuation frequency falls outside the modelling assumptions.

4 A mild-slope equation for broken ice

In order to test the theory developed above, we want to compare with numerical simulations of wave scattering through long finite regions of broken ice of randomly-varying thickness. In pursuit of this goal, we set out in this section to formulate an ODE for wave scattering by means of depth-averaging the original two-dimensional boundary-value problem described in Section 2.

Consider the variational principle $\delta\mathcal{L}=0$ where $\mathcal{L}$ is a functional defined by

[TABLE]

where $D$ denotes the horizontal interval of the fluid/ice domain. Assuming variations which vanish at the boundary of the domain (or satisfy a radiation condition if the domain is infinite) it is readily shown that $\mathcal{L}$ is stationary at $\psi=\phi$ if and only if $\phi$ satisfies the governing equations (2.13), (2.14) and (2.8).

We can generate approximate solutions to the governing equations using the variational principle $\delta L=0$. Specifically, we choose to approximate $\psi\approx\phi$ using the standard mild-slope ansatz that

[TABLE]

where

[TABLE]

represents the depth variation associated with propagating waves and $k=k(x)=k(d(x))$ is the real, positive root of

[TABLE]

being the dispersion relation for $k$ assuming $d(x)$ is locally a constant. The function ${\hat{\psi}}(x)$ is determined by substitution of (4.2) into (4.1) and enforcing $\delta\mathcal{L}=0$ to variations in ${\hat{\psi}}(x)$ which vanish on the boundaries of $D$. After integrating by parts and applying the Leibniz rule, we find that ${\hat{\psi}}(x)$ is determined by solutions of

[TABLE]

where

[TABLE]

We note that

[TABLE]

is zero on $z=-h$ and satisfies $w_{z}=w/(1-d(x))$ on $z=-d(x)$ by virtue of (4.4). Letting

[TABLE]

and substituting into (4.5) we find

[TABLE]

where $r=R(w)$ is defined by

[TABLE]

Henceforth, we will use a dot to denote differentiation with respect to $d$ and note from (4.4) that

[TABLE]

It follows that

[TABLE]

where, after some algebra, it can be found that

[TABLE]

Next, we have

[TABLE]

and note that if we write

[TABLE]

then

[TABLE]

Thus, we arrive at the mild-slope equation for broken ice,

[TABLE]

In the language of Chamberlain and Porter (1995) this version is prefixed with ‘Modified’ since the extra terms proportional to $d^{\prime 2}(x)$ and $d^{\prime\prime}(x)$ were omitted from previous versions of the mild-slope equations (in this case for variable bathmetry and an unloaded free surface). Porter (2020) later showed how a transformation of (4.17) can produce a simplified version of the mild slope equation and we follow that approach below.

We seek to transform ${\hat{\psi}}$ by letting

[TABLE]

where the scaling factor $s(d(x))$ is to be determined. We now have from (4.17) that

[TABLE]

where

[TABLE]

Following Porter (2020) we choose to remove the $d^{\prime\prime}(x)$ term by setting $\tilde{u}_{1}=0$. This can be achieved by noting that $2u_{0}\dot{k}/k=1$ and so

[TABLE]

Using this in (4.20) with $\tilde{u}_{1}=0$ allows us to integrate up to give $u_{0}^{1/2}ks=c$ for some arbitrary constant, $c$ which, without loss of generality, can be set equal to unity. Thus

[TABLE]

defines the scaling factor that eliminates $d^{\prime\prime}(x)$ from the governing equation and it follows that the original equation (4.17) can be recast as

[TABLE]

where the original dependent variable is recovered from

[TABLE]

Now, we can assume the continuous functions $d(x)$ satisfy the mild-slope assumption, i.e.

[TABLE]

which justifies the neglect of the final term in (4.23) so that we end up with the simpler governing equation

[TABLE]

It is notable that this equation was derived by Porter (2020) for scattering of surface waves over variable bathymetry and the only difference here is that $k$ is defined by the dispersion relation (4.4) for variable thickness ice on the surface over a flat bed instead of the dispersion relation for an unloaded surface over a variable bed.

It can be shown that the inclusion of variable depth in addition to variable ice thickness leads to the same equation with $k$ defined by (4.4) but with the constant $h$ replaced by $h(x)$.

5 Results and discussion

The purpose of this section is to compare results from numerical simulations of wave scattering over long finite regions of randomly-varying ice thickness to the theoretical predictions of Section 3. Simulations are calculated using the MSEBI described in Section 4 in which the problem of wave scattering is reduced through depth averaging to solving an ordinary differential equation (ODE) based on the assumption that $d(x)$ is a slowly-varying function. The computation of solutions to the ODE is described in Dafydd and Porter (2024).

For the simulations, $d(x)$ is taken to be a random function defined by $d(x)=d_{0}(1+r(x))$ with $r(x)$ satisfying $\langle r(x)\rangle=0$, $\langle r(x)^{2}\rangle=1$ in $0<x<L$ with $d(x)=d_{0}$ for $x<0$ and $x>L$. The method for defining $r(x)$ is described in Dafydd and Porter (2024). The randomness is characterised by $\sigma$ which is the RMS of vertical variations with respect to $d_{0}$ and $\Lambda$ which is the correlation lengthscale. The attenuation is computed from the eigenvalues of the transfer matrix that arises from the propagation of wave solutions over $0<x<L$ in the same manner as described in Dafydd and Porter (2024). In order to compute approximations to the attenuation coefficient we must fulfil two requirements. The first is that the value of $L$ is large enough and here $L=20\Lambda^{3}/\sigma$ is chosen to anticipate the length required for convergence to be of the form $\langle k_{i}\rangle=C$ for some constant $C$. This choice also benefits from being frequency independent, allowing for multiple values of $k_{0}\Lambda$ per realisation. As in Dafydd and Porter (2024), this form has the benefits of it being frequency independent and convergent at the same point independent of $\sigma$ and $\Lambda$ at higher frequencies. The second is that we average each numerically-derived $k_{i}$ over sufficiently large number, $N$, of independent random realisations of $d(x)$. In Figure 2 the convergence of the average attenuation, $\langle k_{i}\rangle$, with increasing $N$ is shown for values up to 1500 for the MSEBI computations. The three pairs of curves show the differences from different $\sigma$ taken when $k_{0}\Lambda=1$. On this basis we have chosen $N=500$ for subsequent computations.

In Figure 4 we have plotted the dimensionless attenuation coefficient $\Lambda\langle k_{i}\rangle/\sigma^{2}$ against $K\Lambda$ for $\sigma^{2}=0.01,~0.02,~0.04$. The water depth increases from shallow in subplot (a) to deep in subplot (d). In each example, $\Lambda/d_{0}=2$ is fixed and the data for different $\sigma^{2}$ should collapse onto the same curve.

Simulations, averaged over $N=500$ realisations, are computed using the MSEBI method. These results are shown to follow the theoretical predictions of Section 3. The effect of adding the contribution from the infinite sum to the formula for $\langle k_{i}\rangle$ is virtually indistinguishable from the leading order term (3.57). We note that the when the depth is shallow the curves in Figure 4 are similar in character to those presented in Dafydd and Porter (2024).

For small $k_{0}(h-d_{0})$. and therefore small $K(h-d_{0})$ where $K=\omega^{2}/g$, we see that the attenuation is proportional to $\omega^{2}$. For larger depths, i.e. $k_{0}(h-d_{0})\gg 1$, the low frequency behaviour changes and the attenuation tends towards $\omega^{8}$.

For larger values of $\sigma^{2}$, the discrepancy between the theoretical predictions and numerical simulations increases at higher frequencies. This behaviour is anticipated, as the theory is formally valid only for low frequencies and small $\sigma$. Since a similar theoretical result can be obtained from the depth-averaged model in the Appendix, it is likely that the primary source of this discrepancy is the increased value of $\sigma$.

It should be noted that greater agreement with the theory can be achieved through increasing the number of random realisations, $N$, used in the averaging process and increasing the surface length $L$. However, this latter alteration comes at a great computational cost as due to the rapidly-varying nature of the function $d(x)$ in contrast to the slowly-varying $k(x)$ and thus a very small tolerance and step-size is required alongside the use of an implicit ODE solver. In Figure 6, a plot showing the convergence is shown for $k_{0}\Lambda=1.5$ shows that perhaps values of $L/d_{0}\geq 100\Lambda^{3}/\sigma^{2}{d_{0}}^{3}$ maybe be required to achieve high levels of convergence around the peak of the curve. This is computationally unfeasible at the precision we require.

The low frequency attenuation scaling of $\omega^{8}$ for deep water predicted in this paper is outside the range of $\omega^{2}$ to $\omega^{4}$ normally associated with field measurements (see Meylan et al. (2018)). However, the work of Squire et al. (2009) and Meylan et al. (2021) both suggest that at very low frequencies the attenuation data may scale somewhere in the range of $\omega^{8}$ to $\omega^{10}$. Although the current basic model for attenuation based on two-dimensional multiple scattering due to randomness in the thickness of floating ice does not appear to fit data across all frequencies, it may still form a good basis for describing attenuation at low frequencies. Additionally, the work of Herman (2024) suggests that a model’s efficacy does not correspond to its ability to capture the power law of the “apparent attenuation” as these measurements do not adjust for energy loss that is not ice-related and in fact a higher power of $\omega$ may in fact be more desirable than the commonly desired $\omega^{2}$ to $\omega^{4}$. We should be careful not to read too much into other features of our results. However, it is interesting to note that the present model, like that of Dafydd and Porter (2024), exhibits a high-frequency “roll-over” effect at all depths which is a feature of field data reported in numerous work including Wadhams et al. (1988), Liu et al. (1992), Squire et al. (2009), Doble et al. (2015). On the other hand the onset of roll-over in the data is itself disputed with Thomson et al. (2021) providing evidence that this is a statistical effect caused by noisy data.

Within the general framework of the theory developed here, is easy to imagine that any model of multiple scattering through a random environment defined by a spatially-varying local wavenumber will result in a high-frequency roll-over effect. There is therefore an argument that randomness and multiple scattering within a more complex description of floating ice may lead to predicitions of attenuation which follow the field measurements. A modification to the model presented here would imagine that individual floes of ice are separated by small gaps of open water and, instead of moving in constrained heave motion, are allowed to float freely. Ongoing work by the authors have considered such a model to determine a local dispersion relation which assumed identical periodic floes separated by equal gaps in terms of the ice thickness, and the size of the gaps.

6 Conclusions

This paper has considered a linear model for propagation of waves over regions of variable ice in a fluid of finite depth. A simple expression has been analytically derived for the ensemble averaged attenuation of waves propagating through fragmented floating ice over a semi-infinite domain. In this derivation, terms associated with incoherent phase cancellations were removed to preserve only attenuation due to multiple scattering effects. The theory was shown to not only agree with the prior work of Dafydd and Porter (2024) in the shallow water regime, but also compares favourably with numerical simulations based on a finite-depth. The mild-slope equation for broken ice of slowly-varying thickness was derived in the simple two-dimensional case over a constant dept, relevant to the present study, but can easily be extended to three-dimensional scattering over variable bathymetry.

The results have provided a theoretical prediction for the frequency dependence of attenuation through broken ice which has not captured the measurements reported in the literature across all frequencies. However, there are some features of the results which positively correlate to the observed data and suggest that randomness and multiple scattering are a dominant factor in attenuation of waves through regions of broken ice. It is important to note however that this model does not account for many of the complicated aspects of multi-phase, discrete sea ice. Further work such as discretising the blocks of sea ice and adding gaps between them to make them non-continuous are planned, as well as randomising gaps and three-dimensional scattering.

[Funding] L.D. is grateful for the support of an EPSRC (UK) studentship.

[Declaration of interests] The authors report no conflict of interest.

**[Author ORCIDs]

**

L. Dafydd, https://orcid.org/0000-0003-1009-0946

R. Porter, https://orcid.org/0000-0003-2669-0188

Appendix: derivation of attenuation coefficient from the MSE for broken ice

The following section is less rigorous due to being reliant on multiple small parameters, but is included for completeness. Using (4.26) we can trivially perform an asymptotic expansion by defining the ice submergence by the previously defined parameters (2.19-2.21) and then expanding our function in powers of $\sigma$ such that

[TABLE]

in $x>0$ where $X=\sigma^{2}x$ is a slow variable. We can then rewrite (4.26) in the form

[TABLE]

where

[TABLE]

is found from from Taylor expanding $1/k^{2}(d(x))$ in $d(x)$ in powers of $\sigma$. From Dafydd and Porter (2024) we know that (A.2) gives us an attenuation coefficient of the form

[TABLE]

and we can substitute our value for $C_{1}$ in to get

[TABLE]

This is identical to the formula derived for our finite depth model without the small evanescent terms.