Estimating a Crack Growth Parameter for the Sulzberger Ice Shelf from Tsunami Induced Calving

1. Introduction

Ice shelves are floating extensions of land-based glaciers and play a pivotal role in modulating ice discharge into the ocean. Their structural integrity is threatened by the formation and propagation of cracks, which can lead to large-scale calving events and contribute to rising sea level [1,2]. Understanding the vibration response of cracked ice shelves to ocean wave forcing is therefore essential to predict their mechanical behaviour and failure modes. Beginning with the work of [3], it has been suggested that the impact of the repetitive flux of ocean waves on the ice shelf could cause the calving of glacier tongues. This hypothesis has received significant evidence in recent years [4,5]. In [4], they show evidence that seasonal loss of sea ice reduces the natural buffer, allowing ocean swells to flex the vulnerable outer ice-shelf margins and weaken them to the point of calving. This scenario is especially troubling given the decline in Antarctic sea ice [6,7,8].

In this study, we focus on two tsunami triggers of ice shelf calving: the Honshu earthquake and the Hunga-Tonga volcanic eruption. As reported in [9], the Honshu earthquake triggered a tsunami that resulted in the calving of the Sulzberger ice shelf $(76.{538}^{\circ }S,150.{168}^{\circ }W$, east of the Ross ice shelf) on 11 March 2011. As a result of this event, two icebergs formed near Antarctica. On the other hand, On January 15, 2022, around 4:15 a.m., the Hunga Tonga Volcano in the Islands of Tonga in the southwest Pacific Ocean triggered a tsunami in the Pacific Ocean. As a result, the West Sulzberger Ice Shelf (WSIS) experienced a significant calving event after the tsunami event. The calving occurred between January 15 and January 18, 2022 [10,11]. In this event, the WSIS lost approximately 106 ${\mathrm{km}}^2$ of ice, reducing its size to the smallest recorded since 1948 [12]. We chose these events because we can reconstruct to some accuracy the wave forcing and subsequent strain induced on the ice shelf. This, in turn, allows us to estimate the parameters that may control the crack growth in ice shelves. In particular, the evidence from the tsunami induced calving seems to suggest that the slow ringing of the tsunami wave induced crack growth through cyclic loading, implying that this process maybe governed by a Paris type law.

In this paper and building on previous studies of floating elastic beams [13,14,15,16], we model the cracked ice shelf as a thin beam of constant thickness floating in shallow water [17], where the crack is represented as a rotational spring whose stiffness depends on the depth ratio of the crack, following the formulation [16,18,19,20]. The cracked ice shelf domain is divided into three regions, and analytical solutions are constructed in each region using modal decomposition. Continuity conditions at the interfaces and a rotational compatibility condition at the crack yield a system of ten equations for the modal amplitudes. Later, we investigate the dynamic response of a cracked ice shelf subjected to tsunami wave-induced vibrations, with a particular focus on modelling crack growth using Paris’ Law [21,22]. For ice, the first attempt to apply the Paris law was reported in [23], where the authors presented preliminary results on fatigue crack growth on freshwater ice. They found evidence of crack growth that occurred under cyclic loading, although the data was insufficient to determine the validity of the behaviour of ${\displaystyle \frac{da}{dN}}$ and $\Delta K$. Later, in [24], they present preliminary results from fatigue tests in Antarctica. In contrast, several studies in engineering have used Paris’ law to determine the growth of fatigue cracks in various materials [25,26,27]. In [28], they established a procedure for analysing the crack propagation in a ship by using the Paris’ law in combination with the Monte Carlo technique. Also in [29], they study the effect of low temperature on fatigue crack propagation rates in a polar ship using the Paris law in their fatigue crack propagation test. Moreover, we provide a time domain solution and calculate the stress, which we used to compute the stress intensity factor, and fatigue crack growth in modelling Paris’ Law [30]. This approach provides a physics-based framework for simulating progressive fracture in ice shelves and assessing the impact of a tsunami wave on their stability.

The outline of the paper is as follows. In section 2 we investigate in detail the tsunami wave forces as measured near the Sulzberger ice shelf. We show that the waves for both tsunamis were similar. In section 3 we present our simple model for a cracked ice shelf and investigate the response of the ice shelf to the measured tsunami waves in the time domain in section 4. In section 5 we consider the application of Paris law and we give an estimate of the Paris law coefficient which controls the fatigue crack growth rate for an ice shelf, on the assumption the Paris growth law is appropriate.

2. Observational Tsunami Data

We present here measurements of the tsunami wave forcing made by tide gauges. We note that it was not possible to obtain the measurements exact on the Sulzberger ice shelf. For the Honshu Earthquake, the tsunami record used in this study was obtained from the ROBT tide-gauge station (10–20 March 2011) (see Figure 1 (a) ). For the Hunga-Tonga volcanic eruption, we use two tide gauge station datasets showing sea level changes during the volcanic event: the Antarctica Base station $(62.{48}^{\circ }S,59.{66}^{\circ }W)$, located on the northwest of the Antarctic Peninsula, and the Cape Robert station $(77.{03}^{\circ }S,163.{19}^{\circ }E)$, situated near the Ross Ice Shelf [12]. The tsunami wave amplitude for both stations is shown in Figure 1 (b) from 15 to 20 January 2022. We note that the iceberg experiences the tsunami waves as a slowly decaying wave which last for one to two hundred hours. Both tsunami waves are remarkably similar. We note the amplitude is around 10 cm in both cases.

2.1. Methods and Results

Building on the previous observations, we applied a high-pass filter on the tide-gauge record with a cutoff period of 2 hours to smooth the wave amplitude, remove low-frequency tidal and background ocean motion, and isolate the tsunami signal (see Figure 2).

The filtered record was then transformed into the frequency domain using a Fast Fourier Transform (FFT) to convert the signal in the time-domain to its representation in the frequency domain, as in Figure 3 and Figure 4.

The tsunami spectrum exhibits several distinct energy peaks rather than a single dominant frequency. Therefore, to represent the tsunami forcing in our model, a multi-peak Gaussian fitting approach was applied to the filtered tidal-gauge spectrum. Each of the Gaussian functions was given by

$$A\left(f\right)=A_{0}exp\left(-{\displaystyle \frac{{(f-f_0)}^2}{2{\sigma }^2}}\right)$$

(1)

These individual Gaussian components were then superimposed to form a composite spectrum that closely follows the shape and amplitude of the measured tsunami frequency content. Figure 4 shows the fitted Gaussian function for the tsunami amplitude for each tide gauge in both tsunami cases.

This method enables us to use a synthetic spectrum to simulate the broadband characteristics of the real tsunami signal, rather than representing it as a mono-frequency wave.

In our model, we are dealing with a long-wave tsunami and can model the ice shelf as floating on shallow water leading to the following dispersion relation:

$$\omega =\sqrt{gH}k$$

(2)

where $\omega =2\pi f$ is the angular frequency, g is the acceleration of gravity, and H is the depth of the water. Therefore, the wavenumber is equal to

$$k={\displaystyle \frac{\omega }{\sqrt{gH}}}$$

(3)

Now, transforming from the frequency domain to the wavenumber domain, the Gaussian in wavenumber is:

$$A\left(k\right)=A_{0}exp\left(-{\displaystyle \frac{{(k-k_0)}^2}{2{\sigma }_k^2}}\right).$$

(4)

To examine the distribution of tsunami energy across multiple frequency decades, we used Log-log spectral plots. As before, we fitted multiple Gaussian functions to approximate the tsunami frequency spectrum. In logarithmic scaling, the visibility of both dominant low-frequency components and weaker higher-frequency oscillations becomes clearer. In Figure 5, the majority of the tsunami energy is concentrated within the ${10}^{-4}-{10}^{-3}$ Hz band, confirming the long-period nature of the tsunami forcing relevant to ice-shelf vibration.

In Figure 6, we made a compression between the tide-gauge record with the high-pass filter and the inverse Fourier transform of the fitted spectra, to see the precision of the fitted Gaussian tsunami spectrum.

We conclude from our analysis of these wave data for the tsunamis that the forcing in both cases was roughly similar with an amplitude of around 10 cm and with distinct spectral peaks. We note that the origin of these spectral peaks is not clear. We will consider the stress at the crack due to each of these frequencies in turn. To keep the analysis simple we set the peak response to 10 cm.

3. Vibration of a Cracked Ice Shelf Due to Transient Wave Forcing

Our model aims to simulate the motion of a cracked ice shelf to tsunami wave forcing. We begin by simulating the effect of a crack. The methodological analysis to simulate a cracked ice shelf is based on [13] and [16]. For completeness we give the theory here.

To determine the flexure of an ice shelf with a through-thickness crack, we model the shelf as a thin elastic beam floating on the ocean. The crack is assumed to be located at a fixed position $x=L$, and it modifies the boundary conditions of the system (see Figure 7). The displacement of the ice shelf is denoted $u(x,t)$, and the wave forcing is described by the velocity potential $\varphi (x,t)$.

The governing thin-beam approximation leads to

$$D{\displaystyle \frac{{\partial }^{4}u}{\partial x^4}}+{\rho }_{i}h{\displaystyle \frac{{\partial }^{2}u}{\partial t^2}}+{\rho }_{w}gu=-{\rho }_w{\displaystyle \frac{\partial \varphi }{\partial t}},$$

(5)

where $D=\frac{E{h}^3}{12(1-{\nu }^2)}$ is the flexural rigidity of the shelf. Here, $E=11$ GPa is the effective Young’s modulus, $\nu =0.33$ is Poisson’s ratio, ${\rho }_i=922.5\mathrm{kg}{\mathrm{m}}^{-3}$ is the ice density, ${\rho }_w=1024\mathrm{kg}{\mathrm{m}}^{-3}$ is the water density, and $g=9.81\mathrm{m}{\mathrm{s}}^{-2}$ is gravitational acceleration.

By assumption, the water potential $\varphi$ and the vertical displacement u are coupled through the shallow-water equation [17]:

$${\displaystyle \frac{{\partial }^2\varphi }{\partial x^2}}=-{\displaystyle \frac{1}{H}}{\displaystyle \frac{\partial u}{\partial t}},$$

(6)

where H is the (constant) water depth [31]. Figure 7 provides a schematic diagram of the cracked ice shelf.

To proceed, we follow [13] and assume harmonic time dependence. That is, we write.

$$u(x,t)=u\left(x\right)e^{-i\omega t},\varphi (x,t)=\varphi \left(x\right)e^{-i\omega t},$$

where $\omega$ is the angular frequency. This frequency-domain representation is not simply a “simplification,” but a standard Fourier-based approach: it enables us to solve for monochromatic wave responses, from which more general wave packets can be reconstructed by superposition. Substituting these forms yields

$$D{\displaystyle \frac{d^{4}u}{d{x}^4}}-{\omega }^2{\rho }_{i}hu+{\rho }_{w}gu=i\omega {\rho }_w\varphi ,$$

(7)

and

$${\displaystyle \frac{d^2\varphi }{d{x}^2}}={\displaystyle \frac{i\omega }{H}}u.$$

(8)

From these two coupled equations, we derive a sixth-order equation for $\varphi \left(x\right)$:

$$D{\displaystyle \frac{d^6\varphi }{d{x}^6}}-{\omega }^2{\rho }_{i}h{\displaystyle \frac{d^2\varphi }{d{x}^2}}+{\rho }_{w}g{\displaystyle \frac{d^2\varphi }{d{x}^2}}=-{\displaystyle \frac{{\omega }^2{\rho }_w}{H}}\varphi .$$

(9)

The open-water limit (no ice above) reduces to

$${\rho }_{w}g{\displaystyle \frac{d^2\varphi }{d{x}^2}}=-{\displaystyle \frac{{\omega }^2{\rho }_w}{H}}\varphi .$$

(10)

To solve the differential equation for the whole domain, we divide the system into three regions with different properties: (i) the region landward of the crack, (ii) the region seaward of the crack, and (iii) the open water. The governing equations in each region can be summarized as

$$-{\omega }^2\varphi =\left\{\begin{array}{cc}{\displaystyle \frac{H_2}{{\rho }_w}}\left(D{\displaystyle \frac{d^6\varphi }{d{x}^6}}+({\rho }_{w}g-{\omega }^2{\rho }_{i}h){\displaystyle \frac{d^2\varphi }{d{x}^2}}\right),\hfill & x>L,\hfill \\ {\displaystyle \frac{H_1}{{\rho }_w}}\left(D{\displaystyle \frac{d^6\varphi }{d{x}^6}}+({\rho }_{w}g-{\omega }^2{\rho }_{i}h){\displaystyle \frac{d^2\varphi }{d{x}^2}}\right),\hfill & 0<x<L,\hfill \\ Hg{\displaystyle \frac{d^2\varphi }{d{x}^2}},\hfill & x<0,\hfill \end{array}\right.$$

(11)

where L denotes the crack position.

This framework can then be extended to more complicated configurations, such as multiple cracks or spatially varying elastic properties. Here, for clarity, We assume constant ice thickness, Young’s modulus, and density. Thus, the solution can be written as

$$\varphi (x,\omega )=\left\{\begin{array}{cc}{\displaystyle {\sum }_{n=1}^3}{\alpha }_n{e}^{-r_n(x-L)},\hfill & x>L.\hfill \\ {\displaystyle {\sum }_{n=1}^3}{\beta }_n{e}^{s_n(x-L)}+{\displaystyle {\sum }_{n=1}^3}{\gamma }_n{e}^{-s_{n}x},\hfill & 0<x<L.\hfill \\ e^{-ikx}+R_1{e}^{ikx},\hfill & x<0.\hfill \end{array}\right.$$

(12)

In the open water, the wave-numbers are

$$k^2={\displaystyle \frac{{\omega }^2}{gH}}.$$

(13)

and the roots $r_n$ and $s_n$ can be found by the following equations

$${\displaystyle \frac{D}{{\rho }_{w}g}}{r}^6+(1-{\omega }^2{\displaystyle \frac{{\rho }_{i}h}{{\rho }_{w}g}})r^2+{\displaystyle \frac{{\omega }^2}{Hg}}=0.$$

(14)

$${\displaystyle \frac{D}{{\rho }_{w}g}}{s}^6+(1-{\omega }^2{\displaystyle \frac{{\rho }_{i}h}{{\rho }_{w}g}})s^2+{\displaystyle \frac{{\omega }^2}{Hg}}=0.$$

(15)

Here, too, the roots come in pairs. As we mentioned in the previous case, we will choose the roots with positive sign only for both real and imaginary parts for the roots $r_n$ and $s_n$. However, if a root is a complex number, we choose the positive real part.

In the case of the cracked ice shelf, the system is governed by ten boundary and continuity conditions. Four of these are imposed at the free edge ($x=0$), and are expressed as follows:

$$\varphi \left(x^{+}\right)=\varphi \left(x^{-}\right),\mathrm{at}x=0,$$

(16)

$${\varphi }^{\prime }\left(x^{+}\right)={\varphi }^{\prime }\left(x^{-}\right),\mathrm{at}x=0,$$

(17)

$${\varphi }^{\prime \prime \prime \prime }\left(x^{+}\right)=0,\mathrm{at}x=0,$$

(18)

$${\varphi }^{\prime \prime \prime \prime \prime }\left(x^{+}\right)=0,\mathrm{at}x=0.$$

(19)

The remaining six conditions are associated with the crack and the edge at $x=L$, consisting of two continuity conditions at the plate end and four at the cracked section, as described in [19]:

$$\varphi \left(x^{+}\right)=\varphi \left(x^{-}\right),\mathrm{at}x=L,$$

(20)

$${\varphi }^{\prime }\left(x^{+}\right)={\varphi }^{\prime }\left(x^{-}\right),\mathrm{at}x=L,$$

(21)

$${\varphi }^{\prime \prime }\left(x^{+}\right)={\varphi }^{\prime \prime }\left(x^{-}\right),\mathrm{at}x=L,$$

(22)

$${\varphi }^{\prime \prime \prime \prime }\left(x^{+}\right)={\varphi }^{\prime \prime \prime \prime }\left(x^{-}\right),\mathrm{at}x=L,$$

(23)

$${\varphi }^{\prime \prime \prime \prime \prime }\left(x^{+}\right)={\varphi }^{\prime \prime \prime \prime \prime }\left(x^{-}\right),\mathrm{at}x=L.$$

(24)

To account for the influence of the crack on the ice shelf’s slope continuity, an additional compatibility condition is introduced:

$${\varphi }^{\prime \prime \prime }\left(x^{+}\right)+\mu {\varphi }^{\prime \prime \prime \prime }\left(x^{+}\right)={\varphi }^{\prime \prime \prime }\left(x^{-}\right),\mathrm{at}x=L,$$

(25)

where $\mu =5.346hf(a/h)$ represents the rotational spring stiffness modelling the crack. Here, h denotes the ice thickness, and $f(a/h)$ is a compliance function dependent on the crack depth ratio, given by:

$$f(a/h)=1.8624{\left({\displaystyle \frac{a}{h}}\right)}^2-3.95{\left({\displaystyle \frac{a}{h}}\right)}^3+16.375{\left({\displaystyle \frac{a}{h}}\right)}^4-37.226{\left({\displaystyle \frac{a}{h}}\right)}^5+$$

$$76.81{\left({\displaystyle \frac{a}{h}}\right)}^6-126.9{\left({\displaystyle \frac{a}{h}}\right)}^7+172{\left({\displaystyle \frac{a}{h}}\right)}^8-143.97{\left({\displaystyle \frac{a}{h}}\right)}^9+66.56{\left({\displaystyle \frac{a}{h}}\right)}^{10},$$

as reported in [16,19,20,32].

These ten conditions yield a system of linear equations of the form $\mathbf{Mc}=\mathbf{f}$, where $\mathbf{M}$ is the coefficient matrix representing the boundary conditions, and $\mathbf{c}$ is the vector of unknown coefficients.

$$\left[\begin{array}{cccccccccc}0& 0& 0& e^{-s_{1}L}& e^{-s_{2}L}& e^{-s_{3}L}& 1& 1& 1& -1\\ 0& 0& 0& s_1{e}^{-s_{1}L}& s_2{e}^{-s_{2}L}& s_3{e}^{-s_{3}L}& -s_1& -s_2& -s_3& ik\\ 0& 0& 0& s_1^4{e}^{-s_{1}L}& s_2^4{e}^{-s_{2}L}& s_3^4{e}^{-s_{3}L}& s_1^4& s_2^4& s_3^4& 0\\ 0& 0& 0& s_1^5{e}^{s_{1}L}& s_2^5{e}^{s_{2}L}& s_3^5{e}^{s_{3}L}& -s_1^5& -s_2^5& -s_3^5& 0\\ 1& 1& 1& -1& -1& -1& -e^{-s_{1}L}& -e^{-s_{2}L}& -e^{-s_{3}L}& 0\\ -r_1& -r_2& -r_3& -s_1& -s_2& -s_3& s_1{e}^{-s_{1}L}& s_2{e}^{-s_{2}L}& s_3{e}^{-s_{3}L}& 0\\ r_1^2& r_2^2& r_3^2& -s_1^2& -s_2^2& -s_3^2& -s_1^2{e}^{-s_{1}L}& -s_2^2{e}^{-s_{2}L}& -s_3^2{e}^{-s_{3}L}& 0\\ r_1^4& r_2^4& r_3^4& -s_1^4& -s_2^4& -s_3^4& -s_1^4{e}^{-s_{1}L}& -s_2^4{e}^{-s_{2}L}& -s_3^4{e}^{-s_{3}L}& 0\\ -r_1^5& -r_2^5& -r_3^5& -s_1^5& -s_2^5& -s_3^5& s_1^5{e}^{-s_{1}L}& s_2^5{e}^{-s_{2}L}& s_3^5{e}^{-s_{3}L}& 0\\ -r_1^3+\mu r_1^4& -r_2^3+\mu r_2^4& -r_3^3+\mu r_3^4& -s_1^3& -s_2^3& -s_3^3& s_1^3{e}^{-s_{1}L}& s_2^3{e}^{-s_{2}L}& s_3^3{e}^{-s_{3}L}& 0\end{array}\right]\times \left[\begin{array}{c}{\alpha }_1\\ {\alpha }_2\\ {\alpha }_3\\ {\beta }_1\\ {\beta }_2\\ {\beta }_3\\ {\gamma }_1\\ {\gamma }_2\\ {\gamma }_3\\ R_1\end{array}\right]=\left[\begin{array}{c}1\\ ik\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right]$$

Moreover, the equations for the displacement u are

$$-i\omega u(x,\omega )=\left\{\begin{array}{cc}{\displaystyle {\sum }_{n=1}^3}{\alpha }_n{r}_n^2{e}^{-r_n(x-L)},\hfill & x>L.\hfill \\ {\displaystyle {\sum }_{n=1}^3}{\beta }_n{s}_n^2{e}^{s_n(x-L)}+{\displaystyle {\sum }_{n=1}^3}{\gamma }_n{s}_n^2{e}^{-s_{n}x},\hfill & 0<x<L.\hfill \\ i^2{k}^2{e}^{-ikx}+R_1{i}^2{k}^2{e}^{ikx},\hfill & x<0.\hfill \end{array}\right.$$

(26)

Here $L=8000$. Figure 8 shows the displacement u for the cracked ice.

4. Time-Domain Solutions

To simulate the motion in the time domain, we will use the given formula

$$u(x,t)=\Re \left\{{\int }_0^{\infty }\widehat{f}\left(k\right)\varphi (\omega \left(k\right),x)e^{i\omega \left(k\right)t}dk.\right\},$$

(27)

where $k={\displaystyle \frac{\omega }{\sqrt{gH}}}$ is the wave number in the left-hand incident region, and

$$\widehat{f}\left(k\right)=A_{0}exp\left(-{\displaystyle \frac{{(k-k_0)}^2}{2{\sigma }_k^2}}\right).$$

(28)

We take $A_0$ to be 10 cm.

After the simulations, we found that long-period tsunami waves were smoothly transmitted through the Sulzberger Ice Shelf without reflection or noticeable bending at the crack (see Figure 9 and Figure 10). The reason is that the tsunami has low-frequency forcing around $0.1\times {10}^{-3}$ up to $1.8\times {10}^{-3}$ and wavelengths that are extremely large relative to the characteristic dimensions of the ice shelf and the geometry of the crack. Consequently, the incident wave energy is transmitted efficiently through the ice shelf. To further investigate the influence of wave frequency on hydroelastic behaviour, additional simulations were performed using progressively higher forcing frequencies, up to $f=20\times {10}^{-3}\mathrm{Hz}$ (see Figure 11). Increasing the forcing frequency reduces the wavelength of the incident wave, enhancing its interaction with the ice-shelf geometry and crack discontinuity. Under these higher-frequency conditions, the reflection becomes visible, with wave energy reflecting toward the open ocean after interacting with the ice shelf.

These results indicate that tsunami-frequency waves primarily behave as long, shallow-water waves that propagate through the ice shelf with minimal reflection. In contrast, shorter-wavelength forcing produces stronger hydroelastic interactions, localised oscillations, and increased reflected wave energy.

4.1. The Crack Stress in the Time Domain

The bending stress at the crack location in the ice shelf is obtained from the second spatial derivative of the displacement field. The expression is given by:

$$\sigma \left(t\right)=-{\displaystyle \frac{EZ}{1-{\nu }^2}}{u}_{xx}\left(t\right).$$

(29)

where $Z=h/2$ is the distance from the neutral axis to the upper or lower ice surface, and $u_{xx}\left(t\right)$ is the second derivative of the displacement with respect to x at the location of the crack.

We start by taking the second derivative of $u(x,\omega )$ in the frequency domain. The governing expressions for $u^{\prime \prime }(x,\omega )$ are derived from the modal decomposition solution and are given by:

$$-i\omega u^{\prime \prime }(x,\omega )=\left\{\begin{array}{cc}{\displaystyle \sum _{n=1}^3}{\alpha }_n{r}_n^4{e}^{-r_n(x-L)},\hfill & x>L,\hfill \\ {\displaystyle \sum _{n=1}^3}{\beta }_n{s}_n^4{e}^{s_n(x-L)}+{\displaystyle \sum _{n=1}^3}{\gamma }_n{s}_n^4{e}^{-s_{n}x},\hfill & 0<x<L,\hfill \\ i^4{k}^4{e}^{-ikx}+R{i}^4{k}^4{e}^{ikx},\hfill & x<0,\hfill \end{array}\right.$$

(30)

where ${\alpha }_n,{\beta }_n,{\gamma }_n$ are modal coefficients determined from boundary and continuity conditions, $r_n$ and $s_n$ are characteristic roots of the system, and R is the reflection coefficient.

The stress range is then defined by:

$$\Delta \sigma ={\sigma }_{max}-{\sigma }_{min}$$

(31)

Figure 12. Stress distributions along the ice shelf corresponding to four dominant tsunami frequencies.

Figure 12. Stress distributions along the ice shelf corresponding to four dominant tsunami frequencies.

5. Paris’ Law and Crack Growth Phases in Ice Shelves

Predicting crack growth in ice shelves is essential for understanding their stability and the potential impacts on global sea levels. Here, we have used Paris’ law to estimate the fatigue crack growth induced by the tsunami wave and estimate the number of cycles until the ice shelf eventually breaks [33]. The law is given by

$${\displaystyle \frac{da}{dN}}=C{(\Delta K)}^m$$

(32)

where C and m are material fatigue parameters , N is the number of loading cycles, a is the crack depth, and $\Delta K$ is the stress intensity factor range, given by:

$$\Delta K=F(a/h)\Delta \sigma \sqrt{\pi a}$$

(33)

where h is the ice thickness, $\Delta \sigma$ is the tsunami-induced stress range, and $F(a/h)$ is a dimensionless geometry correction factor that accounts for the finite thickness of the ice shelf and the relative crack depth. In this study, the crack was treated as a single-edge crack propagating through the ice thickness, so the correction factor was expressed as a function of the crack-depth ratio $a/h$ [18,21], therefore

$$\Delta K=\Delta \sigma \sqrt{\pi a}\left[{\displaystyle \frac{1.99-\left(\frac{a}{h}\right)\left(1-\frac{a}{h}\right)\left[2.15-3.93\left(\frac{a}{h}\right)+2.70{\left(\frac{a}{h}\right)}^2\right]}{\left(1+2\frac{a}{h}\right){\left(1-\frac{a}{h}\right)}^{3/2}}}\right].$$

(34)

Critical to our work is trying to obtain an estimate of C, which is highly dependent on the material in question.

Table 1. Paris’ law results for the analyzed tsunami frequencies in the Figure 4(a):

f (Hz)	$\mathbf{\Delta }\mathit{\sigma }$(MPa)	$\mathbf{\Delta }\mathit{K}$ (MPa$\sqrt{\mathbf{m}}$)
$0.21\times {10}^{-3}$	2.1007e-04	2.1432e-02
$0.32\times {10}^{-3}$	4.8134e-04	4.9110e-02
$0.50\times {10}^{-3}$	1.1818e-03	1.2058e-01
$0.66\times {10}^{-3}$	1.9883e-03	2.0286e-01
$0.80\times {10}^{-3}$	3.0268e-03	3.0883e-01
$0.98\times {10}^{-3}$	4.5362e-03	4.6289e-01
$1.10\times {10}^{-3}$	5.7200e-03	5.8378e-01
$1.18\times {10}^{-3}$	6.5545e-03	6.6905e-01

Table 1. Paris’ law results for the analyzed tsunami frequencies in the Figure 4(a):

f (Hz)	$\mathbf{\Delta }\mathit{\sigma }$(MPa)	$\mathbf{\Delta }\mathit{K}$ (MPa$\sqrt{\mathbf{m}}$)
$0.21\times {10}^{-3}$	2.1007e-04	2.1432e-02
$0.32\times {10}^{-3}$	4.8134e-04	4.9110e-02
$0.50\times {10}^{-3}$	1.1818e-03	1.2058e-01
$0.66\times {10}^{-3}$	1.9883e-03	2.0286e-01
$0.80\times {10}^{-3}$	3.0268e-03	3.0883e-01
$0.98\times {10}^{-3}$	4.5362e-03	4.6289e-01
$1.10\times {10}^{-3}$	5.7200e-03	5.8378e-01
$1.18\times {10}^{-3}$	6.5545e-03	6.6905e-01

There is no clear evidence of Paris’ law constants for a large-scale ice shelf in any previous studies. For m we adopted a value that is commonly used in fatigue growth, which is $m=3$. For the coefficient C we estimate it using the observed calving timescale of the Sulzberger Ice Shelf following the 2011 and 2022 tsunamis which in both cases was around 3 days.

Our estimate is based on the observation that there was significant energy up to around $1.18\times {10}^{-3}\mathrm{Hz}$, corresponding to a period of $T=847.46\mathrm{s}$. Since the calving event occurred approximately three days after the tsunami [9], the corresponding number of loading cycles is 306. By assuming that the crack propagates from the initial crack depth $a_0=200,\mathrm{m}$ to the critical crack depth $a_c=300,\mathrm{m},$ the required average crack growth per cycle is $3.27\times {10}^{-1}\mathrm{m}/\mathrm{cycle}$. This gives a value of $C=2\times {10}^{-1}{(\mathrm{MPa}/\sqrt{\mathrm{m}})}^{-3}$ As a final remark, we want to make clear that we cannot be sure that Paris law applies to crack growth in ice shelves, and if it does, this can only be seen as giving an order of magnitude estimate.

6. Discussion

We have shown that the tsunami wave forcing, which induced calving events in the Sulzberger ice shelf in 2011 and 2022, was of similar magnitude and spectral range. We simulate this ice shelf with these wave forces and a crack using a simple one dimensional model. We used our results to estimate the value of the critical fatigue growth parameter, C, in the Paris law. We see this as a preliminary study showing how these rare geophysical phenomena can be used to give estimates of critical parameters. The next step is to consider the effect of ocean wave swell on crack growth using our estimated value of C and more sophisticated crack growth and ice shelf models.