Introduction

Viscoelastic material behavior arises in the context of a number of phenomena relevant to the control and evolution of ice-shelf and ice-stream motion and stress fields. Examples of these phenomena include tide-driven grounding-line flexure and migration (e.g. Reference Sayag and WorsterSayag and Worster, 2013; Reference Tsai and GudmundssonTsai and Gudmundsson, 2015), tidally pulsed grounding line ice velocity variations (e.g. Reference GudmundssonGudmundsson, 2011; Reference Rosier, Gudmundsson and GreenRosier and others, 2014), ice-stream stick-slip phenomena (Reference Goldberg, Schoof and SergienkoGoldberg and others, 2014), iceberg calving (e.g. Reference Reeh, Christensen, Mayer and OlesenReeh and others, 2003; Reference ScambosScambos and others, 2009) and, most relevant to the present study, ice-shelf fracture associated with transient surface water loads (e.g. Reference Scambos, Hulbe, Fahnestock and BohlanderScambos and others, 2000; Banwell and MacAyeal, in press). Although treatments of ice-shelf flexure using pure elastic (e.g. Reference SergienkoSergienko, 2005, Reference Sergienko2010, Reference Sergienko2013) or viscous (e.g. Reference Hattersley-SmithHattersley-Smith, 1960; Reference ReehReeh, 1968; Reference Collins and McCraeCollins and McCrae, 1985; Reference Reeh, Christensen, Mayer and OlesenReeh and others, 2003; Reference RibeRibe, 2003, Reference Ribe2012; Reference LaBarbera and MacAyealLaBarbera and MacAyeal, 2011; Reference Slim, Teichman and MahadevanSlim and others, 2012) constitutive laws are familiar to glaciology, viscoelasticity (e.g. Reference MaseMase, 1960; Reference SokolovskySokolovsky, 1969), and especially its nonlinear form (e.g. Reference Findley, Lai and OnaranFindley and others, 1976; Reference Wineman and KolbergWineman and Kolberg, 1995; Reference Vrabie, Ionuţ-Ovidiu and JercaVrabie and others, 2009), which embraces the flow law for ice, is rarely used in the context of understanding ice-shelf phenomena. Often there is strong justification for treating ice-shelf flexure as either exclusively elastic or exclusively viscous: the timescale of the phenomena of interest may be either very short (e.g. Reference Sayag and WorsterSayag and Worster, 2011), in which case elastic flexure treatment is justified, or extremely long (e.g. Reference SchoofSchoof, 2011), in which case viscous flexure treatment is justified. However, the timescales associated with the filling (usually slow) and draining (usually rapid) of surface lakes on ice shelves fall in between the short and long timescales justifying pure elastic and viscous treatments, respectively. Indeed, proposed mechanisms such as that offered by Reference Banwell, MacAyeal and SergienkoBanwell and others (2013) to explain why the sudden disappearance of surface lakes on the Larsen B ice shelf, Antarctica, in February 2002 led to the ice shelf’s disintegration a week later, depend exclusively on viscoelastic effects.

In the present study, we develop a treatment of viscoelastic ice-shelf flexure for the analysis of a single 1 year fill/drain cycle of an idealized supraglacial lake. Our development is distinct in two ways. First, it uses a simple and computationally efficient approximation based on thin-plate theory (e.g. Reference Chou and PanChou and Pan, 1991; Reference Cheng and ZhangCheng and Zhang, 1998; Reference Li, Yang and LuoLi and others, 2009). Secondly, it provides a rudimentary, initial treatment of nonlinear viscoelastic effects associated with Glen’s flow law.

Our treatment is based on well-developed theory and practice within the applied mathematics and theoretical mechanics literature, which, as exemplified by the above-cited literature, is extensive. The intended benefit of the development presented here is to present a simplified formulation described in standard glaciological terms, and to display an example of how the formulation addresses the important problem of supraglacial-lake induced flexure of ice shelves.

Conceptual Model

Where viscoelastic behavior arises in Earth-science applications, it is typically represented by one of several highly idealized conceptual forms (or combinations thereof) consisting of simple mechanical arrangements of springs, dashpots and friction plates. In the present study, we adapt the Maxwell model (commonly cited to have originated in Reference MaxwellMaxwell, 1867), represented by a spring and dashpot in series, because it is useful when elastic and viscous deformations govern the short- and long-term behaviors, respectively, as is the case for many problems in glaciology. To adapt the Maxwell model to the ice-shelf flexure problem, we add another idealized mechanical feature to the spring and dashpot in series: the buoyancy bucket. (An alternative idealization would be to use an additional spring that opposes the load causing deformation of the original spring and dashpot series, such as is done in engineering applications with the use of a Winkler foundation.) The conceptual model we develop here is shown in Figure 1. In this case, three devices, a spring, dashpot and ‘buoyancy bucket’, are arranged in series, and suspended from a rigid anchor point above a body of sea water.

Fig. 1. Idealized spring, dashpot buoyancy bucket system acting as a conceptual metaphor for the response of a floating ice shelf to an imposed surface load. At the initial time, t = 0, the ice shelf is at rest, with the spring and dashpot at their initial, unstrained geometries, and with the bucket only partially submerged. At some time later, t > 0, a mass M is added to the bucket. Ignoring inertial effects, the initial response is for all strain in the system to be caused by extension of the spring, which extends by a distance needed to counterbalance any load that is not compensated by buoyancy associated with the bucket’s position within the water. As t → ∞, the viscous response of the dashpot to the tensile force acting across the spring allows the bucket to sink further, asymptotically approaching a position where buoyancy forces completely compensate the load within the bucket. In this asymptotic final state, strain in the system is entirely associated with displacement of the dashpot piston, as the spring will have relaxed back to its initial, unstrained geometry.

When at rest, and free of any applied loads (represented by weights put inside the bucket), the bucket is suspended at a neutral reference position within the water (to allow both positive and negative perturbations to the load contained within the bucket), which is assumed inviscid. When a load of mass M is added to the bucket, a tensile force of equal magnitude is experienced by both the spring and the dashpot (we disregard inertial effects). This tensile force is equal to the load, Mg, where g is the gravitational acceleration, minus the extra buoyancy force that is caused by the bucket’s additional submersion into the sea water in response to the elastic elongation of the spring. As time progresses, this force will cause the piston in the dashpot to extend in a manner that tends to lower the bucket further into the water, allowing the bucket to generate a greater buoyancy force. As the dashpot piston extends, and the bucket displaces more water, the net force acting across the spring and dashpot will reduce, thus allowing the spring to contract towards its original unextended position. Once sufficient extension of the dashpot piston has allowed the buoyancy force to exactly compensate the imposed load within the bucket, the force acting across the spring and dashpot will vanish, and the system will be in equilibrium.

Thin-Plate Approximation

We adopt a thin-plate treatment of ice-shelf flexure that is applicable to circumstances where the ratio of vertical to horizontal length scales, H and L, respectively, is small (H/L ≪ 1), and where the vertical displacement (assumed constant through the depth of the ice shelf) due to flexure, η, is small compared with H (|η| ≪ H). We also assume that the thickness of the ice shelf does not change significantly when it is deformed by flexure, i.e. that the vertical strain and strain rate are zero to leading order: and , where w is the vertical displacement as a function of z the vertical coordinate, assumed parallel to the thinnest dimension of the ice shelf when it is unfixed and the overdot denotes , where t is time. In this circumstance, we may write w and in terms of η:

(1)

where x and y are the horizontal coordinates. The above expressions for w and provide the basis for the key simplification made possible by the thin-plate approximation. Following the theoretical development by (among many others) Reference Timoshenko and Woinowsky-KriegerTimoshenko and Woinowsky-Krieger (1959), Reference ReehReeh (1968), Reference RibeRibe (2003), Reference SergienkoSergienko (2005, Reference Sergienko2010) and Reference Slim, Teichman and MahadevanSlim and others (2012), we write the leading-order expressions for the horizontal displacements, u and v, and velocities, and , respectively:

(2)

where is a vertical coordinate that is 0 at the central material plane of the ice shelf and at the upper and lower (air and water) surfaces of the ice shelf, and where z is 0 at the ice-shelf base (chosen so that = 0 represents the mid-plane, or neutral surface, of the ice shelf).

The thin-plate approximation also allows the horizontal strains e and strain rates ė to be expressed in a simple form:

(3)

where E and are 2 by 2 tensors having components determined by η only:

(4)

with the definition for being similar.

The thin-plate approximation allows the stress-balance conditions that are three-dimensional in their most general form to be replaced with a two-dimensional treatment of bending moments. The bending-moment components are related to the stresses by the following relation. This is done by defining a vertically integrated form of the stress commonly referred to as the bending-moment tensor M:

(5)

where T_xx , T_yy and T_xy are direct and shear components of the stress tensor T acting in the horizontal directions.

The dynamic pressure, p, defined as the deviation of the pressure field from the hydrostatic pressure due to viscous deformation within the plate, is equal to the zz-component of the deviatoric stress:

(6)

Applying the incompressibility condition,

(7)

to the viscous strain rates (elastic strains are not assumed incompressible), the dynamic pressure may be written

(8)

where we have introduced the viscosity ν assumed constant for the time being. Equations (8) and (3) constitute the simplifying assumptions made possible by assuming the ice shelf to be both thin and undergoing small amplitude deformation.

Application of the Maxwell model

According to the Maxwell model, a single stress field drives both extension of the spring and movement of the dashpot. In our treatment of the ice shelf, we assume a priori that a single bending-moment field determines both elastic and viscous deformation, and that the simple sum of elastic and viscous deformation determines the strain field E and therefore the vertical displacement of the ice shelf η:

(9)

(10)

where subscripts e and f (f for fluid) denote elastic and viscous components of the strain and displacement respectively. The above partition of the total strain and displacement represents the key practical step that allows for viscoelastic treatments of thin-plate flexure, because if the time derivative of the above equations is taken, we obtain

(11)

(12)

We may further deduce from Eqn (12) that the various x and y derivatives of η are given by sums of elastic and viscous components. In particular, we define the elastic, viscous and combined curvatures (second derivatives with respect to x and y) of the vertical displacements using vectors H _e, H _f and H:

(13)

and note that

(14)

(15)

(16)

where

where H is composed using the total flexure η = η _e + η _f in a manner similar to H _e and H _f, and the superscript ^T denotes the transpose of the vector variable on which it appears.

Homogeneous elastic and viscous parameters

We first develop the governing equations under the assumption that all elastic and viscous flow parameters are homogeneous constants that do not vary from place to place within the ice. Applying the assumption that the ice is a Maxwell solid, specifically that a single bending moment, M, determines both the elastic and viscous parts of the deformation, M will determine both H _e and H _f in a unique manner. In particular, we write Eqn (16) using the well-known relations from thin-plate theory that relate η _e and to the bending-moment components:

(17)

where the operators D and V are used to relate components of M to the various derivatives of η. Using the familiar results from elastic and viscous plate flexure, we determine the operators D and V from the following expressions that relate bending-moment components to the elastic and viscous components of the displacement:

(18)

(19)

where E and μ are the Young’s modulus and Poisson’s ratio, respectively, and where ν is the viscosity, which for now we assume to be a constant. With the above expressions, and writing M as a vector with three components as above, we may easily invert the operators D and V exactly:

(20)

and

(21)

Returning to Eqn (17), if we define a temporary variable of convenience Φ:

(22)

Eqn (17) becomes

(23)

The final step in developing governing equations for viscoelastic flexure of an ice shelf is to consider the balance of forces and torques associated with imposed loads on the ice shelf. This consideration gives the familiar equilibrium equation relating bending moment to an applied surface load F(x, y, t) (assumed downward, negative) and to the effect of buoyancy ρ _sw gη x, y, t) (assumed upward, positive when η < 0):

(24)

where ρ _sw is the density of sea water in which the ice shelf floats, and g is the gravitational acceleration at sea level.

We now recognize that we have a closed system of seven equations in seven unknowns (the three components of Φ, the three components of M and the single scalar η):

(25)

(26)

(27)

This system of governing equations applies to the problem of viscoelastic ice-shelf flexure in circumstances where the material properties of the ice shelf are homogeneous, namely that the elastic and viscous parameters are constants.

Two-dimensional axisymmetric viscoelastic flexure

To set up the example of viscoelastic flexure in response to fill/drain cycles of a supraglacial lake with axisymmetric geometry, to be presented in the next section, we record the governing equations expressed in polar coordinates, r and θ, where it is presumed that the lake’s center is at r = 0 and that the lake is perfectly circular, thus rendering derivatives to be zero. Equations (25–27 ) expressed in polar coordinates, under the assumption of axisymmetric geometry, still require vector notation, because the bending moment M has both M_rr and M_θθ components (but the M_rθ = M_θr components are zero, by axisymmetric symmetry):

(28)

(29)

(30)

In the present circumstance,

(31)

(32)

(33)

(34)

and

(35)

Glen Rheology

Modifying the above treatment of viscoelastic ice-shelf flexure to account for the non-Newtonian creep behavior of ice presents a challenge (see, e.g., Reference Findley, Lai and OnaranFindley and others, 1976, ch. 12). According to the above treatment, horizontal strain and strain-rate components vary linearly with ζ. Under the assumption of linear elasticity and viscous flow, this implies that variation of horizontal stress components within the ice shelf (i.e. T_xx , T_yy and T_xy = T_yx ) is also linear in ζ. With Glen’s flow law, the assumed linear variation of horizontal strain rates with ζ then implies that stress varies as , where n is the flow-law exponent.

This may seem like an insurmountable dilemma. Its resolution, however, comes from the thin-plate assumptions we have presented above. According to the assumptions, the bending moment M and η are the principal variables of the problem. We thus proceed with the analysis of viscoelastic creep flexure by, as before, eliminating stress as a variable in favor of the bending moment M.

We begin by expressing Glen’s law using a viscosity ν that is a function of the second invariant of the strain-rate tensor

(36)

where T ^′ is the deviatoric stress tensor, and we have dropped the addition of brackets and subscript f denoting viscous components for notational convenience. The viscosity is a function of the second invariant of the strain rate (e.g. Reference MacAyealMacAyeal, 1989, eqns (39) and (40)):

(37)

where B(ζ) is the flow rate constant, and

(38)

Adhering to the simplification associated with the thin-plate approximation, and taking to be a constant (alternative expressions when B is a function of ζ require evaluation of an integral), namely that all strain-rate components in the above expressions vary linearly with ζ, allows us to derive (see Appendix) the following relation between M and :

(39)

where the effective viscosity, , is

(40)

With the above definition for M given by Eqn (39), the expression for M reduces to the value for constant viscosity given in the previous section when n = 1.

As we will consider the ice-shelf response to fill/drain cycles of surface lakes with axisymmetric symmetry, we record the above development in polar coordinates, r and θ:

(41)

where

(42)

Application to Supraglacial-Lake Induced Flexure

The motivation for deriving the above treatment of viscoelastic response to imposed ice-shelf surface loads is to explore the proposition that fill/drain cycles of standing surface water features such as supraglacial lakes can induce transient stress distributions that contribute to ice-shelf break-up (Reference Banwell, MacAyeal and SergienkoBanwell and others, 2013). To explore this idea, and to demonstrate how viscoelasticity influences ice-shelf response to surface hydrology, we conduct an idealized experiment involving a single fill/drain cycle of a surface lake. The purpose of the experiment is to provide a simple framework for understanding the theoretical development above. We shall relegate a full exploration of parameters, sensitivities and geometric conditions to a separate study presented elsewhere (Banwell and MacAyeal, in press).

In our idealized application, we subject the idealized lake to a fill-and-drain schedule, the first 100 days of which are shown in Figure 2. The schedule requires that the volume fraction of meltwater filling the lake increase linearly over a 60 day period (representing a melt season), with the volume fraction reaching 1 (filling the entire holding capacity of the lake) at 60 days. During day 61, the lake is drained of its water over a 6 hour period, and the lake is assumed to remain dry from that point on for an additional year. The simulation begins with the onset of lake filling at t = 0 (corresponding to a calendar date when melting first begins), and extends for 1 year. The asymmetric timing of the fill/drain schedule is motivated by the idea that the initial volume of the surface feature that is being filled is such that surface meltwater routed to the feature during the 60 day fill period will just fill it to capacity. Possibly this correspondence can be a result of fill/drain cycles occurring during previous years, which match the carrying volume of the lake or crevasse to the water available in a melt season. The 6 hour drainage portion of the schedule is motivated by observations of supraglacial lake drainage in Greenland (e.g. Reference DasDas and others, 2008; Reference Tedesco, Willis, Hoffman, Banwell, Alexander and ArnoldTedesco and others, 2013), and presumes that a small conduit is created by hydrofracture to the ocean below the surface water feature that enlarges over time to complete the drainage in a matter of hours.

Fig. 2. Idealized meltwater fill-and-drain schedule used to simulate the impact of an idealized supraglacial lake on ice-shelf flexure (only first 100 days of a 365 day schedule are shown). When the volume fraction reaches 1, the lake is filled to capacity (filled to full depth) with meltwater. Filling is presumed to take 60 days and drainage is presumed to take 6 hours. For 100 < t < 365 (days), we assume the volume fraction is zero.

The idealized surface lake has a circular basin of radius 500 m that is capable of holding meltwater with a radially uniform depth of 2 m, chosen to be comparable with the deepest lakes on the Larsen B ice shelf observed in 2000 (Reference Banwell, Caballero, Arnold, Glasser, Cathles and MacAyealBanwell and others, 2014). The geometry is taken to be axisymmetric, so the treatment uses the governing equations written in polar coordinates recorded in the previous sections with the origin, r = 0, placed at the lake center. A diagram displaying the idealized geometry of the supraglacial lake is presented in Figure 3.

Fig. 3. Cross-sectional geometry of an idealized axisymmetric supraglacial lake.

The ice thickness is taken to be 200 m, and we apply no-displacement no-bending-moment boundary conditions at r = 10 km. At r = 0 we impose symmetry. Parameters of the simulation are set at arbitrary values representative of ice-shelf conditions: E = 10 GPa, μ = 0.3, n = 3, ρ _sw = 1030 kg m⁻³, g = 9.81 m s⁻¹ and the flow rate constant is . Our choice of Poisson’s ratio μ is based on laboratory measurements (e.g. Reference Jellinek and BrillJellinek and Brill, 1956) and seismic phase speed estimates (e.g. Reference Kirchner and BentleyKirchner and Bentley, 1979); other studies (e.g Reference GudmundssonGudmundsson, 2011; Reference Goldberg, Schoof and SergienkoGoldberg and others, 2014; Reference Rosier, Gudmundsson and GreenRosier and others, 2014) point out that μ may be closer to 0.45 to account for incompressibility of the ice. Indeed, the choice of the Poisson’s ratio, which determines the compressibility in elastic deformation, highlights the fact that our viscoelastic treatment must amalgamate conflicting rheologies: the elastic rheology, which permits compression, or P-waves in seismology, is in conflict with the typical incompressibility condition usually assumed in viscous rheology for glaciological applications. In the computational expression of , we add a moment, 1 Pa m, to regularize the computation when bending moments are near zero or absent. As stated above, these parameters are chosen so as to provide a demonstration useful for visualizing the consequences of viscoelastic behavior. All simulations are done using a finite-element package, COMSOL version 4.

Results

The response of the ice shelf to a 1 year fill/drain cycle of the lake is shown in Figures 4 and 5. At t = 0 there is no vertical displacement, and η(r) is everywhere zero. As the lake comes to full capacity after 60 days of filling, the ice shelf is depressed to ∼58 cm at the center of the lake, r = 0. This depression would continue to increase if the lake were to remain full, as the viscoelastic adjustment to the increasing lake load is slow compared with the 60 day timescale of filling.

Fig. 4. Vertical displacement, η(r, t), at the lake center (r = 0) as a function of time t over a 1 year simulation in which the lake is filled and drained according to the schedule shown in Figure 2.

Fig. 5. Displacement profile η(r, t) of the ice shelf in response to the 1 year fill/drain cycle of the axisymmetric supraglacial lake. From t = 0 to t = 60 days, the lake slowly fills and depresses the ice shelf. The depression at the lake center reaches ∼58 cm. Immediately on draining, the ice shelf rebounds to a depression of ∼36 cm at the lake center in response to the sudden removal of the load; but the rebound is incomplete (i.e. the ice shelf does not return to its initial undeformed configuration), because of viscoelastic deformation during the period of lake filling. As time increases, the ice shelf relaxes toward an undeformed state; but does not complete the adjustment by the end of 1 year. This implies that multiple years of fill/drain cycles can deepen supraglacial lake basins. The forebulge associated with the deformation moves inward toward the center of the lake as the lake recovers from the sudden drainage event.

During day 61 of the simulation, the lake is drained during the initial 6 hours of the day. If there had been no viscoelastic adjustment during the 60 day fill period, and if the depression at the end of the fill period were purely elastic, with no creep, the ice shelf would immediately rebound to its initial, unflexed position held at t = 0. However, as is shown in Figure 5, at the end of day 61, the depression at the lake center, r = 0, is ∼36 cm. This means that during the 60 day fill period, an amount of depression at the lake center <36 cm occurred due to viscoelastic adjustment. The 36 cm depression on day 61 accounts for both this viscoelastic adjustment plus the immediate elastic response to the excess buoyancy associated with the ice shelf being depressed by 36 cm at the end of filling.

The ice shelf remains significantly depressed for 20 days following the sudden drainage, but begins to viscoelastically relax toward its initial condition prior to lake filling, as the excess buoyancy is reduced. With the parameters used during the demonstration, the lake does not return to its initial undeformed state at the end of a 1 year period. This suggests that multiple fill/drain cycles over multiple years can act to deepen lakes beyond the 58 cm depth assumed initially in this idealized experiment, and that multiple years of such cycles may in fact determine the hypsometry of supraglacial lake basins.

The radial and axisymmetric stresses, T_rr and T_θθ , respectively, evaluated at the upper surface of the ice shelf ζ = H/2, and the von Mises stress, evaluated at both surfaces, are shown in Figures 6–8. As in Reference Banwell, MacAyeal and SergienkoBanwell and others (2013) and Reference MacAyeal and SergienkoMacAyeal and Sergienko (2013), tensile stress (both radial and axisymmetric) exists at the ice-shelf base below the lake during lake filling, and the von Mises stress reaches values in excess of 70 kPa on day 60 below the footprint of the lake (r < 500 m). This level of tensile stress implies that a fracture could form (e.g. Reference Albrecht and LevermannAlbrecht and Levermann, 2012) within the lake basin at the ice-shelf base and propagate upward to eventually become a conduit for meltwater drainage to the sea water below. Additionally, during lake filling, a zone of tensile stress, ∼65 kPa in magnitude, also exists on the ice-shelf surface in the pattern of T_rr in an annulus located ∼1.0–1 .5 km from the lake center. After lake drainage, on day 61, the ice shelf’s upper surface becomes tensile within a radius of ∼1 km, again with sufficient amplitude to induce fracture in the bottom of the lake. A zone of tensile stress, ∼65 kPa in magnitude, also exists at the ice-shelf base (cf. the ice-shelf surface in the case of a filling lake), again in the pattern of T_rr in an annulus located 1.0–1 .5 km from the lake center.

Fig. 6. Radial component of the stress, T_rr , at the upper surface of the ice shelf. The value of T_rr at the ice-shelf bottom is −1 times that shown here. Of particular importance is the fact that the stress regime in response to both lake filling and lake drainage is significantly tensile at values approaching 150 kPa in various ranges of radius. At the ice-shelf base, the radial stress is strongly tensile below the lake during the time it is filling. Immediately following drainage, the zone of strong tensile stress at the ice-shelf base moves into an annulus that is located in the region r > ∼1.25 km.

Fig. 7. Azimuthal component of the stress, T_θθ , at the upper surface of the ice shelf. The value of T at the ice-shelf bottom is −1 times that shown here. Of particular importance is the fact that the stress regime in response to both lake filling and lake drainage is significantly tensile at values approaching 150 kPa near the center of the lake at either the ice-shelf surface or the ice-shelf base. At the ice-shelf base, the radial stress is strongly tensile below the lake during the time it is filling. Immediately following drainage, the zone of strong tensile stress at the ice-shelf base moves to the surface of the ice shelf.

Fig. 8. Von Mises stress, , at both the upper and lower surfaces of the ice shelf. A value of T _vM above ∼70 kPa implies that the ice shelf will be subject to fracture damage.

Discussion and Conclusions

The thin-plate treatment of viscoelastic flexure of an ice shelf developed here is formally justified for situations in which the elastic and viscous parameters of the ice shelf can be treated as homogeneous (constants everywhere within the ice). The problems encountered in nature, however, pose important complexities. The most difficult complexity is the fact that creep deformation of ice, according to Glen’s flow law, is a nonlinear function of stress. Elastic deformation of ice, even if the ice is of homogeneous temperature and fabric, is a linear function of stress. These two functional variations of stress with deformation present an incompatibility between the assumptions commonly practised in dealing with thin plates and shallow ice shelves. Further discussion of this incompatibility and potential resolutions is presented by Reference Findley, Lai and OnaranFindley and others (1976, ch. 12). In this study, we have taken a heuristic approach to this incompatibility by assuming that elastic stresses relevant to determining the curvature of the plate deformation and the bending moments are still linear with the vertical coordinate inside the ice shelf, . This allows us to derive an effective viscosity computed from Glen’s flow law that is subsequently used to approximate the creep deformation in response to the elastic stress.

Despite the shortcomings of the thin-plate approximation, the central expectations for viscoelastic flexure of ice shelves that motivated our interest appear to be present in the simple, idealized demonstration presented. As proposed by Reference Banwell, MacAyeal and SergienkoBanwell and others (2013) and Reference MacAyeal and SergienkoMacAyeal and Sergienko (2013), changing meltwater features on ice shelves produce macroscopically varying stress fields that are relevant to damage and fracture within a zone extending several kilometers from the features themselves. This serves to complement the viewpoint expressed by Reference Scambos, Hulbe, Fahnestock and BohlanderScambos and others (2000) that the microscopically varying stress fields associated with surface meltwater’s impact on fracture tip propagation are at play within the process that ultimately breaks up an ice shelf. Simply put, we believe that ice-shelf fracture introduced by large-magnitude flexure stresses acting in the far field relative to the meltwater feature is as important as ice-shelf fracture introduced by stress enhancement in the near field relative to downward-propagating crevasse tips.

Our application to a 1 year fill/drain cycle idealized by a water-fill volume schedule covering tens of days suggests that the ice shelf responds viscoelastically, and that neither an elastic nor a viscous/creep approach alone would capture the essential evolution of the flexure process. Perhaps the most important impact of the viscoelastic deformation associated with a fill/drain cycle is the fact that it allows the rapid drainage of the meltwater feature to produce a flexure stress that is substantial enough to cause fracture of the ice shelf within the lake basin itself. Under an elastic rheology alone, the draining of a lake would put the ice-shelf flexure stress field back to zero. However, under a viscoelastic rheology, deformation during the period of water-volume filling leads to a situation where a ‘permanent’ curvature in the ice shelf is introduced (permanent denoting non-elastic, but still subject to time evolution). This permanent curvature tends to transfer the gravitational load of the surface water to being supported more by buoyancy than by elastic stresses as the lake fills. When the lake drains, the permanent curvature accumulated during the filling phase of the cycle causes excess buoyancy, and this pushes the ice shelf upward with a force capable of producing fractures in the ice shelf. In effect, viscoelasticity allows the sudden draining of the lake to be as traumatic, in terms of possible fracture generation, as the initial loading. Although our results suggest that the application of a 1 year fill/drain cycle only causes fracture within the actual lake basin, over multiple years the lake basin will deepen, allowing larger volumes of water to be accommodated, and therefore facilitating fracture at a distance from the lake basin. As explored further by Banwell and MacAyeal (in press), this process may enable one filling/draining lake to cause other nearby lakes to also drain. As discussed by Reference Banwell, MacAyeal and SergienkoBanwell and others (2013), this is observationally supported by the fact that immediately prior to the collapse of the Larsen B ice shelf, thousands of lakes with an average depth of 0.8 m drained suddenly.

One final process evident from the demonstration of ice-shelf viscoelastic behavior in response to fill/drain cycles is that viscoelasticity may have an important role in determining the surface topography and small-scale roughness patterns of ice shelves that ablate at their surfaces. As the viscoelastic adjustment in the wake of a lake drainage does not completely finish during the intervening months of the annual cycle before the next melt season begins, over repeated melt seasons basins and uplifts associated with viscoelastic flexure will amplify and become important in determining where meltwater loads accumulate.

Perhaps the most novel generalization of the interpretation of the viscoelastic phenomenon displayed in our study is to say that short-term elastic processes ‘leak’ into long-term viscous behavior. Thus, strictly speaking, a proper simulation of long-term behaviors needs to account for the history of short-term processes. As developed further in many of the theoretical mechanics papers on the subject of viscoelasticity (e.g. the Boltzmann relaxation law discussed by Reference Cheng and ZhangCheng and Zhang, 1998), the state variables describing a long-term viscous behavior depend on the history of the short-term forcing. By embracing this important form of causality linking elastic and long-term viscous behaviors, glaciological response to cyclic forcing functions (e.g. seasonal cycles) can more accurately be described.