1 Introduction
A typical feature of both polar regions is that vast areas of ocean are covered by ice. Observations showed that the waves generated in open sea could penetrate far into the ice-covered region (Squire et al. Reference Squire, Dugan, Wadhams, Rottier and Liu1995). The features of the wave will change due to its interaction with the ice sheet. One example is that the reduction of the ice extent and thickness will lead to an overall uptrend of the wave height. The topic of coupled ocean wave and sea ice interactions has increasingly attracted research interests over the last two decades. This is due to the fact the polar ice cover can be greatly affected by the global climate change, which will in turn affect the global environment (Squire Reference Squire2020). It is also due to the requirement of safe navigation of the potential shipping routes and offshore operations in Arctic region (Eliasson et al. Reference Eliasson, Ulfarsson, Valsson and Gardarsson2017).
Field experiments conducted by Robin (Reference Robin1963) and Squire et al. (Reference Squire, Robinson, Langhorne and Haskell1988) have confirmed that the large ice sheet would bend to let the energy pass through in form of flexural-gravity waves. Based on the linear velocity potential theory for fluid flow and thin elastic plate model for ice sheet deflection, a number of solution approaches have been developed for wave scattering in icy waters. When waves propagate from open sea to the ice-covered water or the other way around, the energy will be partially reflected and transmitted due to the change of the dispersion relation. The problem of wave interactions with an homogeneous semi-infinite ice sheet, was solved for the normal incident wave by Fox & Squire (Reference Fox and Squire1990) and for the oblique incident wave by Fox & Squire (Reference Fox and Squire1994), based on the method of matched eigenfunction expansions (MEE). It was found that strong transmission would be more likely for long waves. For the oblique case, when the incident wave angle was larger than a critical value, no wave would transmit into the ice. Matched eigenfunction expansions were also used by Sahoo, Yip & Chwang (Reference Sahoo, Yip and Chwang2001) to the same problem, but with the unknowns found through introducing an orthogonal inner product. In addition, the problem can be solved by the Wiener–Hopf method (WHM), as was given by an unpublished report of Evans & Davies (Reference Evans and Davies1968). Later, some modifications were introduced into the WHM, which enabled the numerical computations to be more easily carried out, e.g. Balmforth & Craster (Reference Balmforth and Craster1999), Chung & Fox (Reference Chung and Fox2002) and Tkacheva (Reference Tkacheva2001, Reference Tkacheva2004). In Balmforth & Craster (Reference Balmforth and Craster1999), the Timoshenko–Mindlin model, which further included the effects of rotary inertia and transverse shear of the plate, was adopted for the ice sheet. The authors then showed that the Kirchhoff–Love model for a thin elastic plate would give similar results to those by the Timoshenko–Mindlin model for a wide range of thickness of the ice sheet. In addition to the MEE and WHM, Linton & Chung (Reference Linton and Chung2003) applied the residue calculus technique (RCT) to the semi-infinite ice sheet, and obtained solutions which were equivalent to those provided by WHM. In some cases, the ice sheet may not extend to infinity. It may have a section of free surface and then continue with ice cover. Chung & Linton (Reference Chung and Linton2005) considered wave propagations through an ice polynya by RCT, and found that perfect transmission could occur at an infinite number of discrete wave frequencies. Based on the Green function method (GFM), Squire & Dixon (Reference Squire and Dixon2001b) studied flexural-gravity wave interactions with an embedded iceberg. It may be noted that the GFM can be also used to solve the wave propagation through an ice floe, as for example in Meylan & Squire (Reference Meylan and Squire1994). Williams & Squire (Reference Williams and Squire2006) considered a more general problem by RCT and WHM, i.e. three connected ice sheets with the first and last ones being semi-infinite. The above individual scatters can be also assembled together to model the wave evolutions through icy waters in a large scale, as has been done in Bennetts & Squire (Reference Bennetts and Squire2012).
Apart from the above two-dimensional (2-D) problems, three-dimensional (3-D) problems have also been considered, which reflect more accurately the physics in many cases. For a circular ice floe, Meylan & Squire (Reference Meylan and Squire1996) solved the problem by the boundary integral equation method, based on either the eigenfunction expansions of the potential and its normal derivative on the circular ice floe or by writing the normal derivative in terms of another integral equation. The problem was then extended to finite water depth by Peter, Meylan & Chung (Reference Peter, Meylan and Chung2004). To simulate the propagation of waves through marginal ice zone (MIZ), ice floe with a complicated shape needs to be considered. This was done by Meylan (Reference Meylan2002) and Wang & Meylan (Reference Wang and Meylan2004) through the boundary element method for fluid flow and a finite element method for the ice sheet. Instead of determining the free modes of vibration or the dry modes of the floe independently, Bennetts & Williams (Reference Bennetts and Williams2010) solved the fluid motion and ice sheet flexure together, i.e. the governing equation for ice sheet flexure was incorporated into the water surface condition for fluid flow through kinematic and dynamic conditions at the interface, as done in most of the other works. In Bennetts & Williams (Reference Bennetts and Williams2010), the wave interactions with a polynya embedded in an infinite ice sheet were also solved. Based on the solution of a single ice floe and the slab clustering method (known as SCM), originally for multiple structures interaction in open sea, Montiel, Squire & Bennetts (Reference Montiel, Squire and Bennetts2016) investigated the wave spectra evolution features in the MIZ which were comprised of tens of thousands of ice floes.
Ice sheets in nature are rarely homogeneous and perfect. One of the imperfections is the crack. Based on the MEE originally developed by Fox & Squire (Reference Fox and Squire1990), Barrett & Squire (Reference Barrett and Squire1996) solved the problem of an infinitely long straight-line crack in an ice sheet on the water surface. Numerical results showed that for short waves most of the energy would be reflected back, which was similar to the case of a wave impinging a semi-infinite ice sheet. It was found that there would be a frequency at which perfect transmission would occur. Based on the Green function for an ice sheet without a crack, an analytical solution for reflection and transmission coefficients could be derived, e.g. by Squire & Dixon (Reference Squire and Dixon2000) for a single crack through an explicit formulae and by Squire & Dixon (Reference Squire and Dixon2001a) for multiple cracks through a matrix equation. Single or multiple parallel infinitely long cracks were also solved by Evans & Porter (Reference Evans and Porter2003) and Porter & Evans (Reference Porter and Evans2006) for a finite water depth. Li, Wu & Ji (Reference Li, Wu and Ji2018b) derived a Green function which satisfied the conditions at the crack. This enabled the expression for the diffraction potential to be written explicitly. The advantage of this method was that it allowed a multipole for a circle to be constructed directly (Li et al. Reference Li, Wu and Ji2018b), and the obtained Green function could be easily extended to multiple cracks and used for a body of arbitrary shape (Li, Wu & Ji Reference Li, Wu and Ji2018a).
For the 3-D problems, the work has been mainly limited to the cracks of some special shapes. Porter & Evans (Reference Porter and Evans2007) considered parallel straight-line cracks. The jumps of the displacement and slope of the crack deflection were expanded into Chebyshev polynomials, which took into account the behaviour at the crack ends. Li, Wu & Shi (Reference Li, Wu and Shi2018c) considered a circular crack. The jumps of the displacement and slope of the crack deflection were expanded into Fourier series and the unknown coefficients were found explicitly. In addition to parallel straight lines and closed circles, the cracks in polar regions may have other shapes. Also, when there is a high concentration of ice floes and the gaps between floes are very narrow, this could be treated as a crack problem. For a crack in a general shape, the problem becomes a major challenge. One of the difficulties is that the free edge conditions at the crack contain high derivatives. When this is applied to an already singular boundary integral equation, it becomes a high-order hypersingular equation. The solution of such an equation is not easy to obtain and therefore the physics in such a case is not easy to capture. In this work, we shall develop a methodology to overcome the difficulty, and then solve the problem of wave interactions with multiple cracks of arbitrary shapes and locations. An individual crack can be either closed or open. In the methodology, we extend the crack in the vertical direction into the fluid domain to form an artificial vertical surface. On this surface, continuity of the velocity and pressure of the fluid will be imposed. The 3-D problem is then decomposed into an infinite number of coupled 2-D problems. The methodology is verified through the analytical solution for the circular crack and then used for a variety of case studies.
The paper is organized as follows. The mathematical problem for flexural-gravity wave propagations in an ice sheet with multiple arbitrary shaped cracks is formulated in § 2, and the conditions satisfied at the crack together with those at the ends of an open crack are described. The solution procedures are constructed in § 3, and the numerical discretization of the hypersingular integral differential equations are given in § 4. Extensive results for cracks with various shapes and distributions are presented and discussed in § 5, and conclusions are made in § 6. In appendix A, the diffraction coefficient corresponding to the far field diffracted wave is derived, and an energy conservation relationship for the scattered cross-section is provided.
2 Mathematical model
We consider the problem of wave interaction with multiple cracks of arbitrary shapes in an ice sheet, as sketched in figure 1. The ice sheet is horizontally infinitely extended, and floating on water of finite depth. To describe the problem, a Cartesian coordinate system 
$O$–
$xyz$ is defined, with the 
$O$–
$xy$ plane being the undisturbed interface of water and ice sheet, and 
$z$-axis pointing vertically upwards. The direction of flexural-gravity incident wave from infinity is assumed to form an angle 
$\unicode[STIX]{x1D6FD}$ with positive 
$x$-axis, and the water depth 
$H$ is assumed to be constant, as shown in figure 1. The shape of crack 
$i$ can be described parametrically by
$$\begin{eqnarray}\unicode[STIX]{x1D730}_{i}(s)=(x(s),y(s)),\quad (-\unicode[STIX]{x1D6FE}_{i}<s<+\unicode[STIX]{x1D6FE}_{i},1\leqslant i\leqslant N),\end{eqnarray}$$ where 
$2\unicode[STIX]{x1D6FE}_{i}$ is the total arclength of crack 
$i$, and 
$s$ is the curvilinear coordinate along the crack.
Figure 1. Coordinate system and sketch of the problem.
Based on the assumption that the fluid is inviscid, incompressible and homogeneous, and its motion is irrotational, the velocity potential 
$\unicode[STIX]{x1D6F7}$ can be introduced to describe the fluid flow. When the amplitude of wave motion is small compared to its length, the linearized velocity potential theory can be further used. For sinusoidal motion in time with radian frequency 
$\unicode[STIX]{x1D714}$, we may write the total velocity potential as
$$\begin{eqnarray}\unicode[STIX]{x1D6F7}(x,y,z,t)=\text{Re}[\unicode[STIX]{x1D719}(x,y,z)\text{e}^{\text{i}\unicode[STIX]{x1D714}t}],\end{eqnarray}$$where
$$\begin{eqnarray}\unicode[STIX]{x1D719}=\unicode[STIX]{x1D719}_{I}+\unicode[STIX]{x1D719}_{D},\end{eqnarray}$$ with 
$\unicode[STIX]{x1D719}_{I}$ and 
$\unicode[STIX]{x1D719}_{D}$ as the incident velocity potential and diffracted velocity potential by the cracks, respectively. The conservation of mass requires that the velocity potential 
$\unicode[STIX]{x1D719}$ should satisfy the following Laplace equation:
$$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}+\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}z^{2}}=0,\end{eqnarray}$$throughout the fluid, where
$$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}=\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}x^{2}}+\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}y^{2}},\end{eqnarray}$$ is the Laplacian in the horizontal plane. Following Squire (Reference Squire2011) and others, the ice sheet is modelled as a thin elastic plate, and its properties (i.e. Young’s modulus 
$E$, Poisson’s ratio 
$\unicode[STIX]{x1D708}$, density 
$\unicode[STIX]{x1D70C}_{i}$ and thickness 
$h$) are assumed to be constant. Assuming that there is no gap between the ice sheet and the water upper surface, we can write the boundary condition at their interface as
$$\begin{eqnarray}(L\unicode[STIX]{x1D6FB}^{4}+\unicode[STIX]{x1D70C}_{w}g-m\unicode[STIX]{x1D714}^{2})\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}z}-\unicode[STIX]{x1D70C}_{w}\unicode[STIX]{x1D714}^{2}\unicode[STIX]{x1D719}=0,\quad (\boldsymbol{r}\neq \unicode[STIX]{x1D730}_{i},z=0),\end{eqnarray}$$ where 
$\boldsymbol{r}=(x,y)$, 
$L=Eh^{3}/[12(1-\unicode[STIX]{x1D708}^{2})]$ is the effective flexural rigidity of the ice sheet, 
$m=\unicode[STIX]{x1D70C}_{i}h$ is the corresponding areal density, 
$\unicode[STIX]{x1D70C}_{w}$ is the density of water and 
$g$ is the acceleration due to gravity. On both sides of the crack, the edges of the ice sheet are assumed to be free to move, which means that the following zero bending moment and modified shear force conditions should be satisfied (Timoshenko & Woinowsky Reference Timoshenko and Woinowsky1959):
$$\begin{eqnarray}{\mathcal{B}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}z}\right)=0\quad \text{and}\quad {\mathcal{S}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}z}\right)=0,\quad (\boldsymbol{r}=\unicode[STIX]{x1D730}_{i},z=0),\end{eqnarray}$$ where the operator 
${\mathcal{B}}$ and 
${\mathcal{S}}$ are, respectively, defined as
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{B}}=\unicode[STIX]{x1D6FB}^{2}-\unicode[STIX]{x1D708}_{0}\left(\sin ^{2}\unicode[STIX]{x1D6E9}\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}x^{2}}+\cos ^{2}\unicode[STIX]{x1D6E9}\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}y^{2}}-\sin 2\unicode[STIX]{x1D6E9}\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}y}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{S}}=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}n}\unicode[STIX]{x1D6FB}^{2}+\unicode[STIX]{x1D708}_{0}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\left[\cos 2\unicode[STIX]{x1D6E9}\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}y}+\frac{\sin 2\unicode[STIX]{x1D6E9}}{2}\left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}y^{2}}-\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}x^{2}}\right)\right], & \displaystyle\end{eqnarray}$$ with 
$\unicode[STIX]{x1D708}_{0}=1-\unicode[STIX]{x1D708}$, and 
$\boldsymbol{n}$ and 
$\boldsymbol{s}$ are the unit vectors along the normal and tangential directions, respectively. Here, 
$\unicode[STIX]{x1D6E9}(s)$ is the angle between positive 
$\boldsymbol{n}$ and the 
$x$-axis, and thus 
$\boldsymbol{n}=(\cos \unicode[STIX]{x1D6E9},\sin \unicode[STIX]{x1D6E9})$ and 
$\boldsymbol{s}=(-\text{sin}\unicode[STIX]{x1D6E9},\cos \unicode[STIX]{x1D6E9})$. By noticing
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}x}{\unicode[STIX]{x2202}s}=-\text{sin}\,\unicode[STIX]{x1D6E9}\quad \text{and}\quad \frac{\unicode[STIX]{x2202}y}{\unicode[STIX]{x2202}s}=\cos \unicode[STIX]{x1D6E9},\end{eqnarray}$$ the operators 
${\mathcal{B}}$ and 
${\mathcal{S}}$ in (2.8) and (2.9) can be also expressed in the curvilinear coordinate system as
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{B}}=\unicode[STIX]{x1D6FB}^{2}-\unicode[STIX]{x1D708}_{0}\left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}s^{2}}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}s}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}n}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{S}}=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}n}\unicode[STIX]{x1D6FB}^{2}+\unicode[STIX]{x1D708}_{0}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}s\unicode[STIX]{x2202}n}-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}s}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\right), & \displaystyle\end{eqnarray}$$ where the partial derivatives with respect to 
$s$ are carried out with respect to the curvilinear coordinate. It should be noticed that (2.12) is in fact an equivalent shear force after the three boundary conditions at a plate edge having been combined into two. The second term in (2.12) is from the twisting moment, which is equivalently replaced by the force perpendicular to the plate. In doing so, when a plate edge has a sharp corner, it leads to a concentrated force at the corner, which has to be set to zero as an additional condition. For an open crack, at the ends of the crack the free edge in fact turns 
$360^{\circ }$. Thus this additional condition should be imposed at the sharp corner of a plate (Timoshenko & Woinowsky Reference Timoshenko and Woinowsky1959)
$$\begin{eqnarray}\left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}s\unicode[STIX]{x2202}n}-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}s}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\right)\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}z}=0,\quad (s\rightarrow \pm \unicode[STIX]{x1D6FE}_{i},z=0).\end{eqnarray}$$The edge conditions described above are on the basis that the ice sheets on both sides of the crack are completely detached from each other, which is suitable when the ratio of width of the gap to the typical length scale of the ice sheet is small (Li et al. Reference Li, Wu and Ji2018b). It might be noticed that the implications of the mathematical model for the ice crack have been reviewed and discussed by Squire (Reference Squire2007). On the flat seabed, the impermeable condition provides
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}z}=0,\quad (z=-H).\end{eqnarray}$$At infinity the radiation condition is imposed, which requires that the diffracted wave should propagate outwards, or
$$\begin{eqnarray}\lim _{r\rightarrow \infty }\sqrt{r}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{D}}{\unicode[STIX]{x2202}r}+\text{i}\unicode[STIX]{x1D705}_{0}\unicode[STIX]{x1D719}_{D}\right)=0,\end{eqnarray}$$ where 
$\unicode[STIX]{x1D705}_{0}$ is the flexural-gravity wavenumber. It may be noticed that 
$\unicode[STIX]{x1D705}_{0}$ is the purely positive real root of the dispersion equation 
$K(\unicode[STIX]{x1D714},\unicode[STIX]{x1D705})=0$ for a flexural-gravity wave in the ice sheet, with
$$\begin{eqnarray}K(\unicode[STIX]{x1D714},\unicode[STIX]{x1D705})\equiv (L\unicode[STIX]{x1D705}^{4}+\unicode[STIX]{x1D70C}_{w}g-m\unicode[STIX]{x1D714}^{2})\unicode[STIX]{x1D705}\tanh (\unicode[STIX]{x1D705}H)-\unicode[STIX]{x1D70C}_{w}\unicode[STIX]{x1D714}^{2}.\end{eqnarray}$$3 Solution procedures
To solve the boundary value problem described in the section above, through using a variable separation method, we may write the diffracted velocity potential 
$\unicode[STIX]{x1D719}_{D}$ in the following form:
$$\begin{eqnarray}\unicode[STIX]{x1D719}_{D}=\mathop{\sum }_{m}\unicode[STIX]{x1D711}_{m}(x,y)\unicode[STIX]{x1D713}_{m}(z).\end{eqnarray}$$ Substituting (3.1) into the governing Laplace equation (2.4), we can obtain 
$(\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D711}_{m})/\unicode[STIX]{x1D711}_{m}=-(\text{d}^{2}\unicode[STIX]{x1D713}_{m}/\text{d}z^{2})/\unicode[STIX]{x1D713}_{m}$. This equation has to be equal to a constant. Letting it be 
$-\unicode[STIX]{x1D705}_{m}^{2}$, we have
$$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D711}_{m}+\unicode[STIX]{x1D705}_{m}^{2}\unicode[STIX]{x1D711}_{m}=0,\end{eqnarray}$$and
$$\begin{eqnarray}\frac{\text{d}^{2}\unicode[STIX]{x1D713}_{m}}{\text{d}z^{2}}-\unicode[STIX]{x1D705}_{m}^{2}\unicode[STIX]{x1D713}_{m}=0,\end{eqnarray}$$ for 
$\unicode[STIX]{x1D711}_{m}$ and 
$\unicode[STIX]{x1D713}_{m}$, respectively. By further applying the boundary conditions in (2.6) and (2.14) to (3.3), we have
$$\begin{eqnarray}\unicode[STIX]{x1D713}_{m}(z)=\frac{\cosh [\unicode[STIX]{x1D705}_{m}(z+H)]}{\cosh (\unicode[STIX]{x1D705}_{m}H)}.\end{eqnarray}$$ Here, 
$\unicode[STIX]{x1D705}_{m}$ are the roots of the dispersion equation (2.16), with 
$\unicode[STIX]{x1D705}_{-2}$ and 
$\unicode[STIX]{x1D705}_{-1}$ as two complex roots with negative imaginary parts and symmetric about the imaginary axis, 
$\unicode[STIX]{x1D705}_{0}$ as the purely positive real root, 
$\unicode[STIX]{x1D705}_{m}$
$(m=1,\ldots ,\infty )$ as an infinite number of purely negative imaginary roots. Taking the summation in (3.1) through all possible eigenvalues and eigenfunctions, we have
$$\begin{eqnarray}\unicode[STIX]{x1D719}_{D}=\mathop{\sum }_{m=-2}^{\infty }\unicode[STIX]{x1D711}_{m}(x,y)\unicode[STIX]{x1D713}_{m}(z).\end{eqnarray}$$ Then the three-dimensional boundary value problems for 
$\unicode[STIX]{x1D719}_{D}$ are transformed into the two-dimensional ones for 
$\unicode[STIX]{x1D711}_{m}$. The incident velocity potential 
$\unicode[STIX]{x1D719}_{I}$ in (2.3) may be written as
$$\begin{eqnarray}\unicode[STIX]{x1D719}_{I}=\unicode[STIX]{x1D711}_{I}(x,y)\unicode[STIX]{x1D713}_{0}(z),\end{eqnarray}$$where
$$\begin{eqnarray}\unicode[STIX]{x1D711}_{I}(x,y)=I\text{e}^{-\text{i}\unicode[STIX]{x1D705}_{0}(x\cos \unicode[STIX]{x1D6FD}+y\sin \unicode[STIX]{x1D6FD})}/\text{tanh}(\unicode[STIX]{x1D705}_{0}H),\end{eqnarray}$$ with 
$I=\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D6EC}/\unicode[STIX]{x1D705}_{0}$ and 
$\unicode[STIX]{x1D6EC}$ as the amplitude of the incident wave.
The two-dimensional Helmholtz equation in (3.2) may be converted into an integral equation through using the Green function
$$\begin{eqnarray}G_{m}(p,q)=\frac{\unicode[STIX]{x03C0}}{2\text{i}}H_{0}^{(2)}(\unicode[STIX]{x1D705}_{m}R),\end{eqnarray}$$ where 
$H_{0}^{(2)}$ is the zeroth-order Hankel function of second kind (Abramowitz & Stegun Reference Abramowitz and Stegun1965), and 
$R$ is the distance between the field point 
$p(x,y)$ and source point 
$q(\unicode[STIX]{x1D709},\unicode[STIX]{x1D702})$ in the two-dimensional horizontal plane. We may note that 
$G_{m}$ in (3.8) satisfies the Helmholtz equation and the radiation condition in (2.15). Applying Green’s identity to 
$G_{m}$ and 
$\unicode[STIX]{x1D711}_{m}$, we have
$$\begin{eqnarray}\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D711}_{m}(p)=\int _{\unicode[STIX]{x1D6E4}}\left[G_{m}(p,q)\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D711}_{m}(q)}{\unicode[STIX]{x2202}n_{q}}-\frac{\unicode[STIX]{x2202}G_{m}(p,q)}{\unicode[STIX]{x2202}n_{q}}\unicode[STIX]{x1D711}_{m}(q)\right]\,\text{d}s_{q},\end{eqnarray}$$ where the integral and normal derivative are carried out with respect to the source point 
$q$, and 
$\unicode[STIX]{x1D6FC}$ is the two-dimensional solid angle at the field point 
$p$. In (3.9), 
$\unicode[STIX]{x1D6E4}$ should include all the lines along the ice sheet edges. However, the closed line at infinity has been excluded due to the fact that both 
$G_{m}$ and 
$\unicode[STIX]{x1D711}_{m}$ satisfy the radiation condition in (2.15), and the integrand becomes zero there. For the ice crack, we may denote the variable on its two sides with superscripts 
$+$ and 
$-$, respectively. Noticing that on the 
$+$ and 
$-$ sides, the 
$\boldsymbol{n}$ in (3.9) are opposite to each other, it can be rewritten as
$$\begin{eqnarray}\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D711}_{m}(p)=\mathop{\sum }_{i=1}^{N}\int _{\unicode[STIX]{x1D6E4}_{i}}\left[G_{m}(p,q)\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D711}}_{m}(q)}{\unicode[STIX]{x2202}n_{q}}-\frac{\unicode[STIX]{x2202}G_{m}(p,q)}{\unicode[STIX]{x2202}n_{q}}\tilde{\unicode[STIX]{x1D711}}_{m}(q)\right]\,\text{d}s_{q},\end{eqnarray}$$ where 
$\unicode[STIX]{x1D6E4}_{i}$ is the route along the crack 
$i$, 
$\boldsymbol{n}$ is now from the 
$-$ side to the 
$+$ side. 
$\tilde{\unicode[STIX]{x1D711}}_{m}$ and 
$\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D711}}_{m}/\unicode[STIX]{x2202}n$ are the corresponding jumps of velocity potential and its normal derivative across the crack, or
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\unicode[STIX]{x1D711}}_{m}=\unicode[STIX]{x1D711}_{m}^{-}-\unicode[STIX]{x1D711}_{m}^{+}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D711}}_{m}}{\unicode[STIX]{x2202}n}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D711}_{m}^{-}}{\unicode[STIX]{x2202}n}-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D711}_{m}^{+}}{\unicode[STIX]{x2202}n}. & \displaystyle\end{eqnarray}$$ Once the values of 
$\tilde{\unicode[STIX]{x1D711}}_{m}$ and 
$\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D711}}_{m}/\unicode[STIX]{x2202}n$ along the cracks are determined, the diffraction potential 
$\unicode[STIX]{x1D719}_{D}$ in any point of fluid can be obtained through (3.5) and (3.10). The main task is then to find these two jumps. To do that, we extend each ice crack vertically from the ice sheet to the seabed to form a vertical surface 
$\unicode[STIX]{x1D6F4}_{i}$, as shown in figure 1. On both sides of the artificial surface 
$\unicode[STIX]{x1D6F4}_{i}$, we denote the total velocity potential 
$\unicode[STIX]{x1D719}$ as 
$\unicode[STIX]{x1D719}^{+}$ and 
$\unicode[STIX]{x1D719}^{-}$, respectively, where the superscripts correspond to those of 
$\unicode[STIX]{x1D711}_{m}$. In a fluid, the pressure and fluid velocity should be continuous across each 
$\unicode[STIX]{x1D6F4}_{i}$, or
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}^{+}=\unicode[STIX]{x1D719}^{-}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}^{+}}{\unicode[STIX]{x2202}n}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}^{-}}{\unicode[STIX]{x2202}n}, & \displaystyle\end{eqnarray}$$ for 
$-H\leqslant z\leqslant 0$. Different from the free surface problem, the eigenfunction in (3.4) is not orthogonal in the usual sense. However, we may introduce the following orthogonal inner product (Sahoo et al. Reference Sahoo, Yip and Chwang2001):
$$\begin{eqnarray}\langle \unicode[STIX]{x1D713}_{m},\unicode[STIX]{x1D713}_{\tilde{m}}\rangle =\int _{-H}^{0}\unicode[STIX]{x1D713}_{m}\unicode[STIX]{x1D713}_{\tilde{m}}\,\text{d}z+\frac{L}{\unicode[STIX]{x1D70C}_{w}\unicode[STIX]{x1D714}^{2}}\left(\frac{\text{d}\unicode[STIX]{x1D713}_{m}}{\text{d}z}\frac{\text{d}^{3}\unicode[STIX]{x1D713}_{\tilde{m}}}{\text{d}z^{3}}+\frac{\text{d}^{3}\unicode[STIX]{x1D713}_{m}}{\text{d}z^{3}}\frac{\text{d}\unicode[STIX]{x1D713}_{\tilde{m}}}{\text{d}z}\right)_{z=0},\end{eqnarray}$$ which provides that 
$\langle \unicode[STIX]{x1D713}_{m},\unicode[STIX]{x1D713}_{\tilde{m}}\rangle =0$ if 
$m\neq \tilde{m}$ and 
$\langle \unicode[STIX]{x1D713}_{m},\unicode[STIX]{x1D713}_{\tilde{m}}\rangle =Q_{m}$ if 
$m=\tilde{m}$, with
$$\begin{eqnarray}Q_{m}=\frac{2\unicode[STIX]{x1D705}_{m}H+\sinh (2\unicode[STIX]{x1D705}_{m}H)}{4\unicode[STIX]{x1D705}_{m}\cosh ^{2}(\unicode[STIX]{x1D705}_{m}H)}+\frac{2L\unicode[STIX]{x1D705}_{m}^{4}}{\unicode[STIX]{x1D70C}_{w}\unicode[STIX]{x1D714}^{2}}\tanh ^{2}(\unicode[STIX]{x1D705}_{m}H).\end{eqnarray}$$ To satisfy the continuity condition of pressure in (3.13), applying the inner product to 
$\unicode[STIX]{x1D719}^{-}-\unicode[STIX]{x1D719}^{+}$ and 
$\unicode[STIX]{x1D713}_{\tilde{m}}$, we have
$$\begin{eqnarray}\displaystyle \langle (\unicode[STIX]{x1D719}^{-}-\unicode[STIX]{x1D719}^{+}),\unicode[STIX]{x1D713}_{\tilde{m}}\rangle & = & \displaystyle \int _{-H}^{0}\left(\unicode[STIX]{x1D719}^{-}-\unicode[STIX]{x1D719}^{+}\right)\unicode[STIX]{x1D713}_{\tilde{m}}\,\text{d}z\nonumber\\ \displaystyle & & \displaystyle +\,\frac{L}{\unicode[STIX]{x1D70C}_{w}\unicode[STIX]{x1D714}^{2}}\left[\frac{\unicode[STIX]{x2202}(\unicode[STIX]{x1D719}^{-}-\unicode[STIX]{x1D719}^{+})}{\unicode[STIX]{x2202}z}\frac{\text{d}^{3}\unicode[STIX]{x1D713}_{\tilde{m}}}{\text{d}z^{3}}+\frac{\unicode[STIX]{x2202}^{3}(\unicode[STIX]{x1D719}^{-}-\unicode[STIX]{x1D719}^{+})}{\unicode[STIX]{x2202}z^{3}}\frac{\text{d}\unicode[STIX]{x1D713}_{\tilde{m}}}{\text{d}z}\right]_{z=0}.\qquad\end{eqnarray}$$Similarly, for the continuity condition of velocity in (3.14), we have
$$\begin{eqnarray}\displaystyle \left\langle \frac{\unicode[STIX]{x2202}(\unicode[STIX]{x1D719}^{-}-\unicode[STIX]{x1D719}^{+})}{\unicode[STIX]{x2202}n},\unicode[STIX]{x1D713}_{\tilde{m}}\right\rangle & = & \displaystyle \int _{-H}^{0}\frac{\unicode[STIX]{x2202}(\unicode[STIX]{x1D719}^{-}-\unicode[STIX]{x1D719}^{+})}{\unicode[STIX]{x2202}n}\unicode[STIX]{x1D713}_{\tilde{m}}\,\text{d}z\nonumber\\ \displaystyle & & \displaystyle +\,\frac{L}{\unicode[STIX]{x1D70C}_{w}\unicode[STIX]{x1D714}^{2}}\left[\frac{\unicode[STIX]{x2202}^{2}(\unicode[STIX]{x1D719}^{-}-\unicode[STIX]{x1D719}^{+})}{\unicode[STIX]{x2202}z\unicode[STIX]{x2202}n}\frac{\text{d}^{3}\unicode[STIX]{x1D713}_{\tilde{m}}}{\text{d}z^{3}}+\frac{\unicode[STIX]{x2202}^{4}(\unicode[STIX]{x1D719}^{-}-\unicode[STIX]{x1D719}^{+})}{\unicode[STIX]{x2202}z^{3}\unicode[STIX]{x2202}n}\frac{\text{d}\unicode[STIX]{x1D713}_{\tilde{m}}}{\text{d}z}\right]_{z=0}.\qquad \quad\end{eqnarray}$$ In the above two equations, 
$z=0$ should be understood as that the crack is approached from ice sheet not from 
$\unicode[STIX]{x1D6F4}_{i}$, and therefore jumps exist. Substituting (3.5) and (3.6) into (2.3), and then the obtained results into (3.17) and (3.18), we can get
$$\begin{eqnarray}\tilde{\unicode[STIX]{x1D711}}_{\tilde{m}}Q_{\tilde{m}}=\frac{L}{\unicode[STIX]{x1D70C}_{w}\unicode[STIX]{x1D714}^{2}}\left[\tilde{D}\frac{\text{d}^{3}\unicode[STIX]{x1D713}_{\tilde{m}}}{\text{d}z^{3}}+\frac{\unicode[STIX]{x2202}^{3}(\unicode[STIX]{x1D719}_{D}^{-}-\unicode[STIX]{x1D719}_{D}^{+})}{\unicode[STIX]{x2202}z^{3}}\frac{\text{d}\unicode[STIX]{x1D713}_{\tilde{m}}}{\text{d}z}\right]_{z=0}\end{eqnarray}$$and
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D711}}_{\tilde{m}}}{\unicode[STIX]{x2202}n}Q_{\tilde{m}}=\frac{L}{\unicode[STIX]{x1D70C}_{w}\unicode[STIX]{x1D714}^{2}}\left[\tilde{S}\frac{\text{d}^{3}\unicode[STIX]{x1D713}_{\tilde{m}}}{\text{d}z^{3}}+\frac{\unicode[STIX]{x2202}^{4}(\unicode[STIX]{x1D719}_{D}^{-}-\unicode[STIX]{x1D719}_{D}^{+})}{\unicode[STIX]{x2202}z^{3}\unicode[STIX]{x2202}n}\frac{\text{d}\unicode[STIX]{x1D713}_{\tilde{m}}}{\text{d}z}\right]_{z=0},\end{eqnarray}$$ where the continuity conditions (3.13) and (3.14) have been used in the first terms on the right-hand sides of (3.17) and (3.18), and 
$\tilde{D}$ and 
$\tilde{S}$ are, respectively, the jumps of displacement and slope across the crack or
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{D}=\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{D}^{-}}{\unicode[STIX]{x2202}z}-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{D}^{+}}{\unicode[STIX]{x2202}z}\right)_{z=0}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{S}=\left(\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}_{D}^{-}}{\unicode[STIX]{x2202}n\unicode[STIX]{x2202}z}-\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}_{D}^{+}}{\unicode[STIX]{x2202}n\unicode[STIX]{x2202}z}\right)_{z=0}. & \displaystyle\end{eqnarray}$$ It may be noticed that in (3.19) and (3.20), the continuity of 
$\unicode[STIX]{x1D711}_{I}$ in (3.6) across the artificial vertical surface 
$\unicode[STIX]{x1D6F4}_{i}$ has been implied. Since 
$\unicode[STIX]{x1D719}_{D}$ should satisfy the Laplace equation in (2.4) or 
$\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}_{D}/\unicode[STIX]{x2202}z^{2}=-\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}_{D}$, we may further write (3.19) and (3.20) as
$$\begin{eqnarray}\tilde{\unicode[STIX]{x1D711}}_{\tilde{m}}Q_{\tilde{m}}=\frac{L}{\unicode[STIX]{x1D70C}_{w}\unicode[STIX]{x1D714}^{2}}\left[\tilde{D}\frac{\text{d}^{3}\unicode[STIX]{x1D713}_{\tilde{m}}}{\text{d}z^{3}}-\frac{\text{d}\unicode[STIX]{x1D713}_{\tilde{m}}}{\text{d}z}\unicode[STIX]{x1D6FB}^{2}\frac{\unicode[STIX]{x2202}\left(\unicode[STIX]{x1D719}_{D}^{-}-\unicode[STIX]{x1D719}_{D}^{+}\right)}{\unicode[STIX]{x2202}z}\right]_{z=0}\end{eqnarray}$$and
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D711}}_{\tilde{m}}}{\unicode[STIX]{x2202}n}Q_{\tilde{m}}=\frac{L}{\unicode[STIX]{x1D70C}_{w}\unicode[STIX]{x1D714}^{2}}\left\{\tilde{S}\frac{\text{d}^{3}\unicode[STIX]{x1D713}_{\tilde{m}}}{\text{d}z^{3}}-\frac{\text{d}\unicode[STIX]{x1D713}_{\tilde{m}}}{\text{d}z}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}n}\left[\unicode[STIX]{x1D6FB}^{2}\frac{\unicode[STIX]{x2202}(\unicode[STIX]{x1D719}_{D}^{-}-\unicode[STIX]{x1D719}_{D}^{+})}{\unicode[STIX]{x2202}z}\right]\right\}_{z=0}.\end{eqnarray}$$ Since both 
$\unicode[STIX]{x1D719}^{-}$ and 
$\unicode[STIX]{x1D719}^{+}$ satisfy the ice sheet edge conditions in (2.7), their difference should also satisfy the same edge conditions. By further noticing that
$$\begin{eqnarray}{\mathcal{B}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{I}^{-}}{\unicode[STIX]{x2202}z}\right)={\mathcal{B}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{I}^{+}}{\unicode[STIX]{x2202}z}\right)\quad \text{and}\quad {\mathcal{S}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{I}^{-}}{\unicode[STIX]{x2202}z}\right)={\mathcal{S}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{I}^{+}}{\unicode[STIX]{x2202}z}\right),\end{eqnarray}$$we can obtain
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}^{2}\frac{\unicode[STIX]{x2202}(\unicode[STIX]{x1D719}_{D}^{-}-\unicode[STIX]{x1D719}_{D}^{+})}{\unicode[STIX]{x2202}z}=\unicode[STIX]{x1D708}_{0}\left(\frac{\unicode[STIX]{x2202}^{2}\tilde{D}}{\unicode[STIX]{x2202}s^{2}}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}s}\tilde{S}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}n}\left[\unicode[STIX]{x1D6FB}^{2}\frac{\unicode[STIX]{x2202}(\unicode[STIX]{x1D719}_{D}^{-}-\unicode[STIX]{x1D719}_{D}^{+})}{\unicode[STIX]{x2202}z}\right]=-\unicode[STIX]{x1D708}_{0}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\left(\frac{\unicode[STIX]{x2202}\tilde{S}}{\unicode[STIX]{x2202}s}-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}s}\frac{\unicode[STIX]{x2202}\tilde{D}}{\unicode[STIX]{x2202}s}\right). & \displaystyle\end{eqnarray}$$Substituting (3.26) and (3.27) into (3.23) and (3.24), we have
$$\begin{eqnarray}\tilde{\unicode[STIX]{x1D711}}_{\tilde{m}}=\frac{L\unicode[STIX]{x1D705}_{\tilde{m}}\tanh (\unicode[STIX]{x1D705}_{\tilde{m}}H)}{Q_{\tilde{m}}\unicode[STIX]{x1D70C}_{w}\unicode[STIX]{x1D714}^{2}}\left[\unicode[STIX]{x1D705}_{\tilde{m}}^{2}\tilde{D}-\unicode[STIX]{x1D708}_{0}\left(\frac{\unicode[STIX]{x2202}^{2}\tilde{D}}{\unicode[STIX]{x2202}s^{2}}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}s}\tilde{S}\right)\right],\end{eqnarray}$$and
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D711}}_{\tilde{m}}}{\unicode[STIX]{x2202}n}=\frac{L\unicode[STIX]{x1D705}_{\tilde{m}}\tanh (\unicode[STIX]{x1D705}_{\tilde{m}}H)}{Q_{\tilde{m}}\unicode[STIX]{x1D70C}_{w}\unicode[STIX]{x1D714}^{2}}\left[\unicode[STIX]{x1D705}_{\tilde{m}}^{2}\tilde{S}+\unicode[STIX]{x1D708}_{0}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\left(\frac{\unicode[STIX]{x2202}\tilde{S}}{\unicode[STIX]{x2202}s}-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}s}\frac{\unicode[STIX]{x2202}\tilde{D}}{\unicode[STIX]{x2202}s}\right)\right].\end{eqnarray}$$Substituting (3.28) and (3.29) into (3.10), we can obtain
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D711}_{m} & = & \displaystyle \unicode[STIX]{x1D712}_{m}\mathop{\sum }_{i=1}^{N}\int _{\unicode[STIX]{x1D6E4}_{i}}\left\{G_{m}\left[\unicode[STIX]{x1D705}_{m}^{2}\tilde{S}+\unicode[STIX]{x1D708}_{0}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\left(\frac{\unicode[STIX]{x2202}\tilde{S}}{\unicode[STIX]{x2202}s}-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}s}\frac{\unicode[STIX]{x2202}\tilde{D}}{\unicode[STIX]{x2202}s}\right)\right]\right.\nonumber\\ \displaystyle & & \displaystyle -\left.\frac{\unicode[STIX]{x2202}G_{m}}{\unicode[STIX]{x2202}n_{q}}\left[\unicode[STIX]{x1D705}_{m}^{2}\tilde{D}-\unicode[STIX]{x1D708}_{0}\left(\frac{\unicode[STIX]{x2202}^{2}\tilde{D}}{\unicode[STIX]{x2202}s^{2}}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}s}\tilde{S}\right)\right]\right\}\,\text{d}s_{q},\end{eqnarray}$$where
$$\begin{eqnarray}\unicode[STIX]{x1D712}_{m}=\frac{L\unicode[STIX]{x1D705}_{m}\tanh (\unicode[STIX]{x1D705}_{m}H)}{\unicode[STIX]{x1D6FC}Q_{m}\unicode[STIX]{x1D70C}_{w}\unicode[STIX]{x1D714}^{2}}.\end{eqnarray}$$We may apply integration by parts to the second term on the right-hand side of (3.30) and use (2.13), which provides
$$\begin{eqnarray}\unicode[STIX]{x1D711}_{m}=\unicode[STIX]{x1D712}_{m}\mathop{\sum }_{i=1}^{N}\int _{\unicode[STIX]{x1D6E4}_{i}}(f_{m}^{D}\tilde{D}+f_{m}^{S}\tilde{S})\,\text{d}s_{q},\end{eqnarray}$$where
$$\begin{eqnarray}\displaystyle & \displaystyle f_{m}^{D}=\frac{\unicode[STIX]{x2202}G_{m}}{\unicode[STIX]{x2202}s_{q}}\unicode[STIX]{x1D708}_{0}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}s}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}-\frac{\unicode[STIX]{x2202}G_{m}}{\unicode[STIX]{x2202}n_{q}}\left(\unicode[STIX]{x1D705}_{m}^{2}-\unicode[STIX]{x1D708}_{0}\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}s^{2}}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle f_{m}^{S}=G_{m}\unicode[STIX]{x1D705}_{m}^{2}+\frac{\unicode[STIX]{x2202}G_{m}}{\unicode[STIX]{x2202}n_{q}}\unicode[STIX]{x1D708}_{0}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}s}-\frac{\unicode[STIX]{x2202}G_{m}}{\unicode[STIX]{x2202}s_{q}}\unicode[STIX]{x1D708}_{0}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}. & \displaystyle\end{eqnarray}$$ Taking normal derivatives in (3.32) at point 
$p$, we further have
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D711}_{m}}{\unicode[STIX]{x2202}n}=\unicode[STIX]{x1D712}_{m}\mathop{\sum }_{i=1}^{N}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}n_{p}}\int _{\unicode[STIX]{x1D6E4}_{i}}(f_{m}^{D}\tilde{D}+f_{m}^{S}\tilde{S})\,\text{d}s_{q}.\end{eqnarray}$$ Now, the integral equations contain the unknown jumps of displacement and slope across the crack. To solve 
$\tilde{D}$ and 
$\tilde{S}$, we may substitute equations (3.32) into (3.5) directly to obtain 
$\unicode[STIX]{x1D719}_{D}$ and (3.35) into the normal derivative of (3.5) to obtain 
$\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{D}/\unicode[STIX]{x2202}n$, and then apply the crack conditions in (2.7), which gives two sets of differential equations, or
$$\begin{eqnarray}\mathop{\sum }_{m=-2}^{\infty }\left[\left(\frac{\text{d}^{3}\unicode[STIX]{x1D713}_{m}}{\text{d}z^{3}}+\unicode[STIX]{x1D708}_{0}\frac{\text{d}\unicode[STIX]{x1D713}_{m}}{\text{d}z}\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}s^{2}}\right)\unicode[STIX]{x1D711}_{m}+\unicode[STIX]{x1D708}_{0}\frac{\text{d}\unicode[STIX]{x1D713}_{m}}{\text{d}z}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}s}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D711}_{m}}{\unicode[STIX]{x2202}n}\right]={\mathcal{B}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{I}}{\unicode[STIX]{x2202}z}\right),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle \mathop{\sum }_{m=-2}^{\infty }\left[\unicode[STIX]{x1D708}_{0}\frac{\text{d}\unicode[STIX]{x1D713}_{m}}{\text{d}z}\left(\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}s^{2}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}s}\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}s^{2}}\right)\unicode[STIX]{x1D711}_{m}+\left(\frac{\text{d}^{3}\unicode[STIX]{x1D713}_{m}}{\text{d}z^{3}}-\unicode[STIX]{x1D708}_{0}\frac{\text{d}\unicode[STIX]{x1D713}_{m}}{\text{d}z}\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}s^{2}}\right)\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D711}_{m}}{\unicode[STIX]{x2202}n}\right]\nonumber\\ \displaystyle & & \displaystyle \quad ={\mathcal{S}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{I}}{\unicode[STIX]{x2202}z}\right),\end{eqnarray}$$ for 
$\boldsymbol{r}\rightarrow \unicode[STIX]{x1D730}_{i}$ and 
$z=0$.
4 Numerical discretization
To find the solution for the boundary integral equations (3.32) and (3.35) as well as the differential equations (3.36) and (3.37) numerically, we may divide the 
$i$th ice crack into 
$M_{i}$ straight-line segments with length, respectively, as 
$\ell _{i}^{1},\ell _{i}^{2},\ldots ,\ell _{i}^{M_{i}}$. On each segment, the jumps of displacement and slope or 
$\tilde{D}$ and 
$\tilde{S}$ are also assumed to be constant and are taken from those at the centre of the segment, where the crack conditions will be enforced. Here, it may be noticed that the solid angle 
$\unicode[STIX]{x1D6FC}$ in the integral equation (3.10) should always be 
$2\unicode[STIX]{x03C0}$, even on the crack, due to that the computational domain exists on both sides of the crack, which is different from that for wave interactions with a body of non-zero thickness. For the first and second partial derivatives with respect to 
$s$, the three-point finite difference scheme for unequal segment will be used (Burden & Faires Reference Burden and Faires2010), or at the 
$i$th segment centre with 
$s=s_{1}$,
$$\begin{eqnarray}\left.\frac{\unicode[STIX]{x2202}f_{i}}{\unicode[STIX]{x2202}s}\right|_{s=s_{1}}=\frac{-h_{2}^{2}\,f_{0}-(h_{1}-h_{2})(h_{1}+h_{2})f_{1}+h_{1}^{2}\,f_{2}}{h_{1}h_{2}(h_{1}+h_{2})}\end{eqnarray}$$and
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}^{2}f_{i}}{\unicode[STIX]{x2202}s^{2}}=2\frac{h_{2}\,f_{0}-(h_{1}+h_{2})f_{1}+h_{1}\,f_{2}}{h_{1}h_{2}(h_{1}+h_{2})},\end{eqnarray}$$ where 
$h_{1}=s_{1}-s_{0}$ and 
$h_{2}=s_{2}-s_{1}$, 
$s_{0}$ and 
$s_{2}$ are the values of 
$s$ at the previous and next segments, respectively. In (4.1) and (4.2), the subscript in 
$f_{i}$ indicates the values at the 
$i$th segment of the crack. Near the edge of an open crack, if there is no required point beyond edge, the one sided three-point equation can be used, or
$$\begin{eqnarray}\displaystyle & \displaystyle \left.\frac{\unicode[STIX]{x2202}f_{i}}{\unicode[STIX]{x2202}s}\right|_{s=s_{0}}=\frac{-h_{2}(2h_{1}+h_{2})f_{0}+(h_{1}+h_{2})^{2}f_{1}-h_{1}^{2}\,f_{2}}{h_{1}h_{2}(h_{1}+h_{2})}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left.\frac{\unicode[STIX]{x2202}f_{i}}{\unicode[STIX]{x2202}s}\right|_{s=s_{2}}=\frac{h_{2}^{2}\,f_{0}-(h_{1}+h_{2})^{2}f_{1}+h_{1}(h_{1}+2h_{2})f_{2}}{h_{1}h_{2}(h_{1}+h_{2})}. & \displaystyle\end{eqnarray}$$ It should be noticed that in (3.35), the normal derivative with respect to 
$\boldsymbol{n}_{p}$ can be moved into the integral. As a result, there is a hypersingular integral which should be computed in the sense of the Hadamard finite part (Martin & Rizzo Reference Martin and Rizzo1989). Substituting (3.32) and (3.35) into (3.36) and (3.37), and using (4.1) and (4.2) for the first and second partial derivatives with respect to 
$s$, we can obtain two sets of linear equations, the total number of which is the same as that for 
$\tilde{D}$ and 
$\tilde{S}$.
When the field point 
$p$ tends to the source point 
$q$, there will be a logarithmic singularity in the Green function 
$G_{m}$, as (Abramowitz & Stegun Reference Abramowitz and Stegun1965)
$$\begin{eqnarray}H_{0}^{(2)}(\unicode[STIX]{x1D705}_{m}R)\rightarrow \frac{2}{\unicode[STIX]{x03C0}\text{i}}\ln (R)+O(1)\quad \text{as }R\rightarrow 0.\end{eqnarray}$$ In such a case, the integrals in (3.32) and (3.35) should be computed in the principle value or Hadamard finite part sense. For the logarithmic singularity, we may separate it out from 
$G_{m}$ and integrate it analytically (Linton & Mciver Reference Linton and Mciver2001). To compute the Hadamard finite part integral, we may denote 
$\boldsymbol{n}_{p}=\left(n_{1}^{p},n_{2}^{p}\right)$ and 
$\boldsymbol{n}_{q}=(n_{1}^{q},n_{2}^{q})$. Then we have
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}^{2}G_{m}}{\unicode[STIX]{x2202}n_{p}\unicode[STIX]{x2202}n_{q}} & = & \displaystyle \frac{\unicode[STIX]{x03C0}\unicode[STIX]{x1D705}_{m}}{2\text{i}R}H_{1}^{(2)}(\unicode[STIX]{x1D705}_{m}R)(\boldsymbol{n}_{p}\boldsymbol{\cdot }\boldsymbol{n}_{q})+\frac{\unicode[STIX]{x03C0}\unicode[STIX]{x1D705}_{m}}{2\text{i}R^{2}}\left[\unicode[STIX]{x1D705}_{m}H_{0}^{(2)}(\unicode[STIX]{x1D705}_{m}R)-\frac{2}{R}H_{1}^{(2)}(\unicode[STIX]{x1D705}_{m}R)\right]\nonumber\\ \displaystyle & & \displaystyle \times \,[(\unicode[STIX]{x1D709}-x)n_{1}^{p}+(\unicode[STIX]{x1D702}-y)n_{2}^{p}][(\unicode[STIX]{x1D709}-x)n_{1}^{q}+(\unicode[STIX]{x1D702}-y)n_{2}^{q}].\end{eqnarray}$$ For 
$p$ and 
$q$ on the same segment, by noticing 
$\boldsymbol{n}_{p}=\boldsymbol{n}_{q}$ and 
$(\unicode[STIX]{x1D709}-x)n_{1}^{p}+(\unicode[STIX]{x1D702}-y)n_{2}^{p}=0$, we may further simplify the above equation as
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}^{2}G_{m}}{\unicode[STIX]{x2202}n_{p}\unicode[STIX]{x2202}n_{q}}=\frac{\unicode[STIX]{x03C0}\unicode[STIX]{x1D705}_{m}}{2\text{i}R}H_{1}^{(2)}(\unicode[STIX]{x1D705}_{m}R).\end{eqnarray}$$ We may also introduce a local coordinate system 
$t\in [-1,+1]$. On each panel of crack 
$i$, we have 
$R=\ell _{i}|t|/2$. Thus
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}^{2}G_{m}}{\unicode[STIX]{x2202}n_{p}\unicode[STIX]{x2202}n_{q}}=\frac{\unicode[STIX]{x03C0}\unicode[STIX]{x1D705}_{m}H_{1}^{(2)}(\unicode[STIX]{x1D705}_{m}\ell _{i}|t|/2)}{\text{i}\ell _{i}|t|}.\end{eqnarray}$$ Similar to (4.5), we have for 
$H_{1}^{(2)}(\unicode[STIX]{x1D705}_{m}R)$ (Abramowitz & Stegun Reference Abramowitz and Stegun1965)
$$\begin{eqnarray}H_{1}^{(2)}(\unicode[STIX]{x1D705}_{m}R)\rightarrow \frac{\unicode[STIX]{x1D705}_{m}R}{\unicode[STIX]{x03C0}\text{i}}\left[\ln (R)-\frac{2}{(\unicode[STIX]{x1D705}_{m}R)^{2}}\right]+O(R)\quad \text{as }R\rightarrow 0.\end{eqnarray}$$Substituting (4.9) into (4.8), we have
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}^{2}G_{m}}{\unicode[STIX]{x2202}n_{p}\unicode[STIX]{x2202}n_{q}}=\frac{4}{\ell _{i}^{2}t^{2}}-\frac{\unicode[STIX]{x1D705}_{m}^{2}}{2}\ln |t|+O(1),\end{eqnarray}$$ for 
$R\rightarrow 0$ or 
$t\rightarrow 0$. Thus, when 
$p$ and 
$q$ are located at the same segment, we may divide the integral of (4.8) into the singular part in (4.10) and the remaining regular part with order 
$O(1)$. The latter one can be computed numerically through the Gauss–Legendre approach. For the former one or the singular part in (4.10), we have
$$\begin{eqnarray}\int _{\ell _{i}}\frac{\unicode[STIX]{x2202}^{2}G_{m}}{\unicode[STIX]{x2202}n_{p}\unicode[STIX]{x2202}n_{q}}\,\text{d}s_{q}=\frac{\ell _{i}}{2}\int _{-1}^{+1}\frac{\unicode[STIX]{x2202}^{2}G_{m}}{\unicode[STIX]{x2202}n_{p}\unicode[STIX]{x2202}n_{q}}\,\text{d}t=\frac{\unicode[STIX]{x1D705}_{m}^{2}\ell _{i}}{2}-\frac{4}{\ell _{i}},\end{eqnarray}$$ where 
$p$ being located at the centre of 
$\ell _{i}$ has been assumed. It might be noticed that the singular term in (4.11) of order two is treated through the Hadamard finite part integral (Martin & Rizzo Reference Martin and Rizzo1989).
5 Results and analysis
To provide meaningful results in physics, the typical values of the parameters of the ice sheet and fluid are taken to be
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle E=5~\text{GPa},\quad \unicode[STIX]{x1D708}=0.3,\quad \unicode[STIX]{x1D70C}_{i}=922.5~\text{kg}~\text{m}^{-3},\quad h=1~\text{m},\\ \displaystyle \unicode[STIX]{x1D70C}_{w}=1025~\text{kg}~\text{m}^{-3},\quad H=100~\text{m},\end{array}\right\}\end{eqnarray}$$ which are the same as those in Squire & Dixon (Reference Squire and Dixon2000), and the considered wavelength 
$\unicode[STIX]{x1D706}_{0}=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D705}_{0}$ will vary from 
$10~\unicode[STIX]{x03C0}\text{m}$ to 
$400~\unicode[STIX]{x03C0}\text{m}$. The length scales of closed and open cracks are chosen to be similar to that in Li et al. (Reference Li, Wu and Shi2018c) and to that in Porter & Evans (Reference Porter and Evans2007), respectively. In the following text, numerical results will be provided in dimensionless form, based on three basic parameters: a characteristic length scale, the density of water 
$\unicode[STIX]{x1D70C}_{w}$ and the acceleration due to gravity 
$g=9.80~\text{m}~\text{s}^{-2}$. To conduct numerical computations, in addition to the discretization of cracks, the infinite series in (3.5) will be truncated at a finite number 
$T-3$ or the first 
$T$ terms will be kept. Higher accuracy can be achieved through smaller segments and more terms being kept in the eigenfunction expansions. It might be noticed that the computational cost very much depends on the number of segments 
$M=M_{1}+\cdots +M_{N}$ used to discretize the integrals (3.32) and (3.35) and the number of terms 
$T$ used to truncate the infinite series (3.36) and (3.37). Generally, the central processing unit (CPU) time required to construct the coefficient matrix is proportional to the square of 
$M$ and linearly to 
$T$. For the case studies here, 
$M_{i}=500$ and 
$T=100$ with 
$N$ up to 4 are found to be enough to provide convergent results for all the considered crack shapes and wavelengths. For 
$N=1$, the developed sequential FORTRAN 90 code WISPICE requires approximately 
$24.8~\text{s}$ to obtain the solution for each wavelength on a computer with a CPU such as i7 6700k. When they become very large, further refinement of the numerical method is needed, such as through a domain decomposition method (Wu & Eatock Taylor Reference Wu and Eatock Taylor1995).
Figure 2. Modulus of the jump of displacement 
$\tilde{\unicode[STIX]{x1D702}}/\unicode[STIX]{x1D6EC}$ and that of slope 
$\tilde{\unicode[STIX]{x1D702}}_{n}/\unicode[STIX]{x1D6EC}$ along the crack with different 
$M_{i}$ and 
$T$ at 
$\unicode[STIX]{x1D705}_{0}=3$ and 
$\unicode[STIX]{x1D6FD}=0$. Solid lines, 
$M_{i}=500$ and 
$T=100$; dashed lines, 
$M_{i}=1000$ and 
$T=200$; open circles, analytical solution from Li et al. (Reference Li, Wu and Shi2018c).
5.1 Convergence study and verification
Computations are first carried out for the wave diffraction by a circular crack with radius 
$a$ in an ice sheet. The case has been investigated by Li et al. (Reference Li, Wu and Shi2018c) and an explicit solution was obtained. Here, 
$a=20~\text{m}$ is taken to be the characteristic length scale, and the same parameters as those in Li et al. (Reference Li, Wu and Shi2018c) have been used, namely 
$\unicode[STIX]{x1D705}_{0}=3$ and 
$\unicode[STIX]{x1D6FD}=0$. Figure 2 shows the modulus of the jump of displacement 
$\tilde{\unicode[STIX]{x1D702}}$ and that of slope 
$\tilde{\unicode[STIX]{x1D702}}_{n}$ across the crack against polar coordinate 
$\unicode[STIX]{x1D703}=\arctan (y/x)$, which are obtained through the kinematic condition on the ice sheet, or
$$\begin{eqnarray}\tilde{\unicode[STIX]{x1D702}}=\frac{\tilde{D}}{\text{i}\unicode[STIX]{x1D714}}\quad \text{and}\quad \tilde{\unicode[STIX]{x1D702}}_{n}=\frac{\tilde{S}}{\text{i}\unicode[STIX]{x1D714}}.\end{eqnarray}$$ Due to symmetry of the problem, only the results for 
$\unicode[STIX]{x1D703}\in [-\unicode[STIX]{x03C0},0]$ are provided. It can be seen from the figure that there is no visible difference between the results obtained by 
$M_{i}=500$, 
$T=100$ and 
$M_{i}=1000$, 
$T=200$, indicating that convergence has been achieved. Here, it may be noticed that the results are obtained through discretizing the crack into segments with equal length. From the figure, it can be also observed that the numerical results are in excellent agreement with the analytical solutions from Li et al. (Reference Li, Wu and Shi2018c), which validates the proposed solution procedure. Further verification of the solution procedure is carried out through the scattered cross-section 
$F$ in (A 4), as shown in figure 3 against wavenumber 
$\unicode[STIX]{x1D705}_{0}$. Numerical computation is carried out through taking 
$M_{i}=500$ and 
$T=100$. It can be seen from this figure that the results obtained by direct integration of 
$U(\unicode[STIX]{x1D703})$ show good agreement with those computed by (A 9), and both these sets of results agree well with the analytical solutions. This again demonstrates the accuracy of the present solution procedure.
Figure 3. Scattered cross-section 
$F$ against wavenumber 
$\unicode[STIX]{x1D705}_{0}$ at 
$\unicode[STIX]{x1D6FD}=0$. Solid line, 
$F$ from (A 4); dashed line, 
$F$ from (A 9); open circles, analytical solution from Li et al. (Reference Li, Wu and Shi2018c).
Figure 4. Modulus of the jump of displacement 
$\tilde{\unicode[STIX]{x1D702}}/\unicode[STIX]{x1D6EC}$ and that of slope 
$\tilde{\unicode[STIX]{x1D702}}_{n}/\unicode[STIX]{x1D6EC}$ along the elliptic crack with different 
$b/a$ at 
$\unicode[STIX]{x1D705}_{0}=2\unicode[STIX]{x03C0}/\sqrt{A}$ and 
$\unicode[STIX]{x1D6FD}=0$ (
$A=\unicode[STIX]{x03C0}ab=400\unicode[STIX]{x03C0}$). Solid lines, 
$b/a=0.5$; dashed lines, 
$b/a=1.0$; dash-dotted lines, 
$b/a=1.5$.
5.2 Wave diffraction by a single closed crack
In this and following sections, the ice thickness 
$h$ will be taken as the characteristic length scale, and 
$M_{i}=500$ and 
$T=100$ will be used for numerical computations, if it is not specifically specified. We first investigate the problem of flexural-gravity wave interactions with a closed crack, and the effect of its shapes on the scattered wave field. For this purpose, closed cracks with elliptic shape are considered, i.e. defined as 
$(x/a)^{2}+(y/b)^{2}=1$, where 
$a$ and 
$b$ are the half-axes in the 
$x$ and 
$y$ directions, respectively. To demonstrate how wave diffraction by the crack will be affected when its shape changes and also to get a better understanding of flow physics for this type of problem, three crack shapes are chosen, namely 
$b/a=0.5,1.0,1.5$ with its confined area kept fixed at 
$A=400\unicode[STIX]{x03C0}$. For discretization, the nodes of the crack are distributed in the following form:
$$\begin{eqnarray}x^{j}=a\cos \left[\frac{(2j-M)\unicode[STIX]{x03C0}}{M}\right]\quad \text{and}\quad y^{j}=b\sin \left[\frac{(2j-M)\unicode[STIX]{x03C0}}{M}\right]\quad \text{with }j=0,1,\ldots ,M-1.\end{eqnarray}$$ The jumps of displacement 
$\tilde{\unicode[STIX]{x1D702}}$ and slope 
$\tilde{\unicode[STIX]{x1D702}}_{n}$ along the crack are, respectively, shown in figures 4(a) and 4(b) against curvilinear coordinate 
$s/\unicode[STIX]{x1D6FE}$, with the incident wavenumber taken to be 
$\unicode[STIX]{x1D705}_{0}=2\unicode[STIX]{x03C0}/\sqrt{A}$. The wavelength corresponding to this 
$\unicode[STIX]{x1D705}_{0}$ is equal to the characteristic scale of the crack or 
$\sqrt{A}$. Here, the incident wave angle is chosen as 
$\unicode[STIX]{x1D6FD}=0$, and thus only half of the 
$s$ range or 
$s/\unicode[STIX]{x1D6FE}\in [-1,0]$ is plotted due to symmetry. It ought to be mentioned that the results given in this figure are convergent, namely 
$M_{i}=500$ and 
$T=100$ are found to be large enough and their larger values give graphically indistinguishable curves. It can be seen from this figure that the jumps of 
$\tilde{\unicode[STIX]{x1D702}}$ and 
$\tilde{\unicode[STIX]{x1D702}}_{n}$ at the crack are very much affected by the shape of the crack, although the areas enclosed by crack are the same. Typically, 
$\tilde{\unicode[STIX]{x1D702}}$ and 
$\tilde{\unicode[STIX]{x1D702}}_{n}$ will oscillate with 
$s/\unicode[STIX]{x1D6FE}$, and the number of oscillations and their peaks and troughs are different when the crack shape varies.
Figure 5. Modulus of the diffraction coefficient 
$U$ as a function of 
$\unicode[STIX]{x1D703}$ and 
$\unicode[STIX]{x1D705}_{0}$ for elliptic cracks with different 
$b/a$ at 
$\unicode[STIX]{x1D6FD}=0$ (
$A=\unicode[STIX]{x03C0}ab=400\unicode[STIX]{x03C0}$). (a) 
$b/a=0.5$; (b) 
$b/a=1.0$; (c) 
$b/a=1.5$.
As given in (A 1), the far field diffracted wave field may be described by a diffraction coefficient 
$U(\unicode[STIX]{x1D703})$, which can be obtained through an integral involving jumps 
$\tilde{D}$ and 
$\tilde{S}$ along the crack, as shown in (A 3). In figure 5, 
$U$ is sketched through a contour plot as a function of polar coordinate 
$\unicode[STIX]{x1D703}$ and wavenumber 
$\unicode[STIX]{x1D705}_{0}$. It can be seen from this figure that there is a strong variation of the far field diffracted waves for cracks with different shapes. It can be also seen from figure 5 that when 
$\unicode[STIX]{x1D705}_{0}\rightarrow 0$, there will be no diffracted wave by the crack, or the ice sheet is equivalent to a continuous one. From a mathematical aspect, when 
$\unicode[STIX]{x1D705}_{0}\rightarrow 0$ or 
$\unicode[STIX]{x1D714}\rightarrow 0$ the ice sheet condition in (2.6) will become 
$\unicode[STIX]{x2202}\unicode[STIX]{x1D719}/\unicode[STIX]{x2202}z=0$. This indicates that the free edge condition in (2.7) at the crack and no concentrated force condition in (2.13) will be automatically satisfied, or the ice sheet will be ‘continuous’. From a physical aspect, when the wavelength 
$\unicode[STIX]{x1D706}_{0}=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D705}_{0}$ tends to infinity, the dynamic effect of wave motion will disappear, and so does the disturbance to the ice sheet, either with a crack or without. When 
$\unicode[STIX]{x1D705}_{0}$ increases, the diffraction coefficient 
$U$ becomes more pronounced, and for a given 
$\unicode[STIX]{x1D705}_{0}$ there is an obvious difference of the distribution of 
$U$ against 
$\unicode[STIX]{x1D703}$ for a crack with a different shape. A common feature is that they all oscillate against 
$\unicode[STIX]{x1D703}$, and a larger 
$\unicode[STIX]{x1D705}_{0}$ will lead to more oscillatory results. This is because more wave cycles will be encountered by a given crack when the wavelength decreases. It can be also seen from figure 5 that for all these three situations and various 
$\unicode[STIX]{x1D705}_{0}$, a large part of the diffracted wave energy appears near the direction of the incident wave. From this figure, we can further observe that when 
$\unicode[STIX]{x1D703}$ is fixed, 
$U$ will oscillate with respect to 
$\unicode[STIX]{x1D705}_{0}$. This can be more clearly seen from figure 6 which shows the 2-D curves of scattered cross-section 
$F$ against 
$\unicode[STIX]{x1D705}_{0}$. For the specific range of 
$\unicode[STIX]{x1D705}_{0}$, three large peaks of 
$F$ can be observed in addition to some local small peaks for different crack shapes. This is one of the typical features for a closed crack. It is well known that for a finite elastic plate in a vacuum, after an initial disturbance of the plate, it will oscillate forever at its natural frequencies if the structural damping effect is ignored. The oscillation amplitude will tend to infinity if the frequency of the continuous exciting force is equal to one of the natural frequencies of the plate. However, for an ice plate inside the crack, when there is fluid in its lower side, its oscillations will generate water waves. Part of energy will transmit through the crack and propagate to infinity, which provides a damping through wave radiation. The peaks at the natural frequencies which will be different from those in the vacuum may be very large, instead of infinity. When the geometry of an ice plate varies, its natural frequencies can be expected to be different (Smith, Meylan & Mcphedran Reference Smith, Meylan and Mcphedran2011). Therefore, for a closed crack with various shapes, its near-natural frequencies will be different or the large peaks will occur at different wavenumbers, as shown in figure 6.
Figure 6. Scattered cross-section 
$F$ against wavenumber 
$\unicode[STIX]{x1D705}_{0}$ for elliptic cracks with different 
$b/a$ at 
$\unicode[STIX]{x1D6FD}=0$ (
$A=\unicode[STIX]{x03C0}ab=400\unicode[STIX]{x03C0}$). Solid line, 
$b/a=0.5$; dashed line, 
$b/a=1.0$; dash-dotted line, 
$b/a=1.5$.
5.3 Wave diffraction by a single open crack
Investigations now move to the question of how differently open cracks will affect the diffracted wave fields. To do this, we consider circular cracks with different open angles, i.e. 
$x=a\cos \unicode[STIX]{x1D703}$ and 
$y=a\sin \unicode[STIX]{x1D703}$ with 
$\unicode[STIX]{x1D703}\in [-\bar{\unicode[STIX]{x1D703}},+\bar{\unicode[STIX]{x1D703}}]$. Here 
$\unicode[STIX]{x1D717}=2(\unicode[STIX]{x03C0}-\bar{\unicode[STIX]{x1D703}})$ is the open angle. Three different cracks are chosen, i.e. 
$\unicode[STIX]{x1D717}=\unicode[STIX]{x03C0}/20$, 
$\unicode[STIX]{x03C0}/3$ and 
$39\unicode[STIX]{x03C0}/20$, with their length fixed as 
$40\unicode[STIX]{x03C0}$. It should be noticed that the corresponding radii of these three cracks are 
$800/39$, 
$24$ and 
$800$, respectively. The incident wavenumber is chosen as 
$\unicode[STIX]{x1D705}_{0}=\unicode[STIX]{x03C0}/20$, or the incident wavelength 
$\unicode[STIX]{x1D706}_{0}=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D705}_{0}=40$ which is the diameter of the case with 
$\unicode[STIX]{x1D717}=0$, and the incident wave angle is chosen as 
$\unicode[STIX]{x1D6FD}=0$. From (A 8) in appendix A, we have the asymptotic behaviours 
$\tilde{D}=O(s^{\unicode[STIX]{x1D6FC}})$ and 
$\tilde{S}=O(s^{\unicode[STIX]{x1D6FC}-1})$ near the tip of the open crack, with 
$\unicode[STIX]{x1D6FC}\geqslant 3/2$, where 
$s$ is the distance measured from the tip. From this, it can be expected that their partial derivatives 
$\unicode[STIX]{x2202}^{2}\tilde{D}/\unicode[STIX]{x2202}s^{2}$ and 
$\unicode[STIX]{x2202}\tilde{S}/\unicode[STIX]{x2202}s$ involved in the integral equations (3.32) and (3.35) could vary very sharply near the tips of the crack. Therefore, to accelerate convergence of the results in the computations, the nodes are distributed in a cosine form along the crack, or
$$\begin{eqnarray}s^{j}=-\unicode[STIX]{x1D6FE}+\unicode[STIX]{x1D6FE}[1-\cos (j\unicode[STIX]{x03C0}/M)]\quad \text{with }j=0,1,\ldots ,M.\end{eqnarray}$$ The jumps of displacement and slope at the crack are shown in figure 7 against curvilinear coordinate 
$s/\unicode[STIX]{x1D6FE}$. As a comparison, the results for a closed circular crack are also provided. It can be seen that the jumps due to the diffracted wave at the cracks with the same length but different curvatures differ from each other significantly. Particularly, the results of larger curvature (
$39/800$ and 
$1/24$) are much more oscillatory against 
$s/\unicode[STIX]{x1D6FE}$ than those of very small curvature (
$1/800$) in the figure. The small curvature case is closer to that of a straight-line crack parallel to the 
$y$-axis. When its length is very large and approaches infinity, the problem tends to a 2-D one, and 
$\tilde{\unicode[STIX]{x1D702}}$ and 
$\tilde{\unicode[STIX]{x1D702}}_{n}$ will be constant along the crack, or their variations will be zero. From figure 7, it is also interesting to see that when the opening of a crack is small and its two ends tend to merge to form a closed crack, the diffracted wave field of the former does not tend to that of the latter, as can be seen from the solid lines and dashed lines. In fact, although the open crack has nearly become a closed crack, the physics is still different. For the latter one, the inside ice sheet is fully separated from the outside ice sheet by the crack, or there is no direct connection between the two ice sheets, and their interactions are through the fluid. However, for an open crack, the ice sheets on both sides of the crack remain fully connected, and there will always be a direct interaction between their motions. When the open angle tend to zero, the connection will become a point at which the jumps 
$\tilde{\unicode[STIX]{x1D702}}$ and 
$\tilde{\unicode[STIX]{x1D702}}_{n}$ should still be zero, as shown in figure 7, which is different from that of the closed one.
Figure 7. Modulus of the jump of displacement 
$\tilde{\unicode[STIX]{x1D702}}/\unicode[STIX]{x1D6EC}$ and that of slope 
$\tilde{\unicode[STIX]{x1D702}}_{n}/\unicode[STIX]{x1D6EC}$ along circular cracks with different open angle 
$\unicode[STIX]{x1D717}$ at 
$\unicode[STIX]{x1D705}_{0}=\unicode[STIX]{x03C0}/20$ and 
$\unicode[STIX]{x1D6FD}=0$ (
$\unicode[STIX]{x1D6FE}=20\unicode[STIX]{x03C0}$). Solid lines, 
$\unicode[STIX]{x1D717}=\unicode[STIX]{x03C0}/20$; dashed lines, 
$\unicode[STIX]{x1D717}=\unicode[STIX]{x03C0}/3$; dash-dotted lines, 
$\unicode[STIX]{x1D717}=39\unicode[STIX]{x03C0}/20$; dotted lines, closed circular crack.
Figure 8. Modulus of the diffraction coefficient 
$U$ as a function of 
$\unicode[STIX]{x1D703}$ and 
$\unicode[STIX]{x1D705}_{0}$ for circular cracks with different open angle 
$\unicode[STIX]{x1D717}$ at 
$\unicode[STIX]{x1D6FD}=0$ (
$\unicode[STIX]{x1D6FE}=20\unicode[STIX]{x03C0}$). (a) 
$\unicode[STIX]{x1D717}=\unicode[STIX]{x03C0}/20$; (b) 
$\unicode[STIX]{x1D717}=\unicode[STIX]{x03C0}/3$; (c) 
$\unicode[STIX]{x1D717}=39\unicode[STIX]{x03C0}/20$; (d) closed circular crack.
The effect of an open crack shape on the diffracted wave in far field is demonstrated in figure 8, through the diffraction coefficient 
$U$ in the contour plot as a function of 
$\unicode[STIX]{x1D703}$ and 
$\unicode[STIX]{x1D705}_{0}$. It can be observed from this figure that 
$U$ tends to be zero as 
$\unicode[STIX]{x1D705}_{0}\rightarrow 0$ for the reason explained in § 5.2. When 
$\unicode[STIX]{x1D705}_{0}$ increases, the diffracted wave become more obvious. From the contour plot, we have that the diffracted wave fields differ from each other significantly for cracks with the same length but different shapes. However, similar to the closed crack case, most of the diffracted wave will propagate in the direction of the incident wave 
$\unicode[STIX]{x1D6FD}$ for all situations. For the crack with a very small curvature, there is also a large part of the diffracted wave in the direction of 
$-\unicode[STIX]{x1D6FD}$, as shown in figure 8(c). This is similar to the waves normally incident to an infinitely straight-line crack of zero curvature, in which waves will transmit in the direction 
$0$ and be reflected back in the direction 
$\unicode[STIX]{x03C0}$. As discussed in the above section, there will be near resonant motions when the crack is closed, which will lead to large amplitude diffracted waves or 
$U(\unicode[STIX]{x1D703})$. When the crack is open, no similar resonant behaviours can be observed, even for a very small open angle. This can be seen from the comparisons of figure 8(a) with 8(d), and more clearly be observed from figure 9 which provides the variations of the scattered cross-section coefficient 
$F$ against 
$\unicode[STIX]{x1D705}_{0}$. From the contour plot in figure 8, we also have that when 
$\unicode[STIX]{x1D705}_{0}$ is fixed, 
$U$ will oscillate with 
$\unicode[STIX]{x1D703}$. As 
$\unicode[STIX]{x1D705}_{0}\rightarrow 0.2$, there will be more cycles of 
$U$ with respect to 
$\unicode[STIX]{x1D703}$ for a moderate open angle in figure 8(b).
Figure 9. Scattered cross-section 
$F$ against wavenumber 
$\unicode[STIX]{x1D705}_{0}$ for circular cracks with different open angle 
$\unicode[STIX]{x1D717}$ at 
$\unicode[STIX]{x1D6FD}=0$ (
$\unicode[STIX]{x1D6FE}=20\unicode[STIX]{x03C0}$). Solid line, 
$\unicode[STIX]{x1D717}=\unicode[STIX]{x03C0}/20$; dashed line, 
$\unicode[STIX]{x1D717}=\unicode[STIX]{x03C0}/3$; dash-dotted line, 
$\unicode[STIX]{x1D717}=39\unicode[STIX]{x03C0}/20$; dotted line, closed circular crack.
5.4 Wave diffraction by multiple cracks
The developed methodology in this paper can be applied to any crack shape and arrangement. Here, we will consider a few cases of wave diffraction by multiple cracks. The first case we consider is two straight-line cracks with the same finite length of 
$2\unicode[STIX]{x1D6FE}_{i}=40\unicode[STIX]{x03C0}$, 
$i=1,2$. One is parallel to the 
$y$-axis with its start and end points at 
$(x_{s}^{1},y_{s}^{1})=(-2.5,-20\unicode[STIX]{x03C0})$ and 
$(x_{e}^{1},y_{e}^{1})=(-2.5,20\unicode[STIX]{x03C0})$, respectively. For the other crack, its start point is fixed at 
$(x_{s}^{2},y_{s}^{2})=(2.5,-20\unicode[STIX]{x03C0})$, but its angle with respect to the 
$x$-axis changes, namely 
$\unicode[STIX]{x1D70F}=\unicode[STIX]{x03C0}/2$, 
$\unicode[STIX]{x03C0}/4$ and 
$0$. Both of these cracks are discretized with nodes distributed in the cosine form in (5.4), which provides a very fine mesh near the tips of crack.
Figure 10. Modulus of the diffraction coefficient 
$U$ as a function of 
$\unicode[STIX]{x1D703}$ and 
$\unicode[STIX]{x1D705}_{0}$ for two straight-line cracks with different arrangements at 
$\unicode[STIX]{x1D6FD}=0$. (a) 
$\unicode[STIX]{x1D70F}=\unicode[STIX]{x03C0}/2$; (b) 
$\unicode[STIX]{x1D70F}=\unicode[STIX]{x03C0}/4$; (c) 
$\unicode[STIX]{x1D70F}=0$; (d) single crack.
Figure 11. Scattered cross-section 
$F$ against wavenumber 
$\unicode[STIX]{x1D705}_{0}$ for two straight-line cracks with different arrangements at 
$\unicode[STIX]{x1D6FD}=0$. Solid line, 
$\unicode[STIX]{x1D70F}=\unicode[STIX]{x03C0}/2$; dashed line, 
$\unicode[STIX]{x1D70F}=\unicode[STIX]{x03C0}/4$; dash-dotted line, 
$\unicode[STIX]{x1D70F}=0$; dotted line, single crack.
Figure 10 provides the contour plot of diffraction coefficient 
$U$ as a function of polar coordinate 
$\unicode[STIX]{x1D703}$ and wavenumber 
$\unicode[STIX]{x1D705}_{0}$, for two straight-line cracks with different arrangements. As a comparison, the results for a single crack or with the second crack above being removed are also provided. Here, the incident wave angle is taken as 
$\unicode[STIX]{x1D6FD}=0$. Different from the case in § 5.3, since the cracks are not arranged symmetrically about 
$y=0$, the obtained results are asymmetric, as shown in the figure. From the figure, it can be seen that the diffracted wave field varies differently for cracks with different arrangements, although they all tend to be zero when 
$\unicode[STIX]{x1D705}_{0}\rightarrow 0$ for the reason discussed in § 5.2. When the two cracks are closely arranged, it is expected that there will be stronger interactions between waves generated by the cracks, as can be observed in figures 10(a) and 10(b). This is clearer in figure 11, which shows the variation of the scattered cross-section 
$F$ against 
$\unicode[STIX]{x1D705}_{0}$. For two cracks parallel to each other or 
$\unicode[STIX]{x1D70F}=\unicode[STIX]{x03C0}/2$, it resembles a 2-D problem. There is much strong variation with 
$\unicode[STIX]{x1D705}_{0}$ in this case. When the right-hand crack is arranged parallel to the 
$x$-axis, or 
$\unicode[STIX]{x1D70F}=0$, the diffracted wave in the far field is very similar to that by a single crack. In figure 11, the dash-dotted line and the dotted line, respectively for these two cases, are nearly overlapping each other. In fact, for a crack parallel to 
$x$-axis, by noticing 
$\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}/\unicode[STIX]{x2202}s=0$, 
$\unicode[STIX]{x2202}/\unicode[STIX]{x2202}s=\unicode[STIX]{x2202}/\unicode[STIX]{x2202}x$ and 
$\unicode[STIX]{x2202}/\unicode[STIX]{x2202}n=-\unicode[STIX]{x2202}/\unicode[STIX]{x2202}y$, the operators 
${\mathcal{B}}$ and 
${\mathcal{S}}$ in (2.11) and (2.12) can be rewritten as
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{B}}=\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}y^{2}}+\unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}x^{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{S}}=-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}\left[\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}y^{2}}+(2-\unicode[STIX]{x1D708})\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}x^{2}}\right]. & \displaystyle\end{eqnarray}$$ For an incident wave with 
$\unicode[STIX]{x1D6FD}=0$ the right-hand sides of (3.36) and (3.37) become
$$\begin{eqnarray}{\mathcal{B}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{I}}{\unicode[STIX]{x2202}z}\right)=\unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}^{3}\unicode[STIX]{x1D719}_{I}}{\unicode[STIX]{x2202}x^{2}\unicode[STIX]{x2202}z}\quad \text{and}\quad {\mathcal{S}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{I}}{\unicode[STIX]{x2202}z}\right)=0,\end{eqnarray}$$ which means that only the second term in (5.5) has a contribution to the results. Then, it can be expected that there will be more diffracted waves for 
$\unicode[STIX]{x1D6FD}=\unicode[STIX]{x03C0}/2$ than those for 
$\unicode[STIX]{x1D6FD}=0$. This feature has also been shown in the numerical results of Porter & Evans (Reference Porter and Evans2007), where they considered only the straight and parallel crack lines. For the case in figure 10(c), an incident wave propagates along the 
$x$-axis. It will arrive with normal incident angle at the left-hand crack first, and generate a diffracted wave which will mainly propagate along and opposite to the positive 
$x$-axis, as demonstrated in § 5.3. The new wave field can be seen as an incident wave for the other crack. When it arrives at the right-hand crack which is parallel to 
$x$-axis, the generated diffracted wave will be very small, as the direction of the main component of the incident wave field is parallel to the 
$x$-axis. This indicates that the results for two cracks will be similar to those for a single crack in figure 10(c).
In figure 12, the jumps of displacement and slope at the left-hand crack are shown for various arrangements of the right-hand crack. The shapes and arrangements of the cracks as well as the incident wave angle are the same as those in figure 10, while the incident wavenumber is taken to be 
$\unicode[STIX]{x1D705}_{0}=0.1$. As can be expected, the jumps of all situations tend to be zero as the point approaches the tips of the crack. Similar to 
$U$ or 
$F$ in figures 10 and 11, from figure 12, it can be observed that the jumps of displacement and slope along the crack for 
$\unicode[STIX]{x1D70F}=0$ are also close to those for a single crack.
Figure 12. Modulus of the jump of displacement 
$\tilde{\unicode[STIX]{x1D702}}/\unicode[STIX]{x1D6EC}$ and that of slope 
$\tilde{\unicode[STIX]{x1D702}}_{n}/\unicode[STIX]{x1D6EC}$ along the left-hand crack for different arrangements of the right-hand crack at 
$\unicode[STIX]{x1D705}_{0}=0.1$ and 
$\unicode[STIX]{x1D6FD}=0$. Solid lines, 
$\unicode[STIX]{x1D70F}=\unicode[STIX]{x03C0}/2$; dashed lines, 
$\unicode[STIX]{x1D70F}=\unicode[STIX]{x03C0}/4$; dash-dotted lines, 
$\unicode[STIX]{x1D70F}=0$; dotted lines, single crack.
Figure 13. Modulus of the jump of displacement 
$\tilde{\unicode[STIX]{x1D702}}/\unicode[STIX]{x1D6EC}$ and that of slope 
$\tilde{\unicode[STIX]{x1D702}}_{n}/\unicode[STIX]{x1D6EC}$ against curvilinear coordinate 
$s/\unicode[STIX]{x1D6FE}_{i}$ for multiple mixed closed and open cracks at 
$\unicode[STIX]{x1D705}_{0}=\unicode[STIX]{x03C0}/20$ and 
$\unicode[STIX]{x1D6FD}=\unicode[STIX]{x03C0}/4$. Solid lines, circular crack at 
$(-40,0)$; dashed lines, circular crack at 
$(20,-20)$; dash-dotted lines, straight-line crack at 
$(-10,20)$; dotted lines, straight-line crack at (50, 0).
Figure 14. Modulus of the diffraction coefficient 
$U$ as a function of 
$\unicode[STIX]{x1D703}$ and 
$\unicode[STIX]{x1D705}_{0}$ at 
$\unicode[STIX]{x1D6FD}=\unicode[STIX]{x03C0}/4$. (a) Multiple mixed closed and open cracks; (b) single closed circular crack.
Figure 15. Scattered cross-section 
$F$ against wavenumber 
$\unicode[STIX]{x1D705}_{0}$ at 
$\unicode[STIX]{x1D6FD}=\unicode[STIX]{x03C0}/4$. Solid line, multiple mixed closed and open cracks; dashed line, single closed circular crack.
The cracks of a uniform ice sheet formed by the excitation of a unidirectional regular wave might have a simple shape. Due to the random nature of the ocean waves, the cracks in polar regions may be generally complex, which can be either open, closed or the mixture of them. This can be partly reflected by the photos taken on site as well as in the laboratory (Kohout et al. Reference Kohout, Williams, Toyota, Lieser and Hutchings2016; Herman, Evers & Reimer Reference Herman, Evers and Reimer2018). Also, when there is a high concentration of ice floes with very narrow gaps between the floes, this could be equivalently treated as a crack problem mathematically, as noted by Squire & Dixon (Reference Squire and Dixon2001a). The shapes of these narrow gaps could be highly complex (e.g. figure 2, Herman et al. (Reference Herman, Evers and Reimer2018)). Therefore, in our final case we consider a more general situation, or the flexural-gravity wave interactions with multiple mixed closed and open cracks, i.e. two closed circular cracks with radius 
$a$ as 
$20$ and centre 
$(x,y)$, respectively, at 
$(-40,0)$ and 
$(20,-20)$, and two straight-line cracks parallel to 
$y$-axis with length 
$2\unicode[STIX]{x1D6FE}_{i}$ as 
$40$ and centre 
$(x,y)$, respectively, at 
$(-10,20)$ and 
$(50,0)$. Here, an oblique incident wave with 
$\unicode[STIX]{x1D6FD}=\unicode[STIX]{x03C0}/4$ is considered.
Figure 13 shows the jumps of displacement 
$\tilde{\unicode[STIX]{x1D702}}$ and slope 
$\tilde{\unicode[STIX]{x1D702}}_{n}$ along the cracks against curvilinear coordinate 
$s/\unicode[STIX]{x1D6FE}_{i}$. Here, the incident wavenumber is taken to be 
$\unicode[STIX]{x1D705}_{0}=\unicode[STIX]{x03C0}/20$, and the corresponding wavelength is equal to the diameter of the circular crack or length of the straight-line crack. Convergence with respect to 
$M_{i}$ and 
$T$ has been verified, and the energy conservation relationship in (A 9) has been used as a partial check. It can be observed from the figure that for the open cracks both 
$\tilde{\unicode[STIX]{x1D702}}$ and 
$\tilde{\unicode[STIX]{x1D702}}_{n}$ will tend to be zero as 
$s\rightarrow \pm \unicode[STIX]{x1D6FE}_{i}$, as expected from the physics. For the specific wavelength or 
$\unicode[STIX]{x1D705}_{0}$, it can be seen that 
$\tilde{\unicode[STIX]{x1D702}}$ and 
$\tilde{\unicode[STIX]{x1D702}}_{n}$ will oscillate with 
$s/\unicode[STIX]{x1D6FE}_{i}$, and more peaks and troughs appear for a closed circular crack than those for a straight- line crack.
In figure 14, contour plots are shown for the diffraction coefficient 
$U$ as a function of 
$\unicode[STIX]{x1D703}$ and 
$\unicode[STIX]{x1D705}_{0}$. It can be observed that similar to the single crack case, the wave energy is mainly diffracted to the direction of the wave incidence. The diffraction coefficient for a single circular closed crack is also shown in figure 14. It is interesting to see that the wavenumber 
$\unicode[STIX]{x1D705}_{0}$ corresponding to the large peaks in figure 14(a) for multiple cracks are similar to those in figure 14(b) for a single crack, although the value of the former case is generally larger than the latter case. This can be more clearly observed from figure 15, which shows the scattered cross-section 
$F$ against 
$\unicode[STIX]{x1D705}_{0}$. Here, three large peaks appear within the range calculated. As discussed in § 5.2, for a closed crack there will be a near resonant wave motions at a set of discrete 
$\unicode[STIX]{x1D705}_{0}$. At these 
$\unicode[STIX]{x1D705}_{0}$, the amplitude of the diffracted wave can be very large. For multiple cracks, there will be interactions between the waves due to different cracks. Thus, the wavenumbers at which large motion may occur may not be exactly the same as those 
$\unicode[STIX]{x1D705}_{0}$ for an individual crack. However, in many cases, these wavenumbers are still close to 
$\unicode[STIX]{x1D705}_{0}$. This is reflected by the results in figure 15.
6 Conclusions
The problem of flexural-gravity wave interactions with multiple cracks in an ice sheet of infinite extent has been solved, in which both the shape and location of the cracks can be arbitrary, and the individual crack can be either open or closed. The problem is first transformed into an infinite number of coupled two-dimensional problems in a horizontal plane. Each two-dimensional problem is converted into an integral equation along the cracks. By writing the potential in a series in the vertical direction and adopting an orthogonal inner product, all the required conditions including those at the crack are satisfied. The integral differential equations are solved using the boundary element method together with the finite difference scheme for the high derivatives along the crack.
For the interaction of wave with a closed crack, the jumps of displacement 
$\tilde{\unicode[STIX]{x1D702}}$ and slope 
$\tilde{\unicode[STIX]{x1D702}}_{n}$ are generally smooth along the crack length. However, there are several cycles of oscillations with peaks and troughs, and the amplitude of the oscillation is not a constant. For a closed crack, the enclosed ice sheet becomes detached from the external one. For an ice plate in a vacuum, it will oscillate forever at one of its natural frequencies after an initial disturbance, if there is no structure damping. When the excitation force has a component at one of the natural frequencies or at resonance, the motion amplitude can be infinite. When the ice plate is on water surface, similar resonant behaviour can occur. However, the motion amplitude will not be infinite because of wave radiation damping. Instead, the motion will show a large amplitude as reflected by the variation of the diffraction coefficient 
$U(\unicode[STIX]{x1D703})$ with 
$\unicode[STIX]{x1D705}_{0}$.
For an open crack, the variations of displacement and slope jumps are less oscillatory. Both jumps tend to zero at the tip of the crack. However, the slope jump changes very rapidly near the tip and its high derivative is singular at the tip. This singularity has to be treated analytically in the integral equation and numerically in the discretization, as has been done successfully in this work. The displacement and slope jumps are not as oscillatory as those of a closed crack. It is interesting to see that even when the two tips of the open crack are very close, the results do not tend to those of a closed one. A single point connection between the internal and external ice sheets make the physics very different.
When there are multiple cracks, both the crack shapes and their arrangements will have effects on the diffracted waves. Interactions between cracks are much stronger in the local field flow than that in the far field flow. This is mainly due to the fact that the evanescent waves of each crack have strong effects only locally. For a straight- line crack, when it is parallel to the local wave direction, its effect on the total flow is small. Its removal does not have much effect on the overall diffraction coefficient.
The developed methodology can be also used to solve flexural-gravity wave diffractions by cracks with corners. In such a case no concentrated force condition in (2.13) should be imposed at the corners. However, the method is limited to the ice sheet with constant physical properties, and the analysis is within the scope of the linear velocity potential theory.
Acknowledgements
This work is supported by the Lloyd’s Register Foundation, to which the authors are most grateful. The Lloyd’s Register Foundation helps to protect life and property by supporting engineering-related education, public engagement and the application of research. This work is also supported by the National Natural Science Foundation of China (grant no. 51709131 and 51879123).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Diffracted wave in far field and the scattered cross-section
After 
$\tilde{D}$ and 
$\tilde{S}$ have been found, similar to Li et al. (Reference Li, Wu and Shi2018c) for wave diffraction by a circular crack, we may write the diffracted velocity potential 
$\unicode[STIX]{x1D719}_{D}$ in the far field as
$$\begin{eqnarray}\lim _{r\rightarrow +\infty }\unicode[STIX]{x1D719}_{D}=\sqrt{\frac{2}{\unicode[STIX]{x03C0}\unicode[STIX]{x1D705}_{0}r}}\text{e}^{-\text{i}(\unicode[STIX]{x1D705}_{0}r-\unicode[STIX]{x03C0}/4)}IU(\unicode[STIX]{x1D703})\unicode[STIX]{x1D713}_{0}(z)/\tanh (\unicode[STIX]{x1D705}_{0}H),\end{eqnarray}$$ where 
$x=r\,\text{cos}\,\unicode[STIX]{x1D703}$ and 
$y=r\,\text{sin}\,\unicode[STIX]{x1D703}$, 
$\unicode[STIX]{x1D709}=r_{0}\,\text{cos}\,\unicode[STIX]{x1D703}_{0}$ and 
$\unicode[STIX]{x1D702}=r_{0}\,\text{sin}\,\unicode[STIX]{x1D703}_{0}$, and the diffraction coefficient 
$U(\unicode[STIX]{x1D703})$ represents the scattered wave energy distribution along the circumferential direction. Substituting (3.32) into (3.5) and noticing
$$\begin{eqnarray}H_{0}^{(2)}(\unicode[STIX]{x1D705}_{m}R)=\sqrt{\frac{2}{\unicode[STIX]{x03C0}\unicode[STIX]{x1D705}_{m}r}}\text{e}^{-\text{i}(\unicode[STIX]{x1D705}_{m}r-\unicode[STIX]{x03C0}/4)}\text{e}^{\text{i}\unicode[STIX]{x1D705}_{m}r_{0}\cos (\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}_{0})}\quad \text{as }r\rightarrow +\infty ,\end{eqnarray}$$we have
$$\begin{eqnarray}\displaystyle U(\unicode[STIX]{x1D703}) & = & \displaystyle \frac{\unicode[STIX]{x03C0}\unicode[STIX]{x1D712}_{0}}{2\text{i}I}\tanh (\unicode[STIX]{x1D705}_{0}H)\mathop{\sum }_{i=1}^{N}\int _{\unicode[STIX]{x1D6E4}_{i}}\left\{\left[\frac{\unicode[STIX]{x2202}G_{\infty }}{\unicode[STIX]{x2202}s_{q}}\unicode[STIX]{x1D708}_{0}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}s}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}-\frac{\unicode[STIX]{x2202}G_{\infty }}{\unicode[STIX]{x2202}n_{q}}\left(\unicode[STIX]{x1D705}_{0}^{2}-\unicode[STIX]{x1D708}_{0}\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}s^{2}}\right)\right]\tilde{D}\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\left(G_{\infty }\unicode[STIX]{x1D705}_{0}^{2}+\frac{\unicode[STIX]{x2202}G_{\infty }}{\unicode[STIX]{x2202}n_{q}}\unicode[STIX]{x1D708}_{0}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}s}-\frac{\unicode[STIX]{x2202}G_{\infty }}{\unicode[STIX]{x2202}s_{q}}\unicode[STIX]{x1D708}_{0}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\right)\tilde{S}\right\}\,\text{d}s_{q},\end{eqnarray}$$ where 
$G_{\infty }=\exp [\text{i}\unicode[STIX]{x1D705}_{0}r_{0}\cos (\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}_{0})]$. Then the scattered cross-section 
$F$ can be defined as follows:
$$\begin{eqnarray}F=\frac{1}{2\unicode[STIX]{x03C0}}\int _{-\unicode[STIX]{x03C0}}^{+\unicode[STIX]{x03C0}}|U(\unicode[STIX]{x1D703})|^{2}\,\text{d}\unicode[STIX]{x1D703},\end{eqnarray}$$which is a measure of the total scattered energy.
It may be noticed that an energy conservation relationship for 
$F$ can be also derived. Applying Green’s identity to 
$\unicode[STIX]{x1D719}$ and its complex conjugate 
$\unicode[STIX]{x1D719}^{\ast }$, we have
$$\begin{eqnarray}\int _{S}\left(\unicode[STIX]{x1D719}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}^{\ast }}{\unicode[STIX]{x2202}n}-\unicode[STIX]{x1D719}^{\ast }\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}n}\right)\,\text{d}s=0,\end{eqnarray}$$ where 
$S$ is comprised of seabed 
$S_{H}$, ice sheet 
$S_{i}$ and a cylindrical surface at infinity 
$S_{\infty }$. Invoking the corresponding boundary conditions together with the Gauss theorem, we have
$$\begin{eqnarray}\displaystyle & & \displaystyle I_{i}+\int _{S_{\infty }}\left(\unicode[STIX]{x1D719}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}^{\ast }}{\unicode[STIX]{x2202}n}-\unicode[STIX]{x1D719}^{\ast }\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}n}\right)\,\text{d}s\nonumber\\ \displaystyle & & \displaystyle \quad +\,\frac{L}{\unicode[STIX]{x1D70C}_{w}\unicode[STIX]{x1D714}^{2}}\int _{\unicode[STIX]{x1D6E4}_{\infty }}\left[\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}^{\ast }}{\unicode[STIX]{x2202}z}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FB}^{2}}{\unicode[STIX]{x2202}n}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}z}-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}z}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FB}^{2}}{\unicode[STIX]{x2202}n}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}^{\ast }}{\unicode[STIX]{x2202}z}\right.\nonumber\\ \displaystyle & & \displaystyle \quad +\left.\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}n\unicode[STIX]{x2202}z}\unicode[STIX]{x1D6FB}^{2}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}^{\ast }}{\unicode[STIX]{x2202}z}-\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}^{\ast }}{\unicode[STIX]{x2202}n\unicode[STIX]{x2202}z}\unicode[STIX]{x1D6FB}^{2}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}z}\right]_{z=0}\,\text{d}s=0,\end{eqnarray}$$where
$$\begin{eqnarray}I_{i}=\frac{2\text{i}\unicode[STIX]{x1D708}_{0}L}{\unicode[STIX]{x1D70C}_{w}\unicode[STIX]{x1D714}^{2}}\mathop{\sum }_{i=1}^{N}\left\{\text{Im}\left[D\left(\frac{\unicode[STIX]{x2202}S^{\ast }}{\unicode[STIX]{x2202}s}-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}s}\frac{\unicode[STIX]{x2202}D^{\ast }}{\unicode[STIX]{x2202}s}\right)+\frac{\unicode[STIX]{x2202}D^{\ast }}{\unicode[STIX]{x2202}s}S\right]_{-\unicode[STIX]{x1D6FE}_{i}}^{+\unicode[STIX]{x1D6FE}_{i}}\right\}_{+}^{-},\end{eqnarray}$$ which is due to the integral over 
$S_{i}$. The crack conditions in (2.7) together with integration by parts with respect to 
$s$ have been used in (A 7). For a closed crack, 
$I_{i}$ in (A 7) is zero explicitly. Based on the argument of the energy density at the crack tips, we have (Norris & Wang Reference Norris and Wang1994)
$$\begin{eqnarray}D=O[(\unicode[STIX]{x1D6FE}_{i}\mp s)^{\unicode[STIX]{x1D6FC}}]\quad \text{and}\quad S=O[(\unicode[STIX]{x1D6FE}_{i}\mp s)^{\unicode[STIX]{x1D6FC}-1}],\quad (s\rightarrow \pm \unicode[STIX]{x1D6FE}_{i},z=0),\end{eqnarray}$$ where 
$\unicode[STIX]{x1D6FC}\geqslant 3/2$. Thus, for an open crack, by noticing the asymptotic expressions of 
$D$ and 
$S$ in (A 8), 
$I_{i}$ in (A 7) is also zero. Then following the similar procedure in Li et al. (Reference Li, Wu and Shi2018c) for a closed circular crack, we have
$$\begin{eqnarray}F=-\text{Re}[U(\unicode[STIX]{x1D6FD})],\end{eqnarray}$$which can be used as a partial check of the accuracy of the numerical results.
