Time Domain Vibration Analysis of an Ice Shelf

Abstract

A method is presented to calculate the vibrations of an ice shelf floating in shallow water under different boundary conditions. One condition is that there is no flux, which reduces all calculations and the other is that there is no pressure at the seaward end of the ice shelf. The effect of these boundary conditions is investigated in detail, and the modes of vibration are also determined. Motion simulations of the system are presented for the potential velocity of the water and the vertical displacement of the ice shelf. These are found through a numerical method, which reduces all calculations to matrix multiplication. The underlying motion is shown to be very complex and difficult to interpret from single-point response measurements. The motion of more realistic ice shelves can be expected to be even more complicated.

1. Introduction

Ice shelves are floating sheets of ice permanently attached to an ice sheet or glaciers. They can exist along the Antarctic coastline and in Arctic regions, covering an area greater than $1.561$ million ${\mathrm{km}}^2$ and occupying $75\%$ of the Antarctic coastline [1]. Various climate projections indicate the rapid melting of ice shelves around Antarctica, and there is evidence that the melting of Antarctic ice shelves has increased dramatically in recent years [1,2]. The collapse of the West Antarctic ice sheet alone could result in a global rise in sea level of up to 3 m [3]. Observation using a satellite radar altimeter for 18 years has shown that losses in the West Antarctic increased by approximately $70\%$ in the past decade [4]. Strong evidence suggests that ocean waves contribute to the calving process of the ice shelf [5,6,7,8]. Seismic instruments have measured the impact of ocean swells on the Ross Ice Shelf in Antarctica, which causes it to vibrate [2,9,10,11]. Furthermore, a tsunami in the Northern Hemisphere has been linked to an Antarctic ice shelf that calved more than 13,000 km away [12]. Improving our understanding of the interactions between ice shelves and the ocean is essential to predict the effects of climate change on coastlines and develop strategies to mitigate them. Studies conducted by Macayeal [13] and Sergienko [14] have revealed that ocean waves can cause significant ice shelf bending, stretching, and compression. Scambos [15] used theory and data to demonstrate that a better understanding of the interactions between ice shelves and waves can help predict the reaction of an ice shelf to a changing environment more accurately.

The most straightforward way to look at an ice shelf is to assume that it is floating on the surface of the sea and oscillating with the same amplitude and period of the incoming waves. The vibration of the ice shelf was first mathematically studied by [16], who proposed that the ice vibrates due to the ocean swell and that this causes the icebergs to calve from the ice shelves. The standard mathematical model of ice shelves is that of a beam or plate floating on an incompressible fluid. This method is used to find stresses and deflections and is based on the early work of Holdsworth and Glenn and Reeh [16,17]. Many subsequent investigations of this vibration have been conducted, and these may be broadly categorized as solving for the modes of vibration or for the wave-induced response. The solution for vibration modes was found in [18,19], assuming a simplified model, and for more complicated cases, by using the finite element method [20,21]. The solution to wave incident forcing was found using the finite element method [22,23,24]. A discussion of the underlying mathematical theory can be found in [25], as well as a numerical solution. More complex cases involving defects in the ice shelf are considered in [26]. Recently, the simplified model for the change in depth was investigated, and simulations were performed in the time domain [27]. We note that the emphasis here is on the time-dependent results, but our methods use the frequency domain results. These frequency domain results could be used to obtain information about the response.

In this paper, we will present the modal vibration of an ice shelf floating, assuming shallow-water approximation. We investigate two boundary conditions at the seaward end, which correspond to no flux and no pressure, which were discussed in detail in [28]. We investigate, in detail, the effect of these boundary conditions and also determine the modes of vibration. We use these modes to simulate the time-dependent motion of the ice shelf. The assumptions underlying our ice shelf model are strong, but this is the first step to simulating ice shelf dynamics in the time domain.

2. Mathematical Model

We begin with the model for the vibration of an ice shelf. We follow the work of [18,19] and assume that we can use the shallow water assumption. This assumption is highly dependent on the frequencies of motion. For ice shelves, the modes of vibration are long periods, and this assumption becomes applicable. We note that it is less applicable to a model of wave forcing, especially for swell waves for which the period, at least compared to the periods of vibration of an ice shelf, is small.

We assume that the ice shelf is free at the open water end and clamped at the landward end, again following the work in [18,19]. The ice shelf is supposed to be of length L, and the shelf and water remain in contact at all points, $x\in [0,L]$, where x is the horizontal co-ordinate. $x=0$ is the landward end, and $x=L$ is the seaward end.

We assume that the motion of the fluid can be described by a potential function, $\Phi$, that satisfies the linear shallow water equation [29]:

$$\frac{{\partial }^2\Phi }{\partial x^2}=-\frac{1}{H}\frac{\partial \eta }{\partial t}.$$

(1)

where H is the uniform depth of water below the ice shelf, and $\eta$ is the displacement of the ice shelf. The standard assumption [18,19] is that there are no-flux conditions at the ends of the cavity, which means that

$$\frac{\partial \Phi }{\partial x}=0,\mathrm{at}x=0,L.$$

(2)

In what follows, we will also consider what may be a more realistic assumption that the pressure at the seaward end is zero, following the work in [28], that is,

$$\begin{array}{cc}\hfill & \frac{\partial \Phi }{\partial x}=0,\mathrm{at}x=0,\hfill \\ \hfill & \Phi =0,\mathrm{at}x=L.\hfill \end{array}$$

(3)

The thin plate approximation is used to model the displacement and flexure of the ice shelf:

$$D\frac{{\partial }^4\eta }{\partial x^4}+{\rho }_{i}h\frac{{\partial }^2\eta }{\partial t^2}+{\rho }_{w}g\eta =-{\rho }_w\frac{\partial \Phi }{\partial t},$$

(4)

where g is the acceleration of gravity $g\approx 9.81$ ${\mathrm{ms}}^{-2}$; the length of the ice shelf here is $L=40$ km, the density of the ice is ${\rho }_i=917\mathrm{kg}{\mathrm{m}}^{-3}$, the density of the water is ${\rho }_w=1024$ kg ${\mathrm{m}}^{-3}$, and $D=\frac{E{h}^3}{12(1-{\nu }^2)}$ is the flexural rigidity of the shelf, where $E=11$ GPa is its effective Young’s modulus and $\nu \approx 0.33$ its Poisson’s ratio. The boundary conditions for the vertical displacement $\eta$ are given as

$$\begin{array}{cc}\hfill & \eta =0,\mathrm{and}\frac{\partial \eta }{\partial x}=0,\mathrm{at}x=0,\hfill \\ \hfill & \frac{{\partial }^2\eta }{\partial x^2}=0,\mathrm{and}\frac{{\partial }^3\eta }{\partial x^3}=0,\mathrm{at}x=L.\hfill \end{array}$$

(5)

The conditions when $x=0$ at the landward end represent no deflection and no slope, and when $x=L$ at the seaward end, the bending moment and shear force will vanish. The equations and a schematic are shown in Figure 1.

The variable $\eta$ can be eliminated by taking the derivative with respect to t in Equation (4) and substituting Equation (1) to obtain an equation containing only $\Phi$:

$$DH\frac{{\partial }^6\Phi }{\partial x^6}+{\rho }_{i}hH\frac{{\partial }^4\Phi }{\partial t^2\partial x^2}+{\rho }_{w}gH\frac{{\partial }^2\Phi }{\partial x^2}-{\rho }_w\frac{{\partial }^2\Phi }{\partial t^2}=0.$$

(6)

Non-dimensional variables are now introduced following the work in [19,30]

$$\widehat{x}=\frac{x}{L_c},\widehat{t}=\frac{t}{t_c},$$

(7)

where

$$L_c=\sqrt[4]{\frac{D}{{\rho }_{w}g}},t_c=\sqrt{\frac{{\rho }_w{L}_c^6}{DH}}.$$

(8)

For clarity, the hats will be dropped from now on.

In the nondimensional variables, Equation (6) becomes

$$\frac{{\partial }^6\Phi }{\partial x^6}+M\frac{{\partial }^4\Phi }{\partial t^2\partial x^2}+\frac{{\partial }^2\Phi }{\partial x^2}-\frac{{\partial }^2\Phi }{\partial t^2}=0.$$

(9)

Equation (9) depends on the single nondimensional parameter

$$M=\frac{{\rho }_{i}hH}{{\rho }_w{L}_c^2}.$$

Note that, for most practical cases, $M\ll 1$. The derivation of M is discussed in detail in [19,30]. If we assume M, we have a nondimensional equation that is independent of any parameters, i.e., it would apply to any ice shelf (obviously, under our assumptions).

Our goal is to find the modes of vibration and use these to simulate the motion of the ice shelf. We begin by determining the time-harmonic solution of Equation (9), expressed as follows:

$$\Phi (x,t)=\mathrm{Re}\left\{X\left(x\right)e^{-i\omega t}\right\}.$$

(10)

Therefore, Equation (9) can be reduced to one involving $X\left(x\right)$ only:

$$\frac{{\partial }^{6}X}{\partial x^6}+(1-M{\omega }^2)\frac{{\partial }^{2}X}{\partial x^2}+{\omega }^{2}X=0,$$

(11)

Thus, the general solution has the form of

$$X\left(x\right)=\sum _{m=1}^6{c}_m{e}^{k_{m}x},$$

(12)

where $k_m$ ($m=1,\dots ,6$) are the roots of k of the characteristic polynomial:

$$k^6+(1-M{\omega }^2)k^2+{\omega }^2=0.$$

(13)

We now have a typical eigenvalue problem to solve, in which we have six unknowns and six homogeneous boundary conditions, and we must find the values of $\omega$ for which there are nontrivial solutions. The conditions are as follows:

$$\frac{\partial X}{\partial x}=0,\frac{{\partial }^{2}X}{\partial x^2}=0\mathrm{and}\frac{{\partial }^{3}X}{\partial x^3}=0\mathrm{at}x=0,$$

(14)

$$\frac{{\partial }^{4}X}{\partial x^4}=0\mathrm{and}\frac{{\partial }^{5}X}{\partial x^5}=0\mathrm{at}x=L.$$

(15)

and either

$$\frac{\partial X}{\partial x}=0,x=L,$$

(16)

for the no-flux or

$$X=0,x=L,$$

(17)

for the no-pressure condition. The no-flux condition assumes that the fluid cannot flow out from the seaward end of the ice shelf. The no-pressure condition assumes that the fluid can flow out but that the pressure remains constant.

3. No-Flux Solution

We begin with the no–flux condition. Enforcing the boundary conditions (14), (15), and (16) at the ends of the ice shelf leads to the system of six linear equations $\mathbf{A}\left(\omega \right)\mathbf{c}=0$:

$$\underset{\mathbf{A}}{\underbrace{\left[\begin{array}{ccccccc}{k}_1& k_2& k_3& k_4& k_5& k_6& \\ k_1^2& k_2^2& k_3^2& k_4^2& k_5^2& k_6^2& \\ k_1^3& k_2^3& k_3^3& k_4^3& k_5^3& k_6^3& \\ k_1{e}^{k_{1}L}& k_2{e}^{k_{2}L}& k_3{e}^{k_{3}L}& k_4{e}^{k_{4}L}& k_5{e}^{k_{5}L}& k_6{e}^{k_{6}L}& \\ k_1^4{e}^{k_{1}L}& k_2^4{e}^{k_{2}L}& k_3^4{e}^{k_{3}L}& k_4^4{e}^{k_{4}L}& k_5^4{e}^{k_{5}L}& k_6^4{e}^{k_{6}L}& \\ k_1^5{e}^{k_{1}L}& k_2^5{e}^{k_{2}L}& k_3^5{e}^{k_{3}L}& k_4^5{e}^{k_{4}L}& k_5^5{e}^{k_{5}L}& k_6^5{e}^{k_{6}L}& \end{array}\right]}}\underset{\mathbf{c}}{\underbrace{\left[\begin{array}{c}{c}_1\\ c_2\\ c_3\\ c_4\\ c_5\\ c_6\end{array}\right]}}=0$$

(18)

Here, in order to avoid the trivial solution, we require the value of $\omega$ to be chosen so that $\mathbf{A}\left(\omega \right)$ is singular. We solve numerically using the minimum of the absolute value of the eigenvalues of $\mathbf{A}$, $\min \left|\right(\mathrm{eig}\left(\mathbf{A}\right)\left)\right|$, and searching for the values of $\omega$ that make this function zero. Once we locate the angular frequency ${\omega }_n$ for which $\mathrm{A}\left({\omega }_n\right)$ is singular, we find the coefficients ${\mathbf{c}}_{\mathbf{m}}(m=1,\dots ,6)$ from the eigenvector associated with the zero eigenvalue.

In the numerical examples that follow, we choose the length of the ice shelf to be $L=40$ km, the thickness of the ice is $h=300$ m, and the depth of the cavity is $H=100$ m following the work in [19]. This leads to a value of $M=0.0162$ and a nondimensional length of $\overline{L}=31.0199$. Figure 2 shows the function $min\left|\right(\mathrm{eig}\left(\mathbf{A}\right)\left)\right|$, with the zeros clearly visible.

Table 1. Non-dimensional angular frequencies, ${\omega }_n$, and dimensional periods, $T_n$, of modes for no-flux case, where the dimensional length of the ice shelf is $L=40$ km when $M=0.0162$ and $M=0$.

	$\mathit{M}=0.0162$		$\mathit{M}=0$
n	${\omega }_n$	$T_n$	${\omega }_n$	$T_n$
1	0.10614	0.67697 h	0.10615	0.67691 h
2	0.21260	0.33799 h	0.21268	0.33787 h
3	0.32035	0.22430 h	0.32061	0.22412 h
4	0.43169	0.16645 h	0.43231	0.16622 h
5	0.55072	0.13047 h	0.55196	0.13018 h
10	1.53493	168.530 s	1.54835	167.069 s
15	3.97533	65.0718 s	4.05164	63.8462 s
20	8.67952	29.8037 s	8.97079	28.8360 s
25	16.2423	15.9264 s	17.0818	15.1437 s
30	27.1612	9.52392 s	29.1544	8.87282 s

shows the nondimensional angular frequencies ${\omega }_n$ and the dimensional periods $T_n$. These values match those previously found in [19]. We also investigate the significance of the term M. When $M=0.0162$, the longest period occurs when $T_1=0.67697$ h with ${\omega }_1=0.10614$. Furthermore, the longest period when $M=0$ occurs at $T_1=0.67691$ h with ${\omega }_1=0.10615$. As the mode number increases, the periods for $M=0.0162$ and $M=0$ decrease significantly while the angular frequencies increase. We can only see a slight difference between $M=0.0162$ and $M=0$ for the lowest frequencies. The effect of non-zero M increases with frequency as expected due to the ${\omega }^2$ term in Equation (13).

The associated mode shape for the potential is given from the zero eigenvector $\mathbf{c}$ as

$$X_n\left(x\right)=\sum _{m=1}^6{c}_m{e}^{k_{m}x}.$$

(19)

Figure 3 shows the first five mode shapes of $X_n\left(x\right)$. For the vertical displacement, the time-harmonic solution is given as

$$\eta (x,t)=Re\left\{V\left(x\right)e^{-i\omega t}\right\}.$$

(20)

To find the mode shape in displacement, we use Equation (1) to construct the vertical displacement of the potential velocity. We have

$$\eta =\frac{H}{i\omega }{\Phi }_{xx}.$$

(21)

So, the mode shape equation for the vertical displacement is given as

$$V_n\left(x\right)=\frac{H}{i{\omega }_n}\sum _{m=1}^6{c}_m{k}_m^2{e}^{k_{m}x}.$$

(22)

Figure 3 shows the shapes of the first five modes of $X_n\left(x\right)$ and $V_n\left(x\right)$ in nondimensional units. We believe that these eigenvectors belong to a generalization of Sturm-Liouville theory, and they can be seen to share many properties. In particular, the number of nodes increases by one each time the mode number increases.

4. No-Pressure Solution

Here, we apply the no-pressure condition at the free end, whereas the clamped end still has the no-flux condition, similar to the work in [28]. Note that although this boundary condition was not considered by [18,19], the results of [28] indicate that it more closely matches the open water boundary. Enforcing the boundary conditions (14), (15), and (17) at the ends of the ice shelf leads to the system of six linear equations $\mathbf{A}\left(\omega \right)\mathbf{c}=0$ as follows:

$$\underset{\mathbf{A}}{\underbrace{\left[\begin{array}{ccccccc}{k}_1& k_2& k_3& k_4& k_5& k_6& \\ k_1^2& k_2^2& k_3^2& k_4^2& k_5^2& k_6^2& \\ k_1^3& k_2^3& k_3^3& k_4^3& k_5^3& k_6^3& \\ e^{k_{1}L}& e^{k_{2}L}& e^{k_{3}L}& e^{k_{4}L}& e^{k_{5}L}& e^{k_{6}L}& \\ k_1^4{e}^{k_{1}L}& k_2^4{e}^{k_{2}L}& k_3^4{e}^{k_{3}L}& k_4^4{e}^{k_{4}L}& k_5^4{e}^{k_{5}L}& k_6^4{e}^{k_{6}L}& \\ k_1^5{e}^{k_{1}L}& k_2^5{e}^{k_{2}L}& k_3^5{e}^{k_{3}L}& k_4^5{e}^{k_{4}L}& k_5^5{e}^{k_{5}L}& k_6^5{e}^{k_{6}L}& \end{array}\right]}}\underset{\mathbf{c}}{\underbrace{\left[\begin{array}{c}{c}_1\\ c_2\\ c_3\\ c_4\\ c_5\\ c_6\end{array}\right]}}=0.$$

(23)

As was the case before, we require $\mathbf{A}$ to be singular to avoid the trivial solution. We find the coefficients ${\mathbf{c}}_{\mathbf{m}}(m=1,\dots ,6)$ from the corresponding zero eigenvector. Figure 4 shows, again, how we search for the critical values of $\omega$, as we did previously.

Table 2. As per Table 1, except we consider the no-pressure case.

	$\mathit{M}=0.0162$		$\mathit{M}=0$
n	${\omega }_n$	$T_n$	${\omega }_n$	$T_n$
1	0.05306	1.35423 h	0.05306	1.35420 h
2	0.15930	0.45107 h	0.15933	0.45098 h
3	0.26627	0.26986 h	0.26642	0.26971 h
4	0.37567	0.19127 h	0.37608	0.19107 h
5	0.49088	0.14638 h	0.49177	0.14612 h
10	1.37544	188.071 s	1.38724	186.472 s
15	3.47204	74.5041 s	3.54283	73.0156 s
20	7.71276	33.5394 s	7.98064	32.4136 s
25	14.7178	17.5761 s	15.4951	16.6944 s
30	24.9871	10.3526 s	26.8477	9.63515 s

shows the dimensional periods, $T_n$, associated with nondimensional angular frequencies, ${\omega }_n$, for the same case considered in Table 1, except for the no-pressure case. When comparing Table 1 and Table 2, we see that the periods are significantly lower for the no-pressure case compared to the no-flux case. This is further shown in Figure 5 and Figure 6, which compare the 20 longest periods for the two cases. We can see that the periods and frequencies seem to lie between each other.

Figure 7 shows the first five mode shapes of $X_n\left(x\right)$ and $V_n\left(x\right)$ for the no-pressure case. Note that they are calculated using Equations (19) and (22), as before.

Missing Mode for the No-Flux Case

We can see from Figure 7 that the lowest mode for the no-pressure case has no nodes. In contrast, the lowest mode for the no-flux case (Figure 3) has one node. This indicates that there is a missing mode for the no-flux case. This missing mode is the nontrivial solution for ${\omega }_0=0$, with the corresponding mode $X_0=1$.

5. Simulating Ice Shelf Vibrations

The main aim of the present work is to use the modes to simulate vibration in the time domain. The potential, $\Phi (x,t)$, is expanded in terms of the modes as

$$\Phi (x,t)=\sum _{n=0\mathrm{or}1}^{\infty }{b}_n{X}_n\left(x\right)cos\left({\omega }_{n}t\right)+\sum _{n=1}^{\infty }\frac{d_n}{{\omega }_n}{X}_n\left(x\right)sin\left({\omega }_{n}t\right),$$

(24)

where $b_n$ depends on the initial value of $\Phi$ at time $t=0$, and $d_n$ depends on the initial time derivative of $\Phi$. Note that for the case of no-flux, the sum of the first term begins at zero.

For the vertical displacement, the general solution is given as

$$\eta (x,t)=\sum _{n=1}^{\infty }{G}_n{V}_n\left(x\right)cos\left({\omega }_{n}t\right)+\sum _{n=1}^{\infty }\frac{S_n}{{\omega }_n}{V}_n\left(x\right)sin\left({\omega }_{n}t\right),$$

(25)

where $G_n$ depends on the initial position at time $t=0$, and $S_n$ depends on the initial velocity.

Note that the sum begins at one for the displacement in both cases. This means that the initial value of the displacement is restricted. Equation (1) must be satisfied so that

$${\int }_0^L-\frac{1}{H}\frac{\partial \eta }{\partial t}dx={\int }_0^L\frac{{\partial }^2\Phi }{\partial x^2}dx=(\frac{\partial \Phi }{\partial x}(L))-(\frac{\partial \Phi }{\partial x}(0)).$$

(26)

Thus, in the case of the no-flux condition, the right-hand side of Equation (26) will vanish. Since this equation applies to each mode, it follows that

$${\int }_0^L{V}_n\left(x\right)dx=0.$$

(27)

Note that this is not the case for the no-pressure condition. This issue is linked to the ${\omega }_0=0$ solution for the potential in the no-flux case. The associated mode, $V_0=0$, cannot be used in the calculation. There is a physical explanation for this: the total fluid under the ice shelf must be conserved (since it is incompressible and cannot flow out) under the no-flux condition. Therefore, when choosing the initial value for displacement for the no-flux conditions, $\eta$ must satisfy the condition

$${\int }_0^L\eta (x,0)dx=0.$$

(28)

In the case of no pressure, there is no such restriction. A similar argument shows that the initial condition for the time-derivative of pressure under the ice shelf must also satisfy a similar condition, as does the time-derivative of displacement for the no-flux case. This point was mentioned in [20] in the context of very large floating structures in an enclosed basin.

Find the Coefficient in the Expansion

Consider the case where we want to determine $b_n$ so that Equation (24) satisfies the initial conditions $\Phi (x,0)=f\left(x\right)$ and ${\Phi }_t(x,0)=0$. This implies that

$$\sum _{n=0}^{\infty }{b}_n{X}_n=f\left(x\right).$$

(29)

By multiplying both sides of Equation (29) by $X_m$ and then integrating the result, we obtain

$${\int }_0^{L}f\left(x\right)X_{m}dx=\sum _{n=0}^{\infty }{b}_n{\int }_0^L{X}_n{X}_{m}dx.$$

(30)

We note here that the modes $X_n$ are not orthogonal with respect to the standard inner product. It may be possible to find a function space in which they are orthogonal, but we were not able to do this. For this reason, we use a numerical method. We approximate the integrals by numerical quadrature as

$$\sum _{i=1}^{N}f\left(x_i\right)X_m\left(x_i\right)w_i=\sum _{n=0}^{\infty }{b}_n\sum _{i=1}^N{X}_n\left(x_i\right)X_m\left(x_i\right)w_i,$$

(31)

where $x_i$ are N-evenly spaced points between 0 and L (note that more complex spacing could be used), and ${\mathbf{w}}_{\mathbf{i}}$ is the associated weight function for the quadrature (we use the trapezoidal rule). We then introduce the following matrix:

$$\mathbf{X}=\left[X_m\left(x_n\right)\right]$$

and vectors $\mathbf{f}=\left[f\left(x_n\right)\right]$, $\mathbf{w}=\left[w_n\right]$, and $\mathbf{b}=\left[b_n\right]$. Then Equation (31) becomes

$$\mathbf{X}\mathrm{diag}\left({\mathbf{w}}_{\mathbf{i}}\right){\mathbf{X}}^T\mathbf{b}=\mathbf{X}\mathrm{diag}\left({\mathbf{w}}_{\mathbf{i}}\right)\mathbf{f},$$

where diag means the diagonal matrix formed from the vector and the superscript T is the matrix transpose.

We can then write the solution for potential (Equation (24)) in the same discrete form (that is, at the points $x_i$) as

$$\Phi (x,t)={\mathbf{X}}^T\mathrm{diag}[cos\left({\omega }_{n}t\right)]\mathbf{b},$$

which gives a very efficient method to calculate the time-dependent motion once the modes have been calculated.

The same method works when calculating the solution for a non-zero initial time derivative of the potential or for the displacement.

6. Results

We present simulations of the ice shelf motion in the time domain. These are best viewed in the movies provided in the Supplementary Materials. We begin by simulating the potential velocity, $\Phi$. Figure 8 shows the simulation of the solution for Equation (24) when we apply the initial conditions $\Phi (x,0)=f\left(x\right)$ and ${\Phi }_t(x,0)=0$. We can see the case of no-flux (blue) and no-pressure (cyan), as well as the initial $f\left(x\right)$ function (red). Here, we take a Gaussian function:

$$f\left(x\right)=exp(-{(x-15)}^2/10).$$

(32)

Figure 9 shows the simulation of the solution for Equation (24) when we apply the initial conditions $\Phi (x,0)=0$ and ${\Phi }_t(x,0)=g\left(x\right)$. Where

$$g\left(x\right)=exp(-{(x-17.5)}^2/10)-exp(-{(x-12.5)}^2/10).$$

(33)

It can be seen in Figure 8 and Figure 9 that the solution matches the boundary conditions for the case of no pressure (cyan).

On the other hand, Figure 10 shows the solution for Equation (25) when we apply the initial conditions $\eta (x,0)=f\left(x\right)$ and ${\eta }_t(x,0)=0$. The vertical displacement, $\eta$, for the case of no-flux (blue), the case of no-pressure (cyan), and the initial function $f\left(x\right)$ (red) are shown. Note that $f\left(x\right)$ was carefully chosen here to satisfy Equation (28):

$$f\left(x\right)=exp(-{(x-17.5)}^2/10)-exp(-{(x-12.5)}^2/10).$$

(34)

In addition, Figure 11 shows the solution for the vertical displacement when we apply the initial conditions $\eta (x,0)=0$ and ${\eta }_t(x,0)=g\left(x\right)$. Here, we take the difference of a Gaussian so that the integral of the function vanishes:

$$g\left(x\right)=exp(-{(x-17.5)}^2/10)-exp(-{(x-12.5)}^2/10).$$

(35)

We can see from Figure 10 and Figure 11 that the simulations of the solution match the boundary conditions (5).

The figures of motion and, more particularly, the accompanying animations of the ice shelf simulations show a complexity that appears to be similar to random processes but is actually based on deterministic principles. The sixth-order equations that model this motion introduce a great deal of complexity in the ice shelf dynamics. Moreover, in reality, our model is a simplified version of a real ice shelf, and we can anticipate it to display even more intricate motion. The intricate deterministic nature of the governing equations implies that a single measurement may not provide a thorough understanding of the wider dynamics. This has implications for the interpretation of point measurements of ice shelf motion.

7. Conclusions

We have explored the time-dependent oscillations of an ice shelf. The solution was constructed by expanding the motion in the mode shapes, which are the eigenfunctions of the system, each of which is associated with an eigenfrequency. In the literature, two energy-conserving boundary conditions have been considered: no-flux, which is more prominent, and no-pressure. The two solutions qualitatively give the same motion, although the no-flux condition can only be used to simulate initial conditions with an integral of zero. We presented time–domain simulations of the solutions for both no-flux and no-pressure cases. The corresponding motion is very complicated, and while deterministic, it has somewhat random or unpredictable behavior due to the sixth-order nature of the equations. This complex motion, arising here in a very idealized model, indicates that actual ice shelf motions will be extremely complex and difficult to interpret from single spatial point measurements of the response. Future work will look to extend this model to three dimensions and include the effect of cracks or other sources of non-uniformity.