2.1. Governing equations

Consider a three-dimensional seawater domain with infinite horizontal extent and constant finite depth, $h$ say. Points in the water are located by Cartesian coordinates $(x,y)$ in the horizontal plane and $z$ in the vertical direction positively oriented upwards. We assume that the free surface at rest coincides with the plane $z=0$ , so the seabed is defined by $z=-h$ . Irregularities in the seabed are not considered here, as $h$ is assumed to be large compared with the wavelength throughout this study.

We seek to model the propagation of a directional ocean wave spectrum through an MIZ composed of thousands to tens of thousands of randomly positioned floating compliant ice floes with circular shape and uniform thickness. We only consider circular ice floes for simplicity, conjecturing that the shape of the floes has a small effect on the mean properties of wave propagation through large random arrays.

The solution method described in § 3 requires clustering of the array of floes into multiple slabs. Without loss of generality, we align the slabs with the $y$ -axis. Let $S$ denote the number of slabs, and let slab $q$ be bounded by ${\it\xi}_{q-1}\leqslant x\leqslant {\it\xi}_{q}$ , $1\leqslant q\leqslant S$ , have width $L_{q}={\it\xi}_{q}-{\it\xi}_{q-1}$ and contain $N_{q}$ floes. A sketch of the geometry is given in figure 1. It should be noted that a floe belongs to a slab if its centre is in the slab bounds.

Ice floe radii $a_{p}^{q}$ , $1\leqslant q\leqslant S$ , $1\leqslant p\leqslant N_{q}$ , are drawn from a prescribed power-law FSD as observed in the field (Toyota et al. Reference Toyota, Takatsuji and Nakayama2006). A parameterization of the FSD will be described in § 5.2. Further, we assume that all floes have constant thickness $D$ and uniform density ${\it\rho}\approx 922.5~\text{kg}~\text{m}^{-3}$ . We do not include floes with different thicknesses and densities in the model in order to limit the number of parameters, although the method we propose can accommodate these extensions.

We consider a multidirectional wave field with small amplitude compared with the wavelength and prescribed angular frequency ${\it\omega}$ . The water is approximated as an inviscid and incompressible fluid with constant density ${\it\rho}_{0}\approx 1025~\text{kg}~\text{m}^{-3}$ and irrotational flow. The linear theory of water waves can then be used to describe the water motion. Assuming time-harmonic conditions, we express the velocity field in the water domain as $(\boldsymbol{{\rm\nabla}},\partial _{z})\text{Re}\{(g/\text{i}{\it\omega}){\it\phi}(x,y,z)\text{e}^{-\text{i}{\it\omega}t}\}$ , where $\boldsymbol{{\rm\nabla}}\equiv (\partial _{x},\partial _{y})$ and $g\approx 9.81~\text{m}~\text{s}^{-2}$ is the acceleration due to gravity. The complex-valued (reduced) potential, ${\it\phi}$ , is governed by Laplace’s equation,

(2.1) $$\begin{eqnarray}{\rm\nabla}^{2}{\it\phi}+\partial _{z}^{2}{\it\phi}=0\quad (-\infty <x,y<\infty ,-h<z<-d),\end{eqnarray}$$

where $z=-d$ describes the upper boundary of the fluid domain, such that $d=0$ when a free surface is present and $d=({\it\rho}/{\it\rho}_{0})D$ (i.e. the Archimedean draught) when the surface is covered by a floe. On the impermeable seabed, we prescribe a no-normal-flow condition

(2.2) $$\begin{eqnarray}\partial _{z}{\it\phi}=0\quad (z=-h).\end{eqnarray}$$

In fluid regions bounded above by a free surface (i.e. $d=0$ ), the potential satisfies the boundary condition

(2.3) $$\begin{eqnarray}\partial _{z}{\it\phi}={\it\alpha}{\it\phi}\quad (z=0),\end{eqnarray}$$

where ${\it\alpha}={\it\omega}^{2}/g$ is a frequency parameter.

Figure 1.

Schematic of the geometry in the horizontal plane $z=0$ .

We prescribe an ambient incident wave field ${\it\phi}^{\mathit{In}}$ travelling in the positive $x$ -direction and defined by a superposition of plane waves with amplitudes that depend continuously on the angle of incidence ${\it\tau}$ (with respect to the $x$ -axis). We express it as

(2.4) $$\begin{eqnarray}{\it\phi}^{\mathit{In}}(x,y,z)={\it\zeta}_{0}(z)\int _{-{\rm\pi}/2}^{{\rm\pi}/2}A^{\mathit{In}}({\it\tau})\text{e}^{\text{i}k_{0}((x-{\it\xi}_{0})\cos {\it\tau}+y\sin {\it\tau})}\,\text{d}{\it\tau},\end{eqnarray}$$

where the incident wave directional spectrum $A^{\mathit{In}}({\it\tau})$ characterizes the angular distribution of ambient wave amplitude at $x={\it\xi}_{0}$ . The function ${\it\zeta}_{0}(z)=\cosh k_{0}(z+h)/\cosh k_{0}h$ describes the vertical motion of the incident wave field. The quantity $k_{0}$ denotes the propagating wavenumber for a wave travelling in the free-surface region and will be defined shortly. Scattering of the ambient wave field by the array of floes gives rise to reflected and transmitted wave components which are expressed as

(2.5a ) $$\begin{eqnarray}{\it\phi}(x,y,z)\approx {\it\phi}^{\mathit{In}}(x,y,z)+{\it\zeta}_{0}(z)\int _{-{\rm\pi}/2}^{{\rm\pi}/2}A^{\mathit{R}}({\it\tau})\text{e}^{\text{i}k_{0}(-(x-{\it\xi}_{0})\cos {\it\tau}+y\sin {\it\tau})}\,\text{d}{\it\tau}\end{eqnarray}$$

as $x\rightarrow -\infty$ and

(2.5b ) $$\begin{eqnarray}{\it\phi}(x,y,z)\approx {\it\zeta}_{0}(z)\int _{-{\rm\pi}/2}^{{\rm\pi}/2}A^{\mathit{T}}({\it\tau})\text{e}^{\text{i}k_{0}((x-{\it\xi}_{S})\cos {\it\tau}+y\sin {\it\tau})}\,\text{d}{\it\tau}\end{eqnarray}$$

as $x\rightarrow \infty$ . The reflected and transmitted wave directional spectra, $A^{\mathit{R}}({\it\tau})$ and $A^{\mathit{T}}({\it\tau})$ , are unknowns of the problem. They characterize the angular distribution of the reflected and transmitted amplitudes at $x={\it\xi}_{0}$ and $x={\it\xi}_{S}$ respectively.

We model the motion experienced by the ice floes using the Kirchhoff–Love theory of thin elastic plates, which assumes that the thickness is small compared with the diameter and vertical deformations are small relative to the thickness. At the water–floe interface the potential then satisfies

(2.6) $$\begin{eqnarray}({\it\beta}\boldsymbol{{\rm\nabla}}^{4}+1-{\it\alpha}d)\partial _{z}{\it\phi}={\it\alpha}{\it\phi}\quad (z=-d),\end{eqnarray}$$

for a floe with thickness $D$ and draught $d=({\it\rho}/{\it\rho}_{0})D$ . The stiffness parameter ${\it\beta}=F/{\it\rho}_{0}g$ is defined in terms of the flexural rigidity of the plate, $F=ED^{3}/12(1-{\it\nu}^{2})$ , where $E\approx 6~\text{G}~\text{Pa}$ is a typical value for the effective Young’s modulus of sea ice (Mellor Reference Mellor and Untersteiner1986) and ${\it\nu}\approx 0.3$ denotes Poisson’s ratio.

We complete the description of the ice floe model by imposing free edge conditions. These are most conveniently expressed using polar coordinates $(r,{\it\theta})$ with the origin at the centre of the floe. For a floe of radius $a$ , we have

(2.7a ) $$\begin{eqnarray}[r^{2}{\rm\nabla}_{r,{\it\theta}}^{2}-(1-{\it\nu})(r\partial _{r}+\partial _{{\it\theta}}^{2})]\partial _{z}{\it\phi}=0\quad (r=a)\end{eqnarray}$$

and

(2.7b ) $$\begin{eqnarray}[r^{3}\partial _{r}{\rm\nabla}_{r,{\it\theta}}^{2}+(1-{\it\nu})(r\partial _{r}-1)\partial _{{\it\theta}}^{2}]\partial _{z}{\it\phi}=0\quad (r=a),\end{eqnarray}$$

where $\boldsymbol{{\rm\nabla}}_{r,{\it\theta}}\equiv (\partial _{r}+1/r,(1/r)\partial _{{\it\theta}})$ . In addition, we assume that the floes do not respond in surge and sway, so that

(2.7c ) $$\begin{eqnarray}\partial _{r}{\it\phi}=0\quad (r=a,\,-d<z<0).\end{eqnarray}$$

2.2. Scattering by a single floe

Each floe scatters the local wave field incident on it, which is the combination of the ambient incident wave field and the wave fields scattered by all other floes. For a given floe, which, as above, is assigned the polar coordinate system $(r,{\it\theta})$ , these local incident and scattered wave potentials are expressed as the truncated eigenfunction expansions

(2.8a ) $$\begin{eqnarray}{\it\phi}_{I}(r,{\it\theta},z)\approx \mathop{\sum }_{m=0}^{M_{1}}{\it\zeta}_{m}(z)\mathop{\sum }_{n=-N}^{N}a_{m,n}\text{J}_{n}(k_{m}r)\text{e}^{\text{i}n{\it\theta}}\quad (r>a)\end{eqnarray}$$

and

(2.8b ) $$\begin{eqnarray}{\it\phi}_{S}(r,{\it\theta},z)\approx \mathop{\sum }_{m=0}^{M_{1}}{\it\zeta}_{m}(z)\mathop{\sum }_{n=-N}^{N}b_{m,n}\text{H}_{n}(k_{m}r)\text{e}^{\text{i}n{\it\theta}}\quad (r>a)\end{eqnarray}$$

respectively, which are solutions to (2.1)–(2.3) in cylindrical coordinates (see, e.g. Peter, Meylan & Chung Reference Peter, Meylan and Chung2003, for a detailed derivation). We have introduced $\text{J}_{n}$ and $\text{H}_{n}$ to denote the Bessel and Hankel functions of the first kind of order $n$ respectively. The amplitudes $a_{m,n}$ and $b_{m,n}$ , $0\leqslant m\leqslant M_{1}$ , $-N\leqslant n\leqslant N$ , are unknowns of the scattering problem. We define the vertical modes as ${\it\zeta}_{m}(z)=\cosh k_{m}(z+h)/\!\cosh k_{m}h$ . The wavenumbers $k_{m}$ , $m\geqslant 0$ , are the solutions $k$ of the dispersion relation for an open water region, that is

(2.9) $$\begin{eqnarray}k\tanh kh={\it\alpha}.\end{eqnarray}$$

We denote the positive real root of (2.9) by $k_{0}$ . It is associated with a wave mode travelling in the horizontal plane. As the water depth $h$ is assumed to be large compared with the wavelength, we have $k_{0}\approx {\it\alpha}$ , so that the wavelength is approximately independent of $h$ . All other $k_{m}$ , $m\geqslant 1$ , are purely imaginary with positive imaginary part and are ordered such that $-\text{i}k_{m}<-\text{i}k_{m+1}$ . They are associated with evanescent vertical wave modes which decay exponentially in the horizontal directions, such that the rate of decay increases for increasing $m$ . In contrast to $k_{0}$ , the values of $k_{m}$ , $m\geqslant 1$ , depend on $h$ . It should be noted that the sums in (2.8) are truncated versions of the corresponding series expansions, with $M_{1}$ and $N$ chosen in order to obtain a sufficient degree of accuracy (see below).

Montiel et al. (Reference Montiel, Bennetts, Squire, Bonnefoy and Ferrant2013) proposed a solution method for the single-floe scattering problem. Using equations (2.1), (2.2), (2.6), (2.7a ) and (2.7b ) they expressed the potential below the floe (i.e. $r<a$ ) as a truncated series of eigenfunctions similar to (2.8a ), but with different wavenumbers characterizing wave modes in the ice-covered water domain. They then used a version of the eigenfunction matching method (EMM) which accommodates the draught of the floe, through (2.7c ), in order to extend the zero-draught EMM of Peter et al. (Reference Peter, Meylan and Chung2003). The EMM produces a mapping between the amplitudes $a_{m,n}$ and $b_{m,n}$ for each angular mode $n$ ; the axisymmetry of the problem decouples the angular modes. In matrix form, the mapping is expressed as

(2.10) $$\begin{eqnarray}\boldsymbol{b}_{n}=\unicode[STIX]{x1D64E}_{n}\boldsymbol{a}_{n},\end{eqnarray}$$

for $-N\leqslant n\leqslant N$ , where $\boldsymbol{a}_{n}$ and $\boldsymbol{b}_{n}$ are column vectors of size $M_{1}+1$ containing the amplitudes $a_{m,n}$ and $b_{m,n}$ respectively for a given $n$ . The matrices $\unicode[STIX]{x1D64E}_{n}$ are squares of size $M_{1}+1$ and form the diffraction transfer matrix (DTM) of the floe when concatenated in a block-diagonal matrix. Extensions of our model to arbitrarily shaped floes could be accommodated following Peter & Meylan (Reference Peter and Meylan2004), who devised a numerical method to compute the DTM of such floes. It should be noted that in this case the DTM loses its block-diagonal property.

We choose the truncation limits $N$ and $M_{1}$ to achieve three-digit accuracy for the scattered energy. In the regime of interest here, $k_{0}a=O(1)$ , this typically requires $M_{1}=O(100)$ vertical modes and $N=O(1)$ angular modes, as demonstrated by Montiel (Reference Montiel2012).

2.3. Multiple scattering and limitations of the direct approach

Our goal is to solve the wave scattering problem deterministically for a large number of floes. Apart from truncations, no approximation will be made with regards to multiple scattering by the array, so that the scattered field due to each floe acts as an incident field on all of the other floes.

The so-called self-consistent approach (or direct matrix method) provides an exact representation of multiple scattering processes and was introduced in the context of ocean wave interactions with floating structures by Kagemoto & Yue (Reference Kagemoto and Yue1986), noting that this approach is standard in many areas concerned with wave scattering by arrays of scatterers (see, e.g., Martin Reference Martin2006). The method is briefly summarized below, and the reader is referred to the investigations cited here for additional details (see also Peter & Meylan (Reference Peter and Meylan2004) in the context of wave interactions with elastic ice floes).

The method describes the wave forcing ${\it\phi}_{I}^{(p)}$ incident upon a floe $p$ , which has its centre located at $(x,y)=(x_{p},y_{p})$ , as the coherent sum of the ambient incident wave ${\it\phi}^{\mathit{In}}$ and the scattered wave fields originating from all of the other floes, i.e.

(2.11) $$\begin{eqnarray}{\it\phi}_{I}^{(p)}={\it\phi}^{\mathit{In}}+\mathop{\sum }_{j,\,j\neq p}{\it\phi}_{S}^{(j)},\end{eqnarray}$$

where ${\it\phi}_{S}^{(j)}$ is the scattered wave potential due to a floe $j$ , with centre located at $(x,y)=(x_{j},y_{j})$ , and the sum over $j$ runs for all floes in the array except $p$ . The incident and scattered wave components are expressed in terms of the eigenfunction expansions (2.8), using the relevant local polar coordinates $(r_{p},{\it\theta}_{p})$ , defined by $(x,y)=(x_{p}+r_{p}\cos {\it\theta}_{p},y_{p}+r_{p}\sin {\it\theta}_{p})$ . Subsequently, the expressions for the scattered waves ${\it\phi}_{S}^{(j)}$ in (2.11) are mapped into the local coordinates $(r_{p},{\it\theta}_{p})$ of floe $p$ . Application of the (reduced) boundary condition (2.10) around floe $p$ then yields the following matrix equation:

(2.12) $$\begin{eqnarray}\boldsymbol{b}_{n}^{(p)}-\mathop{\sum }_{j,\,j\neq p}\mathop{\sum }_{s=-N}^{N}\unicode[STIX]{x1D64E}_{n}^{(p)}\unicode[STIX]{x1D64F}_{n,s}^{(j,p)}\boldsymbol{b}_{s}^{(j)}=\unicode[STIX]{x1D64E}_{n}^{(p)}\boldsymbol{f}_{n}^{(p)},\end{eqnarray}$$

which can be obtained for all floes $p$ . Numerical experiments (see § 4) have shown that only $O(1)$ vertical modes are necessary in (2.8) to resolve wave interactions in an array of floes accurately, assuming that the single-floe solutions are obtained with sufficient accuracy (which requires $M_{1}=O(100)$ vertical modes). The vectors of scattered wave amplitudes $\boldsymbol{b}_{n}^{(p)}$ and matrices $\unicode[STIX]{x1D64E}_{n}^{(p)}$ in (2.12) are then chosen to have size $M_{2}+1$ , where $M_{2}=O(1)$ . They are obtained by truncating the corresponding quantities of size $M_{1}+1$ defined for the single-floe solution.

The resulting system of equations is solved for the scattered wave amplitudes $\boldsymbol{b}_{n}^{(p)}$ for all floes $p$ and angular modes $n$ . Here, the forcing vectors $\boldsymbol{f}_{n}^{(p)}$ , $-N\leqslant n\leqslant N$ , contain the amplitudes of the ambient incident potential (2.4) expressed in the local cylindrical coordinates of floe $p$ . The square matrices $\unicode[STIX]{x1D64E}_{n}^{(p)}$ are analogous to the matrices defined in (2.10) for each floe $p$ . We have also introduced the diagonal matrices $\unicode[STIX]{x1D64F}_{n,s}^{(j,p)}=\text{diag}\{\text{H}_{s-n}(k_{m}\mathscr{R}_{j,p})\text{e}^{\text{i}(s-n){\it\varpi}_{j,p}},0\leqslant m\leqslant M_{2}\}$ of size $M_{2}+1$ , where $(\mathscr{R}_{j,p},{\it\varpi}_{j,p})$ are the polar coordinates of the centre of floe $p$ in the local system associated with floe $j$ . These matrices describe the change of local polar coordinates from floe $j$ to floe $p$ , and their entries are calculated using Graf’s addition theorem (Abramowitz & Stegun Reference Abramowitz and Stegun1970), which couples the angular modes.

The size of system (2.12) grows linearly with the number of floes in the array, and direct inversion will lead to a computational cost increasing with the cubic power of the number of floes. The order-of-scattering method, based on the original paper by Twersky (Reference Twersky1952), has been used to approximate the solution of (2.12) by successive orders of multiple scattering events (see, e.g., Ohkusu Reference Ohkusu, Bishop and Price1974; Mavrakos & Koumoutsakos Reference Mavrakos and Koumoutsakos1987). Mathematically, this is equivalent to solving (2.12) using an iterative scheme, e.g. the Jacobi or Gauss–Seidel method. This method usually leads to performance improvements, although the computational cost is strongly affected by the concentration of floating bodies in the array (Kagemoto & Yue Reference Kagemoto and Yue1986). Numerical experiments conducted by the authors have shown that direct inversion or iterative approaches are limited to arrays of $O(100)$ floes.

3.1. Reflection and transmission by a single slab

Montiel, Squire & Bennetts (Reference Montiel, Squire and Bennetts2015b ) derived a set of relationships (for the special case ${\it\gamma}=0$ ) between the incident and scattered amplitude functions $A_{0;q}^{\pm }({\it\chi})$ on either side of a given slab $q$ as a result of reflection and transmission. Here, we extend these relationships to include evanescent vertical modes (i.e. $m\geqslant 1$ ) and $x$ -decaying wave components (i.e. ${\it\gamma}>0$ ).

We consider the scattering by slab $q$ due to the incident forcing from its left-hand side, ${\it\phi}_{q-1}^{(+)}$ , only. The response to incident forcing from its right-hand side, ${\it\phi}_{q}^{(-)}$ , follows similarly, and the total response to forcing from both sides is calculated via superposition.

The forcing field may be expressed in the local polar coordinates of floe $p$ analogously to (2.8a ). For each angular mode $n$ and vertical mode $m$ , the forcing amplitudes are then given by

(3.4) $$\begin{eqnarray}[\boldsymbol{f}_{n}^{(p)}]_{m}=\text{i}^{n}\int _{{\it\Gamma}_{m}}A_{m;q-1}^{+}({\it\chi})\text{e}^{-\text{i}n{\it\chi}}\text{e}^{\text{i}k_{m}((x_{p}-{\it\xi}_{q-1})\cos {\it\chi}+y_{p}\sin {\it\chi})}\,\text{d}{\it\chi}.\end{eqnarray}$$

The self-consistent method outlined in § 2.3 is used to solve the multiple scattering problem within the slab. It results in a system of equations analogous to (2.12) which yields the mapping

(3.5) $$\begin{eqnarray}\boldsymbol{b}^{[q]}=\unicode[STIX]{x1D63F}^{[q]}\boldsymbol{f}^{[q]},\end{eqnarray}$$

where $\boldsymbol{b}^{[q]}$ ( $\boldsymbol{f}^{[q]}$ ) is a vector of length $(M_{2}+1)(2N+1)N_{q}$ that contains all of the scattered (incident) wave amplitudes contained in $\boldsymbol{b}_{n}^{(p)}$ ( $\boldsymbol{f}_{n}^{(p)}$ ), for $-N\leqslant n\leqslant N$ and $1\leqslant p\leqslant N_{q}$ . The square matrix $\unicode[STIX]{x1D63F}^{[q]}$ has size $(M_{2}+1)(2N+1)N_{q}$ and is the DTM of slab $q$ .

We seek the reflected and transmitted amplitude functions in the form

(3.6a ) $$\begin{eqnarray}A_{m,q-1}^{-}({\it\chi})=\mathop{\sum }_{l=0}^{M_{2}}\int _{{\it\Gamma}_{l}}\mathscr{R}_{m,l;q}^{+}({\it\chi}\,:\,{\it\tau})A_{l;q-1}^{+}({\it\tau})\,\text{d}{\it\tau}\end{eqnarray}$$

and

(3.6b ) $$\begin{eqnarray}A_{m,q}^{+}({\it\chi})=\mathop{\sum }_{l=0}^{M_{2}}\int _{{\it\Gamma}_{l}}\mathscr{T}_{m,l;q}^{+}({\it\chi}\,:\,{\it\tau})A_{l;q-1}^{+}({\it\tau})\,\text{d}{\it\tau}\end{eqnarray}$$

respectively. The functions $\mathscr{R}_{m,l;q}^{+}({\it\chi}\,:\,{\it\tau})$ and $\mathscr{T}_{m,l;q}^{+}({\it\chi}\,:\,{\it\tau})$ are respectively reflection and transmission kernels for each pair of vertical modes $m$ and $l$ .

In order to evaluate the reflection and transmission kernels, we re-express the scattered wave field (2.8b ) due to each floe in the slab as a superposition of plane waves. This is achieved using the following plane wave representation of the outgoing cylindrical harmonics:

(3.7) $$\begin{eqnarray}\text{H}_{n}(k_{m}r)\text{e}^{\text{i}n{\it\theta}}=\left\{\begin{array}{@{}ll@{}}\displaystyle \frac{(-\text{i})^{n}}{{\rm\pi}}\int _{{\it\Gamma}_{m}}\text{e}^{\text{i}n{\it\chi}}\text{e}^{\text{i}k_{m}(x\cos {\it\chi}+y\sin {\it\chi})}\,\text{d}{\it\chi} & (x\geqslant 0),\\ \displaystyle \frac{\text{i}^{n}}{{\rm\pi}}\int _{{\it\Gamma}_{m}}\text{e}^{-\text{i}n{\it\chi}}\text{e}^{\text{i}k_{m}(-x\cos {\it\chi}+y\sin {\it\chi})}\,\text{d}{\it\chi} & (x\leqslant 0),\end{array}\right.\end{eqnarray}$$

where $(x,y)=(r\cos {\it\theta},r\sin {\it\theta})$ . This identity is derived from Sommerfeld’s integral representation of the Hankel function (Sommerfeld Reference Sommerfeld1949) for $m=0$ , while an integral representation of the modified Bessel functions of the second kind $\text{K}_{n}$ has been used for $m\geqslant 1$ (see Linton & Evans Reference Linton and Evans1992, equation (2.12)). To the best of the authors’ knowledge, the two cases have not been unified in this manner before. By substituting (3.7) into (2.8b ) for all floes and modes simultaneously, and using (3.5) and (3.4) in turn, we obtain the following semianalytical expressions for the reflection and transmission kernels after algebraic manipulations:

(3.8a ) $$\begin{eqnarray}\mathscr{R}_{m,l;q}^{+}({\it\chi}\,:\,{\it\tau})=(\boldsymbol{v}_{m;q}^{\mathit{S}-}({\it\chi}))^{\mathit{tr}}\unicode[STIX]{x1D64B}\unicode[STIX]{x1D63F}^{[q]}\unicode[STIX]{x1D64B}^{-1}\boldsymbol{v}_{l;q}^{\mathit{In}+}({\it\tau})\end{eqnarray}$$

and

(3.8b ) $$\begin{eqnarray}\mathscr{T}_{m,l;q}^{+}({\it\chi}\,:\,{\it\tau})=(\boldsymbol{v}_{m;q}^{\mathit{S}+}({\it\chi}))^{\mathit{tr}}\unicode[STIX]{x1D64B}\unicode[STIX]{x1D63F}^{[q]}\unicode[STIX]{x1D64B}^{-1}\boldsymbol{v}_{l;q}^{\mathit{In}+}({\it\tau})+\text{e}^{\text{i}k_{m}L_{q}\cos {\it\chi}}{\it\delta}({\it\chi}-{\it\tau}),\end{eqnarray}$$

where a superscript $tr$ indicates transpose and ${\it\delta}(\cdot )$ denotes the Dirac delta.

In these expressions we have defined the vector $\boldsymbol{v}_{l;q-1}^{\mathit{In}+}({\it\tau})$ of length $N_{q}(2N+1)$ $(M_{2}+1)$ , which provides a change from Cartesian to polar coordinates of the forcing field in the local system of each floe in the slab. Its entries are

(3.9a ) $$\begin{eqnarray}[\boldsymbol{v}_{m;q}^{\mathit{In}+}({\it\tau})]_{\mathit{ind}(l,p,n)}=\text{i}^{n}\text{e}^{-\text{i}n{\it\tau}}\text{e}^{\text{i}k_{l}((x_{p}-{\it\xi}_{q-1})\cos {\it\tau}+y_{p}\sin {\it\tau})}{\it\delta}_{m\,l},\end{eqnarray}$$

where $\mathit{ind}(l,p,n)=lN_{q}(2N+1)+(p-1)(2N+1)+N+n+1$ defines the modal hierarchy (vertical mode, floe number, angular mode) in ordering entries, and ${\it\delta}_{m\,l}$ is the Kronecker delta. In contrast, the vectors $\boldsymbol{v}_{m;q-1}^{\mathit{S}-}({\it\chi})$ and $\boldsymbol{v}_{m;q}^{\mathit{S}+}({\it\chi})$ provide a change from polar to Cartesian coordinates of the scattered field due to each floe, travelling in the leftward and rightward direction respectively. Their entries are

(3.9b ) $$\begin{eqnarray}[\boldsymbol{v}_{m;q}^{\mathit{S}-}({\it\chi})]_{\mathit{ind}(l,p,n)}=\frac{\text{i}^{n}}{{\rm\pi}}\text{e}^{-\text{i}n{\it\chi}}\text{e}^{\text{i}k_{l}((x_{p}-{\it\xi}_{q-1})\cos {\it\chi}-y_{p}\sin {\it\chi})}{\it\delta}_{m\,l}\end{eqnarray}$$

and

(3.9c ) $$\begin{eqnarray}[\boldsymbol{v}_{m;q}^{\mathit{S}+}({\it\chi})]_{\mathit{ind}(l,p,n)}=\frac{(-\text{i})^{n}}{{\rm\pi}}\text{e}^{\text{i}n{\it\chi}}\text{e}^{-\text{i}k_{l}((x_{p}-{\it\xi}_{q})\cos {\it\chi}+y_{p}\sin {\it\chi})}{\it\delta}_{m\,l}.\end{eqnarray}$$

The matrix $\unicode[STIX]{x1D64B}$ is a permutation matrix of size $N_{q}(2N+1)(M_{2}+1)$ used to change the modal hierarchy from that used in the slab DTM $\unicode[STIX]{x1D63F}^{[q]}$ (i.e. floe number, angular mode, vertical mode) to that used in (3.9). It should also be noted that the second term in (3.8b ) represents the contribution from the forcing field to the transmitted field.

3.2. Multiple slabs

At the boundary $x={\it\xi}_{q}$ , the left- and right-travelling amplitude functions take the form

(3.10a ) $$\begin{eqnarray}A_{m,q}^{-}({\it\chi})=\mathop{\sum }_{l=0}^{M_{2}}\int _{{\it\Gamma}_{l}}(\mathscr{R}_{m,l;q+1}^{+}({\it\chi}\,:\,{\it\tau})A_{l;q}^{+}({\it\tau})+\mathscr{T}_{m,l;q+1}^{-}({\it\chi}\,:\,{\it\tau})A_{l;q+1}^{-}({\it\tau}))\,\text{d}{\it\tau}\end{eqnarray}$$

and

(3.10b ) $$\begin{eqnarray}A_{m,q}^{+}({\it\chi})=\mathop{\sum }_{l=0}^{M_{2}}\int _{{\it\Gamma}_{l}}(\mathscr{T}_{m,l;q}^{+}({\it\chi}\,:\,{\it\tau})A_{l;q-1}^{+}({\it\tau})+\mathscr{R}_{m,l;q}^{-}({\it\chi}\,:\,{\it\tau})A_{l;q}^{-}({\it\tau}))\,\text{d}{\it\tau},\end{eqnarray}$$

where the kernel functions $\mathscr{R}_{m,l;q}^{\pm }$ and $\mathscr{T}_{m,l;q}^{\pm }$ are found using the method described in § 3.1. Given that the forcing is provided by (2.4) only, we have $A_{0;0}^{+}({\it\chi})=A^{\mathit{In}}({\it\chi})$ , $A_{0;m}^{+}({\it\chi})=0$ for $1\leqslant m\leqslant M_{2}$ and $A_{S;m}^{-}({\it\chi})=0$ for $0\leqslant m\leqslant M_{2}$ .

A numerical scheme is implemented by discretizing the amplitude and kernel functions using a uniform sampling of the angular parameters ${\it\chi}$ and ${\it\tau}$ . Amplitudes and kernels associated with a travelling vertical mode ( $m=0$ ) are defined on the truncated contour $\widetilde{{\it\Gamma}}_{0}$ introduced at the beginning of § 3 (see figure 2 b). Contour $\widetilde{{\it\Gamma}}_{0}$ is then discretized by selecting $2N_{\mathit{tr}}+1$ samples ${\it\chi}_{i}$ , $-N_{\mathit{tr}}\leqslant i\leqslant N_{\mathit{tr}}$ . Different resolutions are taken for the subinterval $[-{\rm\pi}/2,{\rm\pi}/2]$ of $\widetilde{{\it\Gamma}}_{0}$ and its complex branches. Typically, the resolution chosen for the complex branches is five times as coarse as that of $[-{\rm\pi}/2,{\rm\pi}/2]$ . Likewise, contour $\widetilde{{\it\Gamma}}_{\mathit{e}}$ is discretized using $2N_{\mathit{ev}}+1$ samples.

Combining all the vertical modes, we can define vector versions of the amplitude functions $\boldsymbol{A}_{q}^{\pm }$ containing the value of the corresponding continuous functions $A_{m,q}^{\pm }({\it\chi}_{i})$ at all angular samples ${\it\chi}_{i}$ . Likewise, we obtain matrix versions of the kernel functions after discretization. At this point, a numerical quadrature must be chosen to approximate the integrals involved in (3.10). Although high-order schemes, e.g. Simpson’s rule or Gaussian quadrature, are very accurate for relatively smooth functions, we expect our spectra to be noisy for large random arrays of floes (confirmed by numerical experiments), in which case lower-order quadratures provide more accurate and more efficient estimates. We found that a composite trapezoidal rule gave the best results in terms of convergence. The weighting factors of the trapezoidal rule can be assembled in a diagonal matrix of size $N_{\mathit{ang}}=2N_{\mathit{tr}}+1+M_{2}(2N_{\mathit{ev}}+1)$ , which multiplies the matrix versions of the kernel functions to give the reflection and transmission matrices $\unicode[STIX]{x1D64D}_{q}^{\pm }$ and $\unicode[STIX]{x1D64F}_{q}^{\pm }$ . Consequently, (3.10) is written in the discretized form

(3.11a,b ) $$\begin{eqnarray}\boldsymbol{A}_{q}^{-}=\unicode[STIX]{x1D64D}_{q+1}^{+}\boldsymbol{A}_{q}^{+}+\unicode[STIX]{x1D64F}_{q+1}^{-}\boldsymbol{A}_{q+1}^{-}\quad \text{and}\quad \boldsymbol{A}_{q}^{+}=\unicode[STIX]{x1D64F}_{q}^{+}\boldsymbol{A}_{q-1}^{+}+\unicode[STIX]{x1D64D}_{q}^{-}\boldsymbol{A}_{q}^{-}.\end{eqnarray}$$

The solution to the slab interaction problem is obtained using an efficient iterative technique, which is an extension to the one described by Montiel et al. (Reference Montiel, Squire and Bennetts2015a ) for travelling modes only. Following this approach, at each slab boundary $x={\it\xi}_{q}$ , the unknown amplitude vectors are given by

(3.12a ) $$\begin{eqnarray}\boldsymbol{A}_{q}^{+}=\left(\unicode[STIX]{x1D644}-\unicode[STIX]{x1D64D}_{1,q}^{-}\unicode[STIX]{x1D64D}_{q+1,M}^{+}\right)^{-1}\unicode[STIX]{x1D64F}_{1,q}^{+}\boldsymbol{A}_{0}^{+}\end{eqnarray}$$

and

(3.12b ) $$\begin{eqnarray}\boldsymbol{A}_{q}^{-}=\left(\unicode[STIX]{x1D644}-\unicode[STIX]{x1D64D}_{q+1,M}^{+}\unicode[STIX]{x1D64D}_{1,q}^{-}\right)^{-1}\unicode[STIX]{x1D64D}_{q+1,M}^{+}\unicode[STIX]{x1D64F}_{1,q}^{+}\boldsymbol{A}_{0}^{+},\end{eqnarray}$$

where $\unicode[STIX]{x1D644}$ denotes the identity matrix of order $N_{\mathit{ang}}$ , and $\unicode[STIX]{x1D64D}_{p,q}^{\pm }$ and $\unicode[STIX]{x1D64F}_{p,q}^{\pm }$ are the reflection and transmission matrices respectively for the stack of slabs $p$ to $q$ . These reflection and transmission matrices are computed iteratively, starting from slab 1 alone, which initializes the procedure, to the stack of slabs 1 to $S$ , adding one slab to the stack at each iteration.

We establish a convergence criterion for the numerical integration scheme considered here, based on the energy conservation relation

(3.13) $$\begin{eqnarray}\int _{-{\rm\pi}/2}^{{\rm\pi}/2}|A^{\mathit{R}}({\it\chi})|^{2}\,\text{d}{\it\chi}+\int _{-{\rm\pi}/2}^{{\rm\pi}/2}|A^{\mathit{T}}({\it\chi})|^{2}\,\text{d}{\it\chi}=\int _{-{\rm\pi}/2}^{{\rm\pi}/2}|A^{\mathit{In}}({\it\chi})|^{2}\,\text{d}{\it\chi}\end{eqnarray}$$

being satisfied within a tolerance of $10^{-4}$ , where the reflected and transmitted amplitude functions are given by $A^{\mathit{R}}({\it\chi})=A_{0;0}^{-}({\it\chi})$ and $A^{\mathit{T}}({\it\chi})=A_{0;S}^{+}({\it\chi})$ respectively for $-{\rm\pi}/2\leqslant {\it\chi}\leqslant {\rm\pi}/2$ . It should be noted that we have restricted the integration domain to include travelling wave components only, as they are the only ones to affect the energy balance of the system. We refer to the restriction of the amplitude functions $A_{0;q}^{\pm }({\it\chi})$ , $q=0,\ldots ,S$ , to $-{\rm\pi}/2\leqslant {\it\chi}\leqslant {\rm\pi}/2$ as directional spectra. In particular, $A^{\mathit{R}}({\it\chi})$ and $A^{\mathit{T}}({\it\chi})$ are the reflected and transmitted directional spectra respectively.

We further note that although the energy conservation relation (3.13) is necessary for convergence of our numerical method, it is not sufficient to obtain convergence to the desired solution. In particular, we can always obtain energy conservation for sufficiently high $N_{\mathit{ang}}$ , regardless of the value taken for the truncation parameter ${\it\gamma}$ and the resolution of the complex branches in our numerical approximation of the integration contours $\widetilde{{\it\Gamma}}_{\mathit{e}}$ and $\widetilde{{\it\Gamma}}_{0}$ . Convergence of the solution within the desired tolerance is obtained for a sufficiently large ${\it\gamma}$ , as will be shown in § 4.

5.1. Random sea states

We model the forcing wave field as a random directional sea state, in which wave components travelling at different angles do not interfere coherently, i.e. their phases are uncorrelated. Numerical tests (not shown here) indicate that the reflection properties of an array of floes depend strongly on the directional coherence of the forcing field. Specifically, the reflection coefficient (and therefore the attenuation rate of wave energy) computed for a coherent forcing field is typically higher than that obtained for an incoherent field with the same prescribed incident spectrum. In addition to incoherence, we require the simulated directional sea state to be ergodic in the sense that the wave statistics over the spatial domain are uniform and can be deduced from its properties at a single point in the domain. This property is needed to reduce the variability of the response of the system to a random forcing, as will be discussed below.

A number of methods exist to simulate ergodic directional sea states deterministically (see, e.g., Jefferys Reference Jefferys1987; Miles & Funke Reference Miles and Funke1989). Most methods are based on multiple frequency wave spectra, for which a realization of the random sea state is generated by a double sum over the directional and frequency range of plane waves with random characteristics. As we assume a monochromatic wave forcing, the double sum is replaced by a single sum over $N_{\mathit{tr}}$ directions, giving a free-surface displacement

(5.1) $$\begin{eqnarray}{\it\eta}^{\mathit{In}}(x,y,t)=\mathop{\sum }_{i=1}^{N_{\mathit{tr}}}a_{i}\cos (k_{0}(x\cos {\it\tau}_{i}+y\sin {\it\tau}_{i})-{\it\omega}t+{\it\varepsilon}_{i}).\end{eqnarray}$$

For each wave component travelling at angle ${\it\tau}_{i}$ with respect to the $x$ -axis, the amplitude $a_{i}$ is defined deterministically from an energy spreading function $D({\it\tau})$ , and the phase ${\it\varepsilon}_{i}$ is a random parameter with uniform distribution between $0$ and $2{\rm\pi}$ . For the remainder of the investigation we prescribe a standard cosine-squared energy spreading function, i.e.

(5.2) $$\begin{eqnarray}D({\it\tau})=\frac{2}{{\rm\pi}}\cos ^{2}({\it\tau})\quad (-{\rm\pi}/2\leqslant {\it\tau}\leqslant {\rm\pi}/2),\end{eqnarray}$$

where the constant factor is chosen so that the total energy of the wave field is

(5.3) $$\begin{eqnarray}\int _{-{\rm\pi}/2}^{{\rm\pi}/2}D({\it\tau})\,\text{d}{\it\tau}=1.\end{eqnarray}$$

The energy spreading function (5.2) and the amplitudes of (5.1) are related by $a_{i}=\sqrt{2D({\it\tau}_{i}){\rm\Delta}{\it\tau}}$ for $1\leqslant i\leqslant N_{\mathit{tr}}$ in (5.1), with ${\rm\Delta}{\it\tau}={\rm\pi}/(N_{\mathit{tr}}-1)$ and ${\it\tau}_{i}=(i-1){\rm\Delta}{\it\tau}-{\rm\pi}/2$ .

Although directional incoherence is directly satisfied in (5.1) through the random parameter ${\it\varepsilon}_{i}$ , following Jefferys (Reference Jefferys1987) we show in appendix A that this equation does not simulate an ergodic field, i.e. much spatial variability exists in the mean energy of a generated sea state ( ${\approx}100\,\%$ relative standard error). A simple remedy is to perform averaging over an ensemble of random realizations of the wave field. The method converges slowly, however, and approximately 10 000 realizations of the sea state are necessary to approximate ergodicity with 1 % relative standard error on the incident field mean energy.

We compute the corresponding reflected energy, $R^{2}$ , for the 20-slab grating considered in § 4, with spacing constant ${\it\sigma}=1.5$ . The forcing is defined by the random sea state (5.1) with period $T=9~\text{s}$ and an energy spreading function given by (5.2). Using an ensemble of 1000 realizations of the random forcing field, we calculated the relative standard error on a single estimate of the reflected energy to be approximately $7\,\%$ . Averaging over many realizations of the random wave forcing, the relative error of the average estimate then drops as the inverse square root of the number of realizations, so it is approximately $1\,\%$ for 50 realizations. This contrasts with the 10 000 simulations required to estimate the incident field mean energy with the same tolerance. Similar results were obtained for other wave periods. This analysis suggests that the scattering properties of large arrays have a low sensitivity to random variations in the phase of the wave forcing. For a regular array of ice floes, a relatively small number of random realizations of the sea state then suffices to obtain accurate estimates of the scattering properties of large arrays.

5.2. Simulation of a natural FSD

We model the MIZ as a randomly selected array of floes with different sizes, such that the diameters obey a power-law distribution, which is a standard empirical model of the FSD (see, e.g., Rothrock & Thorndike Reference Rothrock and Thorndike1984; Toyota et al. Reference Toyota, Takatsuji and Nakayama2006, Reference Toyota, Haas and Tamura2011). This is a unique feature of our three-dimensional attenuation model. Further, it allows us to simulate an MIZ with a high concentration, which would otherwise not be possible with a single floe size, e.g. the mean floe size.

We use an approach similar to that of Kohout & Bertino (Reference Dumont, Kohout and Bertino2011) and Williams et al. (Reference Williams, Bennetts, Squire, Dumont and Bertino2013a ) to parameterize the FSD, in which a bounded power-law distribution is used. We define the probability density function $(a/a_{\mathit{min}})^{-{\it\kappa}}$ for $a_{\mathit{min}}\leqslant a\leqslant a_{\mathit{max}}$ , where ${\it\kappa}$ is a constant parameter. Therefore, the probability that a floe has a radius $a<\mathfrak{a}$ is then given by

(5.4) $$\begin{eqnarray}\mathbb{P}(a<\mathfrak{a})=\frac{\mathfrak{a}^{1-{\it\kappa}}-a_{\mathit{min}}^{1-{\it\kappa}}}{a_{\mathit{max}}^{1-{\it\kappa}}-a_{\mathit{min}}^{1-{\it\kappa}}}.\end{eqnarray}$$

The distribution is then discretized so that a small finite number, $N_{b}$ say, of unique floe sizes is considered.

An algorithm to generate a distribution of ice floes in a slab using the bounded power-law FSD described here is described in appendix B. As the MIZ generated in our model is obtained by stacking together a large number of slabs, the same FSD is satisfied for the whole MIZ. The parameters of the algorithm are the dimensions of the ice-covered region in the slab, i.e. the width $L_{x}$ and breadth $L_{y}$ , the number of bins $N_{\mathit{b}}$ , the ice concentration $c$ , the minimum and maximum floe radii $a_{\mathit{min}}$ and $a_{\mathit{max}}$ , and the exponent in the power-law distribution ${\it\kappa}$ . An example of a random array generated using this algorithm is shown in figure 4(a) for parameter values $L_{x}=220~\text{m}$ , $L_{y}=8L_{x}$ , $N_{b}=11$ , $c=0.7$ , $a_{\mathit{min}}=10~\text{m}$ , $a_{\mathit{max}}=100~\text{m}$ and ${\it\kappa}=1.84$ . The last three parameters are the same as those used by Williams et al. (Reference Williams, Bennetts, Squire, Dumont and Bertino2013a ). This range of floe sizes is comparable to the floe sizes typically observed in real MIZs (see, e.g. Toyota et al. Reference Toyota, Takatsuji and Nakayama2006).

To generate a highly concentrated MIZ, our algorithm populates each slab with a large number of small floes. The effect of these small floes on wave interactions is likely to be negligible, however, while it increases the computational cost of the self-consistent method used to calculate the multiple scattering within each slab. Accordingly, we devise below a numerical test to determine the minimum floe size $a_{\mathit{min}}$ contributing to scattering by a large array.

Consider a randomly selected array in a single slab parameterized as before, but with larger breadth $L_{y}=51L_{x}$ . The array generated contains 488 floes. We compute ensemble averages of the reflection coefficient $R$ due to random realizations of the incident sea state with an energy spreading function given by (5.2). The calculations are repeated after removing all of the floes with the smallest floe size successively until $a_{\mathit{min}}=a_{\mathit{max}}$ , in which case the concentration is approximately $5\,\%$ .

Figure 4.

(a) Example of array generated using the random array generator described in appendix B. (b) Ensemble averages of reflection coefficients plotted against the minimum floe radius of the power-law FSD used to generate random realizations of the arrays. Results are given for single- and 10-slab arrays (solid and dashed lines respectively) and wave periods $T=6$ and 9 s (circle and square markers respectively). The wave forcing is a random sea state with normalized cosine square energy spreading.

Figure 4 shows ensemble averages of the reflection coefficient against the minimum floe radius $a_{\mathit{min}}$ , for wave periods $T=6$ and 9 s (solid lines with circle and square markers respectively). Each data point is calculated as the average of 100 simulations, each of which is characterized by a random realization of the array and the incident sea state. For small $a_{\mathit{min}}$ , the reflection coefficient remains roughly constant for both wave periods considered, suggesting that very small floes do not influence the scattering properties of the slab. A change of regime occurs for $a_{\mathit{min}}\approx 28~\text{m}$ at $T=6~\text{s}$ and $a_{\mathit{min}}\approx 55~\text{m}$ at $T=9~\text{s}$ , beyond which the reflection coefficient decreases as $a_{\mathit{min}}$ increases, so that the smallest floes start to contribute to scattering by the slab. These values of $a_{\mathit{min}}$ are used to define the critical minimum floe size $D_{\mathit{crit}}=2a_{\mathit{min}}$ , corresponding to the change of regime described above.

Our estimates of the critical minimum floe size, $D_{\mathit{crit}}\approx 56$ and $110~\text{m}$ , are similar to the corresponding open water wavelengths, ${\it\lambda}_{0}=2{\rm\pi}/k_{0}\approx 56$ and $126~\text{m}$ , for $T=6$  s and 9 s respectively. This suggests that scattering by floes smaller than the forcing wavelength is negligible and that these floes need not be included in the FSD.

We repeat the analysis for arrays composed of 10 slabs. To reduce the computing time, we solve the single-slab problem for 50 random realizations of the array and then perform random permutations of the pool of single-slab arrays to generate random 10-slab arrays (the validity of this approach will be discussed in § 5.3). Ensemble averages of the reflection coefficient for varying $a_{\mathit{min}}$ are plotted in figure 4 for $T=6$ and 9 s (dashed lines with circle and square markers respectively). We observe a two-regime dependence on $a_{\mathit{min}}$ , similar to the single-slab array, with $D_{\mathit{crit}}$ taking the same value for both wave periods, noting that the transition between the regimes is much smoother for $T=6~\text{s}$ . We deduce that the critical minimum floe size does not depend on the size of the array but seems to be an intrinsic property of the FSD, varying with wave period and possibly ice thickness (not studied here).

5.3. Averaging

The multiple-slab interaction technique described in § 3.2 performs efficiently as the computational cost depends linearly on the number of slabs. Computing the reflection and transmission matrices of each slab is more time-consuming, however, as a 2D multiple scattering problem, which is $O(N_{q}^{3})$ expensive, needs to be solved for each slab $q$ containing $N_{q}$ floes. To reduce the number of single-slab solutions to compute, we calculate the reflection and transmission matrices of a fixed number $S_{u}$ of unique slabs and store them, requiring $O(S_{u}N_{\mathit{ang}}^{2})$ memory space. Each realization of a multiple-slab array is then generated from random permutations (allowing repetitions) of the $S_{u}$ unique slabs. Bennetts (Reference Bennetts2011) used this method in a related acoustic problem, where each slab was composed of an infinite regular array of scatterers with different in-row spacings, and found that $S_{u}=50$ was sufficient to take ensemble averages of wave transmission by 100-slab arrays. In our case, the slabs contain only a finite number of scatterers, and the scatterers are of different sizes and not positioned in a regular manner. Therefore, it is unclear whether a value of $S_{u}$ similar to that of Bennetts (Reference Bennetts2011) will be appropriate.

We devised a numerical test to determine a suitable value of $S_{u}$ . Consider a 100-slab array with slab dimensions $L_{x}=220~\text{m}$ and $L_{y}=51L_{x}$ , and an FSD in each slab parameterized as in § 5.2. We compute the reflection coefficient of the array for different values $S_{u}$ at wave periods $T=6$ and 9 s with $a_{\mathit{min}}=37$ and $55~\text{m}$ respectively (using our findings from § 5.2). For each value of $S_{u}$ considered, a single sample of $S_{u}$ unique slabs is generated to perform random permutations, and independent samples are used for different values of $S_{u}$ . Each computed value of the reflection coefficient is the mean of an ensemble of 50 realizations of the array, where each realization is obtained by randomly permuting unique slabs from the same sample of slabs generated. The results of these computations are shown in figure 5, where error bars indicate the standard error of the mean.

Figure 5.

Ensemble average of the reflection coefficient by a 100-slab array computed for different values of $S_{u}$ . The array in each unique slab is generated using the FSD described in § 5.2 with dimensions $L_{x}=220~\text{m}$ and $L_{y}=51L_{x}$ . Results are shown for (a)  $T=6~\text{s}$ and (b)  $T=9~\text{s}$ . Averages are computed from 50 random realizations of the array and wave forcing (with a cosine square spread).

We observe a remarkable consistency in the computed values of the reflection coefficient, with two significant digit accuracy being obtained even for a small number of unique slabs, i.e. $S_{u}=10$ and $2$ for $T=6$  s and 9 s respectively. Specifically, the standard error of the mean values of $R$ for $S_{u}\geqslant 10$ and $2$ (for $T=6$  s and 9 s respectively) is ${<}0.1\,\%$ . It should be noted that each point is computed using an independent set of unique slabs, so that points obtained with different values of $S_{u}$ are uncorrelated. This suggests that a small value of $S_{u}$ is sufficient to simulate scattering by large random arrays of floes. In particular, the heterogeneity introduced in each unique slab (random floe packing) translates to larger scales when the slabs are stacked together. We fix $S_{u}=10$ for the remainder of the investigation, so that the memory space required to store the reflection and transmission matrices is $O(N_{\mathit{ang}}^{2})$ .

We also find that the procedure is very accurate, as each estimate of the reflection coefficient (i.e. for a single random realization of the array and wave forcing) has a relative standard error of the mean of approximately 2.5 and 1 % for $T=6$  s and 9 s respectively. Therefore, after averaging over 50 random realizations, the error drops to approximately 0.35 % and 0.15 % (respectively). We account for some variability in the values of $R$ for $S_{u}\geqslant 10$ by reducing the ensemble size used to average. Specifically, we use an ensemble of 10 random realizations, which gives a relative standard error below $1\,\%$ for both wave periods.

6.1. Attenuation coefficient

The key quantity of existing wave attenuation models in the MIZ is the attenuation coefficient, which defines the rate of exponential attenuation of wave energy in an ice-covered sea. At present, this is the only quantity used to parameterize wave–sea-ice interactions in large-scale IOMs (e.g. Williams et al. (Reference Williams, Bennetts, Squire, Dumont and Bertino2013a ,Reference Williams, Bennetts, Squire, Dumont and Bertino b ) using the scattering/viscous model of Bennetts & Squire (Reference Bennetts and Squire2012b )) and SWMs (Rogers & Orzech (Reference Rogers and Orzech2013) using the viscous models of Liu, Holt & Vachon (Reference Liu, Holt and Vachon1991a ) and Wang & Shen (Reference Wang and Shen2010)).

To extract the attenuation coefficient from our simulations, we fit an exponential curve to the computed data $E^{+}({\it\xi}_{q})$ , $q=0,\ldots ,220$ , i.e.

(6.4) $$\begin{eqnarray}E^{+}(x)\approx E^{+}(0)\text{e}^{-ax},\end{eqnarray}$$

where $a$ is the attenuation coefficient of wave energy and $E^{+}(0)=\exp (a{\it\xi}_{0})$ . We then use linear least-square regression to estimate the expected value $\overline{a}(T)$ for each wave period $T$ (with the overbar denoting the expected value of a random variable). Our approach is similar to experimental measurements, in which the attenuation coefficient is extracted from the actual wave energy profile through a realization of the MIZ.

For each estimated $\overline{a}$ value, we compute the standard error that accounts for the goodness of the least-square fit and the variability of the ensemble averaging process. The statistical method used to estimate the standard error is based on the random effects model and the maximum likelihood method. It is described by Brockwell & Gordon (Reference Brockwell and Gordon2001) in the context of medical science.

Figure 8.

Wave energy attenuation coefficient $\overline{a}$ as a function of wave period $T$ in the range 6–15 s. The error bars represent the standard error of each estimated value of $\overline{a}$ , and account for the goodness of the least-square fit and the variability in the ensemble of simulations.

Figure 8 shows the attenuation coefficient as a function of wave period, in the range $T=6$ –15 s. It should be noted that the order of magnitude of $a$ is $O(10^{-6}{-}10^{-5})$ , which is slightly lower than attenuation coefficients estimated from field observations (see, e.g., Wadhams et al. (Reference Wadhams, Squire, Goodman, Cowan and Moore1988), Meylan et al. (Reference Meylan, Bennetts and Kohout2014), observing that our simulations have not been parameterized to match these data sets in regard to the physical moduli involved). This is not surprising as scattering is the only physical process considered here, while in nature other dissipative processes would contribute to the decay of wave energy in the MIZ. The relatively wide incident directional spectrum used for the simulations (i.e. the cosine-squared spreading function) may also contribute to these low values of the attenuation coefficients.

Qualitatively, the dependence of the attenuation coefficient on the wave period is as expected, as $\overline{a}$ decreases for increasing $T$ . As the wave period becomes smaller, the attenuation coefficient seems to level off. Interestingly, this feature looks similar to the onset of the rollover effect observed in several data sets reported by Wadhams et al. (Reference Wadhams, Squire, Goodman, Cowan and Moore1988). We do not claim, however, that our model is capable of reproducing this effect, which was hypothesized to be a consequence of local wind wave generation, nonlinear wave–wave interactions or dissipative processes, none of which are considered in our model. An extensive sensitivity analysis would be required to interpret this feature in our model, which is beyond the scope of the present case study. We also observe an increase of the attenuation coefficient for $T=15~\text{s}$ , which is probably due to the relatively large variability of the attenuation coefficient between individual simulations when the attenuation rate is very small. It should be noted, however, that a similar effect has been reported by Wadhams et al. (Reference Wadhams, Squire, Goodman, Cowan and Moore1988) in their field data, although these authors did not comment on this feature.

6.2. Rate of spreading and distance to isotropy

As observed in figures 6(b) and 7, our model is capable of reproducing the directional spreading experienced by the wave field as it travels through the MIZ. We utilize the apparent linear relationship between ${\it\sigma}_{1}$ and $x$ observed in figure 6(b) to quantify the directional spreading, by fitting a linear curve to the subset of computed data ${\it\sigma}_{1}({\it\xi}_{q})$ , $q=q_{\mathit{min}},\ldots ,q_{\mathit{max}}$ . Thus,

(6.5) $$\begin{eqnarray}{\it\sigma}_{1}(x)\approx {\it\sigma}_{1}^{0}+sx,\end{eqnarray}$$

where $s$ is the (constant) rate of directional spreading, ${\it\sigma}_{1}^{0}={\it\sigma}_{1}(0)$ , and the lower and upper bounds, $q_{\mathit{min}}$ and $q_{\mathit{max}}$ , are determined manually for each $T$ . A linear least-square regression is used to estimate $\overline{s}(T)$ , as for the attenuation coefficient. The standard errors are also calculated using the statistical method mentioned in § 6.1.

Figure 9.

(a) Rate of directional spreading $\overline{s}$ as a function of wave period $T$ in the range $6$ – $15$  s. (b) Distance to isotropy $\overline{x^{(iso)}}$ plotted over the same range of wave periods. Error bars are computed as in figure 8.

The mean rate of directional spreading, $\overline{s}$ , is plotted as a function of wave period in figure 9(a). For $T\leqslant 10~\text{s}$ , the rate of spreading only depends weakly on the wave period, with $\overline{s}=5{-}6\times 10^{-6}~\text{m}^{-1}$ . The minimum observed at $T=8~\text{s}$ and the maximum at $T=10~\text{s}$ are difficult to interpret physically, although random variability may be a factor as there are relatively large error bars in this regime. For longer wave periods ( $T\geqslant 10~\text{s}$ ), $\overline{s}$ decreases monotonically for increasing wave period, which is sensible as directional spreading is caused by scattering, which diminishes for increasingly long waves.

We now use the more tangible quantity of the distance to isotropy, denoted $x^{(iso)}$ , to describe the spreading. We stated earlier that the wave field becomes isotropic when ${\it\sigma}_{1}={\it\sigma}_{1}^{(iso)}\approx 0.8525$ , so using (6.5) we can infer the distance from the ice edge to reach isotropy to be

(6.6) $$\begin{eqnarray}x^{(iso)}=\frac{{\it\sigma}_{1}^{(iso)}-{\it\sigma}_{1}^{0}}{s}.\end{eqnarray}$$

Care must be taken in estimating $x^{(iso)}$ and its variance for each simulation, as ${\it\sigma}_{1}^{0}$ and $s$ are random variables with an underlying probability distribution, for which we have calculated the expected values and covariance matrix. Equation (6.6) is actually the first-order approximation of the expected value of the ratio distribution. A (better) second-order formula for estimating the expected value of the distance to isotropy is given by

(6.7) $$\begin{eqnarray}\overline{x^{(iso)}}\approx \frac{{\it\sigma}_{1}^{(iso)}-\overline{{\it\sigma}_{1}^{0}}}{\overline{s}}-\frac{\text{Cov}({\it\sigma}_{1}^{0},s)}{\overline{s}^{2}}+\frac{\text{Var}(s)({\it\sigma}_{1}^{(iso)}-\overline{{\it\sigma}_{1}^{0}})}{\overline{s}^{3}},\end{eqnarray}$$

where $\text{Var}$ and $\text{Cov}$ denote the variance and covariance of random variables respectively (see, e.g., Stuart & Ord Reference Stuart and Ord1999). A first-order formula can also be derived for the variance of $x^{(iso)}$ , i.e.

(6.8) $$\begin{eqnarray}\text{Var}(x^{(iso)})\approx \left(\frac{{\it\sigma}_{1}^{(iso)}-\overline{{\it\sigma}_{1}^{0}}}{\overline{s}}\right)^{2}\left[\frac{\text{Var}(s)}{({\it\sigma}_{1}^{(iso)}-\overline{{\it\sigma}_{1}^{0}})^{2}}-2\frac{\text{Cov}({\it\sigma}_{1}^{0},s)}{({\it\sigma}_{1}^{(iso)}-\overline{{\it\sigma}_{1}^{0}})\overline{s}}+\frac{\text{Var}({\it\sigma}_{1}^{0})}{\overline{s}^{2}}\right].\end{eqnarray}$$

Figure 9(b) depicts the expected value of the distance to isotropy $\overline{x^{(iso)}}$ against wave period. The transition between the low-period and high-period regimes is clearly observed here. Isotropy is reached within the first 40 km of the simulated MIZ for $T\leqslant 10~\text{s}$ , and $\overline{x^{(iso)}}$ varies little with $T$ in this regime. On the other hand, $\overline{x^{(iso)}}$ increases abruptly for $T\geqslant 11~\text{s}$ , where the wave field spreads very slowly towards isotropy. For $T>12~\text{s}$ , values of $\overline{x^{(iso)}}$ greater than 500 km are computed, suggesting that long waves experience next to no spreading within the extent of a typical MIZ, i.e. over $O$ (10–100) km. The transition between the two regimes correlates with the prescribed maximum floe diameter of 200 m, which is the open water wavelength of an 11.3 s wave. This finding suggests that waves longer than the maximum floe size do not experience significant directional spreading in the MIZ, which agrees with the observations of Wadhams et al. (Reference Wadhams, Squire, Ewing and Pascal1986).