Kinematic assumptions

We begin by defining $x$ and $y$ as horizontal coordinates, $z$ as the vertical coordinate aligned with the acceleration of gravity $g$ at the Earth surface and zero at the undisturbed long-time sea level, and $t$ as time. The vertical dimensions of the ice shelf are described by surface and basal elevations $S( x,\; \, y,\; \, t)$ and $B( x,\; \, y,\; \, t)$, respectively, and ice thickness $H( x,\; \, y,\; \, t) = S - B$. At time $t = 0$, we assume that the ice shelf is in an undisturbed state with no flexure and is floating such a way as to satisfy Archimedes principle on a local basis (where each vertical ice column, no matter how small in horizontal surface area, displaces a weight of seawater equal to its own weight). We call the state of the ice shelf when floating in this way ‘flotational equilibrium’, which is when:

(1a)$$S( x,\; y,\; t = 0) = \left(1 - {\bar \rho_i\over \rho_{sw}} \right)H( x,\; y,\; t = 0) ,\; $$

(1b)$$B( x,\; y,\; t = 0) = - {\bar \rho_i\over \rho_{sw}} H( x,\; y,\; t = 0) ,\; $$

where $\bar \rho _i( x,\; \, y,\; \, t)$ is the vertical average of the ice density $\rho _i( x,\; \, y,\; \, z,\; \, t)$, and $\rho _{sw}$ is the density of seawater. Using the familiar averaging operator:

(2)$$\bar \cdot = {1\over H} \int_B^S \cdot \ {\rm d}z = {1\over H} \int_{{-H\over 2}}^{{H}/{2}} \cdot \ {\rm d}\zeta ,\; $$

where $\zeta$ is an alternative $z$-coordinate that is zero at the middle surface of the ice shelf, half way between the surface and base,

(3)$$\zeta = z - \left(B + {H\over 2} \right) = z - \left(S - {H\over 2} \right).$$

At $t = 0$, we assume that the initial ice thickness is $H( x,\; \, y,\; \, t = 0) = h( x,\; \, y)$ and that as time progresses, $H$ will be given by:

(4)$$H( x,\; y,\; t) = h( x,\; y) + {\cal H}_s( x,\; y,\; t) + {\cal H}_b( x,\; y,\; t) ,\; $$

where ${\cal H}_s$ and ${\cal H}_b$ are the cumulative amounts of ice thickness gained at the surface and the base of the ice shelf since $t = 0$. (A restatement of mass balance equations comes in the following section.) At times greater than zero, the surface and basal elevations are described by:

(5a)$$S( x,\; y,\; t> 0) = \left(1 - {\bar \rho_i\over \rho_{sw}} \right)h( x,\; y) + {\cal H}_s( x,\; y,\; t) + \eta,\; $$

(5b)$$B( x,\; y,\; t> 0) = - {\bar \rho_i\over \rho_{sw}} h( x,\; y) + {\cal H}_b + \eta ,\; $$

where $\eta$ is the cumulative change in the $z$ elevation of the ice shelf's middle surface (neutral surface) at $t = 0$ for $t> 0$ (when the addition of cumulative ice thicknesses ${\cal H}_s$ and ${\cal H}_b$ at the top and bottom, and the effect loads and bending moments begin to take effect). The ice shelf is defined to be flexed when, except for short term effects of waves and bobbing, $S$ and $B$ do not have the flotational equilibrium values given by Eqns (1a) and (1b). We further write $\eta$ as the sum of a viscous part $\eta _v$, an elastic part $\eta _e$ and an ice-shelf area mean part $\bar {\eta }^{x, y}$,

(6)$$\eta( x,\; y,\; t) = \eta_v( x,\; y,\; t) + \eta_e( x,\; y,\; t) + \bar{\eta}^{x, y}( t) + {\bar \eta}^\circlearrowleft( t) .$$

The terms $\bar {\eta }^{x, y}$ and ${\bar \eta }^\circlearrowleft$ are required to satisfy the constraint that the entire ice shelf is, on average, in flotational equilibrium:

(7)$$\bar{\eta}^{x, y} = \iint_A - { \rho_i\over \rho_{sw}} {{\cal H}}_s\, {\rm d}x\, {\rm d}y + \iint_A \left(1 - { \rho_i\over \rho_{sw}} \right){{\cal H}}_b\, {\rm d}x\, {\rm d}y ,\; $$

and to account for the possibility that surface and basal accumulation will rotate, or tilt, the ice shelf about a horizontal axis while keeping it in flotational equilibrium as a whole:

(8)$$\eqalign{\bar{\eta}^\circlearrowleft = \iint_A - { \rho_i\over \rho_{sw}} \chi^\circlearrowleft( {\cal H}_s) \, {\rm d}x\, {\rm d}y + \iint_A \left(1 - { \rho_i\over \rho_{sw}} \right)\chi^\circlearrowleft ( {\cal H}_b) \, {\rm d}x\, {\rm d}y ,\;} $$

where $A$ is taken to be the area of the ice shelf, and the operator $\chi ^\circlearrowleft ({\cdot})$ represents the averaging operator needed to resolve what parts of ${\cal H}_s$ and ${\cal H}_b$ satisfy the condition of perfect ice shelf tilting. The displacements $\bar {\eta }^{x, y}$ and ${\bar \eta }^\circlearrowleft$ are understood to induce no curvature to the middle surface ice shelf, and thus do not induce flexure stresses. This is why we consider spatially non-uniform surface and basal accumulation to be a process that generates flexure, which is embodied by $\eta _v$ and $\eta _e$.

Stress and strain rate

Having described our definitions of ice-shelf geometry and state of flotation both when flexure effects are present and when they are not, we turn to describing the velocity and strain rate of the ice shelf based on the common assumptions of shallow-shelf and thin-plate approximations described in previous studies (Love, Reference Love1888; MacAyeal, Reference MacAyeal1989). We refrain from describing displacement and strain associated with elastic deformation, but it is similar (MacAyeal and others, Reference MacAyeal, Sergienko and Banwell2015). The key assumption is that the square of the ice-shelf aspect ratio $\epsilon$, defined to be the ratio of typical ice thickness to the typical horizontal extent of the ice-shelf phenomena of interest, is much less than 1. This assumption allows terms in quantities that are ${\cal O}( \epsilon ^2)$ and smaller to be disregarded in formulating the governing equations.

We define $u$ and $v$ to be the $x$ and $y$ horizontal velocities, $w$ to be the vertical velocity, and $\dot {{\bf e}}$ to be the strain-rate tensor (bold faced variables are considered to be vectors or tensors). Considering one component of the horizontal velocity, $u$, we make the following decomposition based on our application of shallow-shelf and thin-plate assumptions:

(9)$$\eqalign{u( x,\; y,\; \zeta,\; t) & = \bar u( x,\; y,\; t) + \dot e_{xz}( x,\; y,\; t) \zeta + {\cal O}( \epsilon^2) + \cdots,\; \cr \hbox{total velocity} \ & = \ z\hbox{-} \ \hbox{and} \ \zeta\hbox{-independent shallow-shelf flow} \cr & + \hbox{linear-in-} \zeta \ \hbox{flexure effect} \cr & + \ \hbox{terms that are small},\; }$$

where $\dot e_{xz}$ is a vertical shear component of strain rate given by:

(10)$$\dot e_{xz} = {\partial u\over \partial z} + {\cal O}( \epsilon^2) + \cdots .$$

We note that the term ${\partial w}/{\partial x}$ does not appear in the above equation, because it is considered to be ${\cal O}( \epsilon ^2)$. To identify vertical shear with just the vertical gradient of the horizontal velocity is common in shallow-shelf, shallow-sheet and thin-plate applications.

Considering only the viscous component of ice-shelf flexure, the vertical strain rates are generated by the rate of change of the viscous component of $\eta$, which is $\eta _v$:

(11)$$\dot e_{xz} = - {\partial \dot \eta_v\over \partial x} ,\; $$

where $\dot \eta _v$ is the time-derivative of $\eta _v$. With this relationship, Eqn (9) may be written:

(12)$$u( x,\; y,\; \zeta,\; t) = \bar u( x,\; y,\; t) - {\partial \dot \eta_v\over \partial x} \zeta + {\cal O}( \epsilon^2) + \cdots .$$

We next consider the horizontal strain rate component $\dot e_{xx}$, which, given Eqn (11) is expressed as a sum:

(13)$$\dot e_{xx} = {\partial \bar u\over \partial x} - {\partial^2 \dot\eta_v\over \partial x^2} \zeta + {\cal O}( \epsilon^2) + \cdots .$$

All the other components of velocity and strain rate are partitioned as sums in the same manner. We note that the second term on the right-hand side of the above equation depends on $\zeta$. This means that this term's average over depth is zero, and it does not contribute directly (by integration over $z$) to the depth-averaged horizontal flow of the ice shelf.

The stress within the ice shelf is determined by gravity (giving pressure $p$ a glaciostatic term) and by deviatoric stress ${\bf T}$ that is related to the strain rates by the constitutive relation for ice. To illustrate just one deviatoric stress term, we consider $T_{xx}$:

(14)$$T_{xx} = 2 \nu {\partial \bar u\over \partial x} - 2 \nu {\partial^2 \dot\eta_v\over \partial x^2} \zeta .$$

where $\nu ( x,\; \, y,\; \, z,\; \, t)$ is a flow-law and ice-parameter (e.g. temperature, ice stiffness parameters) dependent effective viscosity, and we have dropped $< {\cal O}( \epsilon ^2)$-terms from the sum. If we take $\nu ( x,\; \, y,\; \, z)$ to be the sum of its mean over depth and an arbitrary extra function that has no mean (which accounts for variations of ice temperature and stiffness parameters with $z$):

(15)$$\nu( x,\; y,\; \zeta) = \bar \nu( x,\; y) + {\cal S}( \zeta) ,\; $$

where ${\cal S}( \zeta )$ (not to be confused with ice shelf surface elevation $S$) is a function with no mean (i.e. $\bar {{\cal S}}( \zeta ) = 0$). Substitution of this expression for $\nu$ into Eqn (14) gives:

(16)$$T_{xx} = 2 \left(\bar \nu + {\cal S}( \zeta) \right){\partial \bar u\over \partial x} - 2 \left(\bar \nu + {\cal S}( \zeta) \right){\partial^2 \dot\eta_v\over \partial x^2} \zeta .$$

For the depth-average part of the above stress component, we have

(17)$$\eqalign{\bar T_{xx} & = \overline{2 \left(\bar \nu + {\cal S}( \zeta) \right){\partial \bar u\over \partial x} } - \overline{ 2 \left(\bar \nu + {\cal S}( \zeta) \right){\partial^2 \dot\eta_v\over \partial x^2} \zeta } \cr & = 2 \bar\nu {\partial \bar u\over \partial x} - 2 \overline{\left(\zeta {\cal S} \right)} {\partial^2 \dot\eta_v\over \partial x^2} .}$$

In the above expression, we see two viscosity-like parameters, $\bar \nu$ and $\overline {\zeta {\cal S}( \zeta ) }$. Normally, in shallow-shelf flow without flexure, $\bar \nu$ is the effective viscosity that determines the flow. Here, we see as a first-element of coupling that an additional viscosity parameter is involved, namely $\overline {\zeta {\cal S}( \zeta ) }$.

We now consider the bending moment ${\bf M}$. We assume a Maxwell style of viscoelasticity in our treatment of flexure, and so both the elastic and the viscous parts of $\eta$ respond to the same single bending moment, like how a spring and a viscous attenuating piston/cylinder device (dashpot) respond to a force when connected in series. Here we consider how $\eta _v$, the viscous part of $\eta$, depends on $M_{xx}$:

(18)$$\eqalign{M_{xx} & = \int_{-( {H}/{2}) }^{{H}/{2}} \left(T_{xx} - p \right)\zeta\, {\rm d}z \cr & = \int_{-( {H}/{2}) }^{{H}/{2}} \left(2 \left(\bar \nu + {\cal S}( \zeta) \right){\partial \bar u\over \partial x} - 2 \left(\bar \nu + {\cal S}( \zeta) \right){\partial^2 \dot\eta_v\over \partial x^2} \zeta \right)\zeta\, {\rm d}z - \int_{-( {H}/{2}) }^{{H}/{2}} p( x,\; y,\; z,\; t) \zeta\, {\rm d}z \cr & = 2 H \overline{\left({\cal S} \zeta\right)} {\partial \bar u\over \partial x} - 2 {\partial^2 \dot\eta_v\over \partial x^2} \int_{-( {H}/{2}) }^{{H}/{2}} \left(\bar \nu + {\cal S}( \zeta) \right)\zeta^2\, {\rm d}z - \int_{-( {H}/{2}) }^{{H}/{2}} p( x,\; y,\; z,\; t) \zeta\, {\rm d}z \cr & = 2 H \overline{\left({\cal S} \zeta\right)} {\partial \bar u\over \partial x} - \left({ \bar \nu H^3\over 6} + 2 H \overline{\left({\cal S} \zeta^2 \right)} \right){\partial^2 \dot\eta_v\over \partial x^2} - \int_{-( {H}/{2}) }^{{H}/{2}} p( x,\; y,\; z,\; t) \zeta\, {\rm d}z .}$$

We note that the above expression for bending moment depends on three measures of ice viscosity, the depth average viscosity $\bar \nu$ and two measures $\overline {{\cal S}\zeta }$ and $\overline {{\cal S} \zeta ^2}$ that depend on the variation of viscosity with $\zeta$. We thus anticipate that the viscosity that primarily governs shallow-shelf flow, $\bar \nu$, will be different from the viscosity that governs flexure, although this is not the case when ${\cal S} = 0$, i.e. the ice is of constant viscosity.

If viscosity variation with $\zeta$ is presumed to be zero, a simpler expression for $M_{xx}$ emerges. We record this expression for use in describing the demonstration simulations in sections below:

(19)$$M_{xx} = - { \bar \nu H^3\over 6}{\partial^2 \dot\eta_v\over \partial x^2} - \int_{-( {H}/{2}) }^{{H}/{2}} p( x,\; y,\; z,\; t) \zeta\, {\rm d}z .$$

To proceed further with deriving the bending moment $M_{xx}$ we need to determine the pressure $p$. We start with the vertical stress balance:

(20)$${\partial\over \partial x}\left(T_{xz} \right) + {\partial\over \partial y}\left(T_{yz}\right) + {\partial \over \partial \zeta} \left(T_{zz} \right)- {\partial p\over \partial \zeta} - \rho_i( \zeta) g = 0 .$$

As is customary in shallow-shelf and thin-plate approximations we drop the vertical shear terms (assuming they are ${\cal O}( \epsilon ^2)$ or smaller, and substitute for $T_{zz} = 2 \nu ( {\partial w}/{\partial \zeta })$:

(21)$${\partial \over \partial \zeta} \left(2 \nu {\partial w\over \partial \zeta}\right)- {\partial p\over \partial \zeta} - \rho_i g = 0 .$$

We next invoke the incompressibility condition and the definitions of $\dot e_{xz}$ and $\dot e_{yz}$ in terms of the $x$ and $y$ derivatives of $\eta$,

(22)$$\eqalign{{\partial w\over \partial \zeta} & = - \left({\partial u\over \partial x} + {\partial v\over \partial y} \right)\cr & = -\left({\partial \bar u\over \partial x} + {\partial \bar v\over \partial y} \right)- \left({\partial \dot e_{xz}\over \partial x} + {\partial \dot e_{yz}\over \partial y} \right)\zeta \cr & = -\left({\partial \bar u\over \partial x} + {\partial \bar v\over \partial y} \right) + \left({\partial^2 \dot \eta_v\over \partial x^2} + {\partial^2 \dot \eta_v\over \partial y^2} \right)\zeta ,\; }$$

to give:

(23)$$ \eqalign{& - {\partial \over \partial \zeta} \left(2 \nu \left({\partial \bar u\over \partial x} + {\partial \bar v\over \partial y} \right)\right) + {\partial\over \partial \zeta} \left(2 \nu \zeta \left({\partial^2 \dot \eta_v\over \partial x^2} + {\partial^2 \dot \eta_v\over \partial y^2} \right)\right)\cr& - {\partial p \over \partial \zeta} - \rho_i g = 0 .}$$

We define the density, like the viscosity, to be the sum of its depth-average and a term that varies with $\zeta$:

(24)$$\rho_i( x,\; y,\; \zeta,\; t) = \bar \rho_i( x,\; y,\; t) + {\cal R}( x,\; y,\; \zeta,\; t) ,\; $$

where $\bar {{\cal R}} = 0$. Integrating Eqn (23) over the vertical, from a lower limit of $\zeta$ to an upper limit of $\displaystyle {H}/{2}$ where we apply the boundary condition that atmospheric pressure is zero, we obtain:

(25)$$\eqalign{& - 2 \nu \left({\partial \bar u\over \partial x} + {\partial \bar v\over \partial y} \right) + 2 \nu \zeta \left({\partial^2 \dot \eta_v\over \partial x^2} + {\partial^2 \dot \eta_v\over \partial y^2} \right)- p - \bar \rho_i g \left(\zeta - {H\over 2} \right) \cr& - g \int_{\zeta}^{{H}/{2}} {\cal R}( \zeta) \, {\rm d}\zeta = 0 ,\;} $$

which gives:

(26)$$\eqalign{\,p( \zeta) =\; & 2 \nu \left({\partial \bar u\over \partial x} + {\partial \bar v\over \partial y} \right)\cr & + 2 \nu \zeta \left({\partial^2 \dot \eta_v\over \partial x^2} + {\partial^2 \dot \eta_v\over \partial y^2} \right)\cr & - \bar \rho_i g \left(\zeta - {H\over 2} \right)\int_{\zeta}^{{H}/{2}} {\cal R}( \zeta) \, {\rm d}\zeta .}$$

We record for future use the depth average of $p$:

(27)$$\eqalign{\bar p =\;& 2 \bar \nu \left({\partial \bar u\over \partial x} + {\partial \bar v\over \partial y} \right)\cr & + 2 \overline{ \left(\zeta S( \zeta) \right)}\left({\partial^2 \dot \eta_v\over \partial x^2} + {\partial^2 \dot \eta_v\over \partial y^2} \right)\cr & + \bar \rho_i g {H\over 2} \cr & - g \overline{ \left(\int_{\zeta}^{{H}/{2}} {\cal R}( \zeta) \, {\rm d}\zeta \right)} .}$$

For future reference, if the depth variation of viscosity and ice density are zero, i.e. ${\cal S} = 0$ and ${\cal R} = 0$ we obtain:

(28)$$\eqalign{\,p( \zeta) & = - 2 \bar \nu \left({\partial \bar u\over \partial x} + {\partial \bar v\over \partial y} \right)\cr & + 2 \bar \nu \zeta \left({\partial^2 \dot \eta_v\over \partial x^2} + {\partial^2 \dot \eta_v\over \partial y^2} \right)\cr & - \bar \rho_i g \left({H\over 2} - \zeta \right),\; }$$

which gives for the depth-average pressure:

(29)$$\eqalign{\bar p =\;& 2 \bar \nu \left({\partial \bar u\over \partial x} + {\partial \bar v\over \partial y} \right)\cr & + \bar \rho_i g {H\over 2} .}$$

We can now compute the contribution of the pressure to the bending moment using the simplified form in Eqn (28):

(30)$$- \int_{-( {H}/{2}) }^{{H}/{2}} p( x,\; y,\; \zeta,\; t) \zeta\, {\rm d}z = {\bar \nu H^3\over 6} \left({\partial^2 \dot \eta_v\over \partial x^2} + {\partial^2 \dot \eta_v\over \partial y^2} \right) + { \bar\rho_i g H^3\over 12} .$$

This then renders the expression obtained previously (Eqn (19)) into the following form:

(31)$$\eqalign{M_{xx} & = - { \bar \nu H^3\over 6}{\partial^2 \dot\eta_v\over \partial x^2} - {\bar \nu H^3\over 6} \left({\partial^2 \dot \eta_v\over \partial x^2} + {\partial^2 \dot \eta_v\over \partial y^2} \right) + { \bar\rho_i g H^3\over 12} \cr & = - { \bar \nu H^3\over 3}{\partial^2 \dot\eta_v\over \partial x^2} - {\bar \nu H^3\over 6} {\partial^2 \dot \eta_v\over \partial y^2} + { \bar\rho_i g H^3\over 12} .}$$

For what follows, we note that components of ${\bf M} = [ M_{xx} \ M_{yy} \ M_{xy} ] ^T$ are related to the second derivatives of $\dot \eta _v$ by the following linear relation when we disregard $\zeta$ variations in ice viscosity and density, i.e. when ${\cal S}( \zeta ) = 0$ and ${\cal R}( \zeta ) = 0$, respectively (as we shall assume in our demonstration simulations in sections below):

(32)$${\bf M} = {\bf M}_v = - {\bar \nu H^3\over 3} \left[\matrix{ 1 & {1\over 2} & 0 \cr {1\over 2} & 1 & 0 \cr 0 & 0 & {1\over 2} }\right]\cdot \left[\displaystyle\matrix{\displaystyle {\partial^2 \dot{\eta}_v\over \partial x^2} \cr \displaystyle {\partial^2 \dot {\eta}_v\over \partial y^2} \cr \displaystyle {\partial^2 \dot {\eta}_v\over \partial x\partial y} }\right] = {\bf V} \cdot \left[\displaystyle\matrix{\displaystyle {\partial^2 \dot{\eta}_v\over \partial x^2} \cr \displaystyle {\partial^2 \dot {\eta}_v\over \partial y^2} \cr \displaystyle {\partial^2 \dot {\eta}_v\over \partial x\partial y} }\right],\; $$

where we have added the subscript $v$ in the above expression to differentiate the bending moment due to viscous stresses from those due to glaciostatic (subscript $g$) and elastic (subscript $e$). The glaciostatic contribution to the bending moment expressed as

(33)$${\bf M}_g = { \bar \rho_i g H^3\over 12} \left[\matrix{ 1 \cr 1 \cr 0 }\right],\; $$

will be treated as a forcing term to the vertical force balance equation (derived in the section presenting our treatment of flexure below), because it does not depend on $\eta$. We note that the above expression for $M_g$ is not written with terms that depend on $\overline {( \zeta {\cal S}) }$ or ${\cal R}( \zeta )$, however, these terms could be very important in situations where ice shelves have voids or englacial water that disturbs the viscosity and density structure with depth. Such a situation is envisioned for Arctic ice shelves such as the Milne Ice Shelf, where channelized englacial water bodies and voids are detected, and which may contribute to the pervasive rumpling of the ice shelf (Rajewicz, Reference Rajewicz2017).

To avoid repeating the similar, but lengthy, derivation of the elastic contribution for the bending moment ${\bf M}_e$, we simply record it here:

(34)$${\bf M}_e = - {E H^3\over 12\left(1 - \mu^2\right)} \left[\matrix{ 1 & \mu & 0 \cr \mu & 1 & 0 \cr 0 & 0 & \left(1 - \mu \right)}\right]\cdot \left[\displaystyle\matrix{\displaystyle {\partial^2 {\eta}_e\over \partial x^2} \cr \displaystyle {\partial^2 {\eta}_e\over \partial y^2} \cr \displaystyle {\partial^2 {\eta}_e\over \partial x\partial y} }\right] = {\bf D} \cdot \left[\displaystyle\matrix{\displaystyle {\partial^2 {\eta}_e\over \partial x^2} \cr \displaystyle {\partial^2 {\eta}_e\over \partial y^2} \cr \displaystyle {\partial^2 {\eta}_e\over \partial x\partial y} }\right],\; $$

where $E$ is the Young's modulus and $\mu$ is Poisson's ratio, which we assume to be constants. (We do not consider their variation with depth in this study, however, such consideration, if warranted, would follow the above treatment for how ice viscosity and density vary with depth.)

The linear operator expressions are given above for ${\bf V}$ and ${\bf D}$ will be used in the section below on developing the coupled treatment of ice-shelf flow and flexure. They will also be used in the demonstration simulations presented below. We reiterate, however, that the linear operators given above have been simplified under the assumption that contributions from ${\cal S}$ and ${\cal R}$, representing vertical variations in viscosity and density, are zero. Full consideration of variable viscosity and density is relegated to applications that go beyond simple demonstrations.

Ice rheology in the coupled system

As indicated above, the depth-averaged stress appearing in the SSA depends on the depth-averaged effective ice viscosity, and the bending moment depends on both the depth average of, and the first moment of, the effective ice viscosity. This difference can be understood by considering the situation where an anomalously viscous ice layer (e.g. because of temperature or ice impurity) exists at the ice-shelf middle surface. This viscosity anomaly will contribute to a depth average, but not to the first moment. In this study, we do not seek to clarify or demonstrate the effects of ice rheology, e.g, specific variation of effective viscosity with $\zeta$ caused by the flow law and by other effects such as temperature. We thus proceed in the demonstration simulations below to use two different Newtonian viscosities for ice, one that is the depth average $\bar \nu$ relating to the SSA and one that includes an arbitrary (for our demonstrations) additional viscosity variation used for the thin-plate approximation which we call $\nu _f$. In the case of the ice-rumple demonstration presented in the following section, we find that it is convenient to have $\nu _f < \nu$ for the ad hoc way in which we set up the simulation. This could, for example, be explained by ice near the upper and lower surfaces of the ice shelf having lower viscosity than in the centre (i.e. near the middle surface) where the ice is colder.

Pre-stressing of newly accumulated ice layers

In the foregoing derivation, we have focussed on computing bending moments via the integral:

(35)$$\eqalign{& M_{xx} \propto \int_{-( {H}/{2}) }^{{H}/{2}} \cdot \ {\rm d}\zeta =\cr & \int_{-( {h}/{2}) }^{{h}/{2}} \cdot \ {\rm d}\zeta + \int_{{h}/{2}}^{{\cal H}_s} \cdot \ {\rm d}\zeta + \int_{-{\cal H}_b}^{-{( h}/{2}) } \cdot \ {\rm d}\zeta .}$$

Let us consider what happens when, at $t = 0$, the ice thickness is changed from $h$ to $h + {\cal H}_s$ through the instantaneous addition of an accumulated layer of ice ${\cal H}_s$. Let us further suppose that at $t = 0$ the ice shelf is in a state of flexure so that $\eta \ne 0$ and bending moments exist within it. The new addition of ice ${\cal H}_s$ is added by a process (e.g. snow fall and densification) that does not introduce intrinsic bending moments. This means that the original part of the ice shelf contained within the thickness $h$ is pre-stressed, and the new part of the ice shelf contained within the thickness ${\cal H}_s$ is stress-free. As $t> 0$ this disparity between the old and new ice will mean that for the old ice to deform viscoelastically toward its equilibrium shape, it must deform the new ice from its deposited, stress-free shape. This is a bit like pre-stressing in concrete or metallurgy. The proper way to account for this is to track each layer of the ice shelf as it is added. And this can be done via a full Stokes treatment. Since our intention is to reduce complexity, we deal with this subtle problem by assuming that any new ice added to the ice shelf instantaneously inherits the stress state of the ice shelf that it joins.

Vertical force balance due to bending moments

The vertical force balance on ice columns in the ice shelf is used to determine ${{D \eta _v}/{D_t}}$ and ${{D \eta _e}/{D t}}$. In this balance, the divergence of bending moment gradients is equal to the sum of vertical force acting on the ice shelf (Love, Reference Love1888; Slim and others, Reference Slim, Teichman and Mahadevan2012):

(46)$$\rho_i g H {\partial^2 \eta\over \partial t^2} -\left({\partial^2 M_{xx}\over \partial x^2} + 2 {\partial^2 M_{xy}\over \partial x \partial y} + {\partial^2 M_{yy}\over \partial y^2} \right) = f ,\; $$

where ${\bf M} = [ M_{xx} \ M_{yy} \ M_{xy} ] ^T$ is a column vector of bending moments (the order of appearance of its elements has no physical meaning, but is defined for convenience), and $f$ is the sum of forces acting vertically on the plate. We shall further adopt the convention that terms depending on the glaciostatic part of the bending moments ${\bf M}_g$ will be written as a forcing term in the expression for $f$, to be discussed below. This is because ${\bf M}_g$ does not depend on $\eta$. Our first assumption is that the flexure of an ice-shelf is quasistatic, we thus drop the acceleration term on the left-hand side of Eqn (46).

Referring to Eqns (32) and (34), we write:

(47a)$$\eqalign{& {\bf D}^{-1} {\bf M} = \left[\matrix{ \displaystyle {\partial^2 \eta_e\over \partial x^2} \cr \displaystyle{\partial^2 \eta_e\over \partial y^2} \cr \displaystyle{\partial^2 \eta_e\over \partial x\partial y} }\right],\; }$$

(47b)$$\eqalign{& {\bf V}^{-1} {\bf M} = \left[\displaystyle\matrix{\displaystyle {\partial^2 \dot {\eta}_v\over \partial x^2}\cr \displaystyle {\partial^2 \dot {\eta}_v\over \partial y^2} \cr \displaystyle {\partial^2 \dot {\eta}_v\over \partial x\partial y} }\right].}$$

By the definition $\eta = \eta _v + \eta _e$ (and assuming that $\bar {\eta }^{x, y} = \bar \eta ^\circlearrowleft = 0$), we can write:

(48)$$\eqalign{\left[\matrix{ \displaystyle {\partial^2 \dot {\eta}\over \partial x^2} \cr \displaystyle{\partial^2 \dot {\eta}\over \partial y^2} \cr \displaystyle{\partial^2 \dot {\eta}\over \partial x\partial y} }\right]& = \left[\matrix{ \displaystyle {\partial^2 \dot {\eta}_e\over \partial x^2} \cr \displaystyle{\partial^2 \dot {\eta}_e\over \partial y^2} \cr \displaystyle{\partial^2 \dot {\eta}_e\over \partial x\partial y} }\right] + \left[\displaystyle\matrix{\displaystyle {\partial^2 \dot {\eta}_v\over \partial x^2} \cr \displaystyle {\partial^2 \dot {\eta}_v\over \partial y^2} \cr \displaystyle {\partial^2 \dot {\eta}_v\over \partial x\partial y} }\right].}$$

Substitution of the expressions in Eqns (47a) and (47b) into the above equation gives:

(49)$$\left[\displaystyle\matrix{\displaystyle {\partial^2 \dot {\eta}\over \partial x^2} \cr \displaystyle {\partial^2 \dot {\eta}\over \partial y^2} \cr \displaystyle {\partial^2 \dot {\eta}\over \partial x\partial y} }\right] = {\bf D}^{-1} \dot {{\bf M}} + {\bf V}^{-1} {\bf M} .$$

This equation forms the constitutive relation used in our formulation.

In other examples of the application of the viscoelastic constitutive relation in a linear form, such as by Rosier and others (Reference Rosier and Gudmundsson2018), the constitutive relation appears as:

(50)$$\dot e_{ij} = c_1 \overset{\triangle}{\tau}_{ij} + c_2 \tau_{ij} ,\; $$

where $c_1$ and $c_2$ are constants, and the overset triangle denotes the upper convective time derivative (Oldroyd time derivative), having as its leading term the material time derivative indicated by the overdot. We do not use the upper convective time derivative, because its extra terms account for large-amplitude flexure which rotate the neutral surface of the viscoelastic plate significantly out of horizontal. We anticipate that ice-shelf flexure will be of low amplitude and that we may assume that components of stress, strain and strain-rate tensors with indices $xx$, $yy$ and $xy$ can still represent these tensor components on the neutral surface, which will only be slightly tilted out of horizontal.

One final modification is made to the constitutive relation to avoid mixing time derivatives of ${\bf M}$ with undifferentiated ${\bf M}$ in the same equation. We introduce a new variable of convenience, which we shall refer to as a ‘generating potential’, $\mathbf \Phi = [ \Phi _1 \ \Phi _2 \ \Phi _3] ^{{\bf T}}$:

(51)$$\displaystyle \dot {\mathbf \Phi} = {\partial \mathbf \Phi\over \partial t} + {\bf u} \cdot \nabla \mathbf \Phi = {\bf V}^{-1} {\bf M} .$$

This allows Eqn (49) to be written in a manner where all terms are time differentiated:

(52)$$\left[\displaystyle\matrix{\displaystyle {\partial^2 \dot {\eta}\over \partial x^2} \cr \displaystyle {\partial^2 \dot {\eta}\over \partial y^2} \cr \displaystyle {\partial^2 \dot {\eta}\over \partial x\partial y} }\right] = {\bf D}^{-1} \dot {{\bf M}} + \dot {\mathbf \Phi} ,\; $$

or

(53)$$\left({\partial \over \partial t} + {\bf u} \cdot \nabla \right)\left({\mathbf \Phi} + {\bf D}^{-1} {{\bf M}} - \displaystyle\left[\matrix{\displaystyle {\partial^2 {\eta}\over \partial x^2} \cr \displaystyle {\partial^2 {\eta}\over \partial y^2} \cr \displaystyle {\partial^2 {\eta}\over \partial x\partial y} }\right]\right) = 0 .$$

Integration over $t$ gives,

(54)$${\mathbf \Phi} + {\bf D}^{-1} {{\bf M}} - \displaystyle \left[\matrix{\displaystyle {\partial^2 {\eta}\over \partial x^2} \cr \displaystyle {\partial^2 {\eta}\over \partial y^2} \cr \displaystyle {\partial^2 {\eta}\over \partial x\partial y} }\right] = 0 ,\; $$

where we assume the constant of integration to be zero. This assumption is equivalent to assuming that at some initial time (the lower bound of the integration at $t = 0$), there are no viscous contributions to bending moment.

We summarize the full formulation of flexure based on the many assumptions and on our novel approach to working with bending moment and $\eta$ as the principle variables:

(55a)$$-\left({\partial^2 M_{xx}\over \partial x^2} + 2 {\partial^2 M_{xy}\over \partial x \partial y} + {\partial^2 M_{yy}\over \partial y^2} \right) = f ,\; $$

(55b)$${\partial \mathbf \Phi\over \partial t} + {\bf u} \cdot \nabla \mathbf \Phi - {\bf V}^{-1} {\bf M} = \mathbf 0 ,\; $$

(55c)$$\mathbf \Phi + {\bf D}^{-1} {\bf M} - {\bf H} = \mathbf 0 ,\; $$

where ${\bf H}$ is defined as

(56)$${\bf H} = \left[{\partial^2 \eta\over \partial x^2} \ {\partial^2 \eta\over \partial y^2} \ {\partial^2 \eta\over \partial x \partial y} \right]^{{\bf T}} .$$

Normally, in derivations such as those presented here, the effort is to replace the bending moments with terms involving second derivatives of $\eta$. This then yields a single equation that is fourth order in spatial derivatives of $\eta$ when ${\bf u}$ is zero; and fifth order if ${\bf u} \ne 0$. We prefer to work with the system of three equations written above because its largest spatial derivatives are only second order, and this conforms well with the capabilities of Comsol™, which requires all equations of higher than second order spatial derivatives be divided into systems of equations having lower-order spatial derivatives.

We note that our formulation is intended for both short-timescale (seconds to days) and long-timescale processes (years). An argument can be made that, in the case of processes with long timescales, it is unnecessary to retain the elastic part of the viscoelastic constitutive relation we propose here. This may be true. We note, however, that treatment of both elastic and viscous behavior introduces only additional coefficients to the governing partial differential equations (PDEs) expressed above, not new variables; so elastic treatment comes without significant computational cost even if the phenomenon of interest is dominated by the viscous part of flexure. Also, in some applications, such as the differential-ablation phenomena we will discuss in the demonstration simulations below, elastic deformation is introduced to the ice shelf continuously with time (e.g. with each day's spatially variable melting and water movement); and thus elastic effects cannot be so easily separated from viscous effects based on their short-term importance even when the phenomena evolve over many cycles of the viscoelastic relaxation time, commonly referred to as the Maxwell timescale. Finally, the initial elastic response may determine length scales of features that subsequently evolve through viscous deformation. The elastic response, thus, may provide an initial condition for the time-dependent viscous flexure which would otherwise develop differently in the absence of elasticity.

Flexure forcing

The physical meaning of Eqn (55a) is the balance of bending moments and forces $f( x,\; y,\; t)$ caused by a number of physical processes. We represent these processes in the following sum:

(57a)$$f = - \rho_{sw} g \left(\eta_e + \eta_v\right)$$

(57b)$$ + \left(\rho_{sw} - \rho_i \right)g \left({\cal H}_b - \bar{{\cal H}}_b^{x, y} - \bar{{\cal H}}_b^\circlearrowleft \right)$$

(57c)$$- \rho_i g \left({\cal H}_s - \bar{{\cal H}}_s^{x, y} - \bar{{\cal H}}_s^\circlearrowleft \right)$$

(57d)$$- \rho_w g d$$

(57e)$$ + \left(\bar T_{xx} - \bar P_g \right){\partial^2 \eta\over \partial x^2} + 2\bar T_{xy}{\partial^2 \eta\over \partial x \partial y} + \left(\bar T_{yy} - \bar P_g \right){\partial^2 \eta\over \partial y^2}$$

(57f)$$ + {\bar \rho_i g\over 12} \left({\partial^2 H^3\over \partial x^2} + {\partial^2 H^3\over \partial y^2} \right)$$

(57g)$$ + {\cal T}( \theta) $$

(57h)$$ + {\cal A}( p_a) $$

(57i)$$ + O( h_d) .$$

We describe each term in the sum on the right-hand side in order from a to i in Eqns (57a)–(57i):

1. a. The restoring force of buoyancy acting on the initial ice thickness $h$.

2. b. The buoyancy of accumulated mass between the base and the flotational equilibrium value of the basal elevation. This depends on spatial gradients in ${\cal H}_s$ and ${\cal H}_b$ and spatially inhomogeneous basal accumulation. This term arises, for example, in the creation of basal channels (Le Brocq and others, Reference Le Brocp2013; Alley and others, Reference Alley, Scambos, Siegfried and Fricker2016) and inhomogeneous marine ice deposition (Koch and others, Reference Koch, Fitzsimons, Samyn and Tison2015).

3. c. The weight of accumulated mass between the surface and the flotational equilibrium value of the surface elevation. This term arises, for example, in the creation of surface channels (Dow and others, Reference Dow2018, Boghosian, personal communication, 2020).

4. d. The effect of additions or subtractions of surface meltwater (where $\rho _w$ is the density of fresh meltwater). This term arises, for example, where there are filling and draining supraglacial lakes (Banwell and MacAyeal, Reference Banwell and MacAyeal2015; MacAyeal and others, Reference MacAyeal2015; Banwell and others, Reference Banwell, Willis, Macdonald, Goodsell and MacAyeal2019; Robel and Banwell, Reference Robel and Banwell2019) and dolines (Moore, Reference Moore1993; Bindschadler and others, Reference Bindschadler, Scambos, Rott, Skvarca and Vornberger2002).

5. e. The effect of the membrane (horizontal ice-shelf flow) stresses associated with ice-shelf flow having a vertical component when there is acurvature in the middle surface of the ice-shelf (Ribe, Reference Ribe2003; Slim and others, Reference Slim, Teichman and Mahadevan2012). This term arises when there is a second derivative of $\eta$, because a ‘horizontal section’ of the ice shelf which has vertical sides when unflexed will have both its sides tilt slightly ‘face down’ or ‘face up’ when curved. This means that the membrane stress will have a slight vertical coordinate. This term is expected to be important when there are strong compressive membrane stresses in the ice shelf, such as at ice rumples (Matsuoka and others, Reference Matsuoka2015; De Rydt and others, Reference De Rydt, Gudmundsson, Nagler and Wuite2019) and traveling lakes (LaBarbera and MacAyeal, Reference LaBarbera and MacAyeal2011).

6. f. The effect of glaciostatic pressure generated bending moments (derived in a previous section). This term is expected to be important where there are large second derivatives of ice thickness and significant variations of ice density with depth, such as in the case of the Milne Ice Shelf (Rajewicz, Reference Rajewicz2017).

7. g. The effect of thermal stresses ($\theta$ is temperature). This can be important where ice is thin and there is a significant annual or diurnal cycle in surface temperature (Sanderson, Reference Sanderson1978; Bažant, Reference Bažant1992; MacAyeal and others, Reference MacAyeal2015, Reference MacAyeal2017).

8. h. The effects of atmospheric pressure gradients $p_a( x,\; \, y,\; \, t)$ which bear down on the upper surface of the ice shelf such as the inverse barometer effect (Ross, Reference Ross1854).

9. i. The effects of ocean dynamic topography $h_d( x,\; \, y,\; \, t)$ changes due to tides (Rosier and others, Reference Rosier and Gudmundsson2018) and sea swell (Sergienko, Reference Sergienko2010, Reference Sergienko2017).

Flexure is also forced by boundary conditions, such as at grounding lines (Schoof, Reference Schoof2011; Glasser and Gudmundsson, Reference Glasser and Gudmundsson2012; Rosier and others, Reference Rosier and Gudmundsson2018) and ice fronts (Reeh, Reference Reeh1968; Hindmarsh, Reference Hindmarsh2012; Christmann and others, Reference Christmann, Müller and Humbert2019); however, we do not focus on this well known forcing in this study.

In the present study, our focus is on the first, third and fifth terms of the right-hand side of Eqns (57a)–(57i). We thus relegate to future research determination of how relevant the other forcing terms are to ice-shelf behavior, and to what processes they pertain. The second term of Eqns (57a)–(57i) controls the response to accumulation and ablation, and especially their gradients in situations where ablation causes surface channels to incise into the ice shelf, and where debris cover (e.g. the McMIS, Antarctica), and its absence in limited areas, cause pedestalling of high-albedo surfaces (e.g. pedestalled, relict lakes, Macdonald and others, Reference Macdonald2019). The fourth forcing term is of interest when the ice-shelf flow stresses are compressive and give rise to viscous buckling instability that leads to ice rumples along grounding lines where ice flows toward, rather than away from, the grounding line.

To solve the above formulation of viscoelastic flexure, initial and boundary conditions are required. Thus we choose simple, bending moment free boundary conditions for $M$ and $\eta$. For the initial conditions, we will take $\eta$ to be zero or some other function, as the case examined requires. The details of these specifications will depend on the application, as we discuss below for our demonstration simulations.

Differential ablation: pedestalled, relict lakes

In our first idealized demonstration, we investigate the effects of spatially heterogeneous surface ablation in the context of pedestalled, relict lakes we have observed on the McMIS, Antarctica (Macdonald and others, Reference Macdonald2019; Banwell and others, Reference Banwell, Willis, Macdonald, Goodsell and MacAyeal2019) (Fig. 1). Pedestalled, relict lakes are found on debris-covered ice shelves, and each begins as a lake that forms in a debris-covered depression during one melt season. During the first winter, the lake water freezes, becoming superimposed ice that conceals the original debris cover with high-albedo ice. In subsequent years, the superimposed high-albedo ice reduces surface ablation relative to the surrounding debris-covered ice-shelf surface. This causes the surface of the surrounding ice shelf to lose mass at a rate much greater than the spatially limited zone of the original lake, now covered by superimposed ice and free of surface debris. The sharp contrast between the ablation rate at the high-albedo original lake site and the low-albedo debris-covered surroundings (Fig. 1b) leads to circumstances where the second term, ${\cal H}_s$, in the forcing for flexure has a strong spatial gradient at the shoreline of the original lake. This is what is investigated in the second demonstration.

Fig. 1. (a) Aerial photograph (altitude $\sim 500$ m) of Ring Pedestal on the McMIS looking from North to South (Macdonald and others, Reference Macdonald2019) (photo taken in January 2016 by A. Banwell). (b) Photograph at the edge of Peanut Pedestal looking from North to South (photo taken in January 2016 by A. Banwell). The height of the pedestal at its edge is approximately 3 m, however is variable depending on local debris conditions. (c)-(e) True color pan-sharpened Landsat 7 and 8 images showing the initial lakes (c), and development of pedestalled features (d) and (e) over a 7-year period (Macdonald and others, Reference Macdonald2019).

Our demonstration simulation will consider only a cross-section of the ice shelf (along a flowline) containing a relict lake feature. The 1-D version of the ice-shelf flow, mass balance and flexure equations derived above are:

(66a)$${\partial \over \partial x} \left(4 \bar \nu H {\partial u\over \partial x} \right) = -\rho_{sw}g B {\partial B\over \partial x} - \rho_i g H {\partial H\over \partial x} ,\; $$

(66b)$$- {\partial ^2 M_{xx}\over \partial x^2} = - \rho_{sw}g \eta -\rho_i g {\cal H}_s + H\left(4 \nu {\partial u\over \partial x} - {1\over 2} \rho_i g H \right){\partial^2 \eta\over \partial x^2} ,\; $$

(66c)$${\partial \Phi\over \partial t} + u {\partial \Phi\over \partial x} + {3\over \nu_f H^3} M_{xx} = 0 ,\; $$

(66d)$$\Phi - {12 \left(1 - \mu^2\right)\over E H^3} M_{xx} - {\partial^2 \eta\over \partial x^2} = 0 ,\; $$

(66e)$${\partial H\over \partial t} + u {\partial H\over \partial x} = - H {\partial u\over \partial x} .$$

We use a 1-D domain that is 14 km long (motivated by the flowline length of the McMIS depicted in the field study map of Macdonald and others (Reference Macdonald2019)), and specify an inflow boundary condition of $u = 50$ m a$^{-1}$ at the upstream side of the domain. Initial ice thickness is taken to be 50 m. For the ice viscosity, we take $\bar \nu = 10^{16}$ Pa s and $\nu _f = {\bar \nu }/{2}$. This is an arbitrary choice, as explained above, but recognizes that the viscosity applicable in thin-plate flexure will not necessarily equal the depth-averaged viscosity applicable in shallow-shelf flow.

To simulate the effects of variable ablation at the surface, we introduce a debris cover of depth $R( x,\; \, t)$ satisfying the following equation:

(67)$${\partial R\over \partial t} = - \nabla \cdot \left(u R \right),\; $$

with an initial condition given by,

(68)$$R = \left\{\matrix{ 0 \ \hbox{m} \hfill & \ x \in [ 4,\; 5] \ \hbox{km} \hfill \cr 1 \ \hbox{m} \ \hfill & \ \hbox{otherwise} \hfill }\right. .$$

To specify the surface ablation ($\dot A < 0$), we determined where $R < 0.15$ m:

(69)$$\dot A = \left\{\matrix{ 0 \ \hbox{m a}^{-1} & \ R < 0.1 \ \hbox{m} \cr - 1 \ \hbox{m a}^{-1} \ & \ \hbox{otherwise} }\right. .$$

This ablation rate was then used to force the equation for ${\cal H}_s$. We assumed that basal accumulation/ablation was zero everywhere, although there is good evidence that the McMIS is sustained against surface ablation by basal freezing (Koch and others, Reference Koch, Fitzsimons, Samyn and Tison2015).

The results of a 20-year simulation of the pedestal formation process are shown in Figure 2. As seen in Figure 2a, the pedestal has an uplifted core, and is surrounded by a moat that is at a slightly lower elevation than the more distal ice shelf ablation area. The presence of the moat surrounding the pedestal suggests that surface meltwater would likely pool there, and this may be in accord with an aerial photograph of Ring Pedestal on the McMIS in 2008 shown in Figures. 3 and 1d.

Fig. 2. (a) Ice-shelf surface and basal profiles at $t = 20$ a. The basal profile has been moved up by 25 m so that the two profiles are easier to see. The pedestalled, relict lake feature stands about 2 m above the surrounding ice shelf. Its weight is mostly compensated by a larger draft indicated by the downward projection of the ice-shelf base. Some weight compensation is supported by flexure, as indicated by (b) showing the oscillations of $\eta$ near the edges of the pedestal feature.

Fig. 3. Aerial view of Ring Pedestal in 2008, at roughly the time the Landsat 7 image of Fig. 1d. We speculate that the blue ring is the depressed area inside the highest part of the Pedestal, seen in Fig. 2. Immediately outside the blue ring is evidence of the debris cover, which is likely exposed within the wall of the moat that surrounds the pedestal. Beyond the pedestal, the surface debris is covered by a thin layer of seasonal snow. This aerial photograph was taken by Andris Apse, the artist in residence at the New Zealand Antarctic Program's Scott Base, who kindly permitted its use in this study.

Ice rumples and traveling lakes

Ice rumples and traveling lakes are viscoelastic buckling phenomena (Thomas, Reference Thomas1973; Ribe, Reference Ribe2003; LaBarbera and MacAyeal, Reference LaBarbera and MacAyeal2011; Slim and others, Reference Slim, Teichman and Mahadevan2012) evident along margins of ice shelves where both ice flow and strong compressive stress are directed toward the grounding line (Fig. 4). The ice-rumples demonstration presents a case where we assume that the shallow-shelf approximation and the shallow-stream approximation apply on opposite sides of the grounding zone, and the grounding line itself is an internal boundary in the numerical domain. For the shallow-stream approximation, we simply add a basal friction parameter as described by MacAyeal (Reference MacAyeal1989).

Fig. 4. (a) Aerial photograph of the McDonald Ice Rumples on the Brunt Ice Shelf seen from upstream. Image was provided by an on-line facility of U.S. NASA (https://earthobservatory.nasa.gov/). (b) En-echelon, petal-shaped lakes on the GVIIS that propagate along the coast of Alexander Island onto which ice-shelf flow is directed (see LaBarbera and MacAyeal, Reference LaBarbera and MacAyeal2011) (WorldView-2 image, 2011, provided by the Polar Geospatial Center, Imagery ©2011 DigitalGlobe, Inc.). Ice flow is directed toward the coast which is the brown area at the bottom of the photograph. (c) and (d) Landsat 7ETM+ enhanced thematic mapper images taken approximately 1 year apart ((c) in 2008, (d) in 2009, see LaBarbera and MacAyeal (Reference LaBarbera and MacAyeal2011) for details). Red vertical lines provide a fiducial reference to fixed points in the two georegistered images. The leftward movement of the water-filled rumple features (traveling lakes) is apparent between panel (d) and panel (c).

The McDonald Ice Rumples are caused when the BIS episodically pushes up on a seabed shoal (Thomas, Reference Thomas1973; Matsuoka and others, Reference Matsuoka2015; De Rydt and others, Reference De Rydt, Gudmundsson, Nagler and Wuite2019). An aerial photograph of the MIR, showing ice shelf flexure features upstream of the zone of grounding, is shown in Figure 4a. Traveling lakes (LaBarbera and MacAyeal, Reference LaBarbera and MacAyeal2011) are similar to the ice-rumpling features of the MIR, however, they involve continuous ice-shelf flow onto and compression against the shore where the GVIIS abuts Alexander Island (Figs. 4b–d), and they often fill with water (adding to the complexity of forcing) (Banwell and others, Reference Banwell2021). For the traveling lakes, the rumpling process creates sets of en-echelon, petal-shaped depressions that fill with meltwater (Figs. 4b–d) and that appear to propagate with time along the coastline even though the ice-flow is directed against the coastline. Ice rumpling and traveling lakes are related to the same phenomena, i.e. where compressive stress sets up an instability that allows the sinuous structure of $\eta ( x,\; \, y,\; \, t)$ to grow exponentially with time, as demonstrated in the Appendix.

To demonstrate how rumpling is initiated by our coupled flow/flexure formulation, an idealized 1-D simulation of ice rumpling is conducted using the 1-D version of the ice-shelf flow, mass balance and flexure equations:

(70a)$${\partial \over \partial x} \left(4 \bar \nu H {\partial u\over \partial x} \right) = -\rho_{sw}g B {\partial B\over \partial x} - \rho_i g H {\partial H\over \partial x} - \beta^2( x) u ,\; $$

(70b)$$- {\partial ^2 M_{xx}\over \partial x^2} + \rho_{sw}g \eta = H\left(4 \nu {\partial u\over \partial x} - {1\over 2} \rho_i g H \right){\partial^2 \eta\over \partial x^2} ,\; $$

(70c)$${\partial \Phi\over \partial t} + u {\partial \Phi\over \partial x} + {3\over \nu_f H^3} M_{xx} = 0 ,\; $$

(70d)$$\Phi - {12 \left(1 - \mu^2\right)\over E H^3} M_{xx} - {\partial^2 \eta\over \partial x^2} = 0 ,\; $$

(70e)$${\partial H\over \partial t} + u {\partial H\over \partial x} = - 0 \cdot H {\partial u\over \partial x} ,\; $$

where $\beta ^2 = 10^{10}$ Pa s m$^{-1}$ is a basal sliding coefficient (MacAyeal, Reference MacAyeal1989) that is non-zero wherever the base of the ice shelf is in contact with the seabed feature causing the rumples (MacAyeal, Reference MacAyeal1989). Note that we have zeroed out the ice-convergence term on the right-hand side of Eqn (70e) in order to avoid excessive thickening of the ice upstream of the grounding zone. This approximates the 2-D situation where transverse divergence balances longitudinal convergence, as observed for the BIS upstream of the MIR (De Rydt and others, Reference De Rydt, Gudmundsson, Nagler and Wuite2019).

For viscosities $\bar \nu$ and $\nu _f$, we take $10^{14}$ and ${{2\over 3}}\cdot 10^{14}$ Pa s. This choice was made subjectively using a trial-and-error process with a sequence of simulations designed to find a choice of the two viscosities that would yield ice rumple amplitude growth in roughly an 8-year time span (what we believe is the best guess of how MIR grew during the BIS's grounding on the seabed shoal where the rumples are found). The subjective criterion was that we wanted our demonstration to produce rumples that were neither too large nor too small compared with the meters-scale rumples of the MIR. In some simulations, if $\nu _f$ was too high, rumples would not grow; in others, they grew too fast.

The simulation is conducted for 20 years, starting from an undisturbed initial condition ($\eta = 0,\; \, \ h = 50$ m). The numerical domain is 14 km in extent (Fig. 5), with a grounding zone between $x = -2$ km and $x = 0$. An ice front is on the left-hand boundary of the domain ($x = -4$ km) and an inflow boundary is on the right-hand side of the domain. From $t = 0$ a to $t = 8$ a, the inflow condition is that the ice velocity is $-1000$ m a$^{-1}$, after $t = 8$ a, the inflow condition is set to zero. We anticipate that any rumples generated in the domain will grow while the inflow condition causes large compressive stress. Once the inflow becomes zero, we anticipate that the rumples will then decay in amplitude, because there is no work being done at the inflow boundary to drive the gravitational potential energy increase of the system when rumples grow. For boundary conditions on the flexure equations, $\eta$ and ${\partial ^2 \eta }/{\partial x^2}$ are both set to zero for all time at both boundaries (we choose to disregard bending moments at the ice front so as to focus only on the rumples).

Fig. 5. (a) Surface and basal profile at $t = 0$, grounded zone indicated by shading. For the first 8 a of the simulation, flow is from the right to the left, with an input velocity at the right boundary of 1000 m a$^{-1}$ and an input thickness of 50 m. After $t = 8$ a, the inflow is set to 0 m a$^{-1}$ and the rumples are allowed to dissipate. The amplitude of the bell-shaped seabed profile is 10 m. (b) Surface elevation at 7.5 (blue line) and 8.0 (green line) years. (c) Depth-average stress (${1\over H} ( \bar T_xx - \bar p )$) at 8 a. (d) Growth of maximum rumple elevation (crest elevation) as a function of time for the 20-year simulation. The rumples grow in amplitude for the first 8 a while inflow at the right-hand boundary is an energy source for gravitational potential increase associated with rumple growth. After 8 a, when the inflow is set to zero, no energy source for building rumples can overcome the natural tendency for rumples to decay, and so they do.

Initial conditions at $t = 0$ are that $h = 50$ m and $\eta = 0$ m everywhere except in the zone of grounding, where $\eta ( x)$ follows a Gaussian profile (Fig. 5) according to a grounding constraint described below. At the right end of the domain $x = 10$ km in Figure 5, inflow is specified as

(71)$$\eqalign{u & = \left\{\matrix{ -1000 \ \hbox{m a}^{-1} \hfill & t \le 8 \ \hbox{a} \hfill \cr 0 \ \hfill & t > 8 \ \hbox{a} \hfill }\right. ,\; }$$

and $h = 50$ m is specified for all time. The 8-year time span for applying the non-zero inflow condition was based on a purely heuristic reason: to stop rumple growth before their amplitude became unreasonably large (i.e. of the order of the ice thickness). With energy being continually provided to the system prior to the cutoff of inflow, rumple amplitude will grow without bounds, as there is no place for the energy to go other than to build ice thickness and rumples in front of the grounding zone. (Relatively little energy is dissipated in the grounding zone by work done against basal friction because the velocity is so low there.)

To simulate a grounding zone on the sea bed, a bell-shaped seabed profile $g( x)$ is specified for the region of the domain $-2 < x < 0$ km. The maximum height of the bell-shaped profile is 10 m and its height is 0 m at either end of the grounding zone. We adopt a novel, possibly useful in future studies, a technique for constraining the ice to be grounded. We approximate the region within $-2 < x < 0$ km to be ‘floating’ in an extremely dense hydrostatic liquid with density $\rho _r = 10 \rho _{sw}$, and where the ‘liquid surface’ conforms to the bell-shaped seabed profile. In essence, this will yield approximately compliant contact between the ice and bed as $\rho _r \to \infty$. To enforce the compliant contact, we replace the second term in Eqn (70b) ($\rho _{sw} g \eta$) with the following term:

(72)$$\rho_r g \left(\eta - G( x) \right).$$

This enforces $\eta$ to follow the curving surface of $G( x)$ as the ice-shelf flows over the seabed causing the ice to rumple. The advantage of specifying the grounding zone with this novel approximation is that boundary conditions do not have to be specified at the contact between floating and grounded parts of the domain.

As expected from the stability analysis of the Appendix (see also Ribe (Reference Ribe2003) and Slim and others (Reference Slim, Teichman and Mahadevan2012)), the compressive flow against, and over, the top of a grounded region treated with the shallow-stream approximation causes a series of sinusoidal waves in $\eta$, $S$ and $B$ to grow approximately exponentially with time and propagate right-to-left toward the leading edge of the ice rumple (Fig. 5). Ice velocity profiles, stress profiles and the growth of maximum rumple amplitude (maximum $\eta$ in the region $x> 0$) are also presented in Figure 5.

A key question arises from the simulation of 1-D ice rumpling: what process eventually limits the growth of ice-rumple amplitude? In our demonstration simulation, we limited rumple growth by simply cutting off the inflow at the right-hand boundary at $t = 8$ a. This eliminated the energy source (work done by inflowing ice against stress), and so the rumples decayed. We do not know what other physical processes may limit ice-rumple growth, and this remains to be investigated in the future. Certainly, as ice-rumple amplitude grows to a few percent of the ice-thickness, the assumptions (primarily that $\eta$ is small relative to $H$) that allow us to use the reduced-order flow and flexure formulations are violated. Even so, this does not necessarily guide us to an explanation of why ice-rumple growth is limited in the physical world.

We speculate on a number of possible processes that may limit ice-rumple growth, and which must be explored further in the future study. As ice-rumple amplitude increases, troughs in the surface corresponding to over-deepened ice-shelf bottom elevation may themselves ground on the seabed, thus limiting further growth of the rumple amplitude (e.g. deactivating the rumpling process). Another possible process is that as rumple amplitude increases, the horizontal span of the rumpled ice shelf must contract in an accordion-like manner. This perhaps limits the compressive stress within the rumpled ice shelf thereby limiting or removing the impetus for further rumple amplitude growth. Finally, there is a fracture within the rumpled zone that may relieve bending moments and initiate more of a thrust and fold like response to the continued application of compressive stress. This may lead to a situation where the fractured ice no longer can flex and simply pushes forward onto the ice rise as a coarse ice mélange. Investigating these speculative ideas as well as others not yet known will be a priority in our future research.

Two-dimensional simulation of ice rumples and traveling lakes

To demonstrate a numerical simulation of ice rumpling in 2-D, we conduct an idealized experiment using the numerical domain shown in Figure 6. The simulation set-up is identical to that of the 1-D simulation presented above insofar as model parameters and geometry are the same along the $x$-axis center section of the numerical domain (the red dotted line in Figure 6a). One difference is that we chose a value of $\nu _f$ equal to ${1\over 2} \bar \nu$ using the same subjective criteria and a trial-and-error process to obtain meters-scale rumples in the demonstration. Boundary conditions are the same, and the grounding zone that excites the rumpling is treated the same way, except that the bell-shaped seabed profile is taken to be circular in plan view (centered in the center of the circle shown in Fig. 6a). For the 2-D simulation, lower spatial resolution (see the variable finite-element mesh in Figure 6b) was required relative to the 1-D simulations (100 m), because of the computational performance limits of the desktop computer we used.

Fig. 6. (a) Numerical domain for 2-D simulation using same parameters and length scales as the 1-D simulation in Figure 5. The circular zone centered on $x = -1,\; \, \ y = 0$ km represents the area where grounding occurs and where the seabed is a bell-shaped function that has the same variation with radius as the $x$-variation of the bell-shaped function in the 1-D simulation. (b) Finite-element mesh used in the simulation (cubic Lagrange interpolation functions).

In 8 years, the rumples developed to meters-scale amplitude (Fig. 7), and following the cutoff of inflow at $t = 8$ a, the rumple amplitude decayed with time (Fig. 8). We speculate, but reserve investigation for future study, that the lower resolution of the 2-D demonstration relative to the 1-D demonstration excited rumples with a longer wavelength and thus faster growth rate (see the Appendix). Details of what physical process determines the length scale, and thus growth rate, of rumples, and that represented in their numerical simulation, is again left to future study.

Fig. 7. Simulation results (at $t = 8$ a) for the 2-D numerical domain shown in Figure 6. All parameters and simulation choices are the same as in Figure 6, except that this simulation is in two dimensions. (a) Ice-shelf velocity field (m/a). (b) Surface elevation with rumples on the upstream side of the grounding zone (circular area). The inset shows a close up of the curved geometry of the rumples, which wrap around the upstream part of the grounding zone. (c) Depth-average stress in the $x$-direction (${1\over H} ( \bar T_{xx} - \bar p )$, MPa). (d) Longitudinal section of the surface elevation (red dotted line in Figure 6a).

Fig. 8. Growth of maximum rumple elevation (crest elevation) as a function of time for the 20-year simulation. The rumples grow in amplitude for the first 8 a while inflow at the right-hand boundary is an energy source for gravitational potential increase associated with rumple growth. After 8 a, when the inflow is set to zero, the rumples decay. Oscillations in the maximum rumple elevation occur prior to $t = 8$ a, and these are due to rumples being pushed onto the grounded zone, where they quickly decay.

Aside from the difference in rumple amplitude between the 1-D and 2-D simulations, which may involve numerical artefacts in addition to the 2-D geometry of the latter simulation, the 2-D simulation produced features that are compatible with what is seen for the MIR. As shown in Figure 4a, the crests and troughs of the rumples curve around the approximately circular plan-form of the grounding zone (assumed to be the feature just downstream of the rumples). This curvature is developed in the 2-D demonstration shown here (as in the inset of Figure 7b). We refrain from attempting to demonstrate the rumpling style associated with traveling lakes (Figs. 4b–d), leaving it to future research. We do emphasize, however, that there are numerous ice rumpling styles associated with the flow/flexure formulation presented above.