Enigmatic surface rolls of the Ellesmere Ice Shelf

Abstract

The once-contiguous Ellesmere Ice Shelf, first reported in writing by European explorers in 1876, and now almost completely disintegrated, has rolling, wave-like surface topography, the origin of which we investigate using a viscous buckling instability analysis. We show that rolls can develop during a winter season (~ 100 d) if sea-ice pressure (depth-integrated horizontal stress applied to the seaward front of the Ellesmere Ice Shelf) is sufficiently large (1 MPa m) and ice thickness sufficiently low (1–10 m). Roll wavelength initially depends only on sea-ice pressure, but evolves over time depending on amplitude growth rate. This implies that a thinner ice shelf, with its faster amplitude growth rate, will have a shorter wavelength compared to a thicker ice shelf when sea-ice pressure is equal. A drawback of the viscous buckling mechanism is that roll amplitude decays once sea-ice pressure is removed. However, non-Newtonian ice rheology, where effective viscosity, and thus roll change rate, depends on total applied stress may constrain roll decay rate to be much slower than growth rate and allow roll persistence from year to year. Whether the viscous-buckling mechanism we explore here ultimately can be confirmed as the origin of the Ellesmere Ice Shelf rolls remains for future research.

Introduction

The once-contiguous Ellesmere Ice Shelf (Fig. 1) was first encountered by European Arctic explorers in 1876 (Aldrich, 1877), who reported the ice shelf to have had unusual wave-like ‘rolls’ on its surface (Figs. 2 and 3). Since 1876, the Ellesmere Ice Shelf has largely disintegrated, leaving four major remnants: Ward Hunt, Ward Hunt East, Petersen and Milne ice shelves that together covered 504.5 km² in 2015 (and less since the 2020 loss of the front of the Milne Ice Shelf), and nine minor remnants that are individually 0.2–8.9 km² in area (Mueller and others, 2017). As an explanatory note, we use the term ‘ice shelf’ broadly interpreted to distinguish floating ice that is typically greater than tens of meters thickness and deforms under its own weight. Much of what is presently called ice shelf in the Arctic is composed of a mixture of marine and meteoric ice, and originates largely as multiyear landfast sea ice, multiyear landfast sea ice (MLSI) (Dowdeswell and Jeffries, 2017).

Fig. 1. Map of northern Ellesmere Island with the extent of the contiguous Ellesmere Ice Shelf estimated from descriptions by Robert Peary in 1906 (blue, purple and red), remnant ice shelves when they were first mapped from 1959 air photos (purple and red) 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.

Fig. 2. Disraeli Bay and ice rollers on ridge. Sketch of the ice waves in Nares’ expedition report where the account of Pehlam Aldrich, who led a sledging party through the region, is reported (Nares, 1877). With the sketch, Aldrich describes ‘Several low ridges from 30 to 40 feet high, and varying from a few hundred yards to about a mile in length, show up in front of the cliffs. Their general direction being S.E. and N.W., hence on the east coast of the bay they extend at, or nearly at, right angles from the land, while to the south-westward they are nearly parallel with it. I imagine these ridges are composed of hard ice under snow, though I had no means of penetrating it any depth to find whether or no land lay underneath’.

Fig. 3. Remnant of the Ellesmere Ice Shelf (Fig. 1), currently called the Ward Hunt Ice Shelf (Mueller and others, 2017), which displays the characteristic ice-shelf rolls (right side of image, closest to coastline of Ellesmere Island). The image was generated by the data viewer of the 2 m resolution Polar Geospatial Center ArcticDEM dataset (Porter and others, 2018).

The origin of the Ellesmere Ice Shelf, deduced from study of its present-day remnants (Jeffries, 2017), is largely marine, consolidating first from sea ice, and slowly gaining thickness through snow accumulation and basal freezing to the point of becoming integrated into an average ~ 50 m thick ice shelf that deformed in the normal way of ice shelves, i.e. spreading seaward from a grounding line and shearing past lateral boundaries and ice rises (Jeffries, 2017). It had limited glacial sources on land, and fringed the northern coastlines of Ellesmere Island (Mueller and others, 2017), and possibly Greenland (Condron and others, 2020). On its northern side, it faced the highly mobile and broken Arctic Ocean pack ice that compresses toward these coastlines as a result of the Transpolar Drift Stream (Moore and others, 2019). The intriguing surface rolls of the Ellesmere Ice Shelf, still visible on the existing remnants today (Fig. 3), constitute the subject of this study.

We present and evaluate the hypothesis (Holdsworth, 1987) that viscous buckling from sea-ice force (depth-integrated horizontal stress applied to the seaward ice front of the nascent Ellesmere Ice Shelf) creates and preserves the rolls. A schematic diagram of the buckling process from an initially flat ice shelf to one with rolls is shown in Figure 4 (an additional schematic diagram of viscous buckling is Fig. 2 of Brau and others (2013)). Viscous buckling is only one of several disparate processes, including wind-driven snow dune formation and surface hydrology, that have been hypothesized as causal mechanisms over the past few decades (Hattersley-Smith, 1957; Holdsworth, 1987; Jeffries, 2017). We select the viscous-buckling hypothesis for evaluation here because it is yet to be fully analyzed in a manner leading to its verification or denial. Furthermore, the analysis required for this hypothesis can be largely carried out using simple mathematical analysis and models such as recently developed in other contexts (e.g. Brau and others, 2011, 2013). Other hypotheses for roll origin (e.g. involving snow dunes) require detailed knowledge of local wind regimes, melt seasons and surface hydrology processes, for which there is limited information in this region.

Fig. 4. Schematic view (adapted from Fig. 2 of Brau and others (2013)) of an initially flat, floating ice shelf being compressed by sea-ice pressure (P _si). Variables used to measure the deformation in response to the compression are η and ζ, the flexure response and the thickening response, respectively. (A–C) View of the ice-shelf surface rolls developing from zero amplitude (A) to a larger finite amplitude (C) for t > 0. The initial wavelength of the rolls when they are at infinitesimal amplitude, λ _f, is shown in (B). As roll amplitude increases (C), the ice shelf shortens horizontally (measured by increasing $\Xi ( t)$) and the wavelength shortens from its initial value λ(t) < λ_f. In (A–C) the bottom geometry (not shown) is parallel to the surface geometry. (D) and (E) Full thickness cross section along the axis of compression showing ice-shelf thickening response (with rolls artificially absent for clarity). Horizontal dotted lines denote sea level and the increased ice thickness in (E) is denoted by the vertical extent of the shaded region.

The primary driving force of the viscous buckling hypothesis we investigate is a strong sea-ice pressure from the north, which compresses the ice shelf between the sea-ice and the coast of Ellesmere Island. While sea-ice pressure effects have been investigated in the context of iceberg calving (Banwell and others, 2017; Cassotto and others, 2021), here we investigate whether sea ice is actually capable of determining the morphology and flow of an entire ice shelf over significant time periods such as originally suggested by Holdsworth and Traetteberg (1973) from study of ice island T3. We note that while our emphasis is on the origin of the enigmatic rolls of the Ellesmere Ice Shelf, rolls are observed on other ice shelves, such as the McMurdo Ice Shelf, Antarctica (Fig. 5). Results, both supportive and dismissive of the viscous buckling hypothesis, presented here may therefore be applicable to other ice-shelf roll and rumple phenomena.

Fig. 5. (A) Landsat Image Mosaic of Antarctica (LIMA) showing the environs of Cape Armitage and the McMurdo Ice Shelf (Bindschadler and others, 2008). McMurdo Station and Scott Base are indicated by red stars on the southern tip of the cape. Approximate ice-flow directions are shown by curved arrows on the right side of image. Red boxes display the areas pictured in the oblique aerial photos of panels B and C. The scenes the oblique photos are viewing are indicated by the red boxes labeled accordingly, and the view direction for both photos is from the south looking to the north. (B) Oblique aerial photograph (looking north) of rolls along the coastline of Cape Armitage (labeled B in panel A). Rolls in the foreground are in MLSI; rolls in the distance are on the McMurdo Ice Shelf. Roughness on the top of rolls closest to coastline indicates fracture damage. (Photo provided by Anthony Powell, an artist in residence at Scott Base.) (C) Oblique aerial photograph (looking north) of rolls in the open (non-coastal) part of the McMurdo Ice Shelf at the location labeled C in panel A (from the Antarctic Map and Photograph Library, US Geological Survey). The sharp line (running upper center to lower left) between rolls and flat MLSI is the ice front of the McMurdo Ice Shelf. This photo was taken in 1959, and the area of rolls visible has since calved and been transported away as icebergs.

History and background

The first recorded descriptions of the Ellesmere Ice Shelf came from the late-19th century and early-20th century by explorers who made particular note of its unusual, and difficult surface morphology (Aldrich, 1877; Peary, 1907). In his report, expedition leader Captain Sir George Nares recounted a sledge journey led by Lieutenant Pelham Aldrich: ‘After passing Cape Albert Edward, Aldrich refers to the extremely low and level character of the shore, and describes a remarkable formation of what he designates ice “waves”. “Several low ridges from thirty to forty feet high, and varying from a few hundred yards to about a mile in length, show up in front of the cliffs” … The shelving land appears to blend with the ice, which rises in the form of a roller, with a second roller behind it, exactly as water rolls on a beach after a breeze of wind …’ (Nares and Feilden, 1878). About three decades later, Peary (1907) encountered the same terrain and described it as ‘watery hell’ and pronounced ‘the glacial fringe … will form one of the most unique and interesting features of this region to the glacialist … ’.

The first geophysical observations of the Ellesmere Ice Shelf, and of ice islands calved from it (Crary, 1958), conducted in the late-1950s and early-1960s, also noted the unique rolling surface morphology (Hattersley-Smith, 1957; Crary, 1958, 1960). With further research, including radar and borehole drilling, it became clear that the Ellesmere Ice Shelf was an ice shelf that originated primarily by the thickening of land-fast sea ice since ~5500 BP (England and others, 2017), and that only in some areas it was fed by glacier flow from Ellesmere Island. In all early geophysical assessments, and continuing to the present day (e.g. Narod and others, 1988; Mortimer and others, 2012), the rolls on the Ellesmere Ice Shelf (Fig. 3) are identified as exceptional in the understanding of surface morphology of ice shelves throughout the world (Jeffries, 2017).

Past surveys of the rolls on Ellesmere Ice Shelf remnants, MLSI and ice islands (tabular icebergs calved from the Ellesmere Ice Shelf and its remnants) were reviewed by Jeffries (2017). The wavelengths of the rolls were between 87 and 450 m, with a typical wavelength on the order of 210–330 m. Ice-penetrating radar observations reveal challenging aspects of roll geometry to explain. Jeffries (2017) notes that where rolling morphology is observed in the bottom of MLSI, they have the same wavelength, but are offset horizontally by roughly one-quarter of a wavelength from the surface rolls. Additional airborne radar observations have shown cases where the basal topography is flat, or gently undulating with greater wavelength on the order of kilometers (Narod and others, 1988), despite the presence of surface rolls (Holdsworth, 1987). This ice type, often considered to be nascent ice shelf, is thinner (< 20 m) and has rolls ranging from 25 to 127 m in wavelength (but typically 35–80 m; Jeffries, 2017). There is a strong statistical relationship between the MLSI and ice-shelf thickness and the wavelength of the rolls (Jeffries, 2017) which has been cited as a potential clue as to the process or processes responsible for initiating these features (Holdsworth, 1987).

Hypotheses on the origin of rolls

A number of hypotheses for the origin of the rolls on Ellesmere Ice Shelf have been proposed, and none has yet been eliminated through observations and theoretical analysis. The hypotheses listed by Hattersley-Smith (1957), Holdsworth (1987) and Jeffries (2017) all share a common assessment, that surface processes, particularly associated with differential accumulation and ablation and surface hydrology, act to maintain the rolls. The initiation of the rolls, however, is what sets the hypotheses apart.

Hattersley-Smith (1957) discussed several mechanisms for how the rolls originated. His favored mechanism involved elongation of surface melt ponds by consistently roll-parallel prevailing winds. This elongation process, he proposed (and further developed by Crary, 1960), involved wind-driven convection patterns in ponds which would preferentially melt the shore lines in the up-wind and down-wind sides. He further speculated that the present day rolls are ‘fossil snow dunes’, similar to seif dunes in desert terrains, that have been maintained by the elongation of surface melt ponds according to the wind direction.

Following Hattersley-Smith's (1957) initial evaluation of possible origins for rolls, others have been, or can be, added to the list; notably, the addition of basal processes such as internal waves and tides (Narod and others, 1988), and ‘sympathetic’ tides in meltwater ponds (which would be similar to the effect of consistently directed prevailing winds) (MacAyeal and others, 2020).

A class of hypotheses, the one forming the subject of our study, was dismissed by Hattersley-Smith (1957), but has yet to be fully discredited. As described by Crary (1960) and Holdsworth (1987), these hypotheses call on the ice shelf's response to uni-directional compressive stress to initiate the roll orientation, wavelength and amplitude. A plausible origin of the compressive stress is sea-ice pressure perpendicular to the north shore of Ellesmere Island.

Some plausibility is attributed to the compression-stress hypothesis from examples where Antarctic ice shelves display roll-like surface morphology (often referred to as rumples in the Antarctic context). A good example is located in the vicinity of Scott Base and McMurdo Station, on Cape Armitage, Ross Island (Collins and McCrae, 1985). As depicted in Figure 5, rumples are found along the eastern boundary of the cape, where the McMurdo Ice Shelf (25–75 m thick) is pressed into the coastline by the westward flow of the Ross Ice Shelf farther east. As observed in Figure 5a, there are also rolls on the ice shelf well south of Cape Armitage, where the flow is still directed west, but where compressive stresses are likely provided by the more distant coastlines of Brown Peninsula and Black Island. It should be noted, however, that the possibility of thermal origins of rolls has also been explored by Sanderson (1978) in the context of phenomena observed on the George VI Ice Shelf, Antarctica.

Viscous buckling dynamics

To evaluate the sea-ice compression hypothesis for the origin of Ellesmere Ice Shelf's rolls, we construct a simple 2-D (one horizontal and one vertical) model that accounts for flexure and thickness change of a floating ice shelf, treated as a simple viscous fluid. The approach follows MacAyeal and others (2021), wherein the thin plate and shallow shelf approximations are made. We begin by examining sinusoidal, roll-like solutions and determine their temporal growth and decay. Following this, we examine through a series of simple estimation exercises how non-Newtonian ice rheology may allow rolls to persist for significant periods of time (e.g. beyond the period when sea-ice pressure is high).

Governing equations

We consider the deformation of a floating ice shelf that can be termed either an ice shelf or unbroken MLSI. It is fixed by a wall at x = 0 m and ending in a vertical ice front at x = L where compressive forces of sea ice and hydrostatic sea-water pressure are applied in the negative x direction. For simplification, we assume that all properties in one of the two horizontal directions, the one intended to be parallel to the coast of the restraining land mass onto which the ice shelf is compressed, remain uniform. We further assume that there is no ice movement except either to or from the coast of the restraining land mass; thus, making the problem 2-D, with x and z being the horizontal (perpendicular to the restraining land mass, positive toward the sea ice) and vertical spatial coordinates (with z = 0 at sea level and positive upward).

As an initial condition, we take the ice shelf to be uniform in thickness h and everywhere in local hydrostatic equilibrium, where the surface and basal elevations, s and b, respectively, are (1 − γ) h and − γh, respectively, where γ = ρ _i/ρ _sw, and ρ _i and ρ _sw are the densities of ice (917 kg m⁻³) and sea water (1028 kg m⁻³). This represents the imagined initial state of the Ellesmere Ice Shelf as a thin plate of land-fast sea ice having a thickness determined by one or more seasonal cycles. The mid-plane of the ice shelf is designated the neutral plane, and is initially parallel to the ocean surface.

Since our interest is only in how rolls could initiate from compressive forces applied at the ice front, we do not consider accumulation or ablation, taking them to be zero at both the ice-shelf surface and base. This means that all changes to the geometry of the ice shelf are by deformation in response to the applied compression and ice-shelf spreading, and this deformation will be described by three variables η(x,  t), ζ(x,  t) and $\Xi ( t)$ representing vertical deflection of the neutral surface, changes to ice thickness and horizontal (compressive) shortening of the ice shelf, respectively, with x and time t being the horizontal coordinate and time, respectively. (The second horizontal coordinate, y, is assumed to be perpendicular to the compressive forces, and all variables are assumed uniform with y for simplification.) The ice thickness over t is taken to be H = h + ζ(x,  t), with h being a constant initial thickness (not a function of x or t). Following the shallow ice-shelf approximation, we assume that all horizontal strains are independent of z, thus changes to geometry associated with ζ remain in local hydrostatic equilibrium (since the initial state of the ice shelf is also). Changes to geometry that violate local hydrostatic equilibrium, but which maintain hydrostatic equilibrium of the ice shelf as a whole, are described by η(x,  t). Variation of η does not, however, influence ice thickness, hence ζ and η are independent variables describing the geometry of the ice shelf.

The governing equations (derived in Appendix A) are:

(1)$$\rho_{\rm sw}g\eta + P_{\rm si}{\partial^2 \eta\over \partial x^2} + V{\partial^4 \dot \eta\over \partial x^4} = 0,\; $$

(2)$$\rho_{\rm i}g( 1-\gamma) H - {P_{\rm si}\over h} + D\dot H = 0,\; $$

and,

(3)$$\Xi( t) = \int_0^{L-\Xi} \left({1\over 2}\left({\partial \eta( x,\; t) \over \partial x} \right)^2 + {\zeta( t) \over h} \right){\rm d}x$$

where g is the acceleration of gravity, P _si is the net sea-ice force (per unit width of the ice shelf) and the viscous dissipation terms (V and D) are V = (1/3)ν _fH ³, where ν _f is the ice viscosity associated with flexure, and $D = 8 \bar \nu _{\rm eff} H^{-1}$, where $\bar \nu _{\rm eff}$ is the vertical average of the effective viscosity ν_eff(z) (e.g. see MacAyeal and others, 2021). The overdot on variables denotes the time derivative.

For the simple analytical treatment of the next section, we consider the ice to be a Newtonian fluid and take $\nu _{\rm f} = \bar \nu _{\rm eff} = \nu$, where ν is a constant viscosity (the effects of non-Newtonian ice viscosity will be discussed once the solution for Newtonian viscosity is developed). We note that in previous sections as well as in literature cited here, P _si is sometimes referred to as sea-ice pressure, using the term ‘pressure’ in a colloquial manner. In the present study, P _si represents the vertical integral of stress in the sea ice (with units of Pa m or, alternatively, N m⁻¹) per unit width of the ice shelf. This force is assumed to be applied at x = L, the x-coordinate of the ice-shelf frontal boundary.

Equation (1) represents the vertical force balance acting on a viscous plate with constant viscosity ν_f (MacAyeal and others, 2021). The key term in this equation is the second on the left-hand side, P _si(∂²η/∂x ²), where the sea-ice force P _si multiplies the curvature of the neutral axis of the ice shelf. This is the term which rectifies the sea-ice pressure into a vertical force that can deform and displace the ice shelf. Equation (2) represents the evolution of ice thickness in the absence of surface and basal accumulation or ablation. It accounts for normal ice-shelf spreading (first term on the left-hand side) and the effect of back-pressure from the sea ice (second term on the left-hand side). Equation (3) represents the horizontal shortening of the ice shelf as it is compressed by sea-ice pressure. This shortening increases with growth of roll amplitude and with strain-induced thickness change. The overall length of the ice shelf will be L initially, and change to $L-\Xi ( t)$ as deformation evolves (including horizontal extensional deformation when sea-ice pressure is removed). Figure 4 displays the relationship between $\Xi ( t)$ and η(x,  t) and ζ(x,  t).

We note that Eqn (1) requires boundary conditions at x = 0,   L, but that Eqn (2) does not require boundary conditions, because there are no x-derivatives of H. For the present study, where we seek to understand solutions that are periodic over the ice shelf, we shall only consider sinusoidal functions of x to represent solutions of Eqn (1). Specific boundary conditions determine the actual functional variation of η with x, however, since we are interested in only estimating wavelengths and time-behavior of the sinusoidal components of η, we consider the general solutions only.

To best investigate general properties of the system of equations derived above, we adopt dimensionless variables. At the outset, we note that the system derived above already has assumptions about dimensional relationships built-in by our thin plate and shallow ice shelf assumptions (discussed in Appendix A). Specifically, the above equations are valid for:

(4)$$\lambda \gg h\rm\;{ and } \;{\it \eta} \ll {\it h} ,\; $$

where λ is the wavelength of any roll-like solutions that the above system may permit. We mention this because we anticipate that the above system will not permit roll-like solutions when P _si = 0 and λ < h. This issue arises because we have chosen to simplify the problem by not considering the elastic component of deformation and associated energy. By considering visco-elastic deformation, the wavelength of bending has a lower bound (e.g. MacAyeal and others, 2021).

We adopt dimensionless variables (designated by $\tilde \dollar $) to facilitate the solution of Eqns (1) and (2) and to clarify the effects of scales and parameters. We note at the outset that Eqns (1) and (2) are uncoupled, and hence should be rendered dimensionless separately, and thus requires two non-dimensional time variables ${\tilde t_1} \ {\rm and} \; {\tilde t_2}$ for Eqns (1) and (2), respectively. For Eqn (2), we choose:

(5)$$\eqalign{x = h \tilde{x},\; \quad \eta = h \tilde{\eta},\; \ H = h\tilde{H},\; \ t = {V\over \rho_{\rm sw} g h^4} \tilde t_1 \ ,\; \ P_{\rm si} = \rho_{\rm sw} g h^2 \tilde P_{\rm si}} .$$

Substitution into Eqn (1) gives:

(6)$$\left(\rho_{\rm sw} g h \right)\tilde \eta + \left(\rho_{\rm sw} g h^2 \right)\tilde P_{\rm si} \left({h\over h^2} \right)\tilde \eta_{\tilde{x}\tilde x} + \left({\rho_{\rm sw}gh h^4 V\over h^4 V} \right)\dot{\tilde\eta}_{\tilde{x}\tilde{x}\tilde{x}\tilde x} = 0$$

(7)$$\tilde\eta + \tilde P_{\rm si}\tilde \eta_{\tilde{x}\tilde x} + \dot{\tilde \eta}_{\tilde{x}\tilde{x}\tilde{x}\tilde x} = 0$$

where the subscript $\tilde x$ denotes differentiation with respect to $\tilde x$, and the overdot represents time differentiation with respect to $\tilde t_1$. For Eqn (2), we choose:

(8)$$\zeta = h \tilde \zeta \ ,\; \ t = {D\over \rho_{\rm i} g \left(1 - \gamma\right)} \tilde t_2 .$$

Substitution into Eqn (2) gives:

(9)$$\left(\rho_{\rm i} g \left(1 - \gamma \right)h\right)\tilde H - \left({\rho_{\rm sw}g h^2\over h} \right)\tilde P_{\rm si} + \left({\rho_{\rm i} g \left(1 - \gamma\right)D h\over D} \right)\dot{\tilde H} = 0$$

(10)$$\tilde H - \alpha \tilde P_{\rm si} + \dot{\tilde H} = 0$$

where the overdot in this equation represents time differentiation with respect to $\tilde t_2$, and

(11)$$\alpha = {1\over \gamma\left(1 - \gamma \right)} \sim 10 .$$

We note that $\tilde \eta$ and $\tilde \zeta$ evolve with non-dimensional times $\tilde t_1$ and $\tilde t_2$, respectively, that have

(12)$${V\over \rho_{\rm i} g h^4}$$

and

(13)$${D\over \rho_{\rm i} g ( 1 - \gamma) }$$

as their intrinsic timescales, respectively. The ratio of these two timescales is:

(14)$${V ( 1 - \gamma) \over D h^4} = {1 - \gamma\over 24} \approx 0.004,\; $$

when ν _f = ν _eff = ν. This allows us to immediately deduce that growth of $\tilde \eta$ will be much faster than growth of $\tilde \zeta$ when sea-ice pressure is applied. We thus expect at the outset that rolls will primarily manifest themselves as perturbations of $\tilde \eta$.

Results

Solution to homogeneous equations

We note that in the absence of sea-ice pressure, $\tilde P_{\rm si} = 0$, Eqns (7) and (10) are homogeneous (unforced) equations and are written:

(15)$$\dot{\tilde \eta}_{\tilde{x}\tilde{x}\tilde{x}\tilde x} = -\tilde \eta$$

(16)$$\dot{\tilde H} = - \tilde H.$$

Substituting a roll-type (sinusoidal) solution into Eqn (15), with wavenumber $\tilde k = {( 2\pi ) }/{\tilde \lambda }$, wavelength $\tilde \lambda$ and e-folding timescale $\tilde \tau = {1}/{\tilde \sigma }$, i.e. $\tilde \eta \ \propto \ \exp ( \tilde \sigma \tilde t_1 ) \sin ( \tilde k\tilde x )$, we obtain:

(17)$$\tilde k^4 \tilde \sigma = -1$$

which gives

(18)$$\tilde \sigma = -\tilde k^{-4} < 0.$$

We immediately observe that rolls cannot persist in the absence of sea-ice pressure. In fact, the non-dimensional e-folding decay rate $\tilde \tau$ in the absence of sea-ice force is:

(19)$$\tilde \tau = {1\over \tilde \sigma} = - \tilde k^4 = - \left({2\pi\over \tilde \lambda} \right)^4.$$

This means that long wavelengths decay faster than short wavelengths. We can surmise from this that if there is a limited episode of sea-ice pressure which produces a variety of rolls of different wavelengths, the shortest wavelengths will have the greatest chance to persist into the subsequent summer season, when surface melting and related hydrological processes could potentially maintain the rolls.

Turning next to Eqn (16), we can immediately deduce that it will not admit solutions that spontaneously develop sinusoidal, roll-like structure. This is because the equation is first order in $\tilde t_2$ only, i.e. there are no terms which contain derivatives with respect to $\tilde x$. In the absence of sea-ice pressure, $\tilde H$ decays with time, which corresponds to normal, 1-D ice-shelf spreading.

Solutions to inhomogeneous, forced equations

We now seek solutions of Eqn (7) in which $\tilde P_{\rm si}> 0$ (corresponding to compressive force in the ice shelf). Looking first for sinuous rolls (involving only $\tilde \eta$), we again substitute $\tilde \eta \ \propto \ \exp ( \tilde \sigma \tilde t_1 ) \sin ( \tilde k \tilde x)$ into Eqn (7):

(20)$$\tilde \sigma \tilde k^4 - \tilde P_{\rm si}\tilde k^2 + 1 = 0$$

which gives

(21)$$\tilde \sigma = \tilde P_{\rm si}\tilde k^{-2} - \tilde k^{-4}.$$

For the roll solution to grow in time, $\tilde \sigma$ must be >0, and this implies:

(22)$$\tilde \lambda < 2\pi \sqrt{\tilde P_{\rm si}}.$$

This tells us that as $\tilde P_{\rm si}$ grows, wavelength of a growing mode increases. We next determine the fastest growing mode:

(23)$${\partial \tilde \sigma\over \partial \tilde k} = -2\tilde P_{\rm si}\tilde k^{-3} + 4\tilde k^{-5}$$

(24)$$ = 0$$

which implies that

(25)$$\tilde \lambda_{\rm f} = \pi\sqrt{2\tilde P_{\rm si}} \;{\rm which\, is }< \! 2 \pi \sqrt{\tilde P_{\rm si}}.$$

This means that there is a distinct wavelength, $\tilde \lambda _{\rm f}$, that will grow faster than any other. A notable immediate property of $\tilde \lambda _{\rm f}$ is that it depends only on the square root of the sea-ice pressure, and not on ice-thickness, which will be discussed in the next section (Fig. 6A). Substituting $\tilde \lambda _{\rm f}$ into the expression for $\tilde k$ in Eqn (21) gives the value of the fastest e-folding growth rate:

(26)$$\tilde \sigma_{\rm f} = {1\over 4}\tilde P_{\rm si}^2,\; \ \ \rm{which\, gives}$$

(27)$$\tilde \tau_{\rm f} = 4 \tilde P_{\rm si}^{-2}.$$

For future reference, we express the dimensional form of σ and τ:

(28)$$\sigma = {3 P_{\rm si}^2\over 4\rho_{\rm sw}g \nu_{\rm f} H^3} = {1\over \tau}.$$

This shows that the growth rate is inversely proportional to the cube of the ice thickness, and allows us to anticipate that roll creation will be most favored in thin ice shelves.

Fig. 6. (A) Initial wavelength (m) of fastest growing mode as a function of P _si (Pa m) in dimensional units. Initial wavelength is independent of ice thickness H. (B) The e-folding growth timescale (Eqn (28), units of d) as a function of P _si and initial plate thickness h (m). For the results shown here, ice is modeled as a Newtonian fluid with a constant viscosity (10¹³ Pa s) consistent with cold, brine-free ice. Growth rates will increase if viscosity is reduced by effects of marine ice (MacAyeal and Holdsworth, 1986).

We next turn our attention to solutions of Eqn (10) for $\tilde H = 1 + \tilde \zeta ( \tilde x,\; \, \tilde t_2)$. As stated above for the homogeneous form of the equation, this equation is an inhomogeneous first-order ordinary differential equation in $\tilde t_2$ only, and thus solutions with $\tilde x$ dependence will not grow spontaneously with time. The general solution to Eqn (10) is the sum of a solution to the homogeneous equation (where $\tilde P_{\rm si} = 0$) and a particular solution where $\tilde P_{\rm si}> 0$:

(29)$$\tilde H( \tilde x,\; \tilde t_2) = \left(1 - \alpha \tilde P_{\rm si} \right){\rm e}^{-\tilde t_2} + {\alpha}\tilde P_{\rm si} .$$

We note that this solution is convergent in time to a specific ice thickness, where the sea-ice pressure exactly balances the spreading stress:

(30)$$\lim_{\tilde t_2\to\infty} \tilde H = {\alpha}\tilde P_{\rm si}.$$

We again note that the solution for $\tilde H$ (and $\tilde{\zeta}$) does not have $\tilde x$-dependence; hence, we can rule out the possibility of rolls developing from changes in the ice-shelf thickness.

Roll growth rate and wavelength

In Figure 6, we plot the dimensional e-fold growth rate and initial wavelength λ_f(t = 0) for the fastest growing mode of η as functions of the dimensional P _si, and for a viscosity (10¹³ Pa s) consistent with cold, brine-free ice. The value of P _si is varied through 0.1 to 10 MPa m so as to encompass the range of sea-ice pressure sometimes observed on floating ice islands (e.g. Holdsworth and Traetteberg, 1973; Timco, 2010) in the Arctic Ocean. Parameters used to re-dimension the variables are:

(31)$$\eqalign{& h = 1,\; 5,\; 10\, \rm{m},\; \quad {\it g} = 9.81\, \rm{m\,s}^{-2} ,\; \ {\it \rho}_{\rm sw} = 1028 \, \rm{kg\,m}^{-3},\; \ \cr& {\it \nu} = 10^{13} \, \rm{Pa\,s},\; \ {\it V} = {1\over 3} {\it \nu h}^3}.$$

The results show (Fig. 6b) that rolls can develop in relatively thin (i.e. H = 1,  5,  10 m) ice with an e-folding timescale of under ~ 1 year with sea-ice pressures up to 2 MPa m. Rolls grow faster for thinner ice and larger sea-ice pressure. At first, the wavelength of the fastest growing roll (Fig. 6a) is independent of ice thickness, and is ~45 m for a sea-ice pressure of 1 MPa m. Higher sea-ice pressures produce longer wavelength, with wavelength exceeding 100 m for sea-ice pressures of ~5 MPa m. Considering an ice thickness of typical Ellesmere Ice Shelf remnants observed today (~ 50 m), the results show that the e-folding timescale is increased to ~ 2 × 10⁵ d (~550 years). This strongly suggests that rolls must develop in thin ice that is the precursor to the Ellesmere Ice Shelf observed in modern time. The sea-ice pressure required to reduce the e-folding timescale below 100 d for a 50 m ice thickness (not shown in Fig. 6b) is ~50 MPa m, which is above stress scales that ice is capable of sustaining without the risk of fracture (Timco, 2010). Additionally, for thin floating sea ice with thickness of 1 m, the e-folding timescale is increased to ~100 d for P _si reduced to 0.15 MPa m.

The above results show that rolls are unlikely to be created by sea-ice pressure acting on floating ice $\gtrsim$10 m in thickness. The growth rate for ice thicker than ~10 m is so slow that the process would require several centuries of near continuous sea-ice pressure (at ~ 1 MPa m) to build rolls to meters amplitude. A possible mitigating factor to counter this long time scale could be lower ice viscosity, such as is expected when the ice is largely of marine origin (MacAyeal and Holdsworth, 1986). From our analysis, it would take a reduction of ν by 3 orders of magnitude (from ~ 10¹³ to ~ 10¹⁰ Pa s) to achieve an e-folding timescale of ~ 200 d. We do not know if brine content of the marine ice precursor to the Ellesmere Ice Shelf would be sufficient to drop the viscosity by such a large magnitude. We suggest that it is far more likely that the rolls, if created at all by sea-ice pressure, were created in thin (1–10 m) MLSI that initiated the Ellesmere Ice Shelf, with subsequent thickening to ~ 50 m by unrelated processes. This would require that these unrelated processes operate in a manner to preserve the rolls in their initial form. Perhaps, this could be done by bottom ice accretion and/or surface thickening from snowfall.

The fastest growing mode for $\tilde \eta$ (Eqn (7)) is shown by our treatment (Eqn (28)) to have the same initial wavelength, $\tilde \lambda_{\rm f}$, for all values of ice thickness. Ice thickness enters our analysis through the variable V = (1/3)ν _fH ³, the viscous buckling resistance parameter, which appears in only one term of Eqn (7)–that which has the only time derivative of η. Ice thickness and viscosity ν _f can only influence the timescale of the deformation, not the initial spatial scale.

Observations suggest that roll wavelength is variable (Jeffries, 2017) and is correlated with ice thickness; hence, our solutions may appear to be at odds with this key observation. A possible explanation for this discrepancy can be considered. The above analysis (reflected in Fig. 6) only predicts the initial wavelength associated with rolls when they are of infinitesimal amplitude. This initial wavelength will change as the ice shelf continues to shorten in the x direction ($\Xi > 0$) due to continued roll growth (Fig. 4c). Since horizontal shortening, $\Xi ( t) > 0$, is related to both the slope of the rolls (∂η/∂x)² and the fractional ice-thickness growth ζ(t)/h, the developed, evolved rolls will have a wavelength that is inversely proportional to $\Xi$. In particular, at the end of the application of sea-ice pressure, when roll amplitude is maximum, we expect the wavelength to be

(32)$$\lambda( t) = {L-\Xi( t) \over L}\lambda_{\rm f}.$$

This means that at the end of a period of sea-ice compression the wavelength will be shortest in ice that suffers the most overall shortening. Thin plates which develop roll amplitude faster than thick plates, and which fractionally thicken faster than thick plates, will thus have developed, evolved roll wavelength that is shorter. The precise evolution of $\Xi ( t)$ can vary depending on the length of time the sea-ice pressure is applied; hence, it is not possible to generally predict a single correspondence between evolved wavelength and ice thickness. This consideration, however, does address the observation that wavelength tends to be smaller in thinner ice (Jeffries, 2017).

A second possible way to obtain roll wavelength dependence on ice thickness is with elastic deformation, which is influential in choosing the critical initial wavelength for wrinkling of elastic plates and membranes (e.g. Brau and others, 2013). We purposely chose to disregard elastic behavior at the outset of the analysis here, because we were interested in only long-term, permanent deformation, which is what a viscous rheology provides. An assessment of visco-elastic ice rheology is thus called for in future research.

Factors that may control roll persistence

In the previous section, we demonstrated the viscous buckling behavior under the assumption that the ice behaves as a Newtonian fluid, with a single viscosity (10¹³ Pa s). Three important conclusions are evident from the solutions obtained. First, under large compressive sea-ice pressure, sinusoidal flexure of a thin (1–10 m), initially flat and homogeneous floating ice shelf can grow exponentially (while P _si > 0) to yield surface rolls. Second, under the same forcing, ice-thickness can increase uniformly, but does so without developing a roll-like structure. Third, the roll-like features and increased thickness are temporary, and decay exponentially with time once sea-ice pressure is removed.

The question we examine in the present section is: can the non-Newtonian rheology create an asymmetry between the growth and decay phases of roll-like features, thereby allowing rolls to persist in the long term after sea-ice pressure abates? Answering this question is beyond the means of the exact analytic approach used in the treatment of Newtonian rheology in the section above. Therefore, here, we only explore a qualitative answer to this question using various demonstrations and scale estimates loosely based on the analytic solutions derived above. In particular, we will estimate a time series of effective viscosity of the ice shelf through a sequence of roll-building (P _si > 0) followed by roll-decay (P _si = 0). This estimate will show that the growth and decay phases of rolls can operate on different timescales, with the decay phase taking much longer due to reduced effective viscosity when the stress associated with sea ice is absent.

As previously stated, we will use the analytic solutions for η and ζ derived under the assumption of constant Newtonian viscosity to make the estimate and demonstration. This means that there is inconsistency in our approach, and that the estimate and demonstration can only be regarded as a qualitative pointer to results to be obtained in future research, perhaps using numerical methods to cope with the complexity of ice rheology.

Temporal asymmetry of ice viscosity

Rheology of an ice shelf originating through freezing of the ocean surface, augmented by snow accumulation, and of variable age and temperature under varying magnitudes and styles of stress, is complex (e.g. MacAyeal and Holdsworth, 1986). Given our modest goals of only estimating a time series of effective viscosity and demonstrating roll-like structures in ζ, we use a traditional (simplified) treatment of ice rheology,

(33)$$\dot e = A T^n$$

where $\dot e$ is the dominant strain rate, we can write an estimate of the effective viscosity as a function of the dominant stress T:

(34)$$\nu_{\rm eff} = {1\over 2 A T^{n-1}}$$

where A ≈ 2 × 10⁻²⁴ Pa⁻³ s⁻¹ is a value representative of temperate ice (but with zero brine content), and n = 3 (Cuffey and Paterson, 2010). The expression in Eqn (34) represents the heart of the analysis made here: T will be much higher during the application of sea-ice pressure than when the rolls are relaxing, and hence ν _eff will be low when rolls are created and high as they decay. The sea-ice pressure P _si and H ≈ h determine the dominant stress scale T _s in the ice shelf:

(35)$$T_{\rm s} = \left\vert \rho_{\rm i} g ( 1-\gamma) h - {P_{\rm si}\over h} \right \vert \sim {P_{\rm si}\over h} \;\rm{ when } \;{\it P}_{\rm si}> 0.$$

At the thickness scales under consideration (i.e. 1–10 m, which means that the floating ice shelf we consider is initially MLSI rather than fully developed ice shelf) and with P _si in the range 10⁵–10⁶ Pa s, the value of ν _eff is in the range of 2.5 × 10¹¹–2.5 × 10¹⁵ Pa s (which is why in the previous section, we chose to illustrate typical behavior using a viscosity of 10¹³ Pa s).

We next ask whether the stress associated with ice-shelf flexure can be significant compared to P _si, which dominates the effective viscosity for the evolution of H in the absence of flexure. Here, we appeal to the equation relating bending moment and vertical forces, i.e. Eqn (7):

(36)$$\eqalign{& \rho_{\rm sw} g \eta + P_{\rm si}{\partial^2 \eta\over \partial x^2} \sim {\partial^2 M \over \partial x^2} \cr & \rho_{\rm sw} g \vert \eta\vert + P_{\rm si}{\vert \eta\vert \over \lambda^2} \sim {T_{\rm f} h^2\over \lambda^2}}$$

where M is the bending moment which is related to the stress induced by flexure T _f by,

(37)$$M = \int_{-h/2}^{h/2} T_{\rm f} \xi \ {\rm d}\xi \sim T_{\rm f} h^2$$

where ξ is the vertical distance from the neutral (center) plane of the ice shelf. The above considerations give:

(38)$$T_{\rm f} \sim {P_{\rm si} \vert \eta\vert \over h^2} + {\rho_{\rm sw}g \lambda^2 \vert \eta\vert \over h^2} \sim {P_{\rm si}\vert \eta\vert \over h^2} \;\rm{ when } \;{\it P}_{\rm si}> 0 .$$

Using the first of the above scale relationships, a sinusoidal roll profile of η with a 1 m amplitude, with a wavelength of 45 m (the value estimated in the previous section for P _si = 1 MPa m), gives T _f in the range of [10⁻³–1] MPa, for a range of P _si in [0.1–1] MPa m and a range of h in [1–10] m. Combining the estimate of T _s associated with changes in h with T _f gives:

(39)$$T = T_{\rm s} + T_{\rm f} = \left(1 + {\vert \eta( x,\; t) \vert \over h} \right){P_{\rm si}\over h} \;\rm{ when } \;{\it P}_{\rm si}> 0$$

and,

(40)$$T = T_{\rm s} + T_{\rm f} = \rho_{\rm i} g ( 1-\gamma) h + {\rho_{\rm sw}g \lambda^2 \vert \eta\vert \over h^2} \;\rm{ when }\;{\it P}_{\rm si} = 0.$$

We are now in a position to approximate the magnitude of an effective viscosity ν _eff based on the above scaling argument. Taking η = 1 m, h = 10 m and P _si = 1 MPa m, we find

(41)$$\nu_{\rm eff} \approx {h^2\over 2 A P_{\rm si}^2 \left(1 + {2\vert \eta( x,\; t) \vert }/{h} + {\eta( x,\; t) ^2}/{h^2} \right)} \ \approx 2.0 \times 10^{13}\, \rm{Pa\; s} \;\rm{ when } \;{\it P}_{\rm si}> 0$$

and,

(42)$$\nu_{\rm eff} \approx {1\over 2 A \left(\rho_{\rm i} g ( 1-\gamma) h + {\rho_{\rm sw}g \lambda^2 \vert \eta\vert }/{h^2}\right)^2} \ \approx 2.7 \times 10^{14} \, \rm{Pa \;s}\;\rm{ when }\; {\it P}_{\rm si} = 0 .$$

These estimates show that the effective viscosity during roll growth will be ~10 times smaller than that during roll decay. This means that rolls can persist on the surface of the ice shelf once the sea-ice pressure has been removed. If the rolls persist from winter until summer seasons, then perhaps surface hydrological processes may work to preserve them.

Discussion and conclusions

In an effort to identify a possible proximal origin of surface rolls on the Ellesmere Ice Shelf, we have examined one of several hypotheses: that high-magnitude (${\gtrsim }\! 1$ MPa m) compressive sea-ice pressure can induce periodic viscous flexure in an initially flat ice shelf of relatively stiff, brine-free ice (i.e. possibly MLSI). Although our analysis shows that the hypothesis is plausible, formidable difficulties must be contended with to fully accept the hypothesis. Plausibility comes from the result that a periodic flexure profile can grow with an e-folding timescale comparable to the time span of a winter season, when sea-ice pressure is likely to be highest, if the ice thickness is small enough (1–10 m). The greatest difficulty with the mechanism comes from the fact that rolls created by viscous flexure will decay once sea-ice pressure is removed unless some other process serves to maintain them. Specific predictions of our analysis are:

1. (1) rolls grow fastest in thin (1–10 m) ice shelves, with e-folding timescales that are sub-annual or seasonal (winter), and with sea-ice pressure exceeding ~1 MPa m;

2. (2) conversely, rolls in thick (50 m) ice shelves, comparable to present-day thickness of remnants of the Ellesmere Ice Shelf, grow so slowly that hundreds of years would be required to achieve observed amplitudes;

3. (3) the viscous response of ice thickness to compressive sea-ice pressure is uniform in space and, hence, is not initially involved in roll initiation;

4. (4) rolls initiate with a wavelength that is independent of ice thickness;

5. (5) in the finite displacement regime, a shortening ice shelf will yield a shortening evolved wavelength. Wavelength shortening will increase with roll amplitude and, therefore, decrease with ice thickness. The ice thickness variation is a result of slower roll amplitude growth rates in thicker ice (as shown by Eqn (28)).

To address the most formidable difficulty of roll origin by sea-ice pressure, namely that rolls decay when pressure is removed, we undertook a series of scale analysis and estimation exercises to examine how viscous flexure might respond to the non-Newtonian rheology of ice. These qualitative considerations led us to realize that a temporal asymmetry emerges; rolls will develop faster than they decay due to the impact of sea-ice pressure induced stresses on effective ice viscosity. Rolls created rapidly in winter, when sea-ice pressure is high, can thus possibly persist through the summer when sea-ice pressure is low and effective ice viscosity is high. Once the next winter comes and sea-ice pressure is restored, the rolls that have persisted from the previous winter can continue to grow. However, we speculate that other processes, such as spatially inhomogeneous surface ablation in the summer, may be more effective in both preserving and amplifying rolls generated by sea-ice pressure in the winter.

In conclusion, the enigma represented by the unusual rolling surface morphology of the Ellesmere Ice Shelf, both as reported by early explorers from the mid-19th century and visible today in remnants of the Ellesmere Ice Shelf, remains unsolved. The analysis presented here confirms that viscous flexure from sea-ice pressure can provide one of several possible origins for the rolls, but it raises additional doubts about its likelihood. Future research on the viscous flexure hypothesis, such as involving more extensive treatment of ice rheology (including influence of brine content and fracture damage) may quell some of these doubts. The most valuable future effort, however, would be to find an instance of roll origination to observe directly.