2. Field equations

The 3-D elastic model is based on the well-known momentum equations (e.g. Landau and Lifshitz, Reference Landau and Lifshitz1986; Lurie, Reference Lurie2005):

(1)\begin{equation}\left\{ \begin{gathered} \frac{{\partial {\sigma _{xx}}}}{{\partial {x}}} + \frac{{\partial {\sigma _{xy}}}}{{\partial {y}}} + \frac{{\partial {\sigma _{xz}}}}{{\partial {z}}} = \rho \frac{{{\partial ^2}U}}{{\partial {t^2}}}; \hfill \\ \frac{{\partial {\sigma _{yx}}}}{{\partial {x}}} + \frac{{\partial {\sigma _{yy}}}}{{\partial {y}}} + \frac{{\partial {\sigma _{yz}}}}{{\partial {z}}} = \rho \frac{{{\partial ^2}V}}{{\partial {t^2}}}; \hfill \\ \frac{{\partial {\sigma _{zx}}}}{{\partial {x}}} + \frac{{\partial {\sigma _{zy}}}}{{\partial {y}}} + \frac{{\partial {\sigma _{zz}}}}{{\partial {z}}} = \rho \frac{{{\partial ^2}W}}{{\partial {t^2}}} + \rho g; \hfill \\ 0 \lt x \lt L;\,{y_1}(x) \lt y \lt {y_2}(x);{h_b}(x,y) \lt z \lt {h_s}(x,y), \hfill \\ \end{gathered} \right.\end{equation}

where $\left( {XYZ} \right)$ is a rectangular coordinate system with x-axis along the central line (in the direction of wave propagation), and z-axis is pointing vertically upward; $U,{\text{ }}V$ and $W$ are two horizontal and one vertical ice displacements, respectively; $\sigma $ is the stress tensor; and $\rho $ is ice density. The ice shelf is of length L along the central line. The geometry of the ice shelf is assumed to be given by lateral boundary functions ${y_{1,2}}\left( x \right)$ at sides labeled 1 and 2 and functions for the surface and base elevation, ${h_{{\text{s}},{\text{b}}}}\left( {x,y} \right)$, denoted by subscripts s and b, respectively. Thus, the domain on which Eqns (1) are solved is $\Omega=\left\{0 \lt x \lt L,\text{ }y_1\left(x\right) \lt y \lt y_2\left(x\right),\text{ }h_\text{b}\left(x,y\right) \lt z \lt h_\text{s}\left(x,y\right)\right\}$.

Equation (1) can be rewritten in integrodifferential form. This integrodifferential form results from the vertical integration of the momentum Eqn (1) from the current vertical coordinate to the ice surface. In particular, considering the first equation from Eqn (1) and integrating over $z'$ from the current z-coordinate to ${h_{\text{s}}}$ we obtain the following equation:

(2)\begin{equation}\mathop \int \limits_z^{{h_s}} \frac{{\partial {\sigma _{xx}}}}{{\partial x}}{\text{ d}}z' + {\text{ }}\mathop \int \limits_z^{{h_s}} \frac{{\partial {\sigma _{xy}}}}{{\partial y}}{\text{ d}}z' + {\sigma _{xz}}{|_{z = {h_s}}} - {\text{ }}{\sigma _{xz}} = \rho \mathop \int \limits_z^{{h_s}} \frac{{{\partial ^2}U}}{{\partial {t^2}}}\text{d}z'\end{equation}

Next, using the Leibniz integral rule, we replace the first and second terms in Eqn (2) with the following, respectively.

(3.1)\begin{equation}\mathop \int \limits_z^{{h_s}} \frac{{\partial {\sigma _{xx}}}}{{\partial x}}{\text{ d}}z' = \frac{\partial }{{\partial x}}{\text{ }}\mathop \int \limits_z^{{h_s}} {\sigma _{xx}}{\text{d}}z' - {\text{ }}{\left( {{\sigma _{xx}}} \right)_{z = {h_s}}}\frac{{\partial {h_s}}}{{\partial x}};\end{equation}

(3.2)\begin{equation}\mathop \int \limits_z^{{h_s}} \frac{{\partial {\sigma _{xy}}}}{{\partial y}}{\text{ d}}z' = {\text{ }}\frac{\partial }{{\partial y}}{\text{ }}\mathop \int \limits_z^{{h_s}} {\sigma _{xy}}{\text{ d}}z' - {\text{ }}{\left( {{\sigma _{xy}}} \right)_{z = {h_s}}}\frac{{\partial {h_s}}}{{\partial y}}.\end{equation}

Thus, instead of Eqn (2), we obtain the following equation:

(4)\begin{align}&\frac{\partial }{{\partial x}}{\text{ }}\mathop \int \limits_z^{{h_s}} {\sigma _{xx}}{\text{ }}dz' + {\text{ }}\frac{\partial }{{\partial y}}{\text{ }}\mathop \int \limits_z^{{h_s}} {\sigma _{xy}}{\text{ }}dz' - {\text{ }}{\left( {{\sigma _{xx}}} \right)_{z = {h_s}}}\frac{{\partial {h_s}}}{{\partial x}} \nonumber\\ &\quad- {\text{ }}{\left( {{\sigma _{xy}}} \right)_{z = {h_s}}}{\text{ }}\frac{{\partial {h_s}}}{{\partial y}} + {\text{ }}{\left( {{\sigma _{xz}}} \right)_{z = {h_s}}} - {\sigma _{xz}} = \rho \mathop \int \limits_z^{{h_s}} \frac{{{\partial ^2}U}}{{\partial {t^2}}}dz'.\end{align}

Taking in to account that the expression $\frac{{ - {\text{ }}{{\left( {{\sigma _{xx}}} \right)}_{z = {h_s}}}\frac{{\partial {h_s}}}{{\partial x}} - {{\left( {{\sigma _{xy}}} \right)}_{z = {h_s}}}{\text{ }}\frac{{\partial {h_s}}}{{\partial y}} + {{\left( {{\sigma _{xz}}} \right)}_{z = {h_s}}}}}{{\sqrt {{{\left( {\frac{{\partial {h_s}}}{{\partial x}}} \right)}^2} + {{\left( {\frac{{\partial {h_s}}}{{\partial y}}} \right)}^2} + 1} }}$ is the x-component of the force acting on a unit square of the ice surface (e.g. Landau and Lifshitz, Reference Landau and Lifshitz1986) and it is equal to zero accordingly the boundary conditions (ice surface is stress free), we finally obtain the following equation:

(5)\begin{equation}\frac{\partial }{{\partial x}}{\text{ }}\mathop \int \limits_z^{{h_s}} {\sigma _{xx}}{\text{ }}dz' + {\text{ }}\frac{\partial }{{\partial y}}{\text{ }}\mathop \int \limits_z^{{h_s}} {\sigma _{xy}}{\text{ }}dz' - {\sigma _{xz}} = \rho \mathop \int \limits_z^{{h_s}} \frac{{{\partial ^2}U}}{{\partial {t^2}}}dz'.\end{equation}

By performing similar manipulations with the second and third equations from Eqn (1), we obtain the following equations:

(6)\begin{equation}\frac{\partial }{{\partial x}}{\text{ }}\mathop \int \limits_z^{{h_s}} {\sigma _{yx}}{\text{ }}dz' + {\text{ }}\frac{\partial }{{\partial y}}{\text{ }}\mathop \int \limits_z^{{h_s}} {\sigma _{yy}}{\text{ }}dz' - {\sigma _{yz}} = \rho \mathop \int \limits_z^{{h_s}} \frac{{{\partial ^2}V}}{{\partial {t^2}}}dz'\end{equation}

and

(7)\begin{align}&\frac{\partial }{{\partial x}}{\text{ }}\mathop \int \limits_z^{{h_s}} {\sigma _{xx}}{\text{ }}dz' + {\text{ }}\frac{\partial }{{\partial y}}{\text{ }}\mathop \int \limits_z^{{h_s}} {\sigma _{zy}}{\text{ }}dz' - {\sigma _{zz}} = \rho g\left( {{h_s} - z} \right) \nonumber\\ &\quad + \rho \mathop \int \limits_z^{{h_s}} \frac{{{\partial ^2}W}}{{\partial {t^2}}}dz'.\end{align}

Thus, by combining Eqns (5)–(7), we obtain the depth-integrated momentum equations, which are expressed as follows:

(8)\begin{equation}\left\{ \begin{gathered} \frac{\partial }{{\partial {\sigma _x}}}\int_Z^{{h_s}} {{\sigma _{xx}}dz\prime } + \frac{\partial }{{\partial {\sigma _y}}}\int_Z^{{h_s}} {{\sigma _{xy}}dz\prime } - {\sigma _{xz}} = \rho \int_Z^{{h_s}} {\frac{{{\partial ^2}U}}{{\partial {t^2}}}dz\prime } ; \hfill \\ \frac{\partial }{{\partial {\sigma _x}}}\int_Z^{{h_s}} {{\sigma _{yx}}dz\prime } + \frac{\partial }{{\partial {\sigma _y}}}\int_Z^{{h_s}} {{\sigma _{yy}}dz\prime } - {\sigma _{yz}} = \rho \int_Z^{{h_s}} {\frac{{{\partial ^2}V}}{{\partial {t^2}}}dz\prime } ; \hfill \\ \frac{\partial }{{\partial {\sigma _x}}}\int_Z^{{h_s}} {{\sigma _{zx}}dz\prime } + \frac{\partial }{{\partial {\sigma _y}}}\int_Z^{{h_s}} {{\sigma _{zy}}dz\prime } - {\sigma _{zz}} = \rho g({h_s} - z) \\ + \rho \int_Z^{{h_s}} {\frac{{{\partial ^2}W}}{{\partial {t^2}}}dz\prime } ; \hfill \\ 0 \lt x \lt L;\,{y_1}(x) \lt y \lt {y_2}(x);{h_b}(x,y) \lt z \lt {h_s}(x,y). \hfill \\ \end{gathered} \right.\end{equation}

Similar manipulations, for example, yield equations describing ice flow in a two-dimensional ice-flow model (Pattyn, Reference Pattyn2000, Reference Pattyn2002).

Sub-ice water flow is described by the the following wave equation (Holdsworth and Glynn, Reference Holdsworth and Glynn1978):

(9)\begin{equation}\frac{{{\partial ^2}{W_b}}}{{\partial {\text{ }}{t^2}}} = {\text{ }}\frac{1}{{{\rho _w}}}{\text{ }}\frac{\partial }{{\partial {\text{ }}x}}{\text{ }}\left( {{d_0}\frac{{\partial {\text{ }}P'}}{{\partial {\text{ }}x}}} \right) + \frac{1}{{{\rho _w}}}\frac{\partial }{{\partial {\text{ }}y}}\left( {{d_0}\frac{{\partial {\text{ }}P'}}{{\partial {\text{ }}y}}} \right),\end{equation}

where ${\rho _{\text{w}}}$ is the sea-water density; ${d_0}\left( {x,y} \right)$ is the depth of the sub-ice water layer; ${W_{\text{b}}}\left( {x,y,t} \right)$ is the vertical deflection of the ice-shelf base, and ${W_{\text{b}}}\left( {x,y,t} \right) = W\left( {x,y,{h_{\text{b}}}\left( {x,y} \right),t} \right)$; and $P'\left( {x,y,t} \right)$ is the deviation of the sub-ice water pressure from the hydrostatic value.

The boundary conditions for the ice shelf are (i) a stress-free ice surface; (ii) the normal stress exerted by sea water at the ice-shelf free edges and at the ice-shelf base; and (iii) rigidly fixed edges at the grounding line of the ice-shelf. Moreover, the linear combination of the boundary conditions (Konovalov, Reference Konovalov2019) was also applied in the models considered in this study. This linear combination is expressed as

(10)\begin{equation}{\alpha _1}{\text{ }}{F_i}\left( {U,V,W} \right) + {\text{ }}{\alpha _2}{\text{ }}{{\Phi }_i}\left( {U,V,W} \right) = 0,{\text{ }}i = 1,2,3,\end{equation}

where

1. (i) ${F_i}\left( {U,V,W} \right) = 0$ is the typical form of the boundary conditions, that is ${\sigma _{ik}}{\text{ }}{n_k} = {f_i}$, where ${f_i}$ is the given forcing on the boundary and $\vec n$ is the unit vector normal to the surface;

2. (ii) ${{\Phi }_i}\left( {U,V,W} \right) = 0$ is the approximation based on integration of the typical form of the boundary conditions to the momentum equations (Eq (1) or Eq (8));

3. (iii) the coefficients ${\alpha _1}$ and ${\alpha _2}$ satisfy the condition ${\alpha _1} + {\text{ }}{\alpha _2} = 1$.

The boundary conditions for the sea-water layer correspond to the frontal incident wave. They are

1. (i) at $x = 0$: $\frac{{\partial P'}}{{\partial x}} = 0$;

2. (ii) at $y = {y_1}$, $y = {y_2}$: $\frac{{\partial P'}}{{\partial y}} = 0$;

3. (iii) at $x = L$: $P' = {A_0}{\text{ }}{\rho _w}g{\text{ }}{e^{i\omega t}}$, where ${A_0}$ is the amplitude of the incident wave.

The full description of Model 1, based on Eqn (1), is presented in Konovalov (Reference Konovalov2019, Reference Konovalov2020, Reference Konovalov2021a, Reference Konovalov2021b). Model 2 is based on the finite volume method of approximation of momentum equations (1) and is presented in Konovalov (Reference Konovalov2023b, Reference Konovalov2023c). The full description of Model 3, based on Eqn (8), is presented in Konovalov (Reference Konovalov2021c, Reference Konovalov2023a).

3. Model setup

The numerical experiments with forced vibrations were undertaken for a physically idealized ice shelf with the geometry shown in Figure 1. In the undeformed ice shelf, the four edges had coordinates $x = 0,{\text{ }}x = L,{\text{ }}{y_1} = 0,{\text{ }}{y_2} = B$, where $L$ is the plate length along the x-axis and $B$ is the plate width along the y-axis ( $B = {y_2} - {y_1}$, see Eqn (1)).

The ice plate had only one fixed edge (at $x = 0$), while the other edges (at $x = L,{\text{ }}{y_1} = 0,{\text{ }}{y_2} = B$) were free. This is the special case of an ice shelf, which is also known as an ‘ice tongue’ (e.g. Holdsworth and Glynn, Reference Holdsworth and Glynn1978). The intact ice tongue was 2 km in longitudinal extent, 0.1–0.2 km width.

The ‘rolling’ surface morphology (Coffey and others, Reference Coffey2022) was modeled by sinusoidally varying ice thickness:

(11)\begin{equation}H\left( x \right) = {\text{ }}{H_0} + {A_H}{\text{ }}\frac{{{\rho _w}}}{\rho }{\text{ }}\cos \left( {\frac{{2\pi x}}{{\Delta {l_r}}}} \right),\end{equation}

where ${A_H}$ is the amplitude of ice-thickness oscillations, which was considered as a parameter of the models; $\Delta {l_{\text{r}}}$ is the spatial periodicity of the ‘rolls’ ( $\Delta {l_{\text{r}}}$ in the models was equal to 0.5 km); ${H_0}$ in the models was equal to 25 m. Essentially, this ice tongue was considered as a part of the Ward Hunt Ice Shelf.

Taking into account the hydrostatic balance, the elevation of the ice surface ${h_{\text{s}}}\left( {x,y} \right)$ and the elevation of the ice base ${h_{\text{b}}}\left( {x,y} \right)$ are determined by the following equations, respectively:

(12.1)\begin{equation}{h_{\text{s}}}\left( {x,y} \right) = H\left( {1 - \frac{\rho }{{{\rho _{\text{w}}}}}} \right),\end{equation}

(12.2)\begin{equation}{h_{\text{b}}}\left( {x,y} \right) = - H\frac{\rho }{{{\rho _{\text{w}}}}},\end{equation}

where $H$ is the ice thickness (Eqn (11)).

That is, on the Ward Hunt Ice Shelf, both the ice surface and the ice base exhibited sinusoidal changes along the centerline (Fig. 1).

Figure 1. The ice-shelf and the cavity profiles considered in the numerical experiments. 1—ice-shelf surface; 2—ice-shelf base; 3—sea bottom. The amplitude of ice-thickness oscillations ${A_H} = 18{\text{ m}}$. Spatial periodicity ( $\Delta {l_{\text{r}}}$) of the ‘rolls’ is equal to 0.5 km.

The water-layer depth in the case of the intact ice tongue also had sinusoidal variation (Fig. 1).

The periodic structure of the ice tongue is expected to provide Bragg scattering if the double spatial periodicity of the rolls $2{\text{ }}\Delta {l_{\text{r}}}$ is a multiple of the wavelength ${\text{ }}\lambda $. Respectively, the Bragg wavenumbers $k_b^{\left( n \right)}$, at which we expect to observe band gaps in the modeled dispersion spectra, are expressed as

(13)\begin{equation}k_{\text{b}}^{\left( n \right)} = \frac{{\pi {\text{ }}n}}{{{\text{ }}\Delta {l_{\text{r}}}}},{\text{ }}n = 1,2, \ldots \end{equation}

In all performed experiments, the physical properties of ice were defined by the following values: Young’s modulus, E = 9 GPa, and Poisson’s ratio = 0.33 (Schulson, Reference Schulson1999).

4. Numerical experiments

Numerical experiments were carried out using three models with different combinations of parameters ${\alpha _1}$ and ${\alpha _2}$ in Eqn (10). The results presented below were obtained using Model 1 for (i) ${\alpha _1} = 1;{\text{ }}{\alpha _2} = 0$ and (ii) ${\alpha _1} = 0.2;{\text{ }}{\alpha _2} = 0.8$; using Model 2 for (i) ${\alpha _1} = 1;{\text{ }}{\alpha _2} = 0$ and (ii) ${\alpha _1} = 0;{\text{ }}{\alpha _2} = 1$; and using Model 3 for ${\alpha _1} = 1;{\text{ }}{\alpha _2} = 0$.

The supplemental file contains the results obtained from the experiments performed using Models 1 and 3. Here are the results obtained using Model 2.

Figure 2 shows vertical deflections of the ice shelf. In general, using the second type of boundary conditions ( ${{\Phi }_i}\left( {U,V,W} \right)$ in Eqn (10)), the modeling reveal that the deflections are a superposition of a pure bending mode (Lamb-type mode) and a pure torsion mode (Fig. 2b). By determining the distances between the maxima/minima along the centerline deflection profile and then determining the average value, we obtain the wavelength and, accordingly, the wavenumber for a given periodicity $T$ ( $T = \frac{{2\pi }}{\omega }$, where $\omega $ is the frequency of the forcing). Thus, we successively obtain a dispersion curve—the dependence of the wavenumber on the periodicity (or frequency) of the forcing (Fig 3).

Figure 2. The vertical deflections of the ice shelf resulting from the impact of the frontal incident wave were obtained using Model 2 with the period of forcing $T = 4{\text{ s}}$ ( $T = \frac{{2\pi }}{\omega }$, $\omega $ is the frequency of the forcing) in the case of (a) ${\alpha _1} = 1,{\text{ }}{\alpha _2} = 0$ and (b) ${\alpha _1} = 0,{\text{ }}{\alpha _2} = 1$.

Figure 3. Dispersion spectra obtained using Model 2 with ${\alpha _1} = 1,{\text{ }}{\alpha _2} = 0$ for ice shelf geometries differing in the amplitude of ice-thickness oscillations ${A_H}$ (Figure 1): 1— ${A_H} = 0{\text{ m}}$; 2— ${A_H} = 10{\text{ m}}$. The arrowheads on the solid color line (curve 2) indicate the approximate positions of the left and right limits, which approximately define the left and right boundaries of the band gap. The dashed-colored line (in curve 2) indicates the perturbed wavenumber in the band gap. Similar arrowheads were also used in Figure 4 and other figures showing dispersion spectra to indicate of the boundaries of the band gaps.

In Figure 3, dispersion curve 1 was obtained for the ice shelf with constant ice thickness (the case of ${A_H} = 0$). In this case, the dispersion curve consists of sections of monotonic decrease, which are separated by intermode spaces (Figure 3) accompanying the transitions from the $n + 1$ to the $n$ bending Lamb-type mode. In the case of the rolling periodic geometry of the ice shelf (Fig 1), the models yield the dispersion curves that have sections, where the typical relationship with monotonic decrease and intermode spaces is disturbed, as in dispersion curve 2 in Figure 3. In Figure 3, this section of curve 2 is the band gap, which occurs due to Bragg scattering and corresponds to the second Bragg wavenumber $k_{\text{b}}^{\left( 2 \right)} \approx 12.57{\text{ k}}{{\text{m}}^{ - 1}}$. Essentially, these sections with disturbed wavenumber values are defined by comparison with the typical relationship (curve 1 in Fig. 3).

Experiments with ice shelves that have ‘rolling’ surface morphology revealed that there is the threshold value of the amplitude of ice-thickness oscillations ( ${A_H}$ in Eqn (11)), at which the band gaps appear in the dispersion spectra (Fig. 4a) (Konovalov, Reference Konovalov2023a). Essentially, the amplitude of ice-thickness oscillations ( ${A_H}$) determines the depth of cavities at the base of the ice shelf that result from the ‘rolling’ morphology. These cavities are analogous to crevasses at the ice shelf base (Freed-Brown and others, Reference Freed-Brown, Amundson, MacAyeal and Zhang2012).

In Model 2, as in Model 1, the threshold value also depends on the Bragg wavenumber. The first band gap ( $k_{\text{b}}^{\left( 1 \right)} \approx 6.28{\text{ k}}{{\text{m}}^{ - 1}}$), appears in the spectrum at ${A_H} /gt 14{\text{ m}}$ (Fig. 4a), that is, the first threshold value $\left( {{A_H}} \right)_{{\text{th}}}^{\left( 1 \right)} \approx 14\;{\text{m}}$. The second band gap ( $k_{\text{b}}^{\left( 2 \right)} \approx 12.57{\text{ k}}{{\text{m}}^{ - 1}}$), the third band gap ( $k_{\text{b}}^{\left( 3 \right)} \approx 19.04{\text{ k}}{{\text{m}}^{ - 1}}$) and the fourth band gap ( $k_{\text{b}}^{\left( 4 \right)} \approx 25.13{\text{ k}}{{\text{m}}^{ - 1}}$) appear in the spectrum at ${A_{\text{H}}} /gt 1{\text{ m}}$ (Fig 4b and c), that is, the corresponding threshold values $\left( {{A_H}} \right)_{th}^{\left( i \right)} \leqslant 1m;i = 2,3,4$.

Figure 4. Dispersion spectra obtained using Model 2 with ${\alpha _1} = 1,{\text{ }}{\alpha _2} = 0$ for ice-shelf geometries differing in the amplitude of ice-thickness oscillations ${A_H}$ (Figure 1): 1— ${A_H} = 5{\text{ m}}$; 2— ${A_H} = 10{\text{ m}}$; 3— ${A_H} = 12{\text{ m}}$; 4— ${A_H} = 14{\text{ m}}$; 5— ${A_H} = 18{\text{ m}}$; (a) area of the expected first band gap; (b) area of the expected second band gap; (c) area of expected third and fourth band gaps. The arrowheads on the solid color lines indicate the approximate positions of the left and right boundaries of the band gap. The dashed-colored lines indicate the perturbed wavenumber in the band gap.

As in Model 1, the degradation of the amplitude spectrum (in terms of resonances abatement) is also observed at high ${A_H}$ values, exceeding the threshold values ${\left( {{A_H}} \right)_{{\text{th}}}}$ (Fig. 5). In particular, in the amplitude spectrum obtained at the value ${A_H} = 18{\text{ m}}$, which corresponds to the observed fluctuations in ice thickness on the Ward Hunt Ice Shelf, starting from the impact periodicity $T \approx 2.5{\text{ s}}$, there are no resonance peaks in the spectrum (Figure 5).

Figure 5. Amplitude spectra obtained using Model 2 with ${\alpha _1} = 1,{\text{ }}{\alpha _2} = 0$ for ice-shelf geometries differing in the amplitude of ice-thickness oscillations ${A_H}$ (Figure 1): 1— ${A_H} = 5{\text{ m}}$; 2— ${A_H} = 10{\text{ m}}$; 3— ${A_H} = 12{\text{ m}}$; 4— ${A_H} = 14{\text{ m}}$; 5— ${A_H} = 18{\text{ m}}$; (a) area of the expected first band gap (Figure 2a); (b) area of the expected second band gap (Figure 2b).

Figure 6 shows the alignment of the zone of the expected appearance of the first band gap ( $k_{\text{b}}^{\left( 1 \right)} \approx 6.28{\text{ k}}{{\text{m}}^{ - 1}}$) with the resonant peak in the corresponding range of periodicities of the forcing.

In particular, (a) with the amplitude of ice-thickness fluctuations ${A_H}$ equal to $5{\text{ m}}$ (Fig. 6a), the resonance peak is observed at the periodicity ${T_n} \approx 32.54{\text{ s}}$ (i.e. ${T_n} \approx 32.54{\text{ s}}$ is one of the eigenvalues), at which the wavenumber in the dispersion spectrum is about $5.84{\text{ k}}{{\text{m}}^{ - 1}}$ (i.e. ${k_n} \approx 5.84{\text{ k}}{{\text{m}}^{ - 1}}$); (b) with the amplitude of ice-thickness fluctuations ${A_H}$ equal to $10{\text{ m}}$ (Fig. 6b), the resonance peak is observed at the periodicity ${T_n} \approx 44.04{\text{ s}}$, at which the wavenumber in the dispersion spectrum is about $6.16{\text{ k}}{{\text{m}}^{ - 1}}$ (i.e. ${k_n} \approx 6.16{\text{ k}}{{\text{m}}^{ - 1}}$); and (c) with the amplitude of ice-thickness fluctuations ${A_H}$ equal to $12{\text{ m}}$ (Fig. 6c), the resonance peak is observed at the periodicity ${T_n} \approx 55.92{\text{ s}}$, at which the wavenumber in the dispersion spectrum is about $6.39{\text{ k}}{{\text{m}}^{ - 1}}$ (i.e. ${k_n} \approx 6.39{\text{ k}}{{\text{m}}^{ - 1}}$). Respectively, the relative deviation of the corresponding wavenumber ${k_n}$ from the first Bragg wavenumber $k_{\text{b}}^{\left( 1 \right)}$ does not exceed 7%.

Figure 6. Dispersion spectrum and amplitude spectrum, including the area of the expected first band gap, obtained using Model 2 with ${\alpha _1} = 1,{\text{ }}{\alpha _2} = 0$ for ice shelf geometries differing in the amplitude of ice-thickness fluctuations ${A_H}$ (Figure 1): (a) ${A_H} = 5{\text{ m}}$; (b) ${A_H} = 10{\text{ m}}$; (c) ${A_H} = 12{\text{ m}}$.

The distributions of longitudinal stresses ( ${\sigma _{xx}}$) (Figs 7b and 8b) reflect the distributions of vertical deformations along the central line of the ice shelf (Figs 7a and 8a). That is, the maxima $\backslash$minima of longitudinal stresses coincide with maxima $\backslash$minima of vertical deformations in the deformation profile (Figs 7a, b and 8a, b). In particular, for a given eigenmode the distribution of longitudinal stresses has a specific periodic structure, in which the maxima $\backslash$minima are aligned with the antinodes in the mode. Therefore, beyond the band gaps, the longitudinal stress distribution with a specific periodical structure (as in Fig. 7b) is the expected stress distribution in the ice shelf. Vice versa, inside the band gaps the distribution of longitudinal stresses has a quasiperiodic or non-periodic structure (as in Fig. 8b).

Figure 7. (a) Vertical displacement of ice $W$ along the centerline due to the impact of the frontal incident wave. (b) Distribution of longitudinal stress ( ${\sigma _{xx}}$) in a vertical cross-section of the ice shelf along the centerline. (c) Distribution of shear stress ( ${\sigma _{xz}}$) in a vertical cross-section of the ice shelf along the centerline. The amplitude of ice-thickness oscillations ${A_H} = 10{\text{ m}}$, the periodicity of forcing $T = 5\;{\text{s}}$. These distributions were obtained using Model 1 with ${\alpha _1} = 1,{\text{ }}{\alpha _2} = 0$.

Figure 8. (a) Vertical displacement of ice $W$ along the centerline due to the impact of the frontal incident wave. (b) Distribution of longitudinal stress ( ${\sigma _{xx}}$) in a vertical cross-section of the ice shelf along the centerline. (c) Distribution of shear stress ( ${\sigma _{xz}}$) in a vertical cross-section of the ice shelf along the centerline. The amplitude of ice-thickness oscillations ${A_H} = 18{\text{ m}}$, the periodicity of forcing $T = 5\;{\text{s}}$. These distributions were obtained using Model 1 with ${\alpha _1} = 1,{\text{ }}{\alpha _2} = 0$.

The shear stress ( ${\sigma _{xz}}$) is an order of magnitude less than the longitudinal stress (Figs 7c and 8c). Beyond the band gaps, the maximum $\backslash$minimum shear stress is usually observed at the grounding line, where in the models the ice shelf was considered rigidly fixed (Fig. 7c). Within the band gaps, the maximum $\backslash$minimum shear stress is in most cases achieved in a vicinity of ice-shelf terminus, since the incident wave does not penetrate deeply into the ice shelf (Fig. 8c).

5. Discussion and conclusions

1. (1) All suggested models reveal Bragg scattering for an ice shelf with a rolling surface morphology. The modeled Bragg scattering is expressed in the appearance of the anticipated band gaps in the dispersion spectra.

In Konovalov (Reference Konovalov2023a), it was established that there is a threshold value of crevasses depth, at which the first band gap (corresponding to the first Bragg wavenumber) appears in the spectra. For ice shelf with rolling surface morphology, the double amplitude of ice-thickness oscillations ${A_H}$ is a parameter similar to the crevasse depth in a crevasse-ridden ice shelf (Freed-Brown and others, Reference Freed-Brown, Amundson, MacAyeal and Zhang2012). Accordingly, the threshold value of the amplitude $\left( {{A_H}} \right)_{{\text{th}}}^{\left( 1 \right)}$, at which the expected first band gap appears in the spectrum, also exists for the ice shelf with rolling surface morphology.

Analysis of the superposition of the dispersion spectrum and the amplitude spectrum allowed us to establish the following. The amplitude spectra contain resonance peaks corresponding to wavenumbers that are close to the first Bragg value $k_{\text{b}}^{\left( 1 \right)} \approx 6.28{\text{ k}}{{\text{m}}^{ - 1}}$. The appearance of Bragg scattering in a periodic structure, which also moves periodically with the same frequency, implies that the Bragg wavelength ( $\lambda _{\text{b}}^{\left( i \right)}$) and the amplitude of oscillations of the periodic structure ( $a$) satisfy the following condition:

(14)\begin{equation}a \ll \lambda _b^{\left( i \right)}.\end{equation}

In other words, Eqn (14) is the condition for the occurrence of the expected Bragg scattering in a periodically moving structure. Evidently, at resonance, condition (14) is not satisfied. In fact, condition (14) coincides with the condition for applying wave equation (Eqn 9).

Failure to satisfy condition (14) occurs not only when $T = {T_n}$, that is, when the periodicity of the incident wave coincides with the eigenvalue (when $a{\text{ }}\mathop \to \limits_{T \to {T_n}} {\text{ }}\infty $), but also when the periodicity of the forcing falls within the periodicity range containing ${T_n}$ and which is defined by the width of the resonance peak. Therefore, the probability of the alignment of the expected band gap region with the periodicity range, in which the condition (14) is not satisfied, will be higher for a resonant peak with a larger width. The width of the resonance peaks in the amplitude spectra ( $A{\text{ }}vs{\text{ }}T$) increases with increasing periodicity of the forcing (Konovalov, Reference Konovalov2019). Thus, the probability of the alignment is higher for the first Bragg value $k_{\text{b}}^{\left( 1 \right)}$ than for the remaining $k_{\text{b}}^{\left( i \right)},{\text{ }}i = 2,3 \ldots $

These conclusions are confirmed by the results obtained in the present study. That is, the threshold value $\left( {{A_H}} \right)_{{\text{th}}}^{\left( 1 \right)}$ for the appearance of the first band gap is higher than other threshold values $\left( {{A_H}} \right)_{{\text{th}}}^{\left( i \right)},{\text{ }}i = 2,3 \ldots $ In particular, in Model 1 (with ${\alpha _1} = 1,{\text{ }}{\alpha _2} = 0$) $\left( {{A_H}} \right)_{{\text{th}}}^{\left( 1 \right)} \approx 15\;{\text{m}}$, but $\left( {{A_H}} \right)_{{\text{th}}}^{\left( i \right)} \lt 1,{\text{ }}i = 2,3,4$; in Model 2 (with ${\alpha _1} = 1,{\text{ }}{\alpha _2} = 0$) $\left( {{A_H}} \right)_{{\text{th}}}^{\left( 1 \right)} \approx 14\;{\text{m}}$ and $\left( {{A_H}} \right)_{{\text{th}}}^{\left( i \right)} \lt 1,{\text{ }}i = 2,3,4$; in Model 2 (with ${\alpha _1} = 0,{\text{ }}{\alpha _2} = 1$) $\left( {{A_H}} \right)_{{\text{th}}}^{\left( 1 \right)} \approx 15\;{\text{m}}$ and $\left( {{A_H}} \right)_{{\text{th}}}^{\left( i \right)} \approx 2,{\text{ }}i = 2,3,4$.

An increase in the amplitude of ice-thickness oscillations ${A_H}$ (as an analogue of half the depth of ice crevasses in a crevasse-ridden ice shelf) yields shift in the resonance peaks (Konovalov, Reference Konovalov2021a). Thus, a band gap appears in the dispersion spectrum, if the amplitude of ice-thickness oscillations ${A_H}$ becomes higher than the threshold value ${\left( {{A_H}} \right)_{{\text{th}}}}$.

1. (2) The torsional deformation component in the modes creates additional difficulties in treatment of the dispersion spectra, especially in Model 3, in which this component is observed in any case of the ratio ${\alpha _1}/{\alpha _2}$ (see supplemental file). Essentially, the torsional deformation component yields additional intermode spaces (Konovalov, Reference Konovalov2021a) in the dispersion spectra. These additional intermode spaces respectively provide transitions between torsional strain components in the modes and appear as discontinuities in the dispersion curves. Moreover, these discontinuities differ from the discontinuities accompanying the transitions between the components of the Lamb-type bending deformation in the modes, and can be considered, in particular, as a result of Bragg scattering if these discontinuities are located near the Bragg value. Thus, complementary investigation is required to correctly interpret the discontinuity in the dispersion curve. The investigation is based on the combination of dispersion and amplitude spectra (as shown in Fig. 6). Specifically, the discontinuities in the dispersion spectra corresponding to transitions between the torsional strain components in the modes and looking like band gaps, but not corresponding to Bragg scattering of the incident wave, coincide with the resonance peaks in the amplitude spectra (see section 4 in the supplemental file). In other words, these discontinuities, corresponding to transitions between the torsional strain components in the modes, are accompanied by a transition through resonances, while the band gaps corresponding to Bragg scattering are not accompanied by the same transition (see section (1) of this discussion). Thus, this spectral difference allows us to establish the type of discontinuity in the dispersion spectra.

2. (3) In the models considered in this study, the band gap becomes the dominant effect and abates the resonances in the amplitude spectra if the amplitude of ice-thickness oscillations ${A_H}$ exceeds the threshold value ${\left( {{A_H}} \right)_{{\text{th}}}}$. Thus, it can be said that the abatement of the incident wave by ice shelf with rolling surface/base morphology protects the ice shelf from dangerous resonant impact. For example, the range of periodicities, where the first band gap is observed (Fig. 4a), intersects with the range, where infragravity waves were observed: the range of periodicities is 50.250 s (Bromirski and others, Reference Bromirski, Sergienko and MacAyeal2010). That is, the first band gap, in particular, can protect the ice shelf from the impact of infragravity waves. This abatement may explain how multiyear sea ice in the Arctic Ocean along the coast of Ellesmere Island can be sufficiently stable and long-lived to evolve into the ice-shelf, once contiguous along the Ellesmere Island coast, reported by European and American explorers in the late 19th and early 20th centuries.