2. Summary of the model

We consider a two-dimensional scattering problem in which plane-crested monochromatic waves of small amplitude propagate in the positive $x$-direction in $x < 0$ over fluid of constant depth with a surface covered by a continuous layer of fragmented ice of constant thickness. There are no physical mechanisms included in the model for energy dissipation such as fluid viscosity or ice–ice friction. Incident wave energy is partially reflected from, and partially transmitted into, the region $x >0$. This is due to either randomly varying bathymetry or by randomly varying thickness of broken ice (both are illustrated in figure 1) which extends over the interval $0 < x < L$ before returning, in $x > L$, to the same constant values found in $x < 0$. We are interested in monitoring the reflected and transmitted wave energy. In § 4 we set $L = \infty$ so that the randomness extends indefinitely into $x > 0$. In this case all incoming wave energy will be reflected and the focus is determining the attenuation of waves as a function of distance into $x>0$.

Figure 1.

Definition sketch of variable floating broken ice over a variable bed.

Porter (Reference Porter2019) developed a shallow-water (long-wavelength) model for wave scattering over variable bathymetry with no ice cover. This model results from an expansion to second order in a small parameter representing the ratio of vertical to horizontal lengthscales combined with depth averaging and is expressed by

(2.1)\begin{equation} ({\hat {\hat h}}(x) \varOmega'(x))' + K \varOmega(x) = 0, \end{equation}

where $K = \omega ^2/g$, $\omega$ is the angular frequency of the motion, $g$ represents gravitational acceleration and

(2.2)\begin{equation} {\hat {\hat h}}(x) = \frac{h(x)\left(1-\dfrac13 K h(x)\right)}{1 + \dfrac13 v(h) h'^2(x)} \end{equation}

is defined in terms of the fluid depth $h(x)$. Here, $v(h) = 1 + \frac {1}{12} Kh(x)/(1-\frac 13 Kh(x))$ and $v(h) \approx 1$ is a simplification which will be adopted hereafter. The underlying assumptions are expressed by the formal constraint that $K h \ll 1$, although Porter (Reference Porter2019) showed by comparing with exact results for reflected and transmitted wave energy for shoaling beds of finite length, that the model produces accurate predictions up to $Kh \approx 1$.

The dependent variable, $\varOmega$, in (2.1) is related to the time-independent wave elevation $\eta (x)$ obtained under the time-harmonic assumption $\zeta (x,t) = {\rm Re} \{ \eta (x) {\rm e}^{{-\mathrm {i} \omega t}}\}$ by

(2.3)\begin{equation} \eta(x) = \frac{-(\mathrm{i} /\omega)}{\sqrt{1 - \dfrac13 K h(x)}} \left(\varOmega(x) - \frac{\dfrac16 h h'}{1+\dfrac13 h'^2} \varOmega'(x)\right) \end{equation}

and is referred to as the ‘pseudo-potential’ by Toledo & Agnon (Reference Toledo and Agnon2010). It was shown in Porter (Reference Porter2019) that $\varOmega (x)$ and $\varOmega '(x)$ remain continuous at discontinuities in $h'(x)$.

Porter (Reference Porter2019) highlighted the significant improvement in results away from the zero frequency limit that could be achieved when ${\hat {\hat h}}(x) = h(x)$ is replaced by the definition in (2.2), applying in the case of the standard linear shallow-water equation. Thus, the modification in (2.2) includes, in the numerator, the effect of weak dispersion and, in the denominator, a geometric factor indicating a reduction in wave speed over sloping beds. We also remark that (2.1) can also be derived from a linearisation of Boussinesq equations (e.g. Peregrine Reference Peregrine1967) whereby wave amplitudes are assumed sufficiently small compared with $Kh$.

In the Appendix, the model developed by Porter (Reference Porter2019) is extended to include the additional effect of a floating fragmented ice cover. Additional assumptions apply here. Ice is assumed to completely cover the surface of the fluid and is broken into sections which are sufficiently small in horizontal extent and whose thickness varies slowly enough that the submergence of the ice is represented by a continuous function, $d(x)$. Thus, the model is simulating the effect of randomness within the ice cover as rather than from incoming waves approaching the cover. The motion of the ice is constrained in heave (vertical) motion and the expansion to second order of the depth ratio ($\epsilon$ in the Appendix) in the modelling is needed to include the effect of inertia of floating ice. That is, a basic first-order linear shallow-water model neglects vertical accelerations and the effect of ice cover at leading order is manifested only through a reduction in the depth of the fluid from $h(x)$ to $h(x) - d(x)$. Thus, our second-order model extended to incorporate floating ice of submergence $d(x)$ is, see (A38),

(2.4)\begin{equation} ({\hat {\hat d}}(x) \varOmega'(x))' + K \varOmega(x) = 0, \end{equation}

where ${\hat {\hat d}}(x)$ is defined by (A39) and the loaded surface elevation is related to $\varOmega$ by (A40). As before, $\varOmega$ and $\varOmega '$ are continuous even if $d'(x)$ and/or $h'(x)$ is discontinuous.

In $x < 0$ and in $x > L$ we assume $h = h_0$, $d = d_0$ are both constant. Then (2.4) can be solved explicitly and

(2.5)\begin{gather} \varOmega (x) = \mathrm{e}^{\mathrm{i} k_0 x} + R_L \mathrm{e}^{-\mathrm{i} k_0 x}, \quad x < 0, \end{gather}

(2.6)\begin{gather}\varOmega(x) = T_L {\rm e}^{\mathrm{i} k_0 x},\quad x > L, \end{gather}

where $R_L$ and $T_L$ are reflection and transmission coefficients, satisfying $|R_L|^2 + |T_L|^2 = 1$ (energy conservation) and

(2.7)\begin{equation} k_0^2 (h_0 - d_0) = \frac{K}{1 - \dfrac13 K (h_0 + 2 d_0)} \end{equation}

defines the wavenumber, $k_0$, in terms of the frequency, $\omega$. This shallow-water dispersion relation is weakly dispersive, but for sufficiently small frequencies we note that $k_0 \propto \omega$.

3. Description of randomness

We consider wave propagation over a region $0 < x < L$ in which either the bed or the ice thickness randomly varies. We could consider both simultaneously varying, but for clarity consider the two effects separately.

We say that either

(3.1a,b)\begin{equation} d = 0,\quad h(x) = h_0(1+ \sigma r(x)) \end{equation}

or that

(3.2a,b)\begin{equation} h = h_0,\quad d(x) = d_0(1 + \sigma r(x)) \end{equation}

such that $r(x)$ is a random function with mean zero and unit variance. That is,

(3.3a,b)\begin{equation} \langle r \rangle = 0,\quad \langle r^2 \rangle = 1, \end{equation}

implying that $\sigma$ is the root mean square (r.m.s.) of the vertical variations of $h(x)$ or $d(x)$. We ensure that the $r(0) = r'(0) = r(L) = r'(L) = 0$ so that the bed/ice thickness joins the constant values in $x < 0$ and $x > L$ smoothly. The random function $r(x)$ also satisfies the Gaussian correlation relation

(3.4)\begin{equation} \langle r(x) r(x') \rangle = {\rm e}^{-|x-x'|^2/\varLambda^2} \end{equation}

(other models have used an exponential correlation function, but show that it produces only small differences in results). Thus, $\varLambda$ characterises the horizontal length scale of the random bed fluctuations.

4. Analysis of the model

In this section, we assume $L \to \infty$ so that the randomness occupies $x > 0$. The main assumption that is made is that the amplitude of the randomness is small, i.e. $\sigma \ll 1$. We assume $\sigma = O(\epsilon )$ and will expand up to $O(\sigma ^2)$ to be consistent with the $O(\epsilon ^2)$ expansion derived in the Appendix. We note that we can write (2.4) with (A39), (A41) and either (3.1a,b) or (3.2a,b) as

(4.1)\begin{equation} ((1+ \sigma C_1 r(x) - \sigma^2 (C_2r^2(x) + C_3 r'^2(x))) \varOmega')' + k_0^2 \varOmega = 0,\quad x > 0 ,\end{equation}

where terms up to $O(\sigma ^2)$ have been retained, and

(4.2)\begin{equation} \varOmega'' + k_0^2 \varOmega = 0,\quad x < 0 ,\end{equation}

where $k_0$ is defined by (2.7). In (4.1), the coefficients depend on the whether the bed or the thickness of floating ice is represented by the random function $r(x)$. In the case that the bed is varying and the ice is absent, $d_0 = 0$ and

(4.3a–c)\begin{equation} C_1 = \frac{1-{\dfrac23} Kh_0}{1 - {\dfrac13} Kh_0},\quad C_2 = \frac{{\dfrac13} Kh_0}{1 - {\dfrac13} Kh_0},\quad C_3 = {\frac13} h_0^2 \end{equation}

and in the case where the ice is varying and the bed is of constant depth, $h(x) = h_0$ and

(4.4a–c)\begin{align} \left.\begin{gathered} C_1 = \frac{-d_0\left(1+{\dfrac13}K(h_0-4d_0)\right)}{(h_0-d_0)\left(1-{\dfrac13}K(h_0+2d_0)\right)},\quad C_2 = \frac{-{\dfrac23} K d_0^2}{(h_0-d_0)\left(1-{\dfrac13}K(h_0+2d_0)\right)},\\ C_3 = {\frac13} d_0^2. \end{gathered}\right\} \end{align}

The long-wave assumption on which the model is based formally requires $Kd_0 < Kh_0 \ll 1$ and so we do not envisage using the model close to $Kh_0 = 3$ or $K(h_0+2d_0)=3$. The solution to (4.2) is

(4.5)\begin{equation} \varOmega(x) = {\rm e}^{\mathrm{i} k_0 x} + R_\infty {\rm e}^{-\mathrm{i} k_0 x} \end{equation}

and since we anticipate decay of waves into $x \to \infty$ we also impose $\varOmega \to 0$ as $x \to \infty$ and so we must require that $|R_\infty | = 1$; all incident wave energy is reflected.

We make the multiple scales assumption of, e.g. Mei & Li (Reference Mei and Li2004) (but also see other references listed in the introduction) and introduce a slow variable $X = \sigma ^2 x$, writing

(4.6)\begin{equation} \varOmega(x) = \varOmega_0(x,X) + \sigma \varOmega_1(x,X) + \sigma^2 \varOmega_2(x,X) + \cdots. \end{equation}

Accordingly (4.1) becomes

(4.7)\begin{align} & \left[\left(\frac{\partial}{\partial x} + \sigma^2 \frac{\partial}{\partial X} \right) \left( 1 + \sigma C_1 r(x) - \sigma^2 (C_2r^2(x)+ C_3 r'^2(x))) \left(\frac{\partial}{\partial x} + \sigma^2\frac{\partial}{\partial X} \right) \right)+k_0^2 \right] \nonumber\\ &\quad \times(\varOmega_0 + \sigma \varOmega_1 + \sigma^2 \varOmega_2 +\cdots ) = 0,\quad x > 0. \end{align}

The matching conditions at $x=0$ consist of

(4.8)\begin{equation} \varOmega(0^-) = 1 + R_\infty = (\varOmega_0 + \sigma \varOmega_1 + \sigma^2 \varOmega_2+ \cdots)_{x=X=0} \end{equation}

and

(4.9)\begin{equation} \varOmega'(0^-) = ik_0(1 - R_\infty) = \left(\frac{\partial}{\partial x} + \sigma^2\frac{\partial}{\partial X}\right) (\varOmega_0 + \sigma \varOmega_1 + \sigma^2 \varOmega_2+ \cdots)_{x=X=0}. \end{equation}

At leading order, $\varOmega _0$ satisfies the same wave equation (4.2) as in $x < 0$ and its general solution is

(4.10)\begin{equation} \varOmega_0(x,X) = A(X) {\rm e}^{\mathrm{i} k_0 x} + B(X) {\rm e}^{-\mathrm{i} k_0 x}. \end{equation}

This implies that the leading order solution is not explicitly dependent on individual realisations, $r(x)$; $A$ and $B$ will contain information relating to the statistical properties of $r(x)$ however. We require that long-scale variations, $A(X)$ and $B(X)$, tend to zero as $X \to \infty$, whilst $A(0) = 1$ and $B(0) = R_\infty$ are determined from the matching conditions (4.8) and (4.9) at leading order.

Since $|R_\infty | = 1$ there must be no net time-averaged transport of energy flux in $x > 0$ and so we expect that

(4.11)\begin{equation} |A(X)| = |B(X)|. \end{equation}

At $O(\sigma )$ we have

(4.12)\begin{equation} \frac{\partial^2 \varOmega_1}{\partial x^2} + k_0^2 \varOmega_1 ={-}C_1 \frac{\partial}{\partial x} \left(r(x)\frac{\partial \varOmega_0}{\partial x}\right). \end{equation}

Its solution can be determined using the Green's function for the one-dimensional wave equation,

(4.13)\begin{equation} g(x,x') = \frac{{\rm e}^{\mathrm{i} k_0 |x-x'|}}{2 \mathrm{i} k_0}, \end{equation}

satisfying

(4.14)\begin{equation} \frac{\partial^2}{\partial x^2} g + k_0^2 g = \delta(x-x'), \end{equation}

and outgoing as $|x-x'| \to \infty$. The right-hand side of (4.12) is composed of two terms forced by right- and left-propagating waves and the solution $\varOmega _1$, in $x > 0$, is a superposition of solutions derived using $g$ and $\bar {g}$ (where the overbar denotes complex conjugate), respectively, in Green's identity with the two components of $\varOmega _1$ over $x > 0$ and results in

(4.15)\begin{align} \varOmega_1(x,X) &={-} \mathrm{i} k_0 C_1 A(X) \int_0^\infty g(x,x') \frac{\partial}{\partial x'} (r(x') {\rm e}^{\mathrm{i} k_0 x'}) \, \mathrm{d}\kern0.7pt x' \nonumber\\ &\quad +\mathrm{i} k_0 C_1 B(X) \int_0^\infty \bar{g}(x,x') \frac{\partial}{\partial x'} (r(x') {\rm e}^{-\mathrm{i} k_0 x'}) \, \mathrm{d}\kern0.7pt x',\quad x > 0. \end{align}

The use of $\bar {g}$ is non-standard and implies that the component of the first-order solution associated with left-propagating leading-order wave is represented by a distribution of incoming waves. This is required to satisfy the energy balance equation (4.11). Put another way, we require the amplitude, $B(X)$, of the left-going wave to grow as it propagates from right to left, its associated energy being generated from the energy lost to outgoing waves from the right-propagating wave with amplitude $A(X)$.

Integrating by parts once, using $r(0) = 0$ (since the random variations in the bed or the ice continuously joins the constant value set in $x < 0$) gives

(4.16)\begin{align} \varOmega_1(x,X) &={-} \mathrm{i} k_0 C_1 A(X) \int_0^\infty \frac{\partial}{\partial x} g(x,x') r(x') {\rm e}^{\mathrm{i} k_0 x'} \,\mathrm{d}\kern0.7pt x' \nonumber\\ &\quad + \mathrm{i} k_0 C_1 B(X) \int_0^\infty \frac{\partial}{\partial x} \bar{g}(x,x') r(x') {\rm e}^{-\mathrm{i} k_0 x'} \,\mathrm{d}\kern0.7pt x'. \end{align}

Here $\partial _x g = - \partial _{x'} g$ has been used and we note that this function is discontinuous at $x=x'$.

We also remark that $\varOmega _1$ is a random function with zero mean since $\langle \varOmega _1 \rangle = 0$ follows from ensemble averaging (4.16) and using (3.3a,b).

At $O(\sigma ^2)$ we have

(4.17)\begin{align} \frac{\partial^2 \varOmega_2}{\partial x^2} + k_0^2 \varOmega_2 ={-}C_1 \frac{\partial}{\partial x} \left( r(x)\frac{\partial \varOmega_1}{\partial x} \right) - 2 \frac{\partial^2 \varOmega_0}{\partial x \partial X} + \frac{\partial}{\partial x} \left( ( C_2r^2(x)+ C_3 r'^2(x)) \frac{\partial \varOmega_0}{\partial x} \right). \end{align}

We ensemble average the equation using the results from (3.3a,b) and $\langle r'^2 \rangle = 2/\varLambda ^2$ (this can be established using the definition of the derivative as a limit) to give

(4.18)\begin{align} & \frac{\partial^2}{\partial x^2} \langle \varOmega_2 \rangle + k_0^2 \langle \varOmega_2 \rangle ={-} C_1 \frac{\partial}{\partial x} \left\langle r(x) \frac{\partial \varOmega_1}{\partial x} \right\rangle - 2 \mathrm{i} k_0 (A'(X) {\rm e}^{\mathrm{i} k_0 x} - B'(X) {\rm e}^{-\mathrm{i} k_0 x}) \nonumber\\ &\quad -k_0^2 ( C_2 + 2 C_3/\varLambda^2)(A(X) {\rm e}^{\mathrm{i} k_0 x} + B(X) {\rm e}^{-\mathrm{i} k_0 x}). \end{align}

It is instructive to write $\varOmega _1$ from (4.16) in terms of separate wave-like components as

(4.19)\begin{align} & \varOmega_1(x,X) \nonumber\\ &\quad ={-}\frac{C_1 A(X) \mathrm{i} k_0}{2} \left[ {\rm e}^{\mathrm{i} k_0 x} \int_0^x r(x') \, \mathrm{d}\kern0.7pt x' - {\rm e}^{-\mathrm{i} k_0 x} \int_x^\infty r(x') {\rm e}^{2 \mathrm{i} k_0 x'} \,\mathrm{d}\kern0.7pt x' \right] \nonumber\\ &\qquad + \frac{C_1 B(X) \mathrm{i} k_0}{2} \left[ {\rm e}^{-\mathrm{i} k_0 x} \int_0^x r(x') \, \mathrm{d}\kern0.7pt x' - {\rm e}^{\mathrm{i} k_0 x} \int_x^\infty r(x') {\rm e}^{-2 \mathrm{i} k_0 x'} \,\mathrm{d}\kern0.7pt x' \right]. \end{align}

We note that the leading-order right-propagating wave excites both right-propagating waves which accumulate from interactions with the bed to the left of the observation point, $x$, and left-propagating waves which represent the accumulation of upwave reflections from bed interactions to the right of the observation point. Similar comments apply to terms proportional to the leading-order left-propagating wave. The ensemble averaging of the first and third terms of (4.19) in (4.18) lead to a contribution to the attenuation which we describe as ‘fictitious decay’. That is, it is a feature of wave scattering not experienced by individual waves, but which instead originates from phase cancellations from first-order waves when averaged over realisations of $r(x)$. The coefficient multiplying the two ${\rm e}^{\pm \mathrm {i} k_0 x}$ terms under scrutiny is a real integral which depends only on $r(x)$, the geometry and, hence, randomness does not alter the phase of these contributions. This contrasts with the second and fourth terms in (4.19) which correspond to the accumulation of waves that have propagated from the field point $x$ to a point $x'$ and reflected by the bathymetry/broken ice $r(x')$ necessarily encoding randomness into the phase of these contributions. For the purpose of computing the attenuation experienced by individual waves we remove this fictitious decay effect, replacing (4.19) by

(4.20)\begin{align} \varOmega_1(x,X) &= \frac{C_1 A(X) \mathrm{i} k_0}{2} {\rm e}^{-\mathrm{i} k_0 x} \int_x^\infty r(x') {\rm e}^{2 \mathrm{i} k_0 x'}\,\mathrm{d}\kern0.7pt x' \nonumber\\ &\quad -\frac{C_1 B(X) \mathrm{i} k_0}{2}{\rm e}^{\mathrm{i} k_0 x} \int_x^\infty r(x') {\rm e}^{-2 \mathrm{i} k_0 x'} \,\mathrm{d}\kern0.7pt x'. \end{align}

The only term requiring attention now is the first term on the right-hand side of (4.18) where $\varOmega _1$ is given by (4.20). It is straightforward to determine from (4.20) that

(4.21)\begin{align} \left\langle r(x) \frac{\partial \varOmega_1}{\partial x} \right\rangle &={-}\frac{\mathrm{i} k_0}{2} C_1 A(X) {\rm e}^{\mathrm{i} k_0 x} + k_0^2 C_1 A(X) {\rm e}^{\mathrm{i} k_0 x} \nonumber\\ &\quad \times\int_0^\infty {\rm e}^{-\xi^2/\varLambda^2}{\rm e}^{2 \mathrm{i} k_0 \xi} \,\mathrm{d} \xi \nonumber\\ &\quad +\frac{\mathrm{i} k_0}{2} C_1 B(X) {\rm e}^{-\mathrm{i} k_0 x} + k_0^2 C_1 B(X) {\rm e}^{-\mathrm{i} k_0 x} \nonumber\\ &\quad \times\int_0^\infty {\rm e}^{-\xi^2/\varLambda^2} {\rm e}^{-2 \mathrm{i} k_0 \xi} \,\mathrm{d} \xi \end{align}

after using the definition in (3.4) and making a substitution $\xi = x-x'$. As demanded by (4.18), we need to take a further derivative which results in

(4.22)\begin{equation} C_1 \frac{\partial}{\partial x} \left\langle r(x) \frac{\partial \varOmega_1}{\partial x} \right\rangle = \frac{C_1^2 k_0^2}{2} (A(X) F {\rm e}^{\mathrm{i} k_0 x} + B(X) \bar{F} {\rm e}^{-\mathrm{i} k_0 x}), \end{equation}

where

(4.23)\begin{align} F &= 1 +\mathrm{i} k_0 \int_{0}^\infty {\rm e}^{-\xi^2/\varLambda^2} {\rm e}^{2 \mathrm{i} k_0 \xi} \,\mathrm{d} \xi \nonumber\\ &=1 + \frac{\sqrt{\rm \pi}}{2} \mathrm{i} k_0 \varLambda {\rm e}^{-k_0^2 \varLambda^2}(1 + \mathrm{i}\,\mbox{erfi}(k_0 \varLambda)), \end{align}

(see, e.g. Mei & Li Reference Mei and Li2004) and $\mbox {erfi}({\cdot })$ is the imaginary error function.

Armed with (4.22), we return to the governing equation (4.18) for $\langle \varOmega _2 \rangle$ and note that the right-hand side contains secular terms; that is, functions proportional to ${\rm e}^{\pm \mathrm {i} k_0 x}$. These must be removed to avoid unbounded growth in the solution for $\langle \varOmega _2 \rangle$ as $x \to \infty$. In other words, we wish to obtain

(4.24)\begin{equation} \frac{\partial^2}{\partial x^2} \langle \varOmega_2 \rangle + k_0^2 \langle \varOmega_2 \rangle = 0, \end{equation}

requiring $A(X)$ and $B(X)$ to satisfy the solvability conditions

(4.25) \begin{align} 2 \mathrm{i} k_0 A'(X) &={-}k_0^2 A(X) \left(C_1^2 \left( \frac{1}{2} + \frac{\sqrt{\rm \pi}}{4} \mathrm{i} k_0 \varLambda{\rm e}^{-k_0^2\varLambda^2} (1 + \mathrm{i}\,\mbox{erfi}(k_0 \varLambda))\right) \right.\nonumber\\ &\quad \left.\vphantom{\frac{\sqrt{\rm \pi}}{4}} + C_2 + 2 C_3/\varLambda^2\right) \end{align}

and

(4.26) \begin{align} -2 \mathrm{i} k_0 B'(X) &={-}k_0^2 B(X) \left(C_1^2 \left(\frac{1}{2} - \frac{\sqrt{\rm \pi}}{4} \mathrm{i} k_0 \varLambda {\rm e}^{-k_0^2\varLambda^2}(1 - \mathrm{i}\,\mbox{erfi}(k_0 \varLambda)) \right) \right.\nonumber\\ &\quad \left.\vphantom{\frac{\sqrt{\rm \pi}}{4}} + C_2 + 2 C_3/\varLambda^2\right). \end{align}

Solving for $A(X)$ with $A(0) = 1$ gives

(4.27)\begin{equation} A(X) = {\rm e}^{- Q X + \mathrm{i} \kappa X}, \end{equation}

where

(4.28)\begin{equation} Q = \frac{\sqrt{\rm \pi}}{8} C_1^2 k_0^2 \varLambda {\rm e}^{-k_0^2 \varLambda^2} \end{equation}

and

(4.29)\begin{equation} \kappa = C_1^2 \left( \frac{k_0}{4} - \frac{\sqrt{\rm \pi}}{8} k_0^2 \varLambda {\rm e}^{-k_0^2 \varLambda^2} \mbox{erfi}(k_0 \varLambda)\right) + k_0 C_2/2 + k_0 C_3/\varLambda^2. \end{equation}

Meanwhile, solving (4.26) for $B(X)$ with $B(0) = R_\infty$ such that $|R_\infty | = 1$ gives

(4.30)\begin{equation} B(X) = R_\infty {\rm e}^{- Q X - \mathrm{i} \kappa X} \end{equation}

and, thus, (4.16) is satisfied.

Had the first and third terms in (4.19) not been removed and (4.19) not been replaced by (4.20) then, amongst other changes, the expression in (4.28) would have been replaced by $Q = (\sqrt {{\rm \pi} }/{8}) C_1^2 k_0^2 \varLambda (1+{\rm e}^{-k_0^2 \varLambda ^2})$. A similar attenuation factor is determined in the work of Mei et al. (Reference Mei, Stiassnie and Yue2005, § 7.4) and Bennetts et al. (Reference Bennetts, Peter and Chung2015). The additional factor of $+1$, associated with phase cancellation in the ensemble averaging, completely changes the character of attenuation. Bennetts et al. (Reference Bennetts, Peter and Chung2015) highlight the discrepancy between theoretical results and attenuation measured through discrete numerical simulations, most notably in figures 5 and 6 of their paper. Moreover, the expression for $B(X)$ would also change with the factor of $Q$ associated with (4.30) replaced by $Q = (\sqrt {{\rm \pi} }/{8}) C_1^2 k_0^2 \varLambda (-1+{\rm e}^{-k_0^2 \varLambda ^2})$ implying exponential growth towards infinity of the left-propagating wave whilst (4.16) is no longer satisfied.

Returning to (4.10) gives the leading-order solution in $x > 0$ as

(4.31)\begin{align} \varOmega(x) \approx \varOmega_0(x,\sigma^2 x) = {\rm e}^{- \sigma^2 Q x} ({\rm e}^{\mathrm{i} (k_0 + \sigma^2 \kappa) x} + R_\infty {\rm e}^{-{\rm i} (k_0 + \sigma^2 \kappa) x}). \end{align}

Furthermore, since $\langle \varOmega _1 \rangle = 0$, corrections to (4.31) are $O(\sigma ^2)$. From (4.31) the attenuation rate is defined to be

(4.32)\begin{equation} k_i = \sigma^2 Q =\frac{\sqrt{\rm \pi}}{8} k_0^2 \sigma^2 \varLambda C_1^2 {\rm e}^{-k_0^2 \varLambda^2} \end{equation}

with $C_1$ given by (4.3a–c) (or (4.4a–c)), a factor which depends upon $k_0h_0$ (and $d_0/h_0$). In the case of a randomly varying bed with no ice cover and assuming $C_1^2 \approx 1$ since $Kh_0 \ll 1$, the maximum value of $k_i$ will occur at $k_0 \varLambda \approx 1$. This value can be interpreted as being associated with Bragg resonance which occurs close to $k_0 \varLambda = 1$ for periodic beds with periodicity $\varLambda$. Bragg resonance is characterised by coherent multiple reflections. In the case of varying ice $C_1^2\approx d_0^2/(h_0-d_0)^2$ which alters the magnitude of the attenuation, but not the condition $k_0\varLambda \approx 1$ for the maximum.

For $k_0 \varLambda \ll 1$, $k_i \propto k_0^2$, and for $k_0 \varLambda \gg 1$ the attenuation decays exponentially as $k_0 \varLambda$ increases, although we note this limit is outside the long-wavelength assumptions used to develop this model. The latter result holds in this long-wavelength model and contrasts with the conclusions drawn by previous researchers (see, e.g. Devillard et al. (Reference Devillard, Dunlop and Souillard1988) and Mei et al. (Reference Mei, Stiassnie and Yue2005, § 5)) who associate exponential decay in wave attenuation as a finite water depth effect.

These conclusions are based on a long-wave model of wave propagation with randomness described by a continuously varying function. For short wave scattering by floating broken ice, for example, the physics will be different as scattering by discrete ice floes will need to be modelled correctly.

5. Numerical methods and simulations

5.1. Generating a random surface

In order to numerically generate a random function, $r(x)$, with statistical properties (3.1a,b) and (3.4) characterised by the r.m.s. height $1$ and the correlation length $\varLambda$ we implement the weighted moving average method described in Sarris et al. (Reference Sarris, Haslinger, Huthwaite, Nagy and Lowe2021) and originally due to Ogilvy (Reference Ogilvy1988). The function $r(x)$ will be defined at $x = x_i = \textrm {i}\Delta x$ for $i=0,\ldots, V$ where $\Delta x = L/V$; either $\Delta x$ or $V$ can be used as the numerical parameter defining the resolution of the random surface.

We generate the Gaussian weights

(5.1)\begin{equation} w_j = W {\rm e}^{-2(j\Delta x)^2/\varLambda^2} \end{equation}

for $j=-M,\ldots,M$, where $M = \lfloor 4\varLambda /(\Delta x\sqrt {2}) \rceil$ (denoting integer part) is a truncation parameter and $W$ is defined to normalise these values so that

(5.2)\begin{equation} \sum_{j={-}M}^{M} w_j = 1. \end{equation}

Next, we define

(5.3)\begin{equation} \sigma_v^2 = 1/\sum_{j={-}M}^{M}w_j^2 ,\end{equation}

which is used to generate the $2N+1$ uncorrelated random numbers $v_i$, $-N \leq i \leq N$ from a Gaussian distribution with a variance of $\sigma _v$. The height of a random surface at $x=x_i$ is defined by

(5.4)\begin{equation} r_i = \sum_{j={-}M}^{M} w_j v_{j+i+M-N},\quad i=0,\ldots,V ,\end{equation}

requiring $N$ to be defined by $2N = V+2M$. Our theory requires that $r(x) = 0$ at $x=0$, $x=L$ and that these values are approached smoothly from within the interval $x \in (0,L)$. We thus introduce a Tukey smoothing window at either end of the interval of length $\varLambda$ (assumed to be less than $L/2$) via

(5.5) \begin{equation} r(x_i) = \left\{ \begin{array}{l} r_i, \quad V_\varLambda +1 \leq i \leq V-V_\varLambda -1,\\ r_i \left(\dfrac12 - \dfrac12 \cos \left(\dfrac{{\rm i}{\rm \pi}}{V_\varLambda}\right)\right),\quad i=0,\ldots,V_{\varLambda},\\ r_i \left(\dfrac12 - \dfrac12 \cos \left({\rm \pi} \dfrac{V-i}{V_\varLambda}\right)\right),\quad i=V-V_{\varLambda},\ldots,V, \end{array}\right. \end{equation}

where $V_\varLambda = \lfloor \varLambda /\Delta x \rceil$. Numerically, we ensure $V_\varLambda$, which represents the number of points per characteristic length of bed, is sufficiently large.

5.2. Determining decay via a transfer matrix

Simulations of scattering are performed over a region $0 < x < L$ with $L/h_0 \gg 1$. Taking $L$ to be large is done since we wish to compare our results with the theoretical results where $L= \infty$. Thus, we aim to ensure that waves pass over enough of the bed for the effect of randomness to be felt. Attenuation over longer beds can also help suppress multiple scattering effects associated with the junctions at $x=0$ and $x=L$ between constant and random surfaces. However, the method described in the following for determining attenuation is insensitive to multiple scattering effects.

Instead of (2.5), (2.6), let us momentarily express the solution in $x < 0$, $x > L$ more generally as

(5.6)\begin{align} \varOmega(x)=\begin{cases} A_{-}{\rm e}^{\mathrm{i} k_0x} + B_{-}{\rm e}^{-\mathrm{i} k_0x}, & x<0\\ A_+{\rm e}^{\mathrm{i} k_0x} + B_+{\rm e}^{-\mathrm{i} k_0x}, & x>L \end{cases} \end{align}

for complex constants $A_\pm$, $B_\pm$, representing amplitudes of right- and left-propagating waves, respectively, whilst $k_0$ satisfies (2.7).

We encode scattering using either a $2 \times 2$ scattering matrix, $\boldsymbol{\mathsf{S}}$, satisfying

(5.7)\begin{equation} \begin{pmatrix} A_+ \\ B_- \end{pmatrix}=\boldsymbol{\mathsf{S}} \begin{pmatrix} A_- \\ B_+ \end{pmatrix}, \end{equation}

which relates outgoing to incoming waves or a $2 \times 2$ transfer matrix, $\boldsymbol{\mathsf{P}}$, satisfying

(5.8)\begin{equation} \begin{pmatrix} A_+ \\ B_+ \end{pmatrix}= \boldsymbol{\mathsf{P}} \begin{pmatrix} A_- \\ B_- \end{pmatrix}, \end{equation}

which relates waves in $x > L$ to waves in $x < 0$. Energy conservation requires incoming and outgoing wave energy fluxes balance so that $|A_-|^2 + |B_+|^2 = |A_+|^2 + |B_-|^2$ and this implies $\bar{\boldsymbol{\mathsf{S}}}^\textrm {T} \boldsymbol{\mathsf{S}} = {\mathsf I}$ where $\textit {\textsf {I}}$ is the identity and the overbar denotes conjugation; $\boldsymbol{\mathsf{S}}$ is a unitary matrix. Multiplying (5.8) by $(\overline {A}_+, - \overline {B}_+)^\textrm {T}$ results in a similar identity

(5.9)\begin{equation} {\mathsf{E}} {\bar{\mathsf P}}^{\rm T} {\mathsf{E}} {\mathsf{P}} = \boldsymbol{\mathsf{I}},\quad {\mathsf{E}} = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}. \end{equation}

This is sufficient to show that if $\lambda$ is an eigenvalue of $\boldsymbol{\mathsf{P}}$, then so is $\bar {\lambda }$, as is $1/\bar {\lambda }$. The pair of eigenvalues $\lambda _\pm$ of ${\mathsf {P}}$ are therefore either both real, occurring in reciprocal pairs, or complex conjugates lying on the unit circle.

As shown in, for example, Porter & Porter (Reference Porter and Porter2003), the eigenvalues characterise wave propagation across $0 < x < L$: if $\lambda _{\pm }$ are complex conjugates, then there is no attenuation as waves travel from left to right. If, however, $\lambda _{\pm }$ are real, then writing $\lambda _+={\rm e}^{-k_i L}$ and $\lambda _-={\rm e}^{k_i L}$, say, indicate that right- and left-propagating waves are attenuated with the rate $k_i$.

Since the transfer matrix, $\boldsymbol{\mathsf{P}}$, describes the solution over $0 < x < L$ without coupling to the solution in $x < 0$ and $x > L$ its eigenvalues determine decay (or otherwise) without interference from multiple scattering effects associated with waves being reflected at the junctions $x=0$ and $x=L$.

The entries of $\boldsymbol{\mathsf{P}}$ and $\boldsymbol{\mathsf{P}}$ requires us to solve (2.4). We follow Porter (Reference Porter2019), write $x = \xi L$, $p(\xi ) = \varOmega (x) = (1+R)p_1(\xi ) + \mathrm {i} k_0 {{\hat {\hat d}}_0}(1-R)p_2(\xi )$ and $q(\xi ) = {\hat {\hat d}}(x)\varOmega '(x)= (1+R)q_1(\xi ) + \mathrm {i} k_0{{\hat{\hat d}}_0}(1-R)q_2(\xi )$, where ${{\hat{\hat d}}_0} = {\hat {\hat d}}(0)$, and numerically solve the dimensionless coupled first-order system

(5.10a,b)\begin{equation} p_i'(\xi) = (L/{\hat {\hat d}}(L \xi)) q_i(\xi),\quad q'_i(\xi) ={-}KL p_i(\xi),\quad 0 < \xi < 1 \end{equation}

for $i=1,2$ with the initial conditions $p_1(0) = 1$, $q_1(0) = 0$ and $p_2(0) = 0$ and $q_2(0) = 1$. This allows us, after matching to the solution given by (5.6) in $x < 0$ and $x > L$ and with some manipulation of the algebra, to express the solution either using (5.7) with

(5.11)\begin{align} \boldsymbol{\mathsf{S}} &= \begin{pmatrix} \mathrm{i} {{\hat{\hat d}}_0}k_0 p_2(1) - p_1(1) & {\rm e}^{\mathrm{i} k_0 L} \\ \mathrm{i} {{\hat{\hat d}}_0}k_0 q_2(1) - q_1(1) & \mathrm{i} {{\hat{\hat d}}_0}k_0 {\rm e}^{\mathrm{i} k_0 L} \end{pmatrix}^{{-}1} \nonumber\\ &\quad \times\begin{pmatrix} \mathrm{i} {{\hat{\hat d}}_0}k_0 p_2(1) + p_1(1) & {\rm e}^{-\mathrm{i} k_0 L} \\ \mathrm{i} {{\hat{\hat d}}_0}k_0 q_2(1) + q_1(1) & -\mathrm{i} {{\hat{\hat d}}_0}k_0 {\rm e}^{-\mathrm{i} k_0 L} \end{pmatrix} \end{align}

or using (5.8) with

(5.12)\begin{align} \boldsymbol{\mathsf{P}} &= \begin{pmatrix} {\rm e}^{\mathrm{i} k_0L} & {\rm e}^{-\mathrm{i} k_0L} \\ \mathrm{i} {{\hat{\hat d}}_0}k_0 {\rm e}^{\mathrm{i} k_0L} & -\mathrm{i} {{\hat{\hat d}}_0}k_0 {\rm e}^{-\mathrm{i} k_0L} \end{pmatrix}^{{-}1} \nonumber\\ &\quad \times\begin{pmatrix} \mathrm{i}{{\hat{\hat d}}_0}k_0p_2(1)+p_1(1) & -\mathrm{i} {{\hat{\hat d}}_0}k_0p_2(1)+p_1(1) \\ \mathrm{i} {{\hat{\hat d}}_0}k_0q_2(1)+q_1(1) & -\mathrm{i} {{\hat{\hat d}}_0}k_0q_2(1)+q_1(1) \end{pmatrix}. \end{align}

When we set $A_- = 1$ and $B_+=0$, $B_- = R_L$ and $A_+ = T_L$ become the reflection and transmission coefficients to due waves incident from $x < 0$ which are most easily determined from (5.7) with (5.11).

Attenuation, on the other hand, simply requires us to evaluate the pair of eigenvalues of $\boldsymbol{\mathsf{P}}$ from (5.12). The corresponding decay rate is then determined from $k_i = | \ln |\lambda _+||/L$ which, in the case of complex conjugate eigenvalues is zero.

For the ensemble averaging the results we run $N \gg 1$ simulations of different realisations of the bed or the ice thickness and then compute

(5.13a–c)\begin{equation} \langle k_i \rangle = \frac1N \sum_{n=1}^{N} k_i,\quad \langle |R_L| \rangle = \frac1N \sum_{n=1}^{N} |R_L|,\quad \langle |T_L| \rangle = \frac1N \sum_{n=1}^{N} |T_L|, \end{equation}

where the terms under the sum represent the output of each random simulation. Depending on numerical parameters used, computations of the three averages will typically take between 20 and 200 seconds on a standard desktop PC when $N=500$. A standard Runge–Kutta–Fehlberg method is used to solve (5.10a,b).