Ocean wave blocking by periodic surface rolls fortifies Arctic ice shelves

Abstract

The Ward Hunt and Milne ice shelves are the present-day remnants of a much larger ice shelf that once fringed the coast of Ellesmere Island, Canada. These ice shelves possess a unique surface morphology consisting of wave-like rolls that run parallel to the shoreline. Setting aside the question of how these rolls originally developed, we consider the impact of this roll morphology on the stability of the ice shelf. In particular, we examine whether periodic variations in ice-shelf thickness and water depth implied by the rolls prevent the excitation of Lamb waves in the ice shelf. Using a hierarchy of numerical models, we find that there are band gaps in the flexural and extensional modes of the ice shelf, implying the existence of frequency ranges that lack wave motion. We show that an ice shelf with rolls is able to reflect waves in these frequency ranges that are incident upon its ice front, thereby mitigating undue stress and calving. We speculate that the roll morphology provides a “fitness” for survival that explains why rolls are observed in the oldest and thickest multiyear sea ice of the Arctic.

1. Introduction

The most extensive collection of ice shelves remaining in the Arctic is found along the northern coast of Ellesmere Island in the Nunavut region of Canada. These ice shelves, which include the Ward Hunt, Milne, and Peterson ice shelves, are derived from a larger “ancient” ice shelf called the Ellesmere Ice Shelf. This ice shelf formed about 5500 years ago (England and others, 2008) from the buildup of marine ice that thickened through basal accretion and surface accumulation (Dowdeswell and Jeffries, 2017). By the end of the 19th century, this ice shelf covered an area up to 9000  km² in size (England and others, 2017), though episodic calving in the 20th century reduced its extent to less than a tenth of its original size (Fig. 1A). In contrast to ice shelves in Antarctica, the Ellesmere ice shelves have little to no glacial sources, which means that rather than flowing off the land, these ice shelves are comprised of relatively static floating ice known as multiyear landfast sea ice (MLSI). Today, the dominant mechanism of mass loss continues to be iceberg calving, though these ice shelves are modified by a host of in situ processes including basal melting and surface ablation (Braun, 2017; Dowdeswell and Jeffries, 2017). Despite more than a century of losses, these ice shelves continue to constitute the thickest and oldest MLSI in the Northern Hemisphere.

Figure 1. a) Map of northern Ellesmere Island with the past extent of the contiguous Ellesmere Ice Shelf estimated from descriptions by Robert Peary in 1906 (blue), remnant ice shelves when they were first mapped from 1959 air photos (purple) and in 2020 (red) following a calving event on the Milne Ice Shelf. Source: Mueller and others (2017), supplemented by Sentinel-1 imagery from August 2020. b) Surface roll morphology on the Ward Hunt Ice Shelf (west) visible in a digital elevation model (Porter and others, 2018). c) An oblique aerial photograph of the rolls with meltwater visible on the surface, courtesy of Denis Sarrazin, CEN, Université Laval.

The dominant characteristic of these ice shelves is their undulating surface topography (Fig. 1B, C), consisting of a wave-like “roll” pattern with meters-high amplitude and wavelength of hundreds of meters (Hattersley-Smith, 1957). Early explorers to the Ellesmere Ice Shelf (Aldrich, 1877; Peary, 1907) found these rolls notable for encompassing the entire span of the ice-shelf area, being largely uni-directional and oriented parallel to the northern coast of Ellesmere Island. Ice shelves elsewhere in the world display limited roll-like formations or “rumples” in areas of horizontal compressive stress. Unlike the rolls found on the Ellesmere ice shelves, rumples decay in amplitude with distance from a shoreline onto which ice flow is compressively directed (e.g. Collins and McCrae, 1985). In essence, there are no present-day analogs for the roll morphology found on the past and present ice shelves of Ellesmere Island, making roll formation inherently difficult to study. Thus, the exact origin of the rolls remains both a mystery and a subject of speculation (e.g. Hattersley-Smith, 1957; Jeffries, 1992; Coffey and others, 2022).

One approach to unraveling this mystery is to consider the effect of the roll morphology on the overall stability of the ice shelf. In the present work, we explore the possibility that surface rolls could serve as a defense against destabilizing processes related to calving and breakup. While other surface structures such as rifts and crevasses are typically signs of vulnerability in an ice shelf, they can also play a stabilizing role depending on how and where they form (e.g. Glasser and others, 2009; Lai and others, 2020). Because formation and function are closely intertwined, identifying the functional effects of surface rolls may bring us closer to understanding the processes that form them.

Presently, there is only evidence to suggest that the rolls on the Ward Hunt and Milne ice shelves contribute to breakup (Hattersley-Smith, 1957; Mueller and others, 2017). This is demonstrated by the tendency for fractures associated with calving events to form along the troughs between roll crests (Hattersley-Smith, 1957; Mueller and others, 2003). Additional weaknesses induced by the roll morphology result from the extensive meltwater ponds that form in the valleys between rolls (visible in Fig. 1c), leading to flexural stress that can predispose the ice shelf to hydrofracture (Banwell and Macayeal, 2015). Such meltwater effects are frequently considered to be the proximal cause of ice-shelf instability in Antarctica, which is often justified by the observation of copious meltwater on or below the surface of ice shelves just prior to their collapse (Scambos and others, 2000, 2009).

The apparent weaknesses imparted by roll morphology are likely associated with the post-industrial warming of the Arctic, which has led to regionally elevated melt rates (Fisher and others, 2012), thereby thinning these ice shelves over several decades (Mueller and others, 2003; Braun, 2017). However, the Ellesmere Ice Shelf formed thousands of years ago and persisted through a period of time when Arctic melt seasons were relatively short and melt rates were modest (Fisher and others, 2012; Lecavalier and others, 2017). These melt rates were lowest during the Little Ice Age (1600–1900 A.D.), when the ice shelf likely grew and reached its maximum (England and others, 2017; Dowdeswell and Jeffries, 2017), suggesting that the ice shelf was well-suited to the conditions of the period. In this cooler landscape, surface melting may have played an inconsequential role compared to other threats to ice-shelf stability, most notably sea swell waves impinging from the Arctic Ocean. Therefore, we cannot fully understand the functional effects of rolls in ice shelves without first examining their response to impinging sea swell, which may have been the greater threat to ice-shelf stability in past colder conditions. Given their uniform geometry and extensive coverage, we speculate that surface rolls may have been a favorable adaptation that led to the stabilization of the ice shelf during a period of thousands of years.

While studies of recent ice shelf collapse in Antarctica emphasize the role of surface melting and hydrofracture, the additional contribution of sea swell is gaining recognition. Low-frequency sea swell (0.01–0.1 Hz) is routinely observed by seismographs on ice shelves (Cathles and others, 2009; Bromirski and others, 2017; Chen and others, 2018). These waves can possess amplitudes sufficient to induce significant flexural stress and even initiate icequakes within the rifts of an ice shelf (Chen and others, 2019). Furthermore, ice shelf collapse has often occurred at times when augmented sea swell was present in addition to surface meltwater, as witnessed in the collapses of the Larsen B, Wilkins, and Conger-Glenzer ice shelves of Antarctica (Massom and others, 2018; Bromirski and others, 2010; Walker and others, 2022). A combination of distant long-period sea swell (<0.1 Hz) and localized short-period sea swell (0.1–0.4 Hz) were coincident with significant calving on the McMurdo Ice Shelf in 2016 (Banwell and others, 2017).

The aim of this study is to understand how the surface roll morphology affects the interaction of the ice shelf with sea swell in frequency ranges that are relevant to ice shelf fracture. We use a similar approach to Freed-Brown and others (2012), MacAyeal and others (2012), Sergienko (2013), and Konovalov (2019) for the analysis of wave reflection near the ice front by focusing on the eigenmodes of the coupled ice-shelf/ocean system using the elastic physics of the ice shelf as well as the hydrodynamic properties of water underneath the ice shelf. Similar studies have revealed that periodic variations in the flexural rigidity of an ice shelf (Freed-Brown and others, 2012; Konovalov, 2019) or the height of the water column (Davies, 1982; Davies and Heathershaw, 1984; Hara and Mei, 1987; MacAyeal and others, 2012) can block the propagation of both flexural waves and surface-gravity waves respectively, so we expect a similar phenomenon to occur in ice shelves with rolls. This blocking effect is also known as Bragg scattering (Bragg and Bragg, 1913) and may serve as a mechanism for protecting the interior of the ice shelf from wave-induced flexure that could lead to fracture and instability.

The eigenmodes considered here involve periodic, time-dependent elastic deformations of ice that are either vertically symmetric or antisymmetric with respect to the horizontal mid-plane of the ice shelf (Fig. 2). These two types of modes are referred to as S ₀ and A ₀ Lamb waves, respectively, in elastic plate theory (Viktorov, 1967). Lamb waves with more complex vertical structure (labelled A _n and S _n where $n = 1,\; \, 2,\; \, \dots$) are not considered in this study because they operate at frequencies where the wavelength is much shorter than ice thickness (Chen and others, 2018). Glaciologists will identify S ₀ with longitudinal (extensional) waves and A ₀ with transverse (flexural) waves under circumstances where the thickness of the ice shelf is much less than the wavelength of the Lamb wave (i.e. the thin-plate approximation). Flexural-gravity waves occur when an ocean wave penetrates the sub-ice cavity and propagates as a coupled elastic/hydrodynamic wave, while extensional waves are induced by an energy exchange at the ice front that causes full-thickness horizontal compression (Chen and others, 2018). Though flexural-gravity waves have been documented extensively in ice shelves, extensional waves have only been detected recently on the Ross and Pine Island Glacier ice shelves (Chen and others, 2018). There is evidence that these waves may be harmful to ice shelves by facilitating the expansion of existing fractures (Bromirski and others, 2017).

Figure 2. Example of a flexural and extensional wave (left and right, respectively) with the horizontal ice shelf displacement u ₁ shown in grayscale and the velocity potential Φ of the sub-ice ocean layer in rainbow. Instantaneous water movement is indicated by vectors in the water layer (lower panels), and instantaneous elastic deformation of the ice shelf is indicated by the vectors in the upper panels. This illustration is for an ice shelf of length L = 10 km and waves possessing wavenumber k = 0.225 cycles/km. A flexural wave in the ice shelf is coupled to a surface gravity wave in the ocean. Meanwhile, an extensional wave in the ice shelf leads to the thinning of ice in the vertical direction due to the Poisson effect; instead of succumbing to vertical elongation, the fluid flow buttresses the overlying ice shelf, generating peaks on the upper surface of the ice.

In the development below, the influence of roll morphology will be quantified by comparing the eigenmodes and eigenfrequencies associated with rolls to those associated with a homogeneous ice shelf of constant thickness (i.e. without rolls). The comparison will show that the spectrum of eigenfrequencies for the homogeneous ice shelf is continuous, whereas in the presence of rolls the spectrum will have multiple band gaps representing frequency ranges where no eigenmodes exist. Because rolls influence both the elastic properties of the ice and the hydrodynamic properties of the ocean, we shall consider both elastic-only and hydrodynamic-only idealizations of the ice-shelf/ocean in addition to the fully coupled system where both elastic and hydrodynamic effects are operative.

To better illustrate the implications of these band gaps we provide a time-based demonstration where a wave packet impinges on two different ice shelves, one with rolls and one without rolls. For the case of an ice shelf with rolls, the wave packet reflects from the ice front so long as the wave packet is composed of waves with frequency in the band gap. For the case of an ice shelf with no rolls, the wave packet crosses the ice front and continues as a transmitted flexural-gravity wave into the ice shelf interior. Finally in our discussion we speculate on how this wave reflection property may have protected the nascent Ellesmere Ice Shelf as it formed, thereby presenting a possible narrative on how the roll morphology was generated.

2. Model development

In this study we consider an ice shelf with an unconventional geometry, namely, an undulating roll morphology. Because the roll morphology has translational symmetry along the direction of the shoreline, we use a two-dimensional cross-section to represent the ice shelf (Fig. 3). The length L of the Ward Hunt Ice Shelf from ice front to grounding line is around 10 km (Mueller and others, 2017), the wavelength of the rolls ℓ ranges from 90 to 450 m, and the roll height a is between 1–4 m high (Jeffries, 1992). Although these parameters fluctuate across the ice shelf area, we will assume a constant roll height of a = 3.5 m, wavelength ℓ = 333.3 m and ice shelf length L = 10 km. High resolution bathymetric data is unavailable for this region but recent bathymetric maps for the Arctic Ocean area suggest that the depth D in this region is less than 100 m (Jakobsson and others, 2020); we will assume D = 50 m for the purpose of this analysis.

Figure 3. Schematic cross-section of the partial length of a two-dimensional ice shelf, with the shelf front on the left and the grounding line on the right. The surface elevation is S(x), the basal elevation is B(x), the thickness of the ice shelf is H _i(x), the height of the water column is H(x), the height each surface roll is a, the wavelength of each surface roll is ℓ, and the minimum thickness of the ice shelf is T _i. The length scale of the x-axis is on the order of kilometers, while the z-axis is on the order of tens of meters. The sea level is z = 0, the depth is given by D, and the length of the ice shelf is L.

Because the ice shelf is in hydrostatic equilibrium, we expect that the ice shelf is proportionally thicker under the rolls and thinner under the troughs (Hattersley-Smith, 1957). This is supported by the tendency of ice islands to calve from the ice shelf along the troughs, suggesting they are structurally thinner (Hattersley-Smith, 1957). Although the basal morphology of the ice shelf is uncertain due to the shortcomings of radio-echo sounding in saline ice, data from ground-penetrating radar reveals that the base of the ice shelf has a similar undulating topography as its surface (Toeduschuk and Verrall, 1991; Mortimer and others, 2012; Jeffries, 2017). Moreover, it is unlikely that the ice shelf has sufficient rigidity to support a roll pattern on the surface without a hydrostatically balanced reciprocal pattern on the underside (Hattersley-Smith, 1957).

We assume that the weight of each differential section of ice shelf dx is supported by the buoyant force of the section given by Archimedes principle:

(1)$$g \rho_i H_i dx = - g B \rho_w dx$$

where saltwater has density ρ _w = 1026 kg/m³ and ice is assumed to have uniform density of ρ _i = 917  kg/m³ (Timco and Frederking, 1996). Using the fact that the ice shelf thickness can be expressed as H _i(x) = S(x) − B(x) we have that

(2)$$B( x) = {\rho_i \over \rho_i-\rho_w } S( x) $$

which gives the formula for the inferred basal elevation B(x) based on the observed surface elevation S(x). Here we represent the ice shelf surface height S(x) as a cosine function with amplitude a/2 and wavelength ℓ. The maximum ice thickness of the Ward Hunt Ice Shelf is estimated to be around 40–45 m based on ice core data (Jeffries, 1991), seismic methods (Crary, 1958) and freeboard measurements from two ice islands that calved from the ice shelf (Hattersley-Smith, 1963). Subtracting off the size of the surface rolls and the inferred size of the bottom rolls, we can assume that the ice has a minimum thickness of T _i = 10 m in the valleys between roll crests.

2.1 Elastic ice shelf dynamics

Let u(x,  z,  t) represent the displacement field for a continuous medium, where x and z represent horizontal and vertical coordinates, respectively, and t represents time. We can write the conservation of momentum in terms of the stress tensor $\bf \sigma$:

(3)$$\rho_i {\partial^2 {\bf u}\over \partial t^2} = \nabla \cdot {\bf\sigma} $$

where $\nabla = ( \partial /{\partial x} ,\; \, {\partial }/{\partial z} )$. The relationship between the stress tensor σ and the strain tensor ${\boldsymbol \epsilon }$ is given by Hooke's law (Timoshenko and Goodier, 1951):

(4)$$\sigma_{ij} = {E\over 1 + \nu} \left(\epsilon_{ij} + {\nu\over 1-2\nu} \delta_{ij} \epsilon_{kk} \right)$$

where subscripts i,  j,  k represent spatial indices, δ _ij is the Kronecker delta, and ice is assumed to have Young's modulus E = 7 GPa and Poisson ratio ν = 0.33 (Weeks and Assur, 1967; Murat and Lainey, 1982). The strain tensor is directly related to the displacement field by the following definition:

(5)$$\epsilon_{ij} = {1\over 2} \left(u_{i, j} + u_{\,j, i} \right)$$

where subscripts following a comma represent partial derivatives with respect to the designated spatial index. Substituting (4) and (5) into (3) gives us the elastic momentum equation in terms of displacement:

(6)$$\rho_i {\partial^2 {\bf u}\over \partial t^2} = {E ( 1-\nu) \over ( 1 + \nu) ( 1-2\nu) }\nabla ( \nabla \cdot {\bf u} ) - {E\over 2( 1 + \nu) } \nabla \times ( \nabla \times {\bf u}) $$

Though this equation is often written in a slightly different form involving the vector Laplacian (Landau and Lifshitz, 1986; Freed-Brown and others, 2012), the form presented here makes it more clear that the solutions to this equation are plane waves that induce either compression or shear strain. For example, if we take the Helmholtz decomposition of u:

(7)$${\bf u} = \nabla \phi + \nabla \times \psi$$

Then substituting (7) into (6), we arrive at an equation for two separate types of waves, one which propagates longitudinal displacements (parallel to wave vector) with constant speed:

(8)$$c^2_\ell = {E ( 1-\nu) \over \rho_i( 1 + \nu) ( 1-2\nu) }$$

And one which propagates transverse displacements (perpendicular to wave vector) with constant speed:

(9)$$c^2_t = {E \over 2 \rho_i( 1 + \nu) }$$

Comparing these two equations it is seen that longitudinal waves travel more than 4/3 the speed of transverse waves. In the absence of any boundary conditions, the solutions to (6) represent non-dispersive plane waves traveling through an infinite medium.

The nature of these waves changes dramatically when traveling through a solid plate. Such waves are known as Lamb waves, and the application of zero normal stress on the boundaries leads to the characteristic dispersion for symmetric (S _n) modes (Lamb, 1917):

(10)$${\tanh( \beta H_i / 2) \over \tanh( \alpha H_i / 2) } = {4 \alpha \beta k^2\over ( k^2 + \beta^2) ^2}$$

and the reciprocal dispersion for antisymmetric (A _n) modes:

(11)$${\tanh( \beta H_i / 2) \over \tanh( \alpha H_i / 2) } = {( k^2 + \beta^2) ^2\over 4 \alpha \beta k^2}$$

where $k = \displaystyle {2\pi \over \lambda }$ and α, β are given by

(12)$$\alpha^2 = k^2-{\omega^2\over c_\ell^2}$$

(13)$$\beta^2 = k^2-{\omega^2\over c_t^2}$$

where ω = 2π/T is the angular frequency. Equations (10) and (11) have several branches of intractable solutions (Viktorov, 1967) that reveal the complex nature of these waves. Nevertheless, we are only interested in the most fundamental modes (S ₀ and A ₀) in the case where the wavelength of the wave is large relative to the height of the plate (kH _i ≪ 1). In the limit of small wavenumber/frequency, equation (10) can be simplified to give the velocity of extensional (S ₀) waves in a thin plate (derived independently in Landau and Lifshitz, 1986):

(14)$$c^2_\ell = {E\over \rho_i ( 1-\nu^2) }$$

While in the same limit, equation (11) can be reduced to the following relation for flexural (A ₀) waves (also found in Landau and Lifshitz, 1986; Freed-Brown and others, 2012):

(15)$$\omega^2 = {D k^4\over H_i \rho_i}$$

where ω is the frequency, k is the wavenumber, and $D = {EH_i^3}/{12( 1-\nu ^2) }$ is the flexural rigidity. In contrast to extensional waves, flexural waves traveling through a thin plate do not have constant velocity but are instead dispersed according to (15). Note that when the wavelength is small relative to the thickness of the ice shelf (H _i k ≫ 1), the medium is effectively infinite in all directions, in which case the velocity of the flexural waves reverts back to (9). Comparing (14) and (15), one will find that extensional waves travel much faster than flexural waves for small k.

Lastly, we consider conditions on the boundaries of the ice shelf. At the grounding line, we assume for simplicity that the displacement is zero:

(16)$${\bf u}( x = 0,\; z,\; t) = 0$$

When we are looking at purely elastic waves (where ice shelf motion is not coupled to the fluid flow), we assume that the stress on all other boundaries is zero:

(17)$${\bf n} \cdot {\bf\sigma} = 0$$

where n is the vector normal to the boundary. This boundary condition will later be modified once we couple fluid-induced stresses to the boundary of the ice shelf (see below).

2.2 Sub-ice shelf hydrodynamics

We use linear wave theory to model the flow of water in the sub-ice ocean (Stoker, 1957). First we assume that the fluid flow is incompressible, irrotational, and inviscid, so that Laplace's equation is the governing equation for the velocity potential Φ(x,  z,  t):

(18)$$\nabla \cdot {\bf v} = \nabla^2 \Phi = 0$$

Next we consider the conditions on the boundaries of the ocean domain. We use the kinematic boundary condition on the sea floor:

(19)$${\partial \Phi\over \partial {\bf n}} = 0$$

Then we let η(x,  t) represent the water height perturbation at the upper surface. In the absence of the ice shelf, the kinematic condition at the upper surface of the water is expressed as:

(20)$${\partial \Phi\over \partial {\bf n}} = {\partial \eta\over \partial t}$$

We also have the linearized dynamic boundary condition given on the upper surface:

(21)$${\partial \Phi\over \partial t} + g \eta = 0$$

This problem has the well-known solution:

(22)$$\eta( x,\; t) = A \sin( kx - \omega t) $$

(23)$$\Phi( x,\; z,\; t) = {A\omega\over k \sinh( kH) } \cosh( kz) \cos( kx - \omega t) $$

where A is the amplitude of the wave and H is the height of the water column. Equations (22) and (23) will later be used as initial conditions in order to send a wave packet toward the ice shelf. The frequency ω is related to the wavenumber k by the dispersion relationship:

(24)$$\omega^2 = gk\tanh( kH) $$

Now we consider the boundary conditions when the upper surface of the water is the water-ice interface. Because the sub-ice water surface and ice shelf bottom are in constant contact and move synchronously, we write the kinematic boundary condition in terms of the displacement field in the ice u:

(25)$${ \partial \Phi \over \partial {\bf n}} = {\bf n} \cdot {\partial {\bf u}\over \partial t}$$

At the water-ice contact, we require that the pressure is constant across the interface (Sergienko, 2017):

(26)$${\bf n} \cdot {\bf \sigma} = [ \rho_w g ( z + u_2) - P_{w} ] {\bf n} ,\; \quad z = B( x) $$

where B(x) is the spatial elevation of the ice shelf base, and the wave-induced pressure P _w comes from the linearized Bernoulli equation:

(27)$$P_{w} = - \rho_w {\partial \Phi\over \partial t}$$

Since we are concerned with the motion of the ice shelf away from rest, the hydrostatic pressure can be ignored. Also note that the traction vector in (26) is parallel to the normal vector at the boundary, therefore the shear stress on the ice shelf is zero.

The dispersion relation for water-coupled flexural-gravity waves is similar to (24) but with an added ice flexure term:

(28)$$\omega^2 = \left({D\over \rho_w} k^5 + g k \right)\tanh( kH) $$

Characteristic values of the flexural rigidity D for this problem are in the range of 10¹¹ − 10¹⁴ Pa ⋅ m³, therefore flexural-gravity waves may travel several orders of magnitude faster than surface gravity waves in ice-free water.

2.3 Eigenmode analysis

The general solution to both the elastic motion and the fluid flow can be found by using the technique of separation of variables. If we write the ice shelf displacement as

(29)$${\bf u}( x,\; z,\; t) = e^{i \omega t} {\bf u}_0( x,\; z) $$

Then substitution into the equation of motion (6) results in the following eigenvalue problem:

(30)$${\cal L} ( {\bf u}_0 ) = - \omega^2 {\bf u}_0$$

where ${\cal L}$ represents the second-order linear differential operator from equation (6):

(31)$${\cal L}( {\bf v}) = \displaystyle{E\over 2 \rho_i ( 1 + \nu) }\left[{2-2\nu\over 1-2\nu} \nabla \Big(\nabla \cdot {\bf v}\Big)- \nabla \times \Big(\nabla \times {\bf v}\Big)\right]$$

Note that boundary conditions (16) and (17) do not depend on time so they are left unchanged. Equation (30) will have infinitely many eigenvectors u₀ representing spatial modes of elastic motion and eigenfrequencies ω corresponding to the natural frequency of this motion.

The same procedure can be done for the sub-ice fluid flow by assuming a velocity potential of the form

(32)$$\Phi( x,\; z,\; t) = e^{i \omega t} \Phi_0( x,\; z) $$

which leads to Laplace's equation for Φ₀ and the replacement of first order time derivatives ∂/∂t with eigenvalues iω in boundary conditions (25) and (26). Again, the eigenvectors Φ₀ represent the spatial modes of sub-ice wave motion, with ω being their corresponding frequency.

The general solution to the linearized equations of motion can be expressed as a linear combination of eigenmodes. While this eigenvalue problem is difficult to solve analytically with a nonlinear boundary, these equations are readily solved using finite element analysis. Since numerical solutions may vary depending on the discretization used, we compare to analytic theory and test a wide range of discretizations to make sure our results are sufficiently consistent and robust. All calculations were performed using finite-element software COMSOL Multiphysics (version 6.1).

3. Results

In all three models tested (elastic-only, hydrodynamic-only, and coupled elastic/hydrodynamic), we find that adding a spatial roll pattern to the domain significantly alters the waves that are able to propagate. While the eigenmodes of a homogeneous ice shelf occupy a continuous spectrum of frequencies, the eigenmodes of an ice shelf with rolls have several discontinuities, known as band gaps, representing frequencies of waves that are unable to pass through the ice shelf. The lowest order band gap occurs when the wavelength of the wave is twice the wavelength of the rolls $( \displaystyle \lambda = 2\ell )$, and this band gap tends to be the most prominent and most significant for ice-shelf/ocean interactions. When the antinodes of the wave align with the thinnest points in the ice shelf or the sub-ice cavity, the wave is impeded at the largest parts of its phase, making it difficult for the wave to travel. This mechanism explains why we also see band gaps whenever the wavelength of the rolls is some half multiple of the wavelength of the incoming wave $( \displaystyle n \lambda = 2\ell )$, which is known as Bragg's law (Bragg and Bragg, 1913). In our analysis we are mostly concerned with the lowest order band gaps (n = 1 or 2) since higher order band gaps lie outside the range of frequencies that are known to be relevant to ice-shelf/ocean interaction.

For the purely elastic ice shelf model, the eigenmodes are computed using equations (6), (16), and (17) and classified based on whether they are flexural or extensional (Fig. 4). In the absence of any periodic surface structure, the eigenmodes follow the thin plate dispersion relations given in (14) and (15). However if the geometry of the ice shelf is altered by the addition of a roll pattern, the spectrum of eigenmodes will contain several band gaps (Fig. 4). For flexural modes there are two band gaps (0.14–0.3 Hz and 0.7–0.85  Hz), both of which are in frequency ranges of interest. Meanwhile, the extensional modes have a band gap at a broader range (3–6 Hz), which is just above the frequency range of interest.

Figure 4. Spectrum of flexural and extensional wave modes (left and right panels, respectively) found in a fully elastic ice shelf in the absence of an underlying water layer. An ice shelf with constant thickness has a continuous spectrum of eigenmodes, while an ice shelf with rolls has band gaps in several frequency ranges (represented by arrows and shading). The flexural modes (left panel) have two separate band gaps under 1 Hz while the extensional modes (right panel) have a band gap in the 3–6 Hz range. The black lines representing analytic dispersions show that our numerical solutions for the ice shelf without rolls are in agreement with thin plate theory. Both ice shelves have the same average thickness of 27 m. Note that the frequency scales of these two plots are much different, with the extensional modes occupying much higher frequencies than the flexural modes.

Noting that the shape of the sub-ice cavity implied by the rolls may have an impact on the propagation of water waves, we solve for the eigenmodes of these waves separately using (18)–(21) without coupling to the elastic dynamics of the ice shelf (Fig. 5). We find that periodic variations in sub-ice water depth alone are sufficient to induce a band gap, which occurs at a much lower frequency (0.018–0.026 Hz) that is well within the frequency range of interest where potentially damaging sea swell resides. As expected in the absence of rolls, the numerical modes for water waves at constant water depth show ideal agreement with the dispersion relationship (24) given by linear wave theory. Water waves tend to travel much slower than elastic waves in the ice shelf.

Figure 5. Spectrum of eigenmodes of surface gravity waves allowed to propagate freely within the sub-ice cavity in the absence of the effects of ice-shelf elasticity. Both models have the same average depth H = 26 m however the gray markers represent constant depth while the red markers represent depth which varies sinusoidally with 15 m amplitude. The black line represents the analytical dispersion relation for surface gravity waves in a basin with the given height, matching our numerical results. A band gap is indicated by the arrows and light red shading.

Finally, we couple the fluid flow to the elastic motion using boundary conditions (25)–(27). We find that the flexural modes occupy lower frequency ranges than the purely elastic model but higher frequencies than the purely hydrodynamic model (Fig. 6), which is predicted by the dispersion relations (15), (24) and (28). This means that the presence of water will slow down a flexural wave in an ice shelf, while the presence of an ice shelf will speed up a water wave in the ocean. In contrast, the extensional modes are completely unchanged by the presence of water, which suggests that wave coupling to the water flow is relatively weak, since the only contribution to the vertical motion of the ice shelf is the Poisson-induced thinning associated with extension and compression.

Figure 6. Spectrum of eigenmodes in an elastic ice shelf coupled to the sub-ice flow of water. Flexural modes and extensional modes are shown in the left and right panels, respectively. Both ice shelves have the same average thickness of 27 m. The eigenmodes for an ice shelf without rolls are shown in gray, while the eigenmodes for an ice shelf with rolls are shown in color. Analytical dispersion relationships are plotted in black. Periodic variation in both the flexural rigidity of the ice and the water column height implied by the rolls lead to band gaps in both types of wave modes (indicated by arrows and color shading).

When the motion of the ice shelf is coupled to the water underneath, the ice shelf with roll morphology is able to block water waves with relatively long periods (12–36 s), as shown in Fig. 6. This range corresponds to detrimental sea swell waves that fail to be damped by the pack ice cover (Li and others, 2019) and have been previously associated with icequakes (Chen and others, 2019) and ice calving (Banwell and others, 2017). The second order band gap, which occurs for waves with shorter periods (1.5–2.5 s), is likely also important for filtering out some common frequencies of surface gravity waves in the ocean. Again, both the flexural mode spectrum and extensional mode spectrum agree with analytical dispersion relationships (14) and (28) in the case of an ice shelf without rolls.

Lastly, we perform a time-dependent simulation of an incident wave packet interacting with both a rolled ice shelf and a homogeneous ice shelf. When a surface gravity wave with a frequency in a band gap impacts a rolled ice shelf, the ice shelf will respond by flexing heavily at the front and sending all of the energy back out into the ocean (Fig. 7). Though some component of the wave is transmitted, this component decays within a few rolls. Meanwhile, a flat ice shelf will absorb a significant portion of the wave and allow the incident component of the flexural-gravity wave to pass through the ice shelf. This wave can resonate back and forth for a long time due to reflections at the grounding line and at the ice front. Thus, these waves can be particularly damaging since they get “trapped” in the ice shelf and require several passes between the ice front and the grounding line to leak their energy back out into the ocean.

Figure 7. Demonstration of an ocean wave impacting an ice shelf (a) without rolls and (b) with rolls. The frequency of the incident wave (0.06 Hz) was chosen to be within the first flexural mode band gap of the coupled elastic/hydrodynamic system. When the wave initially reaches the ice shelf, there is a transmitted component and a reflected component (depicted by arrows). While the wave is mostly reflected in the rolled ice shelf, there is a significant transmitted component in the homogeneous ice shelf that resonates back and forth between the ice front and grounding line. Only the seaward half of the ice shelf is shown here for ease of display, the scale depicts the length of the horizontal and vertical components. The arrow length is proportional to the wavelength and arrow width is proportional to the amplitude of the depicted wave. Displacement of the water surface and ice shelf is not to scale. (c) Time averaged magnitude of the displacement of the entire ice shelf over the course of 500 s. A homogeneous ice shelf (top panel) is disturbed across its entire length while an ice shelf with rolls (bottom panel) is disturbed only at the ice front and left undisturbed past the first few rolls.

4. Discussion and future work

In this study, we have shown that ice shelves with periodic surface rolls, particularly those found in the Arctic, are able to reflect ocean waves in frequency ranges that are dominated by sea-swell and ordinary gravity-wave forcing, both of which can be harmful to ice shelves. Waves that attempt to impinge upon the ice shelf within these frequency ranges are purely evanescent, meaning that the wave is restricted from propagating past the ice front. Though an ice shelf possessing rolls may flex more at its front edge (Fig. 7), presenting the possibility of fracture at the ice front, we believe that blocking the wave from entering the sub-ice cavity is essential for the long-term stability of the ice shelf. In particular, the flexural mode band gaps, which are found in frequency ranges corresponding to long-period sea swell (0.03–0.08 Hz) and surface-gravity waves (0.4–0.6 Hz), are a crucial mechanism for fortifying Arctic ice shelves.

In contrast, the extensional mode band gaps span higher frequencies (3–6 Hz) than the flexural mode band gaps. In an inviscid medium such as ocean water, these frequencies correspond to capillary waves that are unlikely to have any effect on the ice shelf. For a given frequency, an extensional wave will generally have a much greater wavelength than a flexural wave, therefore the spacing between periodic ice shelf features would need to be much larger in order for these structures to interact with the wave to produce a band gap. For this reason, extensional wave band gaps may be much more relevant to ice shelves with features that are spaced kilometers apart, such as the basal crevasses found on the Larsen C and other Antarctic ice shelves (Luckman and others, 2012; Lai and others, 2020). Ice shelves with features at the right spacing could prevent extensional wave propagation in sea swell frequency ranges, which could be essential to prevent further rifting. This may also explain why some crevasses are relatively stable after they form and some are unstable (Lai and others, 2020). However, we are still uncertain to what extent ocean waves are responsible and capable of generating extensional waves in the ice shelf. To better understand extensional wave propagation, more research needs to be done into the energy transfer process occurring at the ice-shelf/ocean boundary. We also note that it is possible for extensional waves with such large frequencies (3-6 Hz) to be generated through other mechanisms, such as ice rifting and calving or anthropogenic sources.

While the existence of band gaps is quite robust in the presence of some periodic structure, the exact frequency range of the band gap is sensitive to the physical dimensions of the ice shelf-ocean system, such as ice thickness, water depth, and wavelength, amplitude and shape of the periodic structure. While intensifying this structure usually leads to an increase in the size the band gap, the exact dependence of these band gaps on physical parameters is unknown and requires a more sophisticated mathematical treatment which we hope to perform in the future. Because we assumed constant values for these physical parameters when in reality they may vary spatially or they may be unknown, the exact frequency ranges of the band gaps reported here should be considered an estimate rather than exact prediction.

Understanding the dependence of the band gap on physical parameters would also provide valuable insight into the long-term dynamics and processes that shape the ice shelf. For example, knowing how the height, wavelength, and shape of the rolls affect the position of the band gap could tell us whether the process that forms and maintains the ice shelf has adjusted one or more of these parameters over time in order to “tune” out certain frequencies of waves that are most harmful. This would also allow us to predict future changes to the ice shelf, like whether the thinning of the ice shelf or the gradual loss of the roll structure will lead to a permanent loss of the band gap and, by extension, loss of the ice shelf itself.

Another dimension to the problem not considered here is the angle of incidence of the incoming wave. While we have assumed that all waves arrive perpendicular to the ice front, waves impacting the ice shelf at an angle interact differently with the roll morphology, affecting the positioning of a band gap. The exact relationship between lattice spacing, wavelength, and angle of incidence needed to achieve scattering has been described by Bragg's law in the context of crystalline solids (Bragg and Bragg, 1913) and our preliminary results show that this relationship holds for ice shelf structures as well (also studied in Fox and Squire, 1994; Balmforth and Craster, 1999). However, ocean waves traveling through a rolled ice shelf at an angle would likely undergo continuous changes in direction due to the fact that the speed of the flexural-gravity wave depends on the ice thickness. Thus, in order to fully understand how reflection occurs, one would have to simultaneously solve for the nonlinear wave path of the wave due to Snell's law as well as the scattering of the wave due to Bragg's law.

The results presented here strongly imply the presence of band gaps in the Ward Hunt and Milne ice shelves, but verifying their existence and physical relevance requires data and field observations of seismic activity in the Arctic. Deploying seismometers that can measure the polarization and amplitude of incoming flexural-gravity waves would inform us of any peculiar lack of wave energy in certain frequency ranges, which would be evidence of wave scattering. A comparison of this data with background microseism could reveal whether the band gaps correspond to the most prominent waves present in the Arctic Ocean. It is also possible that such measurements will allow us to resolve certain properties of the ice shelf, such as ice thickness, which are currently unkown or disputed.

The work communicated here leads us to speculate that surface rolls may have protected the nascent Ellesmere Ice Shelf as it was being formed. During the process of ice accretion, multiyear sea ice accumulates and the ice shelf progressively grows. At the same time, the sea swell may provide selective pressure, favoring the formation of rolls and breaking up ice that forms in a homogeneous layer of constant thickness. We hypothesize that eventually, surface rolls will dominate over non-roll morphologies, explaining why surface rolls are found on the most extensive ice shelves in the Arctic. In addition to collecting field data, we hope to model ice shelf growth and decay in order to better understand whether ocean wave forcing is sufficient to drive the development of rolls in Arctic-type ice shelves and whether the driving frequency of ocean waves is related to the roll wavelength.

Lastly, there may be several reasons why the roll morphology is unique to this region of the Arctic. These ice shelves face north and are oriented perpendicular to the direction of the transpolar drift, which is responsible for transporting sea ice across the Arctic Ocean and piling it up against the coast of Ellesmere Island. During this process, sea ice can melt and calve; sea ice can ridge and raft; sea ice can tip and tilt. Thus, the roll morphology may simply be a byproduct of the unique processes involved in the formation of the ice shelf from multiyear sea ice. Because Arctic ice shelves are relatively stationary, the mechanism promoting and maintaining the rolls is also able to operate on the ice shelf for long periods of time. Whatever the origin of the rolls may be, the dynamics of sea ice and the geographic orientation of these ice shelves likely create the unique conditions needed to generate this extraordinary structure.