2. Mathematical model

We consider Kirchhoff–Love plate theory with no variation in $y$ (an infinitely wide plate), so the problem is one-dimensional and equivalent to Euler–Bernoulli beam theory, with the flexural rigidity scaled by $1-\nu ^{2}$ , where $\nu$ is Poisson’s ratio. We work on a one-dimensional bounded domain $x \in [-L,L]$ : the left-hand half, $x \in (-L,0)$ , contains a fluid with a free surface, and the right-hand half, $x \in (0,L)$ , contains a fluid overlain by an ice shelf. The thin plate approximation adopted here is consistent with the classical framework for floating-ice flexure (e.g. Timoshenko & Goodier Reference Timoshenko and Goodier1934; Sergienko Reference Sergienko2013). For simplicity, we assume that the ice shelf spans exactly half the domain, though this assumption can be readily extended to more general cases. The fluid and ice shelf are assumed to have constant depth and constant thickness, respectively. We note that in reality, the ice shelf has variable properties in depth (as well as variation in thickness), and all parameters should be seen as representative averages. We also apply the assumption of shallow water, which requires the length scale of any motion to be much greater than the water depth and the ice shelf thickness (since we apply the shallow-water assumption in the open water region). Within the shallow-water framework, the motion of the system is described by the vertical displacement $W(x,t)$ and the fluid velocity potential $\varPhi (x,t)$ . In particular, the vertical displacement $W(x,t)$ represents the displacement of the free surface in the ocean region $(-L,0)$ , and the vertical motion of the ice shelf in the covered region $(0,L)$ . Unlike Meylan (Reference Meylan2002), who considered an elastic plate in an infinite domain, we study an ice shelf in a finite domain. This choice yields a discrete eigenfrequency spectrum, which makes resonances straightforward to identify and resolve; in contrast, the infinite-domain setting leads to a continuous spectrum and scattering resonances that are much harder to isolate numerically. In this framework, $W$ and $\varPhi$ satisfy the following coupled system of partial differential equations derived from shallow-water wave theory and thin plate approximation (Meylan Reference Meylan2002):

(2.1a) \begin{align} D \partial _x^4 W +\rho _I h \partial _t^2 W+\rho _w g W&=- \rho _w \partial _t \varPhi,\quad x\in (0,L), \end{align}

(2.1b) \begin{align} g W &= -\partial _t\varPhi,\qquad\, x\in (-L,0), \end{align}

(2.1c) \begin{align} \partial _x^2 \varPhi &= -\frac {1}{H}\, \partial _t\! W,\quad x\in (-L,0), \end{align}

(2.1d) \begin{align} \partial _x^2 \varPhi &= -\frac {1}{\widetilde {H}}\, \partial _t\! W,\quad x\in (0,L), \end{align}

where $D= {E h^3}/{12(1-\nu ^2)}$ is the ice shelf’s flexural rigidity. Moreover, $\rho _w$ and $\rho _I$ are the densities of water and ice, respectively, which we assume to be $\rho _w=1025$ kg m $^{-3}$ and $\rho _I=922.5$ kg m $^{-3}$ . The thickness of the ice shelf is denoted $h$ , and the depth of the water is denoted $H$ , which implies that the vertical distance between the submerged face of the ice shelf and the seabed (denoted $\widetilde {H}$ ) can be computed as $\widetilde {H}=H-( {\rho _I}/{\rho _w}) h$ . Moreover,  $E$ is the Young’s modulus of the ice shelf, which we take to be $11$ GPa, and the Poisson’s ratio $\nu$ is taken to be 0.33. Acceleration due to gravity is denoted $g=9.81$ m s $^{-2}$ , and the geometry of the problem is illustrated in figure 1.

Figure 1. Schematic of the problem geometry. The domain consists of a bounded ice shelf of length $L$ above a fluid-containing sub-shelf cavity, with a bounded open-water region to the left of the shelf. The coordinate $x=0$ marks the seaward edge of the ice shelf, while $x=L$ denotes its land-connected edge.

Note that (2.1a ) would be equivalent to (2.1b ) if $D$ and $\rho _Ih$ were zero. We specify the following initial conditions at time $t = 0$ :

(2.2a) \begin{align} W(x, 0) = f(x), \end{align}

(2.2b) \begin{align} \frac {\partial\! W}{\partial t}(x, 0) = g(x). \end{align}

At the left-hand boundary ( $x=-L$ ), we assume no flux into the fluid:

(2.3) \begin{align} \frac {\partial \varPhi }{\partial x}\bigg |_{x=-L} = 0. \end{align}

At the interface $x = 0$ between the open water and the ice shelf, two matching conditions and two boundary conditions are applied to ensure appropriate physical behaviour. The matching conditions at $x=0$ , which enforce continuity of pressure and continuity of horizontal flux, are

(2.4a) \begin{align} \varPhi (0^-) &= \varPhi (0^+), \end{align}

(2.4b) \begin{align} H \frac {\partial \varPhi }{\partial x}\bigg |_{x=0^-} &= \widetilde {H} \frac {\partial \varPhi }{\partial x}\bigg |_{x=0^+}, \end{align}

respectively. The ice shelf has a free edge at $x = 0$ , meaning that the bending moment and shear force vanish, implying

(2.4c) \begin{align} \frac {\partial ^4 \varPhi }{\partial x^4}\bigg |_{x=0} = 0, \end{align}

(2.4d) \begin{align} \frac {\partial ^5 \varPhi }{\partial x^5}\bigg |_{x=0} = 0, \end{align}

respectively. At the right-hand boundary $x = L$ , the ice shelf is treated as a clamped elastic plate, and there is assumed to be no horizontal fluid flux in the sub-ice-shelf cavity, yielding the following three boundary conditions:

(2.4e) \begin{align} \frac {\partial \varPhi }{\partial x}\bigg |_{x=L} = 0, \end{align}

(2.4f) \begin{align} \frac {\partial ^2 \varPhi }{\partial x^2}\bigg |_{x=L} = 0, \end{align}

(2.4g) \begin{align} \frac {\partial ^3 \varPhi }{\partial x^3}\bigg |_{x=L} = 0. \end{align}

These initial and boundary conditions fully specify the problem in the time domain. We note that treating the grounded edge of the plate as clamped is a significant simplification of the underlying physical system, and that the present work could be extended by deriving more realistic boundary conditions. We begin by solving for the modes of vibration of the system. Denoting the angular frequency of vibration $\omega$ , the modal displacement and velocity potential of the system are of the form $W(x,t)=\text{Re}(w(x)\,\textrm{e}^{-\mathrm{i} \omega t})$ and the fluid velocity potential $\varPhi (x,t)=\text{Re}(\phi (x)\,\textrm{e}^{-\mathrm{i} \omega t})$ , respectively. Substituting these expressions into (2.1) gives

(2.5a) \begin{align} D \partial _x^4 w - \omega ^2 \rho _I h w +\rho _w g w&= \mathrm{i} \omega \rho _w \phi,\quad x\in (0,L), \end{align}

(2.5b) \begin{align} g w&=\textrm{i} \omega \phi,\quad x\in (-L,0), \end{align}

(2.5c) \begin{align} \partial _x^2 \phi &= \frac {\textrm{i} \omega }{H} w,\quad x\in (-L,L). \end{align}

We use (2.5c ) to eliminate $w$ from (2.5a ) and (2.5b ), which gives rise to the following sixth-order equations for $\phi$ :

(2.6a) \begin{align} D \partial _x^6 \phi -\omega ^2 \rho _I h \partial _x^2 \phi +\rho _w g \partial _x^2 \phi &= -\frac {\omega ^2 \rho _w} {\widetilde {H}} \phi,\quad x\in (0,L), \end{align}

(2.6b) \begin{align} \rho _w g \partial _x^2 \phi &= -\frac {\omega ^2 \rho _w}{H} \phi,\quad x\in (-L,0), \end{align}

which must be solved in tandem with the boundary conditions (2.4). Equation (2.6) can be restated as

(2.7) \begin{align} -\omega ^2 \phi = \begin{cases} gH \partial _x^2 \phi , & x\in (-L,0),\\[3pt] \dfrac {\widetilde {H}}{\rho _w}(D \partial _x^6 \phi + ( \rho _w g -\omega ^2 \rho _I h)\, \partial _x^2 \phi ), & x\in (0,L). \end{cases} \end{align}

This is an eigenvalue problem, in which the eigenvalue $-\omega ^2$ is the negative square of the angular frequency of vibration, and the associated eigenfunctions describe the modes.

3. The discrete spectrum

3.1. Solution for a bounded ice shelf

The general solution to (2.7) for $\omega \gt 0$ is

(3.1a) \begin{align} \phi (x,\omega ) &= \begin{cases} R_1\, \textrm{e}^{\mathrm{i} k x}+R_2\, \textrm{e}^{-\mathrm{i} k x},&x\in (-L,0),\\ \displaystyle {\sum _{n=1}^{3} \beta _n \,\textrm{e}^{-r_n (x-L)}+\sum _{n=1}^{3} \alpha _n\, \textrm{e}^{r_n x}},&x\in (0,L),\\ \end{cases} \end{align}

(3.1b) \begin{align} w(x,\omega ) &= \begin{cases} \displaystyle \frac {-H k^2}{\mathrm{i} \omega } \big ( R_1\, \textrm{e}^{\mathrm{i} k x} + R_2\, \textrm{e}^{-\mathrm{i} k x} \big ), & x \in (-L,0), \\[3pt] \displaystyle \sum _{n=1}^{3} \frac {\widetilde {H}\! r_n^2}{\mathrm{i} \omega } \big ( \beta _n\, \textrm{e}^{-r_n (x-L)} + \alpha _n \, \textrm{e}^{r_n x} \big ), & x \in (0,L), \end{cases} \end{align}

where the wavenumber $k$ is the positive solution to the shallow-water dispersion relation

(3.2) \begin{align} k = \frac {\omega }{\sqrt {g H}}, \end{align}

and $r_n$ are the six roots of the dispersion relation for the ice-covered region:

(3.3) \begin{align} \frac {D}{\rho _w g}r^6+\left(1-\omega ^2 \frac {\rho _I h}{\rho _w g}\right) r^2+\frac {\omega ^2}{\widetilde {H} g}=0. \end{align}

The roots of (3.3) appear in positive and negative pairs. The purely imaginary root was chosen to have a positive imaginary component and is denoted $r_3$ , while the roots with non-zero real parts were chosen to have negative real components and are denoted $r_{1}$ and $r_2$ . The root $r_3$ is the propagating hydroelastic wave, while $r_1$ and $r_2$ are decaying evanescent modes from the interface. Equation (3.1) has eight unknown coefficients, which must be determined from the boundary conditions (2.4). The no-flux condition (2.3) at $x=-L$ implies that

(3.4a) \begin{align} \mathrm{i} k (R_1\, \textrm{e}^{-\mathrm{i} k L}-R_2\, \textrm{e}^{\mathrm{i} k L})=0. \end{align}

The continuity of the potential equation condition at $x=0$ implies that

(3.4b) \begin{align} R_1+R_2-\sum _{n=1}^{3} \alpha _n -\sum _{n=1}^{3} \beta _n\, \textrm{e}^{r_n L}=0, \end{align}

and the continuity of flux condition (2.4b ) implies that

(3.4c) \begin{align} H \mathrm{i} k R_1 - H \mathrm{i} k R_2-\sum _{n=1}^{3} \widetilde {H} r_n \alpha _n+\sum _{n=1}^{3} \widetilde {H} r_n \beta _n\, \textrm{e}^{r_n L}=0. \end{align}

The free-edge conditions (2.4c ) and (2.4d ) imply

(3.4d) \begin{align} \sum _{n=1}^{3} r_n^4 \alpha _n+\sum _{n=1}^{3} r_n^4 \beta _n\, \textrm{e}^{r_n L}&=0, \end{align}

(3.4e) \begin{align} \sum _{n=1}^{3} r_n^5 \alpha _n-\sum _{n=1}^{3} r_n^5 \beta _n\, \textrm{e}^{r_n L} &=0, \\[9pt] \nonumber \end{align}

respectively. Finally, the right-hand boundary at $x = L$ , the three imposed conditions (2.4e ), (2.4f ) and (2.4g ) yield

(3.4f) \begin{align} \sum _{n=1}^{3} r_n \alpha _n\,\mathrm{e}^{r_n L}-\sum _{n=1}^{3} r_n \beta _n &=0, \end{align}

(3.4g) \begin{align} \sum _{n=1}^{3} r_n^2 \alpha _n\,\mathrm{e}^{r_n L}+\sum _{n=1}^{3} r_n^2 \beta _n&=0, \end{align}

(3.4h) \begin{align} \sum _{n=1}^{3} r_n^3 \alpha _n\,\mathrm{e}^{r_n L}-\sum _{n=1}^{3} r_n^3 \beta _n&=0, \end{align}

respectively. Equations (3.4) can be written as a nonlinear eigenvalue problem of the form

(3.5) \begin{align} M(\omega ) \boldsymbol{X} = \mathbf{0}, \end{align}

where

(3.6) \begin{align} M(\omega ) = \begin{pmatrix} \mathrm{i} k \textrm{e}^{-\mathrm{i} k L} & -ik \textrm{e}^{\mathrm{i} k L} & 0 & 0 & 0 & 0 & 0 & 0 \\[3pt] 1 & 1 & -1 & -1 & -1 & -e^{r_1 L} & -e^{r_2 L} & -e^{r_3 L} \\[3pt] H \mathrm{i} k & -H \mathrm{i} k & -\widetilde {H} r_1 & -\widetilde {H} r_2 & -\widetilde {H} r_3 & \widetilde {H} r_1 e^{r_1 L} & \widetilde {H} r_2 e^{r_2 L} & \widetilde {H} r_3 e^{r_3 L} \\[3pt] 0 & 0 & r_1^4 & r_2^4 & r_3^4 & r_1^4 e^{r_1 L} & r_2^4 e^{r_2 L} & r_3^4 e^{r_3 L} \\[3pt] 0 & 0 & r_1^5 & r_2^5 & r_3^5 & -r_1^5 e^{r_1 L} & -r_2^5 e^{r_2 L} & -r_3^5 e^{r_3 L} \\[3pt] 0 & 0 & r_1 e^{r_1 L} & r_2 e^{r_2 L} & r_3 e^{r_3 L} & -r_1 & -r_2 & -r_3 \\[3pt] 0 & 0 & r_1^2 e^{r_1 L} & r_2^2 e^{r_2 L} & r_3^2 e^{r_3 L} & r_1^2 & r_2^2 & r_3^2 \\[3pt] 0 & 0 & r_1^3 e^{r_1 L} & r_2^3 e^{r_2 L} & r_3^3 e^{r_3 L} & -r_1^3 & -r_2^3 & -r_3^3 \end{pmatrix} \end{align}

and

(3.7) \begin{align} \boldsymbol{X} = \begin{pmatrix} R_1 \\ R_2 \\ \alpha _1 \\ \alpha _2 \\ \alpha _3 \\ \beta _1 \\ \beta _2 \\ \beta _3 \end{pmatrix}\!. \end{align}

Non-trivial solutions of (3.5) exist at a discrete set of natural frequencies $\omega _n$ at which $M(\omega _n)$ is singular, where $\mathbf{0}$ denotes the eight-dimensional zero vector. The coefficients within the associated eigenvectors determine the associated mode shapes $w_n(x)$ .

3.2. Numerical procedure for the nonlinear eigenvalue problem

To find solutions of the nonlinear eigenvalue problem (3.5), we first compute $\det (M(\omega ))$ on a discrete grid of values. Candidate eigenfrequencies are identified by locating subintervals within the grid where $\text{Re}(\det (M(\omega )))$ changes sign. These approximate zeros are then refined iteratively. In the refinement step, we solve the generalised eigenvalue problem

(3.8) \begin{align} M(\omega ^{(0)})\,v = \lambda M'(\omega ^{(0)})\,v. \end{align}

We numerically solve for $\lambda$ , which is taken to be the eigenvalue with the smallest absolute value. Then $\omega ^{(1)} = \omega ^{(0)} - \lambda$ provides a refinement of the initial guess $\omega ^{(0)}$ . This method can be used repeatedly until the required level of accuracy is achieved.

For each eigenfrequency $\omega _n$ obtained using this method, the associated null vector $X_n$ of $M(\omega _n)$ determines the coefficients, from which the potential $\phi$ and displacement $w$ are obtained. This iterative procedure is an adaptation of the method developed by Ooi, Nikolovska & Manasseh (Reference Ooi, Nikolovska and Manasseh2008). Further details regarding its implementation can be found in Wilks et al. (Reference Wilks, Montiel, Bennetts and Wakes2025).

3.3. The static mode

The ocean/ice shelf system supports a static mode corresponding to $ \omega _0 = 0$ . Physically, this represents the bending of the ice shelf in response to hydrostatic loading due to a constant sea level $A$ , and there is no associated velocity potential. A similar situation arises under tidal loading (see e.g. Vaughan Reference Vaughan1995). The static mode cannot be obtained in the form (3.1), which is not valid for $k=0$ . We have shown numerically that it does arise in the limit as $\omega \to 0$ , but we have not derived this limit analytically. Instead, we find the static mode by assuming that all time derivatives are zero, and using a modification of (2.1a ):

(3.9) \begin{align} D \partial _x^4 W +\rho _w g W =\rho _w g A, \quad x\in (0,L), \end{align}

and $W=A$ for $x\in (-L,0)$ . The displacement corresponding to the static mode is therefore of the form

(3.10) \begin{align} w_0(x) = \begin{cases} A, & x\in (-L,0), \\ A + C_1 \cos (s(x - L)) + C_2 \sin (s(x - L)) + C_3\,\textrm{e}^{s(x - L)} + C_4\,\textrm{e}^{-s x}, & x\in (0,L), \end{cases} \end{align}

where $ s$ is a complex constant given by

(3.11) \begin{align} s = \textrm{e}^{\textrm{i}\pi /4} \left ( \frac {\rho _w g}{D} \right )^{1/4}\!. \end{align}

To solve for the coefficients $C_1$ , $C_2$ , $C_3$ and $C_4$ , we apply four conditions. At the free edge ( $x = 0$ ), impose zero bending moment, $w''(0) = 0$ , and zero shear force, $w'''(0) = 0$ . At the clamped edge, we impose zero displacement, $w(L) = 0$ , and zero slope, $w'(L) = 0$ . These conditions yield the linear system

(3.12) \begin{align} \begin{bmatrix} 1 & 0 & 1 & \textrm{e}^{-s L} \\ 0 & 1 & 1 & -\textrm{e}^{-s L} \\ {-}\cos (s L) & \sin (s L) & \textrm{e}^{-s L} & 1 \\ {-}\sin (s L) & {-}\cos (s L) & \textrm{e}^{-s L} & -1 \end{bmatrix} \begin{bmatrix} C_1\\C_2\\C_3\\C_4 \end{bmatrix} = \begin{bmatrix} -A \\ 0 \\ 0 \\ 0 \end{bmatrix}\!. \end{align}

The coefficients $C_j$ that solve this system determine the static bending shape of the ice shelf. Figure 2 shows the static mode (labelled Mode 0) for a range of ice shelf thicknesses, among other dynamic modes. The final check of our calculation is that this form matches the one found by the numerical limit and that it also gives a resolution of the identity, i.e. that we recover the initial condition with this mode.

3.4. Orthogonality of the modes of vibration

Following Meylan (Reference Meylan2002) and Hazard & Meylan (Reference Hazard and Meylan2008), we define the inner product as

(3.13) \begin{align} \langle f_1, f_2 \rangle _{\mathcal{V}} = \rho _w g\! \int _{-L}^{L} f_1 \overline {f_2} \, \mathrm{d}x + D\! \int _0^L \frac {\partial ^2\! f_1}{\partial x^2} \frac {\partial ^2 \overline {f_2}}{\partial x^2} \, \mathrm{d}x, \end{align}

where $f_1$ and $f_2$ are arbitrary functions in the associated Hilbert space. The orthogonality of the eigenfunctions, namely

(3.14) \begin{align} \langle w_m, w_n \rangle _{\mathcal{V}} = 0 \quad \text{for} \quad \omega _m \ne \omega _n, \end{align}

plays a key role in enabling a valid modal decomposition of the solution. See Appendix A for the detailed proof of this fact.

4. Time-domain solution

The time evolution is governed by an abstract wave equation

(4.1) \begin{align} \partial _t^2 W + \mathcal{P} W=0, \end{align}

where $\mathcal{P}$ is the compact positive self-adjoint operator governing the transient dynamics of the ocean/ice shelf system. Without ice, we would simply have $\mathcal{P} = {-}\textit{gH} \partial _x^2$ , which yields the classical wave equation with wavespeed $\sqrt {gH}$ . When the ice shelf is present, the form of $\mathcal{P}$ is more complex, and we make no attempt to express it.

The eigenfunctions of $\mathcal{P}$ are those found in § 3, satisfying

(4.2) \begin{align} \mathcal{P} w_n = \omega _n^2 w_n. \end{align}

We assume the solution can be written as a modal expansion

(4.3) \begin{align} W(x,t) = \sum _{n=0}^\infty a_n(t)\, w_n(x), \end{align}

which, upon substitution into (4.1) and an application of orthogonality of $\{w_n\}$ , leads to

(4.4) \begin{align} a_n''(t) = -\omega _n^2\, a_n(t). \end{align}

For $n\gt 0$ , this ordinary differential equation has solutions

(4.5) \begin{align} a_n(t) = C_n \cos (\omega _n t) + D_n\frac {\sin (\omega _n t)}{\omega _n} , \end{align}

where $C_n$ and $D_n$ are as-yet unknown coefficients. In the case $n=0$ , (4.4) has solutions

(4.6) \begin{align} a_0(t)=C_0+D_0t. \end{align}

Thus

(4.7) \begin{align} W(x,t) = (C_0+D_0t)\, w_0(x) + \sum _{n=1}^\infty \left [ C_n \cos (\omega _n t) + D_n\frac {\sin (\omega _n t)}{\omega _n} \right ] w_n(x). \end{align}

We perform simulations using all $N_{\textit{modes}}$ modes whose angular frequencies lie in the interval $[0,1]$ s $^{-1}$ . The modal approximation, obtained from the truncated series representation, is given by

(4.8) \begin{align} W(x,t) = (C_0+D_0t)\,w_0(x) + \sum _{n=1}^{N_{\textit{modes}}} \left [ C_n \cos (\omega _n t) + D_n\frac {\sin (\omega _n t)}{\omega _n} \right ] w_n(x). \end{align}

The coefficients $C_n$ and $D_n$ are determined by the initial conditions. Since we have

(4.9) \begin{align} f(x) &= W(x,0) = \sum _{n=0}^\infty C_n w_n(x), \end{align}

(4.10) \begin{align} g(x) &= \partial _t\! W(x,0) = \sum _{n=0}^\infty D_n w_n(x), \end{align}

an application of orthonormality yields

(4.11) \begin{align} \langle f, w_n \rangle _{\mathcal{V}} = \int _0^L f(x) w_n(x) \, \mathrm{d}x, \end{align}

(4.12) \begin{align} \langle g, w_n \rangle _{\mathcal{V}} = \int _0^L g(x) w_n(x) \, \mathrm{d}x, \end{align}

since the eigenfunctions $w_n(x)$ have been normalised in the previous section, and the modal coefficients are

(4.13a) \begin{align} C_n &= \langle f, w_n \rangle _{\mathcal V}, \end{align}

(4.13b) \begin{align} D_n &= \langle g, w_n \rangle _{\mathcal V}. \end{align}

We evaluate these inner products numerically using midpoint quadrature. Suppose that the domain $(-L,L)$ is partitioned into $N$ subintervals $(x_i-\Delta x/2,x_i+\Delta x/2)$ , where $N$ is an even integer, $x_i$ is the midpoint of the $i$ th interval, and $\Delta x$ is the uniform spacing of the partition. Then we can approximate the coefficients $C_n$ and $D_n$ as

(4.14a) \begin{align} C_n &\approx \rho _w g \sum _{i=1}^N f(x_i)\, w_n(x_i) \Delta x + D \sum _{i = ({N+1})/{2}}^N f(x_i)\, \frac {\partial ^4 w_n(x_i)}{\partial x^4} \, \Delta x, \end{align}

(4.14b) \begin{align} D_n &\approx \rho _w g \sum _{i=1}^N g(x_i)\, w_n(x_i) \Delta x + D \sum _{i = ({N+1})/{2}}^N g(x_i)\, \frac {\partial ^4 w_n(x_i)}{\partial x^4} \, \Delta x. \end{align}

6. Summary

In this work, we developed a mathematical model of a bounded ice shelf coupled to the ocean, combining Kirchhoff–Love plate theory for the ice shelf with linear shallow-water wave theory for the fluid. Matching conditions were applied at the ice–water interface to ensure continuity of the velocity potential, and boundary conditions were imposed at the edges of both the ice shelf and the fluid domain. The resulting nonlinear matrix eigenvalue problem was solved numerically to obtain the natural modes and frequencies of vibration of the coupled system, which form the basis for constructing time-domain solutions via a spectral method. The time-domain simulations illustrated how transient wave packets interact with the ice shelf structure, exciting multiple modes and generating complex interference patterns through boundary reflections. These examples serve primarily to demonstrate the capabilities of the modelling framework and solution method.

The primary aim of this study was to establish a framework for solving initial value problems for a domain with partial ice shelf cover. We note that the model is highly simplified, and consideration of variable thickness, water depth, three-dimensional effects, variation in ice shelf properties with depth, and so on, needs to be incorporated to make qualitative predictions for particular ice shelves. However, any model incorporating such effects would have to include the model presented here as a particular case. Thus this work enables an extension to simulating the sequential break-up of ice shelves, a proposed model of ice shelf break-up in which wave-induced mechanical stresses cause instantaneous fracture events, resulting in disintegration over time. By extending the model to incorporate nonlinear effects (especially break-up), wave scattering, and interaction with floating sea ice or open water, further insights into ice shelf stability and retreat may be achieved. In addition, extending the model to an infinite domain would represent a major change, as it leads to a continuous spectrum of modes and requires a different theoretical framework. A possible improvement of the model is to extend it for very thick ice shelves, where additional wave modes may appear. Another direction is to use a Mindlin plate formulation to capture these effects with more precision.