Interaction between a uniform current and a submerged cylinder in a marginal ice zone

Abstract

The interaction between a uniform current with a circular cylinder submerged in a fluid covered by a semi-infinite ice sheet is considered analytically. The ice sheet is modelled as an elastic thin plate, and the fluid flow is described by the linearised velocity potential theory. The Green function or the velocity potential due to a source is first obtained. As the water surface is divided into two semi-infinite parts with different boundary conditions, the Wiener–Hopf method (WHM) offers significant advantages over alternative approaches and is consequently adopted. To do that, the distribution of the roots of the dispersion equation for fluid fully covered by an ice sheet in the complex plane is first analysed systematically, which does not seem to have been done before. The variations of these roots with the Froude number are investigated, especially their effects or factorisation and decomposition required in the WHM. The result is verified by comparing with that obtained from the matched eigenfunction expansion method. Through differentiating the Green function with respect to the source position, the potentials due to multipoles are obtained, which are employed to construct the velocity potential for the circular cylinder. Extensive results are provided for hydrodynamic forces on the cylinder and wave profiles, and some unique features are discussed. In particular, it is found that the forces can be highly oscillatory with the Froude number when the body is below the ice sheet, whereas such an oscillation does not exist when the body is below the free surface.

1. Introduction

The hydrodynamic problem related to fluid/structure/ice interaction is highly important for the polar region, as well as other cold areas where the water surface may be covered by ice. A better understanding of the physics of the problem is greatly beneficial to environmental protection, engineering operation and safe navigation. The present work focuses on the interaction between an incoming current with a cylinder submerged in the fluid covered by a semi-infinite ice sheet, as a representative case study of an incoming current interacting with an underwater obstacle.

When an ice sheet is of very large horizontal extent, elasticity plays a very important role (Robin 1963; Squire et al. 1988), and in many cases the ice sheet can be modelled as a Kirchhoff–Love plate. When a free surface wave propagates into a semi-infinite ice sheet, or vice versa, there will be a major change of physical properties on the upper surface of the fluid. Wave transmission and reflection will occur. By employing the linearised velocity potential theory for fluid flow, Fox & Squire (1990, 1994) solved the finite water depth problems of normal and oblique incident wave interactions with a semi-infinite ice sheet through the method of matched eigenfunction expansion (MEE), with the free ice edge conditions. Later, Sahoo, Yip & Chwang (2001) extended it to a semi-infinite ice sheet with various edge conditions, such as free, simply supported and clamped. In their work, an inner product with orthogonality was defined to match the solution on the interface. In addition to MEE, this mixed boundary value problem can be also solved by the Wiener–Hopf technique, as in an early work by Evans & Davies (1968), and those by Tkacheva (2001) and Chung & Fox (2002). Linton & Chung (2003) used the residue calculus technique (RCT) for the problem, and confirmed that the solution by MEE is equivalent to that by Wiener–Hopf technique. In addition, the similar problem of infinite water depth was investigated by Chakrabarti (2000), where the Havelock transform (Ursell 1947) was applied to covert the mixed boundary value problem into an integral equation of the Carlemen type over a semi-infinite range. Another interesting hydrodynamic problem in cold regions is wave diffraction by cracks in an ice sheet. Squire & Dixon (2000, 2001) studied the hydroelastic wave propagation in a homogeneous ice sheet with single and multiple infinite length straight cracks floating on fluid with infinite water depth. Evans & Porter (2003) and Porter & Evans (2006) considered the same problems of finite water depth, where the method of vertical mode expansion was used. Later, Porter & Evans (2007) extended it to straight cracks of finite length. In a more recent work, Li, Wu & Ren (2020) further proposed a numerical approach for an ice sheet with multiple cracks of arbitrary shapes. For the problem of two or more ice sheets with different properties and separated by cracks, it can be solved by the Wiener–Hopf method, as done by Marchenko (1993) for a single crack, and by Williams & Squire (2006) for two cracks. In some cases, the effect of the in-plane compressive forces in the ice sheet may need to be included. Das, Sahoo & Meylan (2018) and Barman et al. (2021) considered this effect by incorporating an additional term of the second-order derivatives of deflection into the thin elastic plate model, and identified the phenomenon of wave blocking in such a case.

For wave interaction with a structure, Das & Mandal (2006) investigated the oblique wave scattering by a two-dimensional circular cylinder beneath an ice sheet of infinite extent on the water surface through the Green function method. Later, Li, Wu & Ji (2018c) considered the radiation and diffraction by a circular cylinder submerged below an ice sheet with a crack, where the Green function for an ice sheet with a crack was first derived in an integral form and the solution was obtained by the multipole expansion procedure (Ursell 1949, 1950). Later, Li, Wu & Ji (2018b) extended the work to an ice sheet with multiple cracks. In addition, other similar works about submerged bodies can be also found in Maiti & Mandal (2010) for wave scattering by a thin vertical barrier, as well as by Mondal & Banerjea (2016) for wave diffraction by an inclined porous plate. The interaction between water waves and structures may also occur near the ice edge. In such a case, the ice sheet cannot be treated as infinite. Typically, Sturova (2014) considered the problem of wave radiation by a circular cylinder submerged below a semi-infinite ice sheet. In this work, the Green function was first constructed by the method of MEE, and then the final solution was obtained through the boundary integral equation. Later, this procedure was further employed to solve the radiation of waves by a cylinder submerged in water with ice floe or polynya (Sturova 2015). A similar problem of wave radiation by a submerged elliptic cylinder was studied by Tkacheva (2015) by the Wiener–Hopf technique. In addition to submerged bodies, floating and surface-piercing structures are also common in polar and ocean engineering. Das & Mandal (2009) studied the hydroelastic wave reflection and transmission by a half-immersed circular cylinder confined between two semi-infinite ice sheets by the approach of multipole expansion. Ren, Wu & Thomas (2016) investigated the problem of wave radiation and diffraction by a floating rectangular body in a polynya using the method of MEE. Later, Li, Shi & Wu (2018a) further extended the work to floating bodies with arbitrary shapes by applying the boundary element method.

The studies listed above primarily focus on the periodic waves. In addition, there are also steady waves generated by a body moving forward at constant speed, or a stationary body in a steady incoming current. For the linear free surface problem, Lamb (1924) proposed a first approximation approach for a submerged circular cylinder. Later, Havelock (1936) solved the problem exactly and represented the solution in the form of an infinite series. The linear problem with an ice sheet covering the free surface was considered by Li, Wu & Shi (2019). In their work, the Green function of this steady problem was derived and the multipole expansion procedure was used. Compared with the free surface case, there is a critical Froude number, below which no travelling wave exists away from the body. When the Froude number exceeds the critical value but remains below $1$, a shorter wave is upstream whereas a longer wave is downstream. When the Froude number surpasses $1$, only the shorter wave remains upstream. In this work, we shall consider the problem of a uniform current interaction with a circular cylinder submerged below a semi-infinite ice sheet. In such a case, the derivation of the Green function and multipole expansion becomes far more complex. This problem may be solved by the established Wiener–Hopf technique. However, to locate all the singularities of the dispersion function and decomposition of the complex function are not trivial. Although the study focuses on a circular cylinder, the results do have a much wider application. The analytical formulation also enables us to see many insights into physics related to a submerged body in a current in marginal ice zones.

The paper is arranged as follows. The governing equation and boundary conditions of the velocity potential are presented in § 2. The Green function or the velocity potential due to a single source is derived in § 3.1. The multipoles and the velocity potential due to a submerged circular cylinder are constructed and solved in § 3.2. The formulations of the hydrodynamics forces on a cylinder and the corresponding wave elevation are presented in §§ 3.3 and 3.4, respectively. The numerical results are shown in § 4, followed by the conclusions given in § 5.

2. Governing equation and boundary conditions

The problem of a uniform current interaction with a submerged circular cylinder below a semi-infinite ice sheet is given in figure 1. A Cartesian coordinate system $O$–$xz$ is introduced with the origin at the edge of the ice sheet, and its $x$-axis is along the undisturbed mean water surface, and the $z$-axis points vertically upwards. The homogeneous semi-infinite ice sheet is extended from $x=0$ to $x=+\infty$ with density $\rho _i$ and thickness $h_i$. An incoming current comes from $x=+\infty$ to $x=-\infty$ and will be disturbed by the submerged circular cylinder with radius $a$, whose centre is located at $(x_0,z_0)$.

Figure 1. Coordinate system and sketch of the problem.

The fluid with density $\rho$ and mean water depth $H$ is assumed to be incompressible and inviscid, and its motion is assumed to be irrotational. It is assumed that the body is submerged to submergence levels that causes only small deformations for both free surface and ice sheet. Thus, the linearised velocity potential theory can be used. We may first write the total velocity potential as the summation of the potential due to the current and the potential due to the disturbance by the cylinder, or

(2.1)\begin{equation} \varPhi\left(x,z\right)={-}Ux+\lim_{\epsilon\rightarrow0^{-}}\phi\left(x,z\right){\rm e}^{-\epsilon t}, \end{equation}

where $U$ denotes the speed of the incoming current, and $U>0$. Here $\epsilon <0$ is introduced in (2.1), similar to that in Lighthill (1978, pp. 265–268 and pp. 364–366), which means that a disturbance grows from $t=-\infty$ to the present. This helps the Fourier transform to be performed with respect to $x$ as well as its inverse transform, when the solution of the problem is sought. The disturbed velocity potential $\phi (x,z)$ is governed by the Laplace equation in the entire fluid domain as

(2.2)\begin{equation} \nabla^2\phi= \frac{\partial^2\phi}{\partial x^2} + \frac{\partial^2\phi}{\partial z^2}=0,\quad -\infty < x <{+}\infty, \ -H \leq z \leq 0. \end{equation}

From Wehausen & Laitone (1960), the boundary condition on the free surface, or $x<0$, gives

(2.3)\begin{equation} \left(\frac{\partial}{\partial t} - U\frac{\partial}{\partial x}\right)^2\left(\phi {\rm e}^{-\epsilon t}\right)+ g\frac{\partial \left(\phi {\rm e}^{-\epsilon t}\right)}{\partial z} = 0,\quad -\infty < x < 0, \ z=0, \end{equation}

where $g$ is the acceleration due to gravity. From Li et al. (2019), the boundary conditions on the ice sheet or $x>0$ gives

(2.4) \begin{align} &\left(L\frac{\partial^4}{\partial x^4} + m_i\frac{\partial^2}{\partial t^2} +\rho g\right)\frac{\partial \left(\phi {\rm e}^{-\epsilon t}\right)}{\partial z}\nonumber\\ &\quad +\rho\left(\frac{\partial }{\partial t} - U\frac{\partial }{\partial x}\right)^2\left(\phi {\rm e}^{-\epsilon t}\right) = 0,\quad 0 < x <{+}\infty,\ z = 0, \end{align}

where $L=Eh_i^3/[12(1-\nu ^2)]$ represents the flexural rigidity and $m_i=\rho _i h_i$ denotes the mass per unit length of the ice sheet, $E$ and $\nu$ denote the Young's modulus and Poisson ratio, respectively. Equations (2.3) and (2.4) can be further simplified as

(2.5)\begin{gather} \left(\epsilon + U\frac{\partial}{\partial x}\right)^2\phi + g\frac{\partial \phi}{\partial z} = 0,\quad -\infty < x < 0, \ z=0, \end{gather}

(2.6)\begin{gather} \left(L\frac{\partial^4}{\partial x^4} + m_i\epsilon^2 +\rho g\right)\frac{\partial \phi}{\partial z}+\rho\left(\epsilon + U\frac{\partial }{\partial x}\right)^2\phi = 0,\quad 0 < x <{+}\infty,\ z = 0. \end{gather}

The impermeable condition on the surface of the circular cylinder $S_B$ can be written as

(2.7)\begin{equation} \frac{\partial\phi}{\partial n} = Un_x,\quad \text{on}\ S_B, \end{equation}

where $\boldsymbol {n}=(n_x,n_z)$ denotes the unit normal vector of $S_B$. Similarly, the impermeable condition on the seabed can be expressed as

(2.8)\begin{equation} \frac{\partial\phi}{\partial z} = 0,\quad z={-}H. \end{equation}

At the edge of the ice sheet, the free edge is assumed, and its conditions can be written as

(2.9)\begin{equation} \frac{\partial^2\eta}{\partial x^2}=0,\quad \text{and}\quad\frac{\partial^3\eta}{\partial x^3}=0,\quad x = 0^{+},\ z = 0, \end{equation}

where $\eta$ denotes the wave elevation and can be determined from Li et al. (2019)

(2.10)\begin{equation} \eta = \left\{ \begin{array}{@{}ll} \dfrac{U}{g}\dfrac{\partial \phi}{\partial x}, & x<0,\\ \dfrac{L}{\rho g U}\dfrac{\partial^4\phi}{\partial x^3 \partial z}+\dfrac{U}{g}\dfrac{\partial \phi}{\partial x}, & x > 0, \end{array} \right.,\quad z \rightarrow 0^-. \end{equation}

For a non-zero $\epsilon$, $\phi$ will tend to zero at $x\rightarrow \pm \infty$. When $\epsilon =0$ is taken, $\phi$ becomes an oscillatory function at infinity. In such a case, the radiation condition at far field can be expressed as

(2.11)\begin{equation} \frac{\partial \phi}{\partial x} = w_{{\pm}}(x,z)\quad \text{as}\ x\rightarrow \pm\infty, \end{equation}

where $w_{\pm }(x,z)$ represents a wavy function oscillatory with $x$ at $x\rightarrow \pm \infty$. The group velocities of the waves at upstream and downstream are larger and smaller than $U$, respectively.

3. Solution procedure

3.1. Velocity potential due to a single source: the Green function

We may first define the non-dimensional variables based on the density of fluid $\rho$ together with $g$ and $H$. Typically, $\phi ^{\prime }=\phi /H\sqrt {gH}$, $x^{\prime }=x/H$, $z^{\prime }=z/H$ and $\epsilon ^{\prime }=\epsilon \sqrt {H/g}$. For convenience, the primes are omitted in the following. In such a case, (2.5) and (2.6) can be re-expressed as

(3.1)\begin{gather} \left(\epsilon + F\frac{\partial}{\partial x}\right)^2\phi + \frac{\partial \phi}{\partial z} = 0,\quad -\infty < x < 0, \ z=0, \end{gather}

(3.2)\begin{gather} \left(D\frac{\partial^4}{\partial x^4} + M\epsilon^2 +1\right)\frac{\partial \phi}{\partial z}+\left(\epsilon + F\frac{\partial }{\partial x}\right)^2\phi = 0,\quad 0 < x <{+}\infty,\ z = 0, \end{gather}

where $F=U/\sqrt {gH}$ denotes the depth-based Froude number, $D=L/\rho gH^4$ and $M=m_i/\rho H$.

The Green function $G(x,z;x_0,z_0)$ is first introduced, which is the velocity potential at a field point $(x,y)$ induced by a single source at $(x_0,z_0)$. Here $G$ satisfies the following equation

(3.3)\begin{equation} \nabla^2G=2{\rm \pi} \delta\left(x-x_0\right)\delta\left(z-z_0\right),\quad -\infty < x <{+}\infty, \ -1 \leq z \leq 0, \end{equation}

and the boundary conditions in (2.8)–(3.2), where $\delta (x)$ is the Dirac delta function. To find $G$, we may apply the Fourier transform to $G$

(3.4)\begin{equation} \hat{G}(\alpha,z)=\int_{-\infty}^{+\infty}G(x,z){\rm e}^{{\rm i}\alpha x}\,{{\rm d}\kern 0.06em x}. \end{equation}

Here $\alpha$ is extended to the complex plane. As the conditions in (3.1) and (3.2) are divided by the sign of $x$, we may use the Wiener–Hopf technique (Noble 1958) and further introduce the following Fourier transforms

(3.5)\begin{equation} \left. \begin{aligned} \hat{G}_{-}(\alpha,z) & =\int_{-\infty}^{0}G(x,z){\rm e}^{{\rm i}\alpha x}\,{{\rm d}\kern 0.06em x},\\ \hat{G}_{+}(\alpha,z) & =\int_{0}^{+\infty}G(x,z){\rm e}^{{\rm i}\alpha x}\,{{\rm d}\kern 0.06em x}. \end{aligned} \right\} \end{equation}

Because of the presence of $\epsilon$ in (3.1) and (3.2), $G$ decays as $x\rightarrow \pm \infty$ at a rate of ${\rm e}^{-\epsilon _0|x|}$. Here $\epsilon _0>0$ and tends to zero when $\epsilon \rightarrow 0^-$. Following Noble (1958), $\hat {G}_+(\alpha,z)$ is analytical in the region ${\rm Im}\{\alpha \}>-\epsilon _0$ and $\hat {G}_-(\alpha,z)$ is analytical in the region ${\rm Im}\{\alpha \}<\epsilon _0$, respectively. From (3.4) and (3.5), we have

(3.6)\begin{equation} \hat{G}(\alpha,z)=\hat{G}_{+}(\alpha,z) + \hat{G}_{-}(\alpha,z). \end{equation}

Applying (3.4) to (3.3), we obtain

(3.7)\begin{equation} \frac{\partial^2\hat{G}}{\partial z^2} - \alpha^2\hat{G}=2{\rm \pi} {\rm e}^{{\rm i}\alpha x_0}\delta(z-z_0). \end{equation}

Based on the boundary condition in (2.8), (3.7) can be solved as

(3.8)\begin{equation} \hat{G}=A(\alpha)C(z,\alpha) + \frac{2{\rm \pi} {\rm e}^{{\rm i}\alpha x_0}}{\alpha} \left\{ \begin{array}{@{}ll} S(z,\alpha)C(z_0,\alpha), & z > z_0,\\ C(z,\alpha)S(z_0,\alpha), & z < z_0 , \end{array} \right. \end{equation}

where $A(\alpha )$ is an unknown coefficient and

(3.9)\begin{equation} \left. \begin{aligned} S(z,\alpha) & = \sinh\alpha(z+1),\\ C(z,\alpha) & = \cosh\alpha(z+1). \end{aligned} \right\} \end{equation}

Substituting (3.5) and (3.8) into (3.1) and (3.2) and keeping only the leading term of $\epsilon$, we have

(3.10)\begin{align} &D_{+}(\alpha)+D_{-}(\alpha)=A(\alpha)K_1^{\epsilon}(\alpha,F)\nonumber\\ &\qquad +\frac{2{\rm \pi} {\rm e}^{{\rm i}\alpha x_0}}{\alpha}\left(\alpha \cosh\alpha - F^2\alpha^2 \sinh\alpha -2{\rm i}\epsilon F\alpha \sinh\alpha \right)C(z_0,\alpha), \end{align}

(3.11)\begin{align} &F_{+}(\alpha)+F_{-}(\alpha)=A(\alpha)K_2^{\epsilon}(\alpha,F)\nonumber\\ &\qquad +\frac{2{\rm \pi} {\rm e}^{{\rm i}\alpha x_0}}{\alpha}\left[\left(D\alpha^4+1\right)\alpha \cosh\alpha - F^2\alpha^2 \sinh\alpha - 2{\rm i}\epsilon F\alpha \sinh\alpha\right]C(z_0,\alpha), \end{align}

where

(3.12)\begin{equation} K_{i}^{\epsilon}(\alpha,F)=K_{i}(\alpha,\alpha F + {\rm i}\epsilon)\approx K_{i}(\alpha, \alpha F)-{\rm i}\epsilon\text{sgn}(\alpha),\quad i = 1,2, \end{equation}

and

(3.13a)\begin{gather} K_{1}(\alpha,\alpha F)=\alpha \sinh\alpha - \left(\alpha F\right)^2\cosh\alpha, \end{gather}

(3.13b)\begin{gather} K_{2}(\alpha,\alpha F)=\left(D\alpha^4+1\right)\alpha \sinh\alpha - \left(\alpha F \right)^2\cosh\alpha. \end{gather}

In (3.10) and (3.11) $D_{\pm }(\alpha )$ and $F_{\pm }(\alpha )$ are the Fourier transforms of (3.1) and (3.2), respectively, defined in (3.5). It should be noted that the main effect of $\epsilon$ term in (3.12) is that the location of $K_i^{\epsilon }(\alpha,F)=0$ ($i=1, 2$) will be changed slightly from those corresponding to $K_i(\alpha,\alpha F)=0$. When $\alpha$ is a fully complex number, such a change may be trivial. However, when $\alpha$ is a real number, such a change becomes significant, as it will affect the path at the singularities when the inverse Fourier transform is performed, and affect the decomposition in the Wiener–Hopf method. For this reason, the sign function ${\rm sgn}(\alpha )$ is used in (3.12) as the term is significant only when $\alpha$ is real. It should be noted that $\epsilon$ will also slightly change the double root of $K_i(\alpha,\alpha F)=0$ at $\alpha =0$. However, such a change reflected in the integration path of the inverse transform at $\alpha =0$ is actually equivalent to adding a constant term to the Green function (Li et al. 2019), and will not affect the results. Thus, the effect of $\epsilon$ on the roots at $\alpha =0$ is neglected.

From the boundary conditions in (3.1) and (3.2), we have $D_-(\alpha )=F_+(\alpha )=0$. To obtain $D_+(\alpha )$ and $F_-(\alpha )$, we may eliminate $A(\alpha )$ from (3.10) and (3.11), and remove the trivial $\epsilon$ terms, which provides

(3.14)\begin{equation} F_{-}(\alpha) = D_{+}(\alpha)K^{\epsilon}(\alpha,F)-\frac{2{\rm \pi} DF^2\alpha^6}{K_1^{\epsilon}(\alpha, F)}C(z_0,\alpha){\rm e}^{{\rm i}\alpha x_0}, \end{equation}

where

(3.15)\begin{equation} K^{\epsilon}(\alpha,F)=\frac{K_2^{\epsilon}(\alpha, F)}{K_1^{\epsilon}(\alpha, F)}. \end{equation}

Following the procedure of Wiener–Hopf method, in the complex plane of $\alpha$, (3.15) may be factorised as

(3.16)\begin{equation} K^{\epsilon}(\alpha, F)=K^{\epsilon}_{-}(\alpha, F)K^{\epsilon}_{+}(\alpha, F), \end{equation}

where $K_-^{\epsilon }(\alpha,F)$ are $K_+^{\epsilon }(\alpha,F)$ are analytical in their own regions in the complex plane. To do that, we need to find the distribution of the roots of $K_i^{\epsilon }(\alpha,F)=0$ ($i=1,2$).

As shown in McCue & Stump (2000) for $K_1(\alpha,\alpha F)=0$ and Appendix A for $K_2(\alpha,\alpha F)=0$, $K_1(\alpha,\alpha F)=0$ and $K_2(\alpha,\alpha F)=0$ both have an infinite number of roots at $\alpha =\pm \mathcal {k}_m$ ($m=0,1,2,\ldots$) and $\alpha =\pm \kappa _m$ ($m=-1,0,1,\ldots$), respectively. The properties of the roots are related to the Froude number $F$. Typically, for $K_1(\alpha,\alpha F)=0$, $\mathcal {k}_0$ is a positive real root when $F<1$, and $\mathcal {k}_0$ is a purely positive imaginary root between $0$ and ${\rm \pi} {\rm i}/2$ when $F>1$. Here $\mathcal {k}_m$ ($m\geq 1$) are all purely positive imaginary roots, and $\mathcal {k}_m$ is between $m{\rm \pi} {\rm i}$ and $(m{\rm \pi} +{\rm \pi} /2){\rm i}$. For $K_2(\alpha, \alpha F)=0$, there is a critical Froude number $F_c$ (Li et al. 2019). When $0< F< F_c$, $\kappa _{-1}$ and $\kappa _0$ are two complex roots with positive imaginary part, which satisfy $\kappa _0=-\bar {\kappa }_{-1}$ and ${\rm Re}\{\kappa _{-1}\}>0$. When $F_c< F<1$, $\kappa _{-1}$ and $\kappa _0$ become two positive real roots with $\kappa _{-1}>\kappa _0$. When $F>1$, $\kappa _{-1}$ remains to be a positive real root, but $\kappa _{0}$ becomes a purely positive imaginary root between $0$ and ${\rm \pi} {\rm i}/2$. Similar to $\mathcal {k}_m$ ($m\geq 1$), $\kappa _m$ ($m\geq 1$) are all purely negative imaginary roots between $m{\rm \pi} {\rm i}$ and $(m{\rm \pi} +{\rm \pi} /2){\rm i}$.

The root of $K_1^{\epsilon }(\alpha,F)=0$ in (3.12) corresponding to $\pm \mathcal {k}_0$ can be obtained as $\pm \mathcal {k}_0-{\rm i}\epsilon _0^{\prime }$ ($\epsilon _0^{\prime }\rightarrow 0$). Similarly, the roots of $K_2^{\epsilon }(\alpha,F)=0$ in (3.12) corresponding to $\pm \kappa _m$ ($m=-1, 0$) can be expressed as $\pm \kappa _m-{\rm i}\varepsilon _m^{\prime }$ ($\varepsilon _m^{\prime }\rightarrow 0$) respectively. When $\mathcal {k}_0$ and $\kappa _m$ ($m=-1, 0$) are real, we have

(3.17)\begin{equation} \left. \begin{aligned} \epsilon_0^{\prime} & ={-}\epsilon\text{sgn}(\mathcal{k}_0)\left/\frac{\partial K_1(\mathcal{k}_0,\mathcal{k}_0F)}{\partial \alpha} \right. ={-}\epsilon \text{sgn}(-\mathcal{k}_0)\left/\frac{\partial K_1(-\mathcal{k}_0,-\mathcal{k}_0F)}{\partial \alpha},\right.\\ \varepsilon^{\prime}_m & ={-}\epsilon\text{sgn}(\kappa_m)\left/\frac{\partial K_2(\kappa_m,\kappa_m F)}{\partial \alpha} \right.={-}\epsilon\text{sgn}(-\kappa_m)\left/\frac{\partial K_2(-\kappa_m, -\kappa_m F)}{\partial \alpha}\right. . \end{aligned} \right\} \end{equation}

From (3.13a), it can be shown that ${\partial K_1(\mathcal {k}_0,\mathcal {k}_0 F)}/{\partial \alpha }<0$. From (3.17), $\epsilon _0^{\prime }<0$, which means $\textrm {Im}\{\pm \mathcal {k}_0-\textrm {i}\epsilon _0^{\prime }\}>0$. Similarly, from (3.13b), ${\partial K_2(\kappa _{-1},\kappa _{-1}F)}/{\partial \alpha }>0$ and ${\partial K_2(\kappa _0,\kappa _0F)}/{\partial \alpha }<0$ and, therefore, $\varepsilon _{-1}^{\prime }>0$, $\varepsilon _{0}^{\prime }<0$, which gives $\textrm {Im}\{\pm \kappa _{-1}-\textrm {i}\varepsilon _{-1}^{\prime }\}<0$ and $\textrm {Im}\{\pm \kappa _{0}-\textrm {i}\varepsilon _{0}^{\prime }\}>0$. When $\mathcal {k}_0$ and $\kappa _m$ ($m=-1, 0$) are not purely real, $\epsilon _0^{\prime }=0$ and $\varepsilon _m^{\prime }=0$. Thus, we may write

(3.18a)\begin{gather} \left. \begin{array}{ll@{}} \epsilon_0^{\prime} <0, & \text{Im}\{\mathcal{k}_0\}=0\ (0< F<1),\\ \epsilon_0^{\prime} =0, & \text{Im}\{\mathcal{k}_0\}\neq0\ (F>1), \end{array} \right\} \end{gather}

(3.18b)\begin{gather} \left. \begin{array}{ll@{}} \varepsilon_{{-}1}^{\prime} >0, & \text{Im}\{\kappa_{{-}1}\}=0\ (F>F_c),\\ \varepsilon_{{-}1}^{\prime} =0, & \text{Im}\{\kappa_{{-}1}\}\neq0\ (0< F< F_c), \end{array} \right\} \end{gather}

(3.18c)\begin{gather} \left. \begin{array}{ll@{}} \varepsilon_{0}^{\prime} >0, & \text{Im}\{\kappa_{0}\}=0\ (F_c< F<1),\\ \varepsilon_{0}^{\prime} =0, & \text{Im}\{\kappa_{0}\}\neq0\ (0< F< F_c \ \text{and} \ F>1). \end{array} \right\} \end{gather}

In such a case, based on the Weierstrass factorisation, we have

(3.19a)$$\begin{gather} K_{1}^{\epsilon}(\alpha,F)=(1-F)^2\alpha^2\left(1-\frac{\alpha}{\mathcal{k}_0 - {\rm i}\epsilon_0^{\prime}}\right)\left(1+\frac{\alpha}{\mathcal{k}_0 + {\rm i}\epsilon_0^{\prime}}\right)\prod_{m=1}^{+\infty}\left(1-\frac{\alpha^2}{\mathcal{k}_m^2}\right), \end{gather}$$

(3.19b) $$\begin{gather}K_{2}^{\epsilon}(\alpha,F)=(1-F)^2\alpha^2\Bigg(1-\frac{\alpha}{\kappa_{{-}1} - {\rm i}\varepsilon_{{-}1}^{\prime}}\Bigg)\Bigg(1+\frac{\alpha}{\kappa_{{-}1} + {\rm i}\varepsilon_{{-}1}^{\prime}}\Bigg)\nonumber\\ \qquad\qquad\qquad\qquad\ \ \times\left(1-\frac{\alpha}{\kappa_{0} - {\rm i}\varepsilon_{0}^{\prime}}\right)\left(1+\frac{\alpha}{\kappa_{0} + {\rm i}\varepsilon_{0}^{\prime}}\right)\prod_{m=1}^{+\infty}\left(1-\frac{\alpha^2}{\kappa_m^2}\right). \end{gather}$$

From (3.13), it can be shown that $\mathcal {k}_m=\textrm {i}(m+\tfrac {1}{2}){\rm \pi} +O(m^{-1})$ and $\kappa _m=\textrm {i}m{\rm \pi} +O(m^{-3})$ when $m\rightarrow +\infty$, which mean that (3.19) is convergent. We may define

(3.20)\begin{equation} \mathcal{K}_{i}^{\epsilon}(\alpha,F)=K_{i}^\epsilon(\alpha,F)/\alpha^2,\quad i = 1,2, \end{equation}

where $\mathcal {K}_i^{\epsilon }(\alpha,F)$ ($i=1, 2$) is bounded at $\alpha =0$. In (3.19), when $F\rightarrow 1$, $(1-F^2 )\rightarrow 0$, $\mathcal {k}_0\rightarrow 0$ and $\kappa _0\rightarrow 0$, and it can be found from (3.13) that ${\lim }_{F\rightarrow 1}({1-F^2})/{\mathcal {k}_0^2}={\lim }_{F\rightarrow 1}({1-F^2})/{\kappa _0^2}=\tfrac {1}{3}$. Therefore, the Weierstrass factorisations in (3.19) are still valid. Substituting (3.19) into (3.15), we have

(3.21) \begin{align} K^{\epsilon}(\alpha,F)&=\frac{\mathcal{k}_0\gamma_{{-}1}^+(\alpha,F)\gamma_{{-}1}^-(\alpha,F)\gamma_{0}^+(\alpha,F)\gamma_{0}^-(\alpha,F)\left[(\kappa_{{-}1}\!+\!{\rm i}\varepsilon_{{-}1}^{\prime})+\alpha\right]\left[(\kappa_{0}\!-\!{\rm i}\varepsilon_{0}^{\prime})-\alpha\right]}{\kappa_{{-}1}\kappa_{0}\beta_0^+(\alpha,F)\beta_{0}^-(\alpha,F)\left[(\mathcal{k}_0 \!-\! {\rm i}\epsilon_0^{\prime})\!-\!\alpha\right]}\nonumber\\ &\quad \times\prod_{m=1}^{+\infty}\frac{\mathcal{k}_m^2\left(\kappa_m^2-\alpha^2\right)}{\kappa_m^2\left(\mathcal{k}_m^2-\alpha^2\right)}, \end{align}

where

(3.22a)\begin{gather} \beta_{0}^+(\alpha,F)=\left\{ \begin{array}{@{}ll} 1, & \text{Im}\{\mathcal{k}_0\}=0,\\ \dfrac{\mathcal{k}_0+\alpha}{\mathcal{k}_0}, & \text{Im}\{\mathcal{k}_0\}>0, \end{array} \right.\quad \beta_{0}^-(\alpha,F)=\!\left\{ \begin{array}{@{}ll} \dfrac{(\mathcal{k}_0+{\rm i}\epsilon_0^{\prime})+\alpha}{\mathcal{k}_0}, & \text{Im}\{\mathcal{k}_0\}\!=\!0,\\ 1, & \text{Im}\{\mathcal{k}_0\}>0, \end{array} \right. \end{gather}

(3.22b) \begin{gather}\left. \begin{array}{l} \gamma_{{-}1}^+(\alpha,F)=\left\{ \begin{array}{@{}ll} \dfrac{(\kappa_{{-}1}-{\rm i}\varepsilon_{{-}1}^{\prime})-\alpha}{\kappa_{{-}1}}, & \text{Im}\{\kappa_{{-}1}\}=0,\\ 1, & \text{Im}\{\kappa_{{-}1}\}>0, \end{array} \right. \nonumber\\ \gamma_{{-}1}^-(\alpha,F) =\left\{ \begin{array}{@{}ll} 1, & \text{Im}\{\kappa_{{-}1}\}=0,\\ \dfrac{\kappa_{{-}1}-\alpha}{\kappa_{{-}1}}, & \text{Im}\{\kappa_{{-}1}\}>0. \end{array} \right. \end{array}\right\} \end{gather}

(3.22c)\begin{gather} \gamma_{0}^+(\alpha,F)=\left\{ \begin{array}{@{}ll} 1, & \text{Im}\{\kappa_0\}=0,\\ \dfrac{\kappa_{0}+\alpha}{\kappa_{0}}, & \text{Im}\{\kappa_{0}\}>0, \end{array} \right.\quad \gamma_{0}^-(\alpha,F)=\left\{ \begin{array}{@{}ll} \dfrac{(\kappa_{0}+{\rm i}\varepsilon_0^{\prime})+\alpha}{\kappa_0}, & \text{Im}\{\kappa_0\}=0,\\ 1, & \text{Im}\{\kappa_0\}>0. \end{array} \right. \end{gather}

From (3.12), (3.13) and (3.15), it can be shown that $K^{{\epsilon }}(\alpha,F)\sim O(|\alpha |^3)$ when $|\alpha |\rightarrow +\infty$. Using (3.16), $K_\pm ^{\epsilon }(\alpha,F)$ can be defined as

(3.23a)\begin{gather} K_+^{\epsilon}(\alpha,F)=\frac{\gamma_{{-}1}^+(\alpha,F)\gamma_{0}^+(\alpha,F)}{\beta_0^+(\alpha,F)}\frac{\left[(\kappa_{{-}1}+{\rm i}\varepsilon_{{-}1}^{\prime})+\alpha\right]}{\kappa_{{-}1}}\prod_{m=1}^{+\infty}\frac{\mathcal{k}_m(\kappa_m+\alpha)}{\kappa_m(\mathcal{k}_m+\alpha)}, \end{gather}

(3.23b)\begin{gather} K_-^{\epsilon}(\alpha,F)=\frac{\gamma_{{-}1}^-(\alpha,F)\gamma_{0}^-(\alpha,F)}{\beta_0^-(\alpha,F)}\frac{\mathcal{k}_0\left[(\kappa_{0}-{\rm i}\varepsilon_{0}^{\prime})+\alpha\right]}{\kappa_0\left[(\mathcal{k}_0-{\rm i}\epsilon_0^{\prime})-\alpha\right]}\prod_{m=1}^{+\infty}\frac{\mathcal{k}_m(\kappa_m-\alpha)}{\kappa_m(\mathcal{k}_m-\alpha)}, \end{gather}

which is to ensure that $K_+^{\epsilon }(\alpha,F)$ and $K_-^{\epsilon }(\alpha,F)$ are analytical in the upper and lower half complex planes, respectively. As the singularities may move across the real axis when $F$ changes. $\beta _{0}^{\pm }(\alpha,F)$ and $\gamma _{-1}^\pm (\alpha,F)$, $\gamma _{0}^\pm (\alpha,F)$ are introduced to ensure $K_\pm ^{\epsilon }(\alpha,F)$ remain analytical in their corresponding regions at different ranges of $F$. In particular, when $\textrm {Im}\{\mathcal {k}_0\}=0$ ($0< F<1$) and $\textrm {Im}\{\kappa _{0}\}=0$ ($F_c< F<1$), the terms $({(\mathcal {k}_0+\textrm {i}\epsilon _0^{\prime })+\alpha })/{\mathcal {k}_0}$ and $({(\kappa _0+\textrm {i}\varepsilon _0^{\prime })+\alpha })/{\kappa _0}$ will appear in $K_-^{\epsilon }(\alpha,F)$, otherwise it will appear in $K_+^{\epsilon }(\alpha,F)$. When $\textrm {Im}\{\kappa _{-1}\}=0$ ($F>F_c$), the term $({(\kappa _{-1}-\textrm {i}\varepsilon _{-1}^{\prime })-\alpha })/{\kappa _{-1}}$ will be in $K_+^{\epsilon }(\alpha,F)$, otherwise it will appear in $K_-^{\epsilon }(\alpha,F)$.

In (3.19), we may define $\tau = \min \{ |\textrm {Im}\{\mathcal {k}_0 - \textrm {i}\epsilon _0^{\prime }\}|, |\textrm {Im}\{\kappa _{-1}-\textrm {i}\varepsilon _{-1}^{\prime }\}|, |\textrm {Im}\{\kappa _{0} - \textrm {i}\varepsilon _{0}^{\prime }\}|,$ $|\textrm {Im}\{\kappa _{1}\}|,|\textrm {Im}\{\mathcal {k}_1\}|\}$, which makes $|\textrm {Im}\{\mathcal {k}_m\}|\geq \tau$ and $|\textrm {Im}\{\kappa _m\}|\geq \tau$ for all $m$. Then, $K_+^{\epsilon }(\alpha,F)$ and $K_-^{\epsilon }(\alpha,F)$ in (3.23) are analytical in the complex planes $S_+$ with $\textrm {Im}\{\alpha \}>-\tau$ and $S_-$ with $\textrm {Im}\{\alpha \}<\tau$, respectively, as well as have no zero in $S_+$ and $S_-$, respectively. As shown in (B7) of Appendix B, $K_+^{\epsilon }(\alpha,F) \sim O(|\alpha |^{5/2})$ and $K_-^{\epsilon }(\alpha,F) \sim O(|\alpha |^{1/2})$ as $|\alpha |\rightarrow +\infty$, Thus, the asymptotic behaviour of $K_+^{\epsilon }(\alpha,F)K_-^{\epsilon }(\alpha,F)$ is consistent with that of $K^{\epsilon }(\alpha,F)$. Substituting (3.16) and (3.20) into (3.14), and dividing $K_-^{\epsilon }(\alpha,F)$ on both sides, we have

(3.24)\begin{equation} \frac{F_{-}(\alpha)}{K_{-}^{\epsilon}(\alpha,F)}=D_+(\alpha)K_{+}^{\epsilon}(\alpha,F)-\frac{2{\rm \pi} D F^2\alpha^4C(z_0,\alpha){\rm e}^{{\rm i}\alpha x_0}}{\mathcal{K}_1^{\epsilon}(\alpha,F)K_{-}^\epsilon(\alpha,F)}. \end{equation}

The second term on the right-hand side needs to be further decomposed. Noting (3.15) and (3.16), we may write

(3.25)\begin{align} \frac{DF^2\alpha^2C(z_0,\alpha){\rm e}^{{\rm i}\alpha x_0}}{\mathcal{K}_1^{\epsilon}(\alpha,F)K_{-}^{\epsilon}(\alpha,F)}=\frac{DF^2\alpha^2K_+^{\epsilon}(\alpha,F)C(z_0,\alpha){\rm e}^{{\rm i}\alpha x_0}}{\mathcal{K}_2^{\epsilon}(\alpha,F)}=M_{+}(\alpha,x_0,z_0)+M_{-}(\alpha,x_0,z_0), \end{align}

$M_{\pm }(\alpha,x_0,z_0 )$ are analytical in $S_\pm$, respectively. Following the procedure in Noble (1958), when $x_0>0$, $M_-(\alpha,x_0,z_0)$ can be obtained by the Cauchy integral in the upper half-plane as

(3.26)\begin{align} M_{-}(\alpha,x_0,z_0)&={-}\frac{DF^2}{2{\rm \pi}{\rm i}}\int_{-\infty -{\rm i}\sigma}^{+\infty -{\rm i}\sigma}\frac{\zeta^2K_{+}^{\epsilon}(\zeta,F)C(z_0,\zeta){\rm e}^{{\rm i}\zeta x_0}}{\mathcal{K}_2^{\epsilon}(\zeta,F)(\zeta - \alpha)}\,{\rm d}\zeta\nonumber\\ &= DF^2\sum_{\zeta\in R_+}\frac{\zeta^2K_{+}^{\epsilon}(\zeta,F)C(z_0,\zeta){\rm e}^{{\rm i}\zeta x_0}}{\mathcal{K}_2^{\epsilon\prime}(\zeta,F)(\zeta - \alpha)}, \end{align}

where $0<\sigma <\tau$, $R_+$ is the set containing all the roots of $K_2^{\epsilon }(\alpha,F)=0$ on the complex plane $S_+$, which can be expressed as

(3.27)\begin{equation} R_+{=} \left\{ \begin{array}{@{}ll} \left\{\kappa_m|m={-}1,0,1,\ldots\right\}, & 0< F< F_c,\\ \left\{\kappa_m|m=0,1,2,\ldots\right\}\cup\left\{{\pm}\kappa_{0}-{\rm i}\varepsilon_0^{\prime}\right\}, & F_c< F<1,\\ \left\{\kappa_m|m=0,1,2,\ldots\right\}, & F>1. \end{array} \right. \end{equation}

Subsequently, $M_+(\alpha,x_0,z_0)$ can be obtained from

(3.28)\begin{equation} M_{+}(\alpha,x_0,z_0)=\frac{DF^2\alpha^2C(z_0,\alpha){\rm e}^{{\rm i}\alpha x_0}}{\mathcal{K}_1^{\epsilon}(\alpha,F)K_{-}^{\epsilon}(\alpha, F)} - M_{-}(\alpha,x_0,z_0). \end{equation}

When $x_0<0$, the Cauchy integral can be applied in the lower half-plane to obtain $M_{+}(\alpha,x_0,z_0)$ or $S_-$. Substituting (3.25) into (3.24), we have

(3.29)\begin{equation} \frac{F_{-}(\alpha)}{K_{-}^{\epsilon}(\alpha,F)}+2{\rm \pi}\alpha^2M_{-}(\alpha,x_0,z_0)=D_{+}(\alpha)K_{+}^{\epsilon}(\alpha,F)-2{\rm \pi} \alpha^2M_{+}(\alpha,x_0,z_0). \end{equation}

The functions on the left- and right-hand sides of (3.29) are analytical in the complex domains $\textrm {Im}\{\alpha \}<\tau$ and $\textrm {Im}\{\alpha \}>-\tau$, respectively. As these domains overlap within $-\tau <\textrm {Im}\{\alpha \}<\tau$, based on the theorem of analytical continuation, these two functions should be identical and analytical in the entire complex plane. From the Liouville's theorem, such a function should be a polynomial $2{\rm \pi} Q(\alpha )$, which provides

(3.30a)\begin{gather} \frac{F_{-}(\alpha)}{K_{-}^{\epsilon}(\alpha,F)}+2{\rm \pi}\alpha^2M_{-}(\alpha,x_0,z_0)=2{\rm \pi} Q(\alpha), \end{gather}

(3.30b)\begin{gather} D_{+}(\alpha)K_{+}^{\epsilon}(\alpha,F)-2{\rm \pi} \alpha^2M_{+}(\alpha,x_0,z_0) = 2{\rm \pi} Q(\alpha). \end{gather}

Using (3.11) and (3.30a), $A(\alpha )$ can be expressed as

(3.31)\begin{align} A(\alpha)&=\frac{2{\rm \pi} Q(\alpha)K_{-}^{\epsilon}(\alpha,F)}{K_2^{\epsilon}(\alpha F)}-\frac{2{\rm \pi} \alpha^2M_{-}(\alpha,x_0,z_0)K_{-}^{\epsilon}(\alpha,F)}{K_2^{\epsilon}(\alpha,F)}\nonumber\\ &\quad -\frac{2{\rm \pi} {\rm e}^{{\rm i}\alpha x_0}}{\alpha}\frac{\left[(D\alpha^4 + 1)\alpha \cosh\alpha - F^2\alpha^2 \sinh\alpha\right]C(z_0,\alpha)}{K_2^{\epsilon}(\alpha,F)}. \end{align}

Substituting (3.31) into (3.8), the Green function can be obtained through the inverse Fourier transform

(3.32)\begin{equation} G=\frac{1}{2{\rm \pi}}\int_{-\infty}^{+\infty}\hat{G}{\rm e}^{-{\rm i}\alpha x}\,{\rm d}\alpha. \end{equation}

Letting $\epsilon \rightarrow 0^-$, it provides

(3.33)\begin{equation} G = G_1 + G_2 + G_{ice}, \end{equation}

where

(3.34a)\begin{gather} G_1 = \int_{-\infty}^{+\infty}\frac{Q(\alpha)K_{-}(\alpha,F)C(z,\alpha){\rm e}^{-{\rm i}\alpha x}}{\alpha^2\mathcal{K}_2(\alpha,F)}\,{\rm d}\alpha, \end{gather}

(3.34b)\begin{gather} G_2={-}\int_{-\infty}^{+\infty}\frac{M_{-}(\alpha,x_0,z_0)K_{-}(\alpha,F)C(z,\alpha) {\rm e}^{-{\rm i}\alpha x}}{\mathcal{K}_2(\alpha,F)}\,{\rm d}\alpha, \end{gather}

(3.34c)\begin{gather} G_{ice} = \ln r + \ln r^{\prime}-\int_{-\infty}^{+\infty}\frac{{\rm e}^{-|\alpha|}\left(D\alpha^4 + 1 + F^2|\alpha|\right)C(z,\alpha)C(z_0,\alpha) {\rm e}^{-{\rm i}\alpha(x-x_0)}}{\alpha^2\mathcal{K}_2(\alpha,F)}\,{\rm d}\alpha, \end{gather}

with $r=\sqrt {(x-x_0)^2 + (z-z_0)^2}$ and $r^{\prime }=\sqrt {(x-x_0)^2 + (z+z_0+2)^2}$, $\mathcal {K}_i(\alpha,F)=K_i(\alpha,\alpha F)/\alpha ^2$ ($i=1, 2$). In (3.34c) $G_{ice}$ corresponds to the third term in (3.31) and represents the Green function of the fluid fully covered by an ice sheet, which is obtained by applying a procedure similar to that in Li et al. (2019). It should be noted that the singularity at $\alpha =0$ in the integrands of $G$ has been treated through the Hadamard regulation and Cauchy principal value. This is equivalent to adding a constant term to $G$. It will not affect the physics, which involves only the spatial derivatives of $G$. When $\epsilon \rightarrow 0^-$, $K_2^{\epsilon }(\alpha,F)\rightarrow K_2(\alpha,\alpha F)$ and we may consider how some roots approach the real axis of $\alpha$. For $K_-(\alpha,F)$, when $F<1$, $\pm \mathcal {k}_0-\textrm {i}\epsilon _0^{\prime }$ approach the real axis of $\alpha$ from above. For $\mathcal {K}_2(\alpha,F)$, when $F>F_c$, $\pm \kappa _{-1}-\textrm {i}\varepsilon _{-1}^{\prime }$ approach from below, and when $F_c< F<1$, $\pm \kappa _{0}-\textrm {i}\varepsilon _{0}^{\prime }$ approach from above. In such a case, the integration path in $G$ from $-\infty$ to $+\infty$ should pass under the poles at $\pm \mathcal {k}_0$ when $0< F<1$, and should pass under the poles at $\pm \kappa _0$ and pass over the poles at $\pm \kappa _{-1}$ when $F_c< F<1$, and should pass over the poles at $\pm \kappa _{-1}$ only when $F>1$.

The Green function can be also obtained by using (3.10) and (3.30b), which gives

(3.35)\begin{equation} G = G_3 + G_4 + G_{water}, \end{equation}

where

(3.36a)\begin{gather} G_3 = \int_{-\infty}^{+\infty}\frac{Q(\alpha)C(z,\alpha){\rm e}^{-{\rm i}\alpha x}}{\alpha^2K_{+}(\alpha,F)\mathcal{K}_1(\alpha,F)}\,{\rm d}\alpha, \end{gather}

(3.36b)\begin{gather} G_4 =\int_{-\infty}^{+\infty}\frac{M_{+}(\alpha,x_0,z_0)C(z,\alpha){\rm e}^{-{\rm i}\alpha x}}{K_{+}(\alpha,F)\mathcal{K}_1(\alpha,F)}\,{\rm d}\alpha, \end{gather}

(3.36c)\begin{gather} G_{water} = \ln r + \ln r^{\prime}-\int_{-\infty}^{+\infty}\frac{{\rm e}^{-|\alpha|}\left(1 + F^2|\alpha|\right)C(z,\alpha)C(z_0,\alpha) {\rm e}^{-{\rm i}\alpha(x-x_0)}}{\alpha^2\mathcal{K}_1(\alpha,F)}\,{\rm d}\alpha, \end{gather}

where $G_{water}$ denotes the Green function for free surface flow, and it should be noted that $G_1=G_3$. The integration path is the same as that discussed previously for (3.34).

The ice sheet deflection $\xi$ at $x>0$ can be obtained by applying (3.33) and (3.34) to (2.10), which provides

(3.37)\begin{equation} \xi(x,x_0,z_0)=\xi_1(x,x_0,z_0) + \xi_2(x,x_0,z_0) +\xi_{ice}(x,x_0,z_0),\quad x>0, \end{equation}

where

(3.38a)\begin{gather} \xi_1 = \frac{{\rm i}}{F}\int_{-\infty}^{+\infty}\frac{Q(\alpha)K_{-}(\alpha,F)\left[D\alpha^3S(z,\alpha) - F^2 C(z,\alpha)\right]}{\alpha\mathcal{K}_2(\alpha,F)} {\rm e}^{-{\rm i}\alpha x}\,{\rm d}\alpha,\quad z\rightarrow 0^-, \end{gather}

(3.38b)\begin{gather} \xi_2 =\frac{{\rm i}}{F} \int_{-\infty}^{+\infty}\frac{M_{-}(\alpha,x_0,z_0)K_{-}(\alpha,F)\sinh\alpha}{\mathcal{K}_2(\alpha,F)} {\rm e}^{-{\rm i}\alpha x}\,{\rm d}\alpha, \end{gather}

(3.38c)\begin{gather} \xi_{ice} ={\rm i}F\int_{-\infty}^{+\infty}\frac{C(z_0,\alpha)}{\alpha \mathcal{K}_2(\alpha,F)} {\rm e}^{-{\rm i}\alpha(x-x_0)}\,{\rm d}\alpha, \end{gather}

$z\rightarrow 0^-$ is used in $\xi _1$ is to ensure the convergence of the integral. Similarly, substituting (3.35) and (3.36) into (2.10), the free surface wave elevation at $x<0$ can be written as

(3.39)\begin{equation} \xi(x,x_0,z_0)=\xi_3(x,x_0,z_0) + \xi_4(x,x_0,z_0) +\xi_{water}(x,x_0,z_0),\quad x<0, \end{equation}

where

(3.40a)\begin{gather} \xi_3 ={-}{\rm i}F\int_{-\infty}^{+\infty}\frac{Q(\alpha)C(z,\alpha)}{\alpha K_{+}(\alpha,F)\mathcal{K}_1(\alpha,F)} {\rm e}^{-{\rm i}\alpha x}\,{\rm d}\alpha,\quad z \rightarrow 0^-, \end{gather}

(3.40b)\begin{gather} \xi_4 ={-}\frac{{\rm i}}{F} \int_{-\infty}^{+\infty}\frac{M_{+}(\alpha,x_0,z_0)\sinh\alpha}{K_{+}(\alpha,F)\mathcal{K}_1(\alpha,F)} {\rm e}^{-{\rm i}\alpha x}\,{\rm d}\alpha, \end{gather}

(3.40c)\begin{gather} \xi_{water} ={\rm i}F\int_{-\infty}^{+\infty}\frac{C(z_0,\alpha)}{\alpha \mathcal{K}_1(\alpha,F)} {\rm e}^{-{\rm i}\alpha(x-x_0)}\,{\rm d}\alpha. \end{gather}

Using (3.13b) in (3.38a), and (3.13a) in (3.40a), and removing the zero integral terms with no singularity, $\xi _1$ and $\xi _3$ can be written in the following form

(3.41a)\begin{gather} \xi_1 ={-}\frac{{\rm i}}{F}\int_{-\infty}^{+\infty}\frac{Q(\alpha)K_{-}(\alpha,F)\sinh\alpha}{\alpha^2\mathcal{K}_2(\alpha,F)} {\rm e}^{-{\rm i}\alpha x}\,{\rm d}\alpha+\frac{\rm \pi}{F}Q(0),\quad x>0, \end{gather}

(3.41b)\begin{gather} \xi_3 ={-}\frac{{\rm i}}{F}\int_{-\infty}^{+\infty}\frac{Q(\alpha)\sinh\alpha}{\alpha^2K_{+}(\alpha,F)\mathcal{K}_1(\alpha,F)}{\rm e}^{-{\rm i}\alpha x}\,{\rm d}\alpha-\frac{\rm \pi}{F}Q(0),\quad x<0, \end{gather}

where the term $\pm ({{\rm \pi} }/{F})Q(0)$ results from the difference of the residues at $\alpha =0$ of the two integrands in (3.38a) and (3.40a). In (3.41a), $K_-(\alpha,F)\sim O(|\alpha |^{1/2})$ and ${\mathcal {K}_2(\alpha,F)}/{\sinh \alpha }\sim O(|\alpha |^3)$ as $|\alpha |\rightarrow +\infty$. If there are terms of $\alpha ^n$ ($n\geq 4$) in $Q(\alpha )$, the integral will be divergent. Thus, based on the discussion above, $Q(\alpha )$ can at most be a cubic polynomial as

(3.42)\begin{equation} Q(\alpha) = b + c\alpha + {\rm d}\alpha^2 + f\alpha^3, \end{equation}

where $b$, $c$, $d$ and $f$ are four unknown coefficients. However, there are only two edge conditions, which means two of them are undetermined. Here, common in this kind of problem, when the flow leaves the edge of the plate, we impose the Kutta condition, which is achieved by assuming the free surface and ice sheet have the same elevation and same slope at $x=0$, or

(3.43)\begin{equation} \left. \begin{aligned} \xi\left(0^+,x_0,z_0\right) & = \xi\left(0^-,x_0,z_0\right),\\ \frac{\partial \xi}{\partial x}\left(0^+,x_0,z_0\right) & = \frac{\partial \xi}{\partial x}\left(0^-,x_0,z_0\right). \end{aligned} \right\} \end{equation}

Substituting (3.37)–(3.41) into (3.43) and using (3.42), we have $\xi (0^+,x_0,z_0)-\xi (0^-,x_0,z_0)=2{\rm \pi} b/F=0$, which gives $b=0$. We also have $({\partial \xi }/{\partial x})(0^+,x_0,z_0)-({\partial \xi }/{\partial x})(0^-,x_0,z_0)=f[I^{\prime \prime \prime }(0^-)-I^{\prime \prime \prime }(0^+)]$, where

(3.44)\begin{equation} I(x)={-}\frac{{\rm i}}{F}\int_{-\infty}^{+\infty}\frac{K_{-}(\alpha,F)\sinh \alpha}{\alpha \mathcal{K}_2(\alpha,F)}{\rm e}^{-{\rm i}\alpha x}\,{\rm d}\alpha={-}\frac{{\rm i}}{F}\int_{-\infty}^{+\infty}\frac{\sinh\alpha}{\alpha\mathcal{K}_1(\alpha,F)K_+(\alpha,F)} {\rm e}^{-{\rm i}\alpha x}\,{\rm d}\alpha. \end{equation}

Invoking $\sinh \alpha =({\mathcal {K}_2(\alpha,F)-\mathcal {K}_1(\alpha,F)})/{D\alpha ^3}$, (3.44) can be expressed as

(3.45)\begin{equation} I(x) ={-}\frac{{\rm i}}{DF}\int_{-\infty}^{+\infty}\frac{1}{\alpha^4}\left[K_{-}(\alpha,F)-\frac{1}{K_{+}(\alpha,F)}\right] {\rm e}^{-{\rm i}\alpha x}\,{\rm d}\alpha. \end{equation}

Noting that $K_\pm (\alpha,F)$ are analytical in $S_\pm$, respectively, we may consider the integral in $S_-$ when $x>0$ and in $S_+$ when $x<0$. This gives

(3.46)\begin{equation} I(x) = \left\{ \begin{array}{ll} \dfrac{{\rm i}}{DF}\displaystyle\int_{-\infty}^{+\infty}\dfrac{{\rm e}^{-{\rm i}\alpha x}}{\alpha^4K_+(\alpha,F)}\,{\rm d}\alpha, & x>0,\\ \dfrac{{\rm i}}{DF}\displaystyle\int_{-\infty}^{+\infty}\dfrac{K_{-}(\alpha,F){\rm e}^{-{\rm i}\alpha x}}{\alpha^4}\,{\rm d}\alpha, & x<0, \end{array} \right. \end{equation}

where the integral paths in (3.46) for $x>0$ and $x<0$ should pass under and over the pole at $\alpha =0$, respectively. It can be shown from (3.46) that $I^{\prime \prime \prime }(0^-)\rightarrow \infty$ while $I^{\prime \prime \prime }(0^+)$ is finite, which provides $f=0$. By further substituting (3.37), (3.38b), (3.38c) and (3.41a) into the edge conditions in (2.9), and using (3.44), we have

(3.47)\begin{equation} \left[ \begin{array}{cc} I^{\prime\prime}(0^+) & {\rm i}I^{\prime\prime\prime}(0^+)\\ I^{\prime\prime\prime}(0^+) & {\rm i}I^{(4)}(0^+) \end{array} \right] \left[ \begin{array}{c} c\\ d \end{array} \right]= \left[ \begin{array}{c} \dfrac{\partial^2}{\partial x^2}\left[\xi_2(0^+,x_0,z_0)+\xi_{ice}(0^+,x_0,z_0)\right]\\ \dfrac{\partial^3}{\partial x^3}\left[\xi_2(0^+,x_0,z_0)+\xi_{ice}(0^+,x_0,z_0)\right] \end{array} \right]. \end{equation}

Differentiating (3.46) with respect to $x$, letting $x\rightarrow 0^+$ and then applying the theorem of residue in the upper half-plane of $\alpha$, we obtain

(3.48)\begin{equation} \left. \begin{aligned} I^{\prime\prime}(0^+) & =\frac{2{\rm \pi}}{DF}K_+^{\prime}(0,F),\\ I^{\prime\prime\prime}(0^+) & ={-}\frac{2{\rm \pi} {\rm i}}{DF},\\ I^{(4)}(0^+) & = 0. \end{aligned} \right\} \end{equation}

The right-hand side of (3.47) can be also treated in a similar way as in (3.44)–(3.46) and (3.48). This gives,

(3.49)\begin{equation} \left. \begin{aligned} \frac{\partial^2}{\partial x^2}\left[\xi_2(0^+,x_0,z_0)+\xi_{ice}(0^+,x_0,z_0)\right] & ={-}\frac{2{\rm \pi}}{DF}M_{-}\left(0,x_0,z_0\right),\\ \frac{\partial^3}{\partial x^3}\left[\xi_2(0^+,x_0,z_0)+\xi_{ice}(0^+,x_0,z_0)\right] & = 0. \end{aligned} \right\} \end{equation}

Substituting (3.48) and (3.49) into (3.47), we have

(3.50)\begin{equation} \left. \begin{aligned} c & = 0,\\ d & = M_{-}\left(0,x_0,z_0\right). \end{aligned} \right\} \end{equation}

Using this with (3.34) and (3.36), and noticing $M_+(0,x_0,z_0)=-M_-(0,x_0,z_0)$, the Green function may be written as

(3.51a)\begin{gather} G = G_{ice} - \int_{-\infty}^{+\infty}\frac{N_{-}\left(\alpha,x_0,z_0\right)K_{-}(\alpha,F)C(z,\alpha) {\rm e}^{-{\rm i}\alpha x}}{\mathcal{K}_2(\alpha,F)}\,{\rm d}\alpha, \end{gather}

(3.51b)\begin{gather} G = G_{water} + \int_{-\infty}^{+\infty}\frac{N_{+}\left(\alpha,x_0,z_0\right)C(z,\alpha) {\rm e}^{-{\rm i}\alpha x}}{K_{+}(\alpha,F)\mathcal{K}_1(\alpha,F)}\,{\rm d}\alpha, \end{gather}

where $N_{\pm }(\alpha,x_0,z_0 )=M_{\pm }(\alpha,x_0,z_0 )-M_\pm (0,x_0,z_0 )$.

3.2. Multipole expansion for a submerged circular cylinder

Once the Green function has been determined, the potentials due to multipoles or higher-order singularities can be found by differentiating the Green function in (3.33) and (3.34) with respect to the position of the source ($x_0,z_0$). Define $x-x_0=r\sin \theta$ and $z-z_0=r\cos \theta$, and apply the following operator (Wu 1998)

(3.52)\begin{equation} \left(D_{{\pm}}\right)^n={-}\frac{1}{2^{n-1}(n-1)!}\left(\frac{\partial }{\partial z_0}\pm{\rm i}\frac{\partial}{\partial x_0}\right)^n. \end{equation}

Note that the Green function here is a real function since the problem is steady. In fact, we may use $K_{-}(-\alpha,F)=\overline {K_{-}}(\alpha,F)$ and $N_{-}(-\alpha,x_0,z_0)=\overline {N_{-}}(\alpha,x_0,z_0)$, (3.51a) and (3.51b) become

(3.53a)\begin{gather} G = G_{ice}-2\text{Re}\left\{\int_{0}^{+\infty}\frac{N_{-}(\alpha,x_0,z_0)K_{-}(\alpha,F)C(z,\alpha) {\rm e}^{-{\rm i}\alpha x}}{\mathcal{K}_2(\alpha,F)}\,{\rm d}\alpha\right\}, \end{gather}

(3.53b)\begin{gather} G = G_{water}+2\text{Re}\left\{\int_{0}^{+\infty}\frac{N_{+}(\alpha,x_0,z_0)C(z,\alpha) {\rm e}^{-{\rm i}\alpha x}}{K_{+}(\alpha,F)\mathcal{K}_1(\alpha,F)}\,{\rm d}\alpha\right\}, \end{gather}

where

(3.54a)\begin{gather} G_{ice} = \ln r + \ln r^{\prime}-2\text{Re}\left\{\int_{0}^{+\infty}\frac{{\rm e}^{-\alpha}\left(D\alpha^4 + 1 + F^2\alpha\right)C(z,\alpha)C(z_0,\alpha){\rm e}^{-{\rm i}\alpha(x-x_0)}}{\alpha^2\mathcal{K}_2(\alpha,F)}\,{\rm d}\alpha\right\}, \end{gather}

(3.54b)\begin{gather} G_{water} = \ln r + \ln r^{\prime}-2\text{Re}\left\{\int_{0}^{+\infty}\frac{{\rm e}^{-\alpha}\left(1 + F^2\alpha\right)C(z,\alpha)C(z_0,\alpha){\rm e}^{-{\rm i}\alpha(x-x_0)}}{\alpha^2\mathcal{K}_1(\alpha,F)}\,{\rm d}\alpha\right\}. \end{gather}

In such a case, $(D_+)^n$ and $(D_-)^n$ lead to a pair of conjugate functions. Thus, we may apply only $(D_+)^n$ here. Using (3.26), (3.51a) and

(3.55) \begin{equation} \left. \begin{gathered} \left(D_{+}\right)^n \ln r = \frac{{\rm e}^{{\rm i}n\theta}}{r^n},\\ \left(D_{+}\right)^n \big[C(z_0,\alpha){\rm e}^{{\pm}{\rm i}\alpha x_0}\big] ={-}\frac{({\mp} 1)^n}{(n-1)!}\alpha^n \exp({\pm {\rm i}\alpha x_0 \mp \alpha(z_0+1)}), \end{gathered} \right\} \end{equation}

we have

(3.56)\begin{equation} G_n = G_{n,1} + G_{n,ice}, \end{equation}

where

(3.57a)\begin{equation} G_{n,1} = \frac{1}{(n-1)!}\int_{-\infty}^{+\infty}\frac{\hat{N}_n(\alpha,x_0,z_0)K_{-}(\alpha,F)C(z,\alpha) {\rm e}^{-{\rm i}\alpha x}}{\mathcal{K}_2(\alpha,F)}\,{\rm d}\alpha, \end{equation}

(3.57b) \begin{align} G_{n,ice}&=\frac{{\rm e}^{{\rm i}n\theta}}{r^n}+\frac{({-}1)^n}{(n-1)!}\int_{0}^{+\infty}\alpha^{n-1} \exp({-\alpha(z+z_0+2)-{\rm i}\alpha(x-x_0)})\,{\rm d}\alpha\nonumber\\ &\quad+\frac{({-}1)^n}{(n-1)!}\int_{-\infty}^{+\infty}\nonumber\\ &\quad\times\frac{\alpha^{n-2} {\rm e}^{-|\alpha|}\left[(D\alpha^4\!+\!1)\!+\!F^2|\alpha|\right]C(z,\alpha) \exp({-\alpha(z_0\!+\!1)\!-{\rm i}\alpha(x\!-x_0)})}{\mathcal{K}_2(\alpha,F)}\,{\rm d}\alpha, \end{align}

(3.58)\begin{align} \hat{N}_n(\alpha,x_0,z_0)&={-}(n-1)!\left(D_+\right)^nN_{-}(\alpha,x_0,z_0)\nonumber\\ &=DF^2({-}1)^n\alpha\sum_{\zeta\in R_+}\frac{\zeta^{n+1}K_+(\zeta,F) \exp({{\rm i}\zeta x_0 -\zeta (z_0+1)})}{\mathcal{K}_2^{\prime}(\zeta,F)(\alpha - \zeta)}. \end{align}

Then, the velocity potential due to a submerged cylinder can be expressed in a multipole expansion form as

(3.59)\begin{equation} \phi={\rm Re}\left\{\sum_{n=1}^{+\infty}a^nA_nG_n\right\}, \end{equation}

where $A_n$ are unknown coefficient. Substituting (3.56) and (3.57) into (3.59) and using (Abramowitz & Stegun 1968)

(3.60)\begin{equation} \left. \begin{aligned} \exp({\alpha(z+1)-{\rm i}\alpha x}) & = \exp({\alpha(z_0+1)-{\rm i}\alpha x_0})\sum_{l=0}^{+\infty}\frac{\alpha^l r^l}{l!} {\rm e}^{-{\rm i}l\theta},\\ \exp({-\alpha(z+1)-{\rm i}\alpha x}) & = \exp({-\alpha(z_0+1)-{\rm i}\alpha x_0})\sum_{l=0}^{+\infty}\frac{({-}1)^l\alpha^l r^l}{l!} {\rm e}^{-{\rm i}l\theta}, \end{aligned} \right\} \end{equation}

the velocity potential can be expressed in the polar coordinate $(r,\theta )$ as

(3.61) \begin{equation} \phi = {\rm Re}\left\{\sum_{n=1}^{+\infty}\left(\frac{a}{r}\right)^n+\sum_{n=1}^{+\infty}\sum_{l=0}^{+\infty}\frac{a^nA_n}{(n-1)!}\frac{r^l}{l!}\big[J_+(n,l) {\rm e}^{{\rm i}l\theta} + J_{-}(n,l){\rm e}^{-{\rm i}l\theta}\big]\right\}, \end{equation}

where

(3.62a)\begin{align} J_{+}(n,l)&=\frac{({-}1)^l}{2}\int_{-\infty}^{+\infty}\frac{\alpha^l\hat{N}_n(\alpha,x_0,z_0)K_{-}(\alpha,F) \exp({-\alpha(z_0+1)-{\rm i}\alpha x_0})}{\mathcal{K}_2(\alpha,F)}\,{\rm d}\alpha\nonumber\\ &\quad+\frac{({-}1)^{n+l}(n+l-1)!}{\left[2(z_0+1)\right]^{n+l}}\nonumber\\ &\quad +\frac{({-}1)^{n+l}}{2}\int_{-\infty}^{+\infty}\frac{\alpha^{l+n-2}{\rm e}^{-|\alpha|}\left(D\alpha^4 + 1 +F^2|\alpha|\right) {\rm e}^{{-}2\alpha(z_0+1)}}{\mathcal{K}_2(\alpha,F)}\,{\rm d}\alpha, \end{align}

(3.62b)\begin{align} J_{-}(n,l)&=\frac{1}{2}\int_{-\infty}^{+\infty}\frac{\alpha^l\hat{N}_n(\alpha,x_0,z_0)K_{-}(\alpha,F) \exp({\alpha(z_0+1)-{\rm i}\alpha x_0})}{\mathcal{K}_2(\alpha,F)}\,{\rm d}\alpha\nonumber\\ &\quad +\frac{({-}1)^{n}}{2}\int_{-\infty}^{+\infty}\frac{\alpha^{l+n-2} {\rm e}^{-|\alpha|}\left(D\alpha^4 + 1 +F^2|\alpha|\right)}{\mathcal{K}_2(\alpha,F)}\,{\rm d}\alpha. \end{align}

Substituting (3.61) into (2.7) and noting $n_x=-\sin \theta =-(({\textrm {e}^{\textrm {i}\theta } - \textrm {e}^{-\textrm {i}\theta }})/{2\textrm {i}})$, we obtain the following system of linear equations to solve $A_n$,

(3.63)\begin{align} -l\frac{A_l}{a}+\sum_{n=1}^{+\infty}\frac{a^{l+n-1}J_{+}(n,l)}{(n-1)!(l-1)!}A_n+\sum_{n=1}^{+\infty}\frac{a^{l+n-1}\bar{J}_{-}(n,l)}{(n-1)!(l-1)!}\bar{A}_n={-}{\rm i}F\delta_{l1},\quad l = 1,2,\ldots, \end{align}

where $\delta _{l1}$ denotes the Kronecker delta function. Equation (3.63) can be solved by applying the conjugate to obtain another set of equations, or by separating the real and imaginary parts, respectively.

3.3. Hydrodynamic forces on the submerged circular cylinder

The hydrodynamic forces can be determined by the integration of the hydrodynamic pressure over the surface of the cylinder, which gives

(3.64)\begin{equation} -{\rm i}F_R + F_L ={-}a\int_{-{\rm \pi}}^{\rm \pi}\left(p{\rm e}^{{\rm i}\theta}\right)_{r=a}\,{\rm d}\theta, \end{equation}

where $F_R$ and $F_L$ represent the drag and lift forces, respectively, $p$ denotes the non-dimensionalised pressure difference between the hydrodynamic pressure and atmospheric pressure. As explained in Wu (1991), the local disturbance of the flow field near the cylinder may not be small and the nonlinear terms in the Bernoulli equation should be retained, which gives

(3.65)\begin{equation} p={-}\tfrac{1}{2}\big[\boldsymbol{\nabla}\left(\phi -Fx\right)\boldsymbol{\cdot}\boldsymbol{\nabla}\left(\phi -Fx\right) - F^2\big]. \end{equation}

Through using a similar procedure in Li et al. (2019), (3.64) gives

(3.66)\begin{equation} -{\rm i}F_R + F_L = \frac{2{\rm \pi}}{a}\sum_{n=1}^{+\infty}n(n+1)A_n\bar{A}_{n+1}. \end{equation}

3.4. Wave elevation and ice sheet deflection

The ice sheet deflection can be obtained by substituting (3.56), (3.57) and (3.59) into (2.10). The ice sheet deflection $\eta$ at $x>0$ can be expressed as

(3.67)\begin{equation} \eta ={-}{\rm Im}\left\{\sum_{n=1}^{+\infty}\frac{a^nA_n}{(n-1)!}\left[\eta_1(n)+\eta_{ice}(n)\right]\right\},\quad x>0, \end{equation}

where

(3.68)\begin{equation} \left. \begin{aligned} \eta_1(n) & ={-}\frac{1}{F}\int_{-\infty}^{+\infty}\frac{\hat{N}_n(\alpha, x_0, z_0)K_{-}(\alpha, F)\sinh\alpha}{\mathcal{K}_2(\alpha,F)} {\rm e}^{-{\rm i}\alpha x}\,{\rm d}\alpha,\\ \eta_{ice}(n) & ={-}({-}1)^nF\int_{-\infty}^{+\infty}\frac{\alpha^{n-1}}{\mathcal{K}_2(\alpha,F)} \exp({-\alpha(z_0+1)-{\rm i}\alpha (x-x_0)})\,{\rm d}\alpha. \end{aligned} \right\} \end{equation}

Using a similar procedure, the free surface wave elevation at $x<0$ can be obtained as

(3.69)\begin{equation} \eta ={-}{\rm Im}\left\{\sum_{n=1}^{+\infty}\frac{a^nA_n}{(n-1)!}\left[\eta_2(n)+\eta_{water}(n)\right]\right\},\quad x<0, \end{equation}

where

(3.70)\begin{equation} \left. \begin{aligned} \eta_2(n) & = \frac{1}{F}\int_{-\infty}^{+\infty}\frac{\check{N}_n(\alpha, x_0, z_0)\sinh\alpha}{K_{+}(\alpha, F)\mathcal{K}_1(\alpha,F)} {\rm e}^{-{\rm i}\alpha x}\,{\rm d}\alpha,\\ \eta_{water}(n) & ={-}({-}1)^nF\int_{-\infty}^{+\infty}\frac{\alpha^{n-1}}{\mathcal{K}_1(\alpha,F)} \exp({-\alpha(z_0+1)-{\rm i}\alpha (x-x_0)})\,{\rm d}\alpha, \end{aligned} \right\} \end{equation}

and

(3.71)\begin{align} \check{N}_n(\alpha,x_0,z_0)&={-}(n-1)!\left(D_+\right)^nN_{+}(\alpha,x_0,z_0)\nonumber\\ &\quad =\frac{DF^2({-}1)^n\alpha^{n+2}\exp({{\rm i}\alpha x_0 - \alpha(z_0+1)})}{\mathcal{K}_1(\alpha,F)K_{-}(\alpha,F)}-\hat{N}(\alpha,x_0,z_0). \end{align}

4. Numerical results and discussion

The dimensional physical parameters used in the present calculation are chosen to be the same as those in Li et al. (2019) or

(4.1a–f)$$\begin{gather} E = 5\,{\rm GPa},\quad \nu = 0.3, \quad \rho_i=922.5\, {\rm kg}\,{\rm m}^{{-}3}, \quad \rho=1025\,{\rm kg}\,{\rm m}^{{-}3}, \nonumber\\ H = 40\, {\rm m},\quad g = 9.81\,{\rm m}\,{\rm s}^{{-}2}. \end{gather}$$

The computations are conducted based on the parameters in (4.1a–f), unless otherwise stated. All the numerical results are presented in a non-dimensionalised form as used in § 3. The numerical results are calculated by truncating the infinite series in (3.19) and (3.26) at $m=200$, as well as the series of multipole terms in (3.59) at $n=6$, which has been verified to ensure convergence of the results.

4.1. Analysis of the distribution of the roots of $K_2(\alpha,\alpha F)=0$

The distribution of the roots corresponding to $K_2(\alpha,\alpha F)=0$ has been proved analytically in Appendix A. Here, we give a graphic description in figure 2. When $0< F< F_c$, there are four fully complex roots that are conjugate to each other and exhibit symmetrical distribution in the complex plane, as illustrated in figure 2(a,c). As $F$ approaches $F_c$, the four complex roots move towards the real axis of $\alpha$, and eventually become two real double roots at $F=F_c$ (as shown in figure 2a–c). In such a case, the corresponding dispersion equation in (3.13b) should satisfy $K_2(\pm \kappa _c,\pm \kappa _cF_c)={\partial K_2(\pm \kappa _c,\pm \kappa _c F_c)}/{\partial \alpha }=0$. The Green function at $F=F_c$ will be infinite. Direct numerical solution at this point is not practical. One way to treat this is to modify the equation as in Liu & Yue (1993) and Yang, Wu & Ren (2022). Here we perform the calculation only at $F$ sufficiently close to $F_c$. When $F>F_c$, the two real double roots become four different real roots, as illustrated in figure 2(d,e). As $F$ continues to increase and tends to $1$, two of the real roots gradually approach $0$. When $F>1$, these two roots become purely imaginary and are located between $0$ and $\pm ({{\rm \pi} }/{2})\textrm {i}$, respectively.

Figure 2. The contours of $\log [K_2(\alpha,\alpha F)]$ at different at different Froude number, the roots of $K_2(\alpha,\alpha F)=0$ are marked by red solid circles: (a) $F=0.2$; (b) $F=0.6$; (c) $F=0.78$; (d) $F=0.79$; (e) $F=0.9$; (f) $F=0.99$; (g) $F=1.01$; (h) $F=1.1$; (i) $F=1.6$. Here $h_i=1/40$, $D=1.78\times 10^{-2}$ and $F_c=0.7868$.

4.2. Verification of the Green function

The Green function derived in § 3 using the Wiener–Hopf technique is validated through a comparison with that obtained by the method of MEE, where the detailed derivation is provided in Appendix C. A comparison of the wave elevation induced by a single source at three typical Froude numbers is given in figure 3. It can be seen that there is no visible difference between the results by these two approaches, which confirms that the Green function derived by Wiener–Hopf is consistent with that by MEE.

Figure 3. A comparison of the wave elevations by the Wiener–Hopf technique and the method of MEE: (a) $F=0.3$; (b) $F=0.9$; (c) $F=1.2$. Here $h_i=1/40$, $D=1.78\times 10^{-2}$, $F_c=0.7868$, $x_0 = 1$ and $z_0 = -0.5$.

Figure 3 also shows that when $F< F_c$ (figure 3a), there is no travelling wave in the ice sheet and there is a wave of $\mathcal {k}_0$ in the free surface region. When $F_c< F<1$ (figure 3b), in the ice sheet region, there is a travelling wave of $\kappa _{-1}$ in the upstream. The downstream $\kappa _{0}$ wave is not evident because the region $0< x< x_0=1$ is small. In the free surface region, there is still the $\mathcal {k}_0$ wave. When $F>1$ (figure 3c), the only travelling wave is $\kappa _{-1}$ at upstream. All these are consistent with the previous analysis. In addition, there is a mean elevation at $x\rightarrow \pm \infty$, which is induced by the pole at $\alpha =0$ in the integrands of (3.38c) and (3.40c), and its value can be evaluated as $\pm {\rm \pi}F/(1-F^2)$ by using the theorem of residue. This phenomenon is due to the fact that $G$ is due to a source which generates a net flow into the fluid. For a dipole obtained by taking derivatives with respect to $x_0$ and $z_0$, this term becomes a constant, which no longer generates a mean free surface elevation and mean current. For higher-order multipole, this term disappears.

4.3. Hydrodynamic forces on a submerged circular cylinder

We now consider the scenarios of a circular cylinder submerged beneath the water surface. The submergence of the cylinder $z_0/a=-4$ is used here. The resistance $F_R$ and lift $F_L$ on the cylinder against $F$ at $6$ different values of $D$ ($h_i/a$) are given in figure 4, where $x_0/a=8$ and particularly $D=0$ corresponds to the results of the free surface case. For the resistance, we may apply a procedure similar to that in appendix C of Yang, Wu & Ren (2021) (although the symmetry conclusion of the Green function in appendix B is wrong, due to the mistake that the real part of the complex function is not taken, the procedure in appendix C is correct) to obtain the far-field formula

(4.2)\begin{equation} F_R = I_0 + I_{+\infty} + I_{-\infty}, \end{equation}

where

(4.3) \begin{align} \left. \begin{aligned} I_0 & = \frac{1}{2}[(\eta^2)_{x=0^+}-(\eta^2)_{x=0^-}]\\ I_{+\infty} & =\left\{\begin{gathered} \frac{1}{2}\int_{{-}1}^{0}\left(\frac{\partial \phi}{\partial x} \frac{\partial \phi}{\partial x}-\phi\frac{\partial^2\phi}{\partial x^2}\right)_{x={+}\infty}\,{\rm d}z -\frac{1}{2}\left(\frac{D}{F}\frac{\partial^4\phi}{\partial x^3\partial z}+ F\frac{\partial \phi}{\partial x}\right)^2_{x={+}\infty,z=0}\\ -\frac{1}{2F}\left[\frac{\partial\phi}{\partial z}\left(F\phi+\frac{D}{F}\frac{\partial^3\phi}{\partial x^2\partial z}\right)\right]_{x={+}\infty,z=0}+\frac{D}{2F^2}\left[\left(\frac{\partial^2\phi} {\partial x\partial z}\right)^2_{x={+}\infty,z=0} - \left(\frac{\partial \phi}{\partial z} \frac{\partial^3\phi}{\partial x^2\partial z}\right)\right] \end{gathered} \right\}\\ I_{-\infty}&={-}\frac{1}{2}\int_{{-}1}^{0}\left(\frac{\partial \phi}{\partial x}\frac{\partial \phi}{\partial x}-\phi\frac{\partial^2\phi}{\partial x^2}\right)_{x={-}\infty}\,{\rm d}z\nonumber\\ &\quad +\frac{1}{2}F^2\left(\frac{\partial \phi}{\partial x}\right)^2_{x={-}\infty,z=0}+\frac{1}{2}\left(\phi\frac{\partial \phi}{\partial z}\right)_{x={-}\infty,z=0}. \end{aligned} \right\} \end{align}

Here $I_0=0$ because the Kutta condition is satisfied at $x=0$; and $I_{\pm \infty }$ can be determined by using the far-field expression of $\phi$. From (3.57) and (3.59), we obtain

(4.4)\begin{equation} \phi(x,z)=\left\{ \begin{array}{ll} H(1-F){\rm Re}\left\{\chi_{-}(\mathcal{k}_0){\rm e}^{-{\rm i}\mathcal{k}_0 x}+\chi_{-}(-\mathcal{k}_0) {\rm e}^{{\rm i}\mathcal{k}_0x}\right\}C(z,\mathcal{k_0}), & x\rightarrow -\infty,\\ H(F-F_c){\rm Re}\left\{\chi_{+}(\kappa_{{-}1}){\rm e}^{-{\rm i}\kappa_{{-}1} x}+\chi_{+}(-\kappa_{{-}1}) {\rm e}^{{\rm i}\kappa_{{-}1}x}\right\}C(z,\kappa_{{-}1}), & x\rightarrow +\infty, \end{array} \right. \end{equation}

where $H(x)$ denotes the Heaviside step function, and

(4.5) \begin{equation} \left. \begin{aligned} \chi_{-}(\alpha) & =2{\rm \pi}{\rm i}\sum_{n=1}^{+\infty}\frac{a^nA_n}{(n-1)!}\frac{\hat{N}_n(\alpha,x_0,z_0)}{\mathcal{K}_1^{\prime}(\alpha,F)K_+(\alpha,F)}\\ \chi_{+}(\alpha) & ={-}2{\rm \pi}{\rm i}\sum_{n=1}^{+\infty}\frac{a^nA_n}{(n-1)!}\\ &\quad \times\left\{ \begin{aligned} & \frac{\hat{N}_n(\alpha,x_0,z_0)K_{-}(\alpha,F)}{\mathcal{K}_2^{\prime}(\alpha,F)}\\ & \quad +\frac{({-}1)^n\alpha^{n-2}{\rm e}^{-|\alpha|}\left(D\alpha^4 + 1 + F^2|\alpha|\right) \exp({-\alpha(z_0+1)+{\rm i}\alpha x_0})}{\mathcal{K}_2^{\prime}(\alpha,F)} \end{aligned} \right\} \end{aligned} \right\}. \end{equation}

Substituting (4.4) into (4.3) and using (4.2), we have $F_R$ as

(4.6)\begin{equation} F_R=\left\{ \begin{aligned} & H(F-F_c)\frac{\left|\chi_{+}(\kappa_{{-}1})+\bar{\chi}_{+}(-\kappa_{{-}1})\right|^2\sinh \kappa_{{-}1}}{4F^2}\frac{\partial K_2(\kappa_{{-}1},\kappa_{{-}1}F)}{\partial \alpha}\\ & \quad - H(1-F)\frac{\left|\chi_{-}(\mathcal{k}_0)+\bar{\chi}_{-}(-\mathcal{k}_0)\right|^2\sinh \mathcal{k}_0}{4F^2}\frac{\partial K_1(\mathcal{k}_0,\mathcal{k}_0F)}{\partial \alpha} \end{aligned} \right\}. \end{equation}

As presented in § 3.1, ${\partial K_1(\mathcal {k}_0,\mathcal {k}_0 F)}/{\partial \alpha }<0$ and ${\partial K_2(\kappa _{-1},\kappa _{-1} F)}/{\partial \alpha }>0$ and, thus, in (4.6) $I_{\pm \infty } >0$. The far-field results can be used to verify the near-field ones, as given in figure 4 for $D=1.78\times 10^{-2}$ ($h_i/a=0.2$). In addition, it can be seen that with the decrease of $D$, the influence of the ice sheet gradually becomes weak, which leads the results to become consistent with that in the free surface ocean. When $D=1.78\times 10^{-5}$ ($h_i/a=0.02$), the curves are nearly identical to those of $D=0$.

Figure 4. Forces on a circular cylinder at different values of $D$ and $h_i$: (a) resistance and (b) lift. Here $H/a=8$, $z_0/a=-4$ and $x_0/a=8$.

In figure 4, it can be seen that both $F_R$ and $F_L$ are quite small but non-zero when $F$ is small. This behaviour is different from the case fully covered by an ice sheet (Li et al. 2019), where $F_R=0$ when $F< F_c$. Here $\phi$ is a wavy function as $x\rightarrow -\infty$ when $F<1$, which makes $I_{-\infty } \neq 0$. Consequently, $F_R$ is non-zero. When $F$ tends to $F_c$, the imaginary parts of $\kappa _0$ and $\kappa _{-1}$ become very small. The downstream $\kappa _0$ wave will decay very slowly and it will then affect the free surface region significantly, where there will be a travelling wave. As a result, as $F\rightarrow F_c$, a rapid rise can be found in $F_R$ as well as $F_L$. In fact, a peak can be seen when $F$ is very close to $F_c$. After that, a drop can be observed from the curves of $F_R$ and $F_L$, and $F_L$ changes from positive (attraction to the free surface) to negative (repulsion). When $F\rightarrow 1$, noticeable abrupt changes can be observed from $F_R$ and $F_L$, which is due to the fact that the downstream $\kappa _0$ wave in the ice sheet region and $\kappa _0$ wave in the free surface disappear.

The hydrodynamic forces on the circular cylinder are also affected by the longitudinal position $x_0$ of the cylinder. The resistance and lift on a cylinder at $x_0>0$ are given in figure 5. When $F< F_c$, if $x_0$ moves towards $x=+\infty$, $F_R$ gradually decreases and $F_L$ gradually increases. At $x_0/a=80$, $F_R$ and $F_L$ are nearly visually identical to those of the case fully covered by an ice sheet (Li et al. 2019). Such an asymptotic behaviour of $F_R$ and $F_L$ can be also observed when $F>1$, and is already achieved when $x_0/a=8$. By contrast, when $F_c<1< F$, the behaviour of $F_R$ and $F_L$ vs $F$ at different $x_0$ are quite different. As $x_0$ gradually increases, a highly and persistent oscillatory behaviour can be seen in the curves of $F_R$ and $F_L$. The reason behind such phenomena can be explained from the Green function. If we let $x_0\rightarrow +\infty$ in (3.26), and note from (3.27) that only $\kappa _0$ is real at $F_c< F<1$, which provides

(4.7)\begin{align} M_{-}(\alpha,x_0,z_0)&\rightarrow H\left((F-F_c)(1-F)\right)\frac{DF^2\kappa_0^2C(z_0,\kappa_0)}{\mathcal{K}_2^{\prime}(\kappa_0,F)}\nonumber\\ &\quad \times\left[\frac{K_{+}(\kappa_{0},F){\rm e}^{{\rm i}\kappa_0x_0}}{\alpha - \kappa_0}-\frac{K_+(-\kappa_0,F) {\rm e}^{-{\rm i}\kappa_0x_0}}{\alpha+\kappa_0}\right],\quad x_0\rightarrow +\infty. \end{align}

Substituting (4.7) into (3.33) and (3.34), we have the $G$ as

(4.8) \begin{align} G&\rightarrow G_{ice} - H\left((F-F_c)(1-F)\right)\frac{DF^2\kappa_0^2C(z_0,\kappa_0)}{\mathcal{K}_2^{\prime}(\kappa_0,F)}\nonumber\\ &\quad \times\int_{-\infty}^{+\infty}\left[\frac{K_{+}(\kappa_{0},F){\rm e}^{{\rm i}\kappa_0x_0}}{\alpha \!-\! \kappa_0}-\frac{K_+(-\kappa_0,F) {\rm e}^{-{\rm i}\kappa_0x_0}}{\alpha+\kappa_0}\right]\frac{\alpha K_{-}(\alpha,F)C(z, \alpha) {\rm e}^{-{\rm i}\alpha x}}{\mathcal{K}_2(\alpha,F)}\,{\rm d}\alpha,\nonumber\\ &\quad\quad\quad x_0 \rightarrow +\infty. \end{align}

From (4.8), it can be found that $G\rightarrow G_{ice}$ when $0< F< F_c$ and $F>1$, which makes $F_R$ and $F_L$ in figure 5 tend to those in fluid fully covered by an ice sheet. However, when $F_c< F<1$, an additional term with oscillatory components $\exp ({\pm \textrm {i}\kappa _0x_0})$ always exist in $G$. In fact, $\kappa _0$ is an implicit function of $F$, which can be seen from the dispersion equation in (3.13b). In such a case, oscillatory behaviours are expected in the curves of $F_R$ and $F_L$ vs $F$, as shown in figure 5. The larger $x_0$ is, the more oscillatory the result will be.

Figure 5. Forces on a circular cylinder at different values of $x_0$ ($x_0>0$): (a) resistance and (b) lift. Here $H/a=8$, $z_0/a=-4$, $h_i/a=0.2$ and $D=1.78\times 10^{-2}$.

Investigations are also conducted for the resistance and lift on a circular cylinder at $x_0<0$, as presented in figure 6. The properties of $F_R$ and $F_L$ are very different from those in figure 5 at $x_0>0$. In particular, $F_R$ and $F_L$ vary smoothly in the entire range of $F$ except near $F=F_c$ and $F=1$. In addition, when $x_0/a\leq -4$, there is hardly any visible difference with the results in the fully free surface. In fact, we have

(4.9)\begin{equation} M_+(\alpha,x_0,z_0)\rightarrow 0,\quad x_0 \rightarrow -\infty. \end{equation}

Substituting (4.9) into (3.35) and (3.36), it provides

(4.10)\begin{equation} G\rightarrow G_{water},\quad x_0\rightarrow -\infty. \end{equation}

Hence, $F_R$ and $F_L$ always tend to those in the free surface problem when $x_0\rightarrow -\infty$.

Figure 6. Forces on a circular cylinder at different values of $x_0$ ($x_0<0$): (a) resistance and (b) lift. Here $H/a=8$, $z_0/a=-4$, $h_i/a=0.2$ and $D=1.78\times 10^{-2}$.

The above result can also be understood from the physical point of view. When $x_0>0$, there will be a travelling wave $\kappa _0$ behind the cylinder when $F_c< F<1$. When this wave arrives at the ice edge, it will be reflected and come to the cylinder. It will be further reflected by the cylinder and go back to the cylinder. This forward and backward process leads to the oscillatory behaviour. When $x_0<0$, there is no free surface travelling wave far ahead of the cylinder and there will be no reflection by the ice sheet edge. Thus, when $|x_0|$ is sufficiently large, the result tends to that corresponding to the free surface.

According to Tuck (1965), the free surface nonlinear effects for a circular cylinder may be ignored for the case considered in figures 4–6. Here, to acquire some more in-depth understanding of the nature of the hydrodynamic forces at different Froude numbers, we may further use some smaller submergence to highlight the features, as shown in figure 7. It can be seen that the hydrodynamic forces on the cylinder increase rapidly generally, as $|z_0|/a$ reduces.

Figure 7. Forces on a circular cylinder at different values of $z_0$: (a) resistance and (b) lift. Here $H/a=8$, $x_0/a=8$, $h_i/a=0.2$ and $D=1.78\times 10^{-2}$.

4.4. Wave profiles generated by a circular cylinder

The wave profiles $\eta (x)$ with $x_0>0$ at $0< F< F_c$ are plotted in figure 8. It can be seen that there is a regular wave with wavenumber $\kappa _0$ in the free surface region, whereas in the region covered by an ice sheet, the wave amplitude decays very quickly. Such a phenomenon is due to the fact that travelling wave cannot exist in the ice sheet when $F< F_c$. With the increase of $F$, although the amplitude of the entire wave gradually increases, compared with results in fluid fully covered by an ice sheet (Li et al. 2019), a very large wave trough near $x=x_0$ occurred before $F_c$ is not observed here. The wave profiles at $F_c< F<1$ is given in figure 9. When $F$ is near $F_c$, two distinct waves become evident. One with wavenumber $\kappa _0$ in the free surface region and the other with wavenumber $\kappa _{-1}$ ahead of the circular cylinder, whereby the amplitude of the wave at $x < 0$ surpasses marginally that at $x > 0$. As $F$ continues to increase, the amplitude of the wave at $x>x_0$ progressively decreases, whereas that at $x<0$ increases. When $F\rightarrow 1^-$, the wave component $\mathcal {k}_0\rightarrow 0$, its wavelength tends to infinity and its wave amplitude becomes very large. In figure 10, the wave profile at $F>1$ is presented. In such a case, the wave components $\mathcal {k}_0$ and $\kappa _0$ disappear from the downstream region, whereas the wave component $\kappa _{-1}$ remains in the upstream region. When $F\rightarrow 1^+$, a significant wave elevation is observed near $x=8a$, which corresponds to the longitudinal position of the centre of the circular cylinder. As $F$ continues to increase, this large wave elevation becomes smaller. Similar phenomenon is also observed in Li et al. (2019).

Figure 8. Wave profile for $0< F< F_c$. Here $H/a=8$, $z_0/a=-4$, $x_0/a=8$, $h_i/a=0.2$ and $D=1.78\times 10^{-2}$.

Figure 9. Wave profile for $F_c< F<1$. Here $H/a=8$, $z_0/a=-4$, $x_0/a=8$, $h_i/a=0.2$ and $D=1.78\times 10^{-2}$.

Figure 10. Wave profile for $F>1$. Here $H/a=8$, $z_0/a=-4$, $x_0/a=8$, $h_i/a=0.2$ and $D=1.78\times 10^{-2}$.

There are many similarities between the waves here and those in figure 3 due to a source. There is, however, one distinctive difference. Here we do not have a marked mean surface elevation at infinity. The reason for this is that for a source there is net flow into the fluid. For a cylinder, the leading term will be a dipole and there will be no net flow into the fluid. Thus, there is no mean free surface elevation.

5. Conclusions

The interaction of a uniform current with a circular cylinder submerged in a fluid covered by a semi-infinite ice sheet has been investigated analytically. The solution procedure is based on the linearised velocity theory for fluid and the Kirchhoff–Love plate theory for ice sheet. The Green function of the problem has first been derived through the Wiener–Hopf technique, which has been verified with the result by the method of MEE and has been found to be consistent. Based on this, the velocity potentials for multipoles have been obtained directly by differentiating the position of the source directly, which have then been used to construct the velocity potential due to a submerged circular cylinder.

When deriving the Green function, understanding of the distribution of the roots of the dispersion equation for a fluid fully covered by a homogeneous ice sheet is crucial. The root distribution here is quite different from that in the periodic problems of wave radiation and diffraction, and related to the Froude number. In particular, for the depth-based Froude number $F$, when $0< F< F_c$, where $F_c$ is the critical Froude number, there are four symmetrical complex roots $\pm \kappa _{-1}$ and $\pm \kappa _{0}$ with $\bar {\kappa }_{-1}=-\kappa _{0}$ and an infinite number of purely imaginary roots $\kappa _m$ ($m=1,2,\ldots$). When $F_c< F<1$, $\kappa _{-1}$ and $\kappa _{0}$ become two purely positive real roots with $\kappa _{-1}>\kappa _{0}$. When $F>1$, $\kappa _{0}$ becomes a purely positive imaginary root but $\kappa _{-1}$ remains to be real.

From the solution, in addition to satisfying the free edge conditions, the obtained Green function also satisfies the Kutta condition at the ice edge, which ensures the continuity of both the wave elevation and slope at $x=0$. Such a phenomenon is quite different from the problem of an incoming wave, where jumps on wave elevation and slope are allowed at the ice edge. The hydrodynamic forces on a submerged circular cylinder under different flexural rigidity ($D$) of the ice sheet shows that the result is consistent with those in the free surface problem when $D\rightarrow 0$. Different from the conclusion that the resistance is zero when $F< F_c$ when water surface is fully covered by the ice sheet (Li et al. 2019). Here, $F_R$ is never zero because the free surface wave exists even when $F< F_c$.

On the curves of resistance $F_R$ and lift $F_L$ vs $F$, a peak is observed when $F$ is close to $F_c$, whereas a sudden jump is found at $F=1$, which is similar to those in the case fully covered by an ice sheet. Here $F_R$ and $F_L$ are also affected by the longitudinal position $x_0$ of the cylinder. In particular, when $x_0\rightarrow +\infty$, $F_R$ and $F_L$ tend to those in fluid fully covered by an ice sheet only when $0< F< F_c$ and $F>1$. When $F_c< F<1$, a highly oscillatory behaviour is observed, which is actually induced by successive wave reflection due to the cylinder and the ice edge. Such reflection does not exist in other ranges of $F$. By contrast, $F_R$ and $F_L$ vary smoothly when $x_0<0$, and they always tend to the results in free surface problem when $x_0\rightarrow -\infty$.

The travelling free surface wave always exists when $0< F<1$. There exists no travelling wave related to the ice sheet when $0< F< F_c$. When $F_c< F<1$, there will be two travelling waves, or $\kappa _0$ and $\kappa _{-1}$ waves, related to the ice sheet. When the body is submerged below the ice sheet, $\kappa _{-1}$ wave will be before the body whereas the $\kappa _0$ wave will be behind the body. The latter will reach the ice sheet edge, transmit into the free surface, and also be reflected back to the cylinder. When the body is submerged below the free surface, only $\kappa _{-1}$ wave will travel to far upstream of the ice sheet region. When $F>1$, the only travelling wave is $\kappa _{-1}$ wave which will propagate into far upstream.

It is worth mentioning that the derived Green function can also be used to construct the boundary integral equation for bodies with arbitrary shapes. In addition, the formulation here can be extended easily to other types of edge conditions. Moreover, the solution procedure can be further applied to ice sheets with imperfections, such as cracks.