## Introduction

The once-contiguous Ellesmere Ice Shelf (Fig. 1) was first encountered by European Arctic explorers in 1876 (Aldrich, Reference Aldrich1877), 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^{2} 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^{2} in area (Mueller and others, Reference Mueller, Copland and Jeffries2017). 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, Reference Dowdeswell and Jeffries2017).

The origin of the Ellesmere Ice Shelf, deduced from study of its present-day remnants (Jeffries, Reference Jeffries2017), 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, Reference Jeffries2017). It had limited glacial sources on land, and fringed the northern coastlines of Ellesmere Island (Mueller and others, Reference Mueller, Copland and Jeffries2017), and possibly Greenland (Condron and others, Reference Condron, Joyce and Bradley2020). 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, Reference Moore, Schweiger, Zhang and Steele2019). 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, Reference Holdsworth1987) 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 (Reference Brau, Damman, Diamant and Witten2013)). 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, Reference Hattersley-Smith1957; Holdsworth, Reference Holdsworth1987; Jeffries, Reference Jeffries2017). 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, Reference Brau2011, Reference Brau, Damman, Diamant and Witten2013). 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.

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, Reference Banwell2017; Cassotto and others, Reference Cassotto, Burton, Amundson, Fahnestock and Truffer2021), 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 (Reference Holdsworth and Traetteberg1973) 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.

## 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, Reference Aldrich1877; Peary, Reference Peary1907). 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, Reference Nares and Feilden1878). About three decades later, Peary (Reference Peary1907) 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, Reference Crary1958), conducted in the late-1950s and early-1960s, also noted the unique rolling surface morphology (Hattersley-Smith, Reference Hattersley-Smith1957; Crary, Reference Crary1958, Reference Crary1960). 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, Reference England, Evans and Lakeman2017), 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, Reference Narod, Clarke and Prager1988; Mortimer and others, Reference Mortimer, Copland and Mueller2012), 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, Reference Jeffries2017).

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 (Reference Jeffries2017). 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 (Reference Jeffries2017) 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, Reference Narod, Clarke and Prager1988), despite the presence of surface rolls (Holdsworth, Reference Holdsworth1987). 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, Reference Jeffries2017) which has been cited as a potential clue as to the process or processes responsible for initiating these features (Holdsworth, Reference Holdsworth1987).

## 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 (Reference Hattersley-Smith1957), Holdsworth (Reference Holdsworth1987) and Jeffries (Reference Jeffries2017) 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 (Reference Hattersley-Smith1957) 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, Reference Crary1960), 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, Reference Narod, Clarke and Prager1988), and ‘sympathetic’ tides in meltwater ponds (which would be similar to the effect of consistently directed prevailing winds) (MacAyeal and others, Reference MacAyeal, Willis, Banwell, Macdonald and Goodsell2020).

A class of hypotheses, the one forming the subject of our study, was dismissed by Hattersley-Smith (Reference Hattersley-Smith1957), but has yet to be fully discredited. As described by Crary (Reference Crary1960) and Holdsworth (Reference Holdsworth1987), 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, Reference Collins and McCrae1985). 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 (Reference Sanderson1978) 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 (Reference MacAyeal2021), 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^{−3}) and sea water (1028 kg m^{−3}). 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:

and,

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)*ν* _{f}*H* ^{3}, 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, Reference MacAyeal2021). 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^{−1}) 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, Reference MacAyeal2021). The key term in this equation is the second on the left-hand side, *P* _{si}(∂^{2}η/∂*x* ^{2}), 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:

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, Reference MacAyeal2021).

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:

Substitution into Eqn (1) gives:

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:

Substitution into Eqn (2) gives:

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

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

and

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

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:

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:

which gives

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:

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):

which gives

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

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

which implies that

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:

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

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.

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$:

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

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^{13} 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, Reference Holdsworth and Traetteberg1973; Timco, Reference Timco2010) in the Arctic Ocean. Parameters used to re-dimension the variables are:

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^{5} 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, Reference Timco2010). 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, Reference MacAyeal and Holdsworth1986). From our analysis, it would take a reduction of *ν* by 3 orders of magnitude (from ~ 10^{13} to ~ 10^{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)*ν* _{f}*H* ^{3}, 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, Reference Jeffries2017) 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*)^{2} 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

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, Reference Jeffries2017).

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, Reference Brau, Damman, Diamant and Witten2013). 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^{13} 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, Reference MacAyeal and Holdsworth1986). 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,

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*:

where *A* ≈ 2 × 10^{−24} Pa^{−3} s^{−1} is a value representative of temperate ice (but with zero brine content), and *n* = 3 (Cuffey and Paterson, Reference Cuffey and Paterson2010). 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:

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^{5}–10^{6} Pa s, the value of *ν* _{eff} is in the range of 2.5 × 10^{11}–2.5 × 10^{15} Pa s (which is why in the previous section, we chose to illustrate typical behavior using a viscosity of 10^{13} 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):

where *M* is the bending moment which is related to the stress induced by flexure *T* _{f} by,

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

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^{−3}–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:

and,

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

and,

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) 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) 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) 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) rolls initiate with a wavelength that is independent of ice thickness;

(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.

## Acknowledgements

The research reported here was initiated by Coffey as part of an undergraduate research internship at the University of Chicago supported by the Jeff Metcalf Undergraduate Research Internship Fund and advised by MacAyeal. Support to MacAyeal and Banwell was provided by the US National Science Foundation (NSF-OPP) under grants 1841467 awarded to the University of Chicago, and 1841607 to the University of Colorado Boulder, respectively. This report was prepared by Sergienko under award NA18OAR4320123 from the National Oceanic and Atmospheric Administration, USDepartment of Commerce. The statements, findings, conclusions and recommendations are those of the author(s) and do not necessarily reflect the views of the National Oceanic and Atmospheric Administration, or the US Department of Commerce. Copland and Mueller acknowledge the Natural Sciences and Engineering Research Council of Canada, and ArcticNet Network of Centres of Excellence Canada, for funding. We thank Anthony Powell for allowing us to use his photograph displayed in Figure 5b. Photograph in Figure 5c is from the Antarctic Map and Photograph Library, US Geological Survey. We thank Gerald Holdsworth for both being an inspiration for this study and for helpful advice on the manuscript. We thank the scientific editor, Alan Rempel, and two anonymous referees for providing constructive advice, editing suggestions and for alerting us to a previously unrealized implication of our analysis.

## Appendix A: Derivation of the governing equations

Equations (1) and (2) are readily derived as 1-D simplifications of the treatment in MacAyeal and others (Reference MacAyeal2021). They can also be derived using a Lagrange/Hamilton approach much more easily using *η* and *ζ* as generalized variables (MacAyeal and Thomas, Reference MacAyeal and Thomas1982; Schoof, Reference Schoof2006; Bassis, Reference Bassis2010; Goldberg, Reference Goldberg2011). We present this derivation here for the benefit of interested readers.

We begin by describing the kinematic relationship between the variables *η* and *ζ* and the overall change in ice-plate geometry, particularly its horizontal span. This geometry is shown schematically in Figure 4. As a simplifying assumption, we take the ice front to be vertical throughout the deformation, and describe its position *x* _{if}(*t*) using the variable $\Xi ( t)$, which is the distance the ice-shelf span has contracted (shortened) in response to sea-ice pressure:

where *X* is the material coordinate describing horizontal distance along the plate at *t* = 0, θ(*X*) is the local slope of the neutral axis, and *H* = *h* + ζ (we use *H* and ζ interchangeably as they differ only by a constant),

We note that for θ ≪ 1,

In addition, for ζ ≪ 1 we may approximate ζ/*H* as ζ/*h*. With these approximations, the relation between $\Xi$ and the system's coordinates *η* and *ζ* is:

Deformation of the ice shelf involves the action of conservative forces (gravitational potential energy and work done by hydrostatic pressure and applied sea-ice pressure) and dissipative forces (deviatoric stress induced strain rate associated with creep and flow). The Lagrangian and its subsequently derived Euler–Lagrange equations that dictate how *η* and *ζ* evolve with time will be constructed in two parts, a traditional conservative part involving potential (and kinetic, but not applicable in this study by assumption) energy and a dissipative part involving non-conserved energy associated with work done against internal deviatoric stress as dictated by ice rheology. The addition of a dissipative function to the Euler–Lagrange equations follows the development of Rayleigh as described in Minguzzi (Reference Minguzzi2015). The conservative part of the Lagrangian involves the gravitational potential energy density function Δ*Q*,

and the work Δ*W* done to change the potential energy in the absence of dissipation:

In each of these functions, the Δ( · ) notation signifies that potential energy is relative to the undeformed initial state and that work is relative to the initial state where the impetus for deformation, compression by the sea ice, is absent, i.e. $\Xi = 0$. Rayleigh dissipation functions (Minguzzi, Reference Minguzzi2015), ${\cal D}_\eta$ and ${\cal D}_\zeta$, added to the Euler–Lagrange equations for *η* and *ζ*, respectively, are taken to be,

where the overdot $\dot {( ) }$ signifies the time derivative ∂/∂*t*. We note that when *V* = (1/3)ν_{f}*H* ^{3}, where ν_{f} is the ice viscosity associated with flexure, the dissipation for ice-plate flexure is identical to that derived in classical thin viscous plate theory (MacAyeal and others, Reference MacAyeal2021). Likewise, when $D = 8 \bar \nu _{\rm eff} H^{-1}$, where $\bar \nu _{\rm eff}$ is the vertical average of the effective viscosity *ν* _{eff}(*z*) (which is stress dependent when the flow law for ice is considered), the dissipation for ice-shelf spreading becomes identical to that for an unconfined floating ice tongue.

Setting up the Lagrangian ${\cal L}$, we determine a potential energy density function *F*:

where the prime notation denotes the partial derivative with respect to *x*. With the formulations for Δ*W* and Δ*Q* given above, the Euler–Lagrange equations associated with *η* and *ζ*, constructed from the Lagrangian and the dissipation functions (Minguzzi, Reference Minguzzi2015), are found to be:

and,

We note that in constructing the Euler–Lagrange equations, we use the form that is applicable when there are *x*-derivatives of *η* (none are involved for ζ = *H* − *h*) as outlined in Appendix B.

## Appendix B: Potential energy associated with *η*(*x*, *t*) ≠ 0 and *ζ*(*x*, *t*) ≠ 0

The Lagrangian involves the potential energy of conservative forces, which means that the Lagrangian does not depend on the path taken from the initial ice-plate geometry to its deformed geometry, it only depends on the initial and final states.

When the ice shelf is compressed, space beyond the initial position of the ice front (assumed vertical), between $x = L-\Xi$ and *x* = *L*, is replaced with sea water:

For the body of the deformed ice shelf, accounting for ice replacing air at the top and ice replacing water at the bottom:

We next note that:

This now allows simplification:

Thus, for the total change in potential energy Δ*Q*, we have:

We next consider the work Δ*W* done by the sea-ice pressure force and the hydrostatic pressure of sea water acting on the vertical ice front of the ice shelf. In addition, the hydrostatic pressure force in the ice shelf must be worked against. These forces work over the distance $\Xi$. Evaluating the net force acting on the ice front, and composing a work term for use in the Lagrangian, we have: