We are looking for a model of ice-sheet flexure that can be coupled with a subglacial hydrological model for an arbitrary number and size of lakes, which suggests starting with a simpler model. Because we expect the amount of uplift caused by subglacial lakes to be much smaller than the thickness of the overlying ice, we model flexure of the ice sheet using the thin plate (Kirchoff) approximation. [Our assumptions are similar to those of MacAyeal et al. (2015), though the details of the derivation differ somewhat; also cf. Evatt and Fowler (2007) for use of a viscous beam model in a jökulhaup problem.] The moment equilibrium equation for a plate can be expressed over a horizontal plane with coordinates (x, y) as (e.g., Boresi and Schmidt, 2003).

where the bending moments M_ij are taken per unit length of the j coordinate line (and assumed symmetric, so M_xy = M_yx) and P is the load per unit area. We use the subscript notation for partial derivatives, e.g., ∂_x = ∂/∂x. The bending moments are defined in terms of the stresses σ_ij in the plate as:

where the integration is from bottom to top of a plate of thickness h, and z is the (positive upward) vertical coordinate.

In order to express the bending moments in terms of vertical displacement, it is necessary to choose an ice rheology that describes the relationship between stress and strain. When considering all time scales, the simplest rheology that includes both instantaneous elastic response and long-time viscous behavior is the Maxwell rheology. The strain rate in a Maxwell viscoelastic material is the sum of two components: an elastic strain rate depending on the stress rate and a viscous strain rate depending on the stress. (Related quantities, e.g., displacement and velocity, can also be expressed as the sum of viscous and elastic components.) For plane stress, the Maxwell rheology gives the strain rates ϵ˙ij in terms of the deviatoric stresses σ_ij and stress rates σ˙ij as (e.g., Turcotte and Schubert, 2002):

for elastic (Young's) modulus E and viscosity η. Assuming that ice is incompressible is equivalent to taking the Poisson ratio ν to be 0.5, so that the pressure p is eliminated; note that this is equivalent to assuming that both the viscous and elastic strain components depend on the deviatoric stresses only.

We also use the relations:

between strains and curvature of the plate (e.g., Turcotte and Schubert, 2002) in order to express the bending moments (2–4) in terms of the vertical displacement w(x, y).

For Antarctic subglacial lakes, which form from basal meltwater, the time scale of filling and drainage can be months to decades, so that we expect the viscous component of w to dominate. However, the reverse will likely be the case for the uplift processes in Greenland, where basal hydrology can be strongly affected by rapid surface meltwater drainage (Dow et al., 2015). (By comparison, the intermediate timescales of ice response to Antarctic supraglacial lakes MacAyeal et al., 2015, make both components significant.) As we are primarily interested in the subglacial Antarctic case here, we will begin by considering the viscous obstacle problem, and later briefly comment on the elastic obstacle problem. In the viscous case (limit as E → ∞), the strain rate and deviatoric stress are then related by 2ηϵ˙ij=σij.

Integrating (2–4) using (8–10) and viscous rheology leads to:

Defining D ≡ ηh³/6 as the viscous analog to elastic rigidity, (1) becomes:

which can be further rearranged as:

Note that we have not assumed D to be constant in this derivation, but in that case the bracketed terms cancel and we are left with the simpler bilaplacian form of the bending equation.

Here we consider subglacial lake formation and associated ice flexure in the context of an entire hydrological system (at least as modeled) rather than independently. In a state-of-the-art 2D subglacial hydrological model with both channelized and distributed water flow (e.g., Werder et al., 2013), lakes form where water flow into an area (such as an overdeepening) exceeds the downstream hydrological capacity of the system. The resulting high water pressure can eventually contribute to lake drainage by enhancing downstream flux and the formation of channels, but pressure above overburden can persist in the lake for several years before this occurs (Dow et al., 2016). It is reasonable to expect that sustained overpressure should cause upward flexure of the ice sheet that would decrease water pressure but possibly cause the lake to spread. We therefore focus on areas where water pressure exceeds overburden and their immediate surroundings. Our intent is first to examine the extent of ice-sheet uplift around a subglacial lake at a given time in the above context. In section 4 we discuss how the results inform models of the hydrological system.

For flexure associated with subglacial lake formation, the load P on the ice sheet above the lake is subglacial water pressure p_w minus the overburden pressure q = ρ_igh, for an ice sheet with density ρ_i and thickness h, and acceleration due to gravity g; that is, the load is the negative of the effective pressure q − p_w typically used in glacier hydrology. It follows that upward flexure occurs when the effective pressure is negative and vice versa, with the condition that w ≥ 0 because zero displacement corresponds to the ice resting on its bed. Because the ice sheet must first be fully supported by the subglacial hydrological system, only water pressure above overburden contributes to upward flexure. It is useful to define the upward component of P:

so that P = p⁺ − q. We note that although p⁺ can be discontinuous, this does not imply a discontinuity in the water pressure p_w, which should be assumed continuous throughout the hydrological model domain.

Before applying the bending equation to an ice sheet, we nondimensionalize it over a circular domain for convenience. For (15) we use the scale R for horizontal distance, the scale W_t for uplift rate, the scale V for viscosity, and the scale H for ice thickness. Because we expect the subglacial water pressure and overburden to be of similar magnitude, we take ρ_igH as the scale for both, using ρ_i = 920 kg/m³ and g = 9.81 m/s². The nondimensional equation is then:

where r is the nondimensional radius,

and

spatial derivatives are in terms of the scaled coordinates, and all lowercase variables are dimensionless regardless of previous usage. We note that D_nc comprises the terms accounting for the effect of spatially variable h and/or η.

We can define the lake as the area where p⁺ > 0, so that a circular lake has nondimensional radius r_L where p⁺ becomes zero. (Here we are focusing on upward flexure of the ice sheet; see e.g., Evatt et al., 2006; Evatt and Fowler, 2007 for lake drainage.) The simplest assumption is to solve (17) applying the “clamped” boundary conditions w = 0 and normal derivative ∂_nw = 0 (i.e., the solution is “flat” across the lake boundary) at r = r_L (cf. the laccolith problem in Turcotte and Schubert, 2002) to find the uplift rate above the lake only. Coupling an ice flexure model with a hydrological model (for any lake shape) will be easier and less computationally expensive if this assumption is justifiable.

However, the stiffness of the ice sheet may result in the uplift area being somewhat larger than the area of positive upward load, even though p⁺ = 0 for r > r_L. It is necessary to solve an obstacle problem to obtain the largest radius such that the ice inside this radius is uplifting and the ice outside remains in contact with the bed. [An obstacle problem requires solving for the free boundary between a region in contact with an obstacle, here the bed, and a region above the obstacle, here the uplift area of the ice sheet. See, e.g., Evans (1998) for a formal mathematical treatment.] Without changing p⁺, we iterate to find the maximum r such that the solution ∂_tw of (17) is nonnegative everywhere within the circle of (nondimensional) radius r (with clamped boundary conditions applied at radius r in each iteration). We will call this maximum value of r the uplift radius r_U. That is, we apply the (nondimensional) difference between water pressure and overburden as the load above the lake and require that any ice outside the lake and lifting away from the bed support its own weight. In the iteration, any value of r greater than the eventual solution r_U results in some region where ∂_tw < 0, that is, in the nonphysical situation where the ice is sinking into the bed (Figure 1). The result of solving the obstacle problem will be the radius of the uplift area surrounding the lake, as well as the uplift rate everywhere within this area.

We have so far set up the viscous obstacle problem for uplift of the ice-sheet area surrounding a subglacial lake at one given time. We first solve this problem over a broad range of parameters, before returning in the Discussion to consider the likely effects over time of flexural uplift on the (modeled) hydrological system.

The flexure equations are solved by the finite element method, using the GetFEM++ library (Y. Renard and J. Pommier, http://getfem.org) via its MATLAB interface. GetFEM++ allows the use of higher order elements that can directly handle clamped boundary conditions; we use the third-order, continuously differentiable Hsieh-Clough-Tocher triangular element (Ciarlet, 1978). We use a mesh with over 9,000 nodes, generated by the DistMesh MATLAB package (Persson and Strang, 2004), which provides a highly uniform mesh over a disk while allowing us to require a node at the origin. We apply a simple bisection method to solve for the largest r such that ∂_tw ≥ 0 everywhere within the circle of radius r.

We consider the obstacle problem (17) over a reasonable range of the nondimensional variables η, h, p, and r_L, with Γ given by the chosen scales for the problem. For simplicity, we assume spatially constant ice thickness h and viscosity η (i.e., D_nc = 0), and consider spatially variable h in section 3.3. The nondimensional results are easiest to evaluate if the scalings are taken from a solution of the dimensional problem. In this case, we use H = 1,000 m, R = 5,000 m, V = 10¹⁸ Pa s, and W_t = 0.10 m a⁻¹.

With variable bed topography, we expect real lakes (or at least lakes in subglacial hydrological models, e.g., Dow et al., 2016) to have maximum pressure at some point inside the lake, with pressure decreasing more gradually to overburden at the lake boundary. (For our purposes in this study, zero pressure above overburden defines the lake boundary.) Many different pressure distributions meeting these conditions are possible. For a smooth profile, we take the shape of the pressure above overburden to be a cubic function:

so that p^* = 1 at the lake center and p^* = 0 at and outside the lake boundary r_L, and its derivative is zero at both the center and the boundary. The full pressure above overburden is then p^* multiplied by the maximum nondimensional value p_max at the lake center, so that p+-h=p*·pmax over the lake.

We find that neither Γ nor η affects the solution for r_U, though the resulting solution for ∂_tw at radius r_U is directly proportional to Γ and inversely proportional to η. Furthermore, r_U scales with r_L; that is, the ratio of uplift radius to lake radius is independent of the size of the lake.

We are left to consider the effects of the nondimensional pressure p and thickness h on the solution of the obstacle problem. The range of ice thicknesses for which subglacial lakes have been observed is well known; here we will use 500–3,000 m. However, water pressure in subglacial lakes is not well constrained. In current subglacial hydrological models that are not yet coupled to ice flexure models, the modeled water pressure in lakes can become quite high as it is not eased by uplift. We therefore run the model across a very broad range of pressures above overburden, from 5 to 1,000 kPa. Figures 2, 3 show model results for low and high pressures, plotted separately for ease of viewing. Figure 4 presents the results in terms of ice thickness to emphasize the nonlinearity with respect to this variable.

Figures 2, 3 show that r_U/r_L for a given ice thickness is a weakly nonlinear function of the pressure above overburden, particularly for lower pressures. For example, the H = 1, 000 m results in Figure 2 can be linearly fitted with adjusted coefficient of determination R̄2=0.969 (not to be confused with the radius scale R), although the shape of the fit and the individual residuals are much better with a cubic (R̄2=0.998) fit. In contrast, r_U/r_L for a given pressure above overburden (Figure 4) is a more strongly nonlinear function of ice thickness. For example, the 100 kPa pressure above overburden results can be linearly fitted with an R̄2 of only 0.814, but a cubic fit (consistent with the dependence of ice stiffness on thickness) produces a much better agreement (R̄2=0.998).

Because subglacial lakes vary in shape, solving the obstacle problem for noncircular lakes is clearly important. We can make a start toward determining the effect of lake shape by considering elliptical lakes. For simplicity, we assume spatially constant ice thickness h and viscosity. If we repeat the derivation of (17) using scales of A for x and B for y (instead of R for both), we arrive at:

where a_L is the lake's minor semiaxis length, Γ_A is (19) with A replacing R, and as before all lower case variables are nondimensional. The shape of the pressure above overburden is a version of (20) recast for an ellipse:

Similar to the circular case, we iterate to find the maximum a such that ∂_tw ≥ 0 everywhere within the ellipse with minor semiaxis length a and major semiaxis length b, calling the solution a_U. By setting up the obstacle problem to be solved for only one of the nondimensional semiaxis lengths, we have assumed that the uplift area will have the same shape as the lake, with the ratio of the semiaxis lengths fixed at a/b. This assumption, which seems reasonable for an initial study, allows the obstacle problem to be solved by a tractable one-parameter iteration.

Having already explored the effects of the other parameters in the a = b circular case, we focus on changing the shape of the lake by varying a/b. We note that a_U/a_L shows the same scale independence as r_U/r_L for (17), so we can consider only a/b without concern for the size of the lake. The results (Figure 5) show that lake shape does indeed matter. With the assumption that the uplift area has the same shape as the lake, a circular lake produces the largest relative uplift area, and any narrowing by reducing a/b results in progressively smaller relative uplift areas. However, the maximum uplift rate ∂_tw_max is nearly inversely proportional to a/b, slightly more than doubling for a/b = 0.5.

Due to the effects of topography and basal friction, it is reasonable to expect the ice above a lake to vary in thickness. For nonconstant ice thickness h, we solve (17) with the D_nc terms included. To vary the thickness, we choose a cubic function similar to our pressure distribution (20):

where h₀ is the relative thickness at the lake center, h^* = 1 at the lake boundary, and the derivative is zero at both the center and the boundary. Ice thickness h above the lake is then h^* scaled by the maximum (constant) value of thickness outside the lake.

As expected, r_U/r_L does increase for thinner ice and decrease for thicker ice over the lake (Figure 6). For example, r_U/r_L changes by ~ 0.006 (H = 500 m) and ~ 0.001 (H = 3, 000 m) for ± 10% thickness change at the center of the lake (using the cubic profile (20) with 100 kPa maximum overpressure). The maximum uplift rate increases approximately linearly with thinning, increasing for example by ~ 10.5% for 10% thinning of H = 1, 000 m and 100 kPa maximum overpressure. The effect of spatial variations in the ice thickness over the lake by up to ± 10% is very similar to that caused by ± 10% changes in the ice thickness in the “constant thickness” case (section 3.1), with differences in r_U/r_L on the order of 1 × 10⁻⁴. Given this result and the nonlinear dependence of r_U/r_L on ice thickness, moderate thinning or thickening over the lake can noticeably affect the size of the uplift area.

As noted earlier, different plausible profiles for water pressure across a lake can be supposed. In order to better determine the effect of the overpressure distribution, we introduce several new functions to complement the cubic function (20) used in section 3.1:

and

These functions are chosen to be one at the lake center and zero at and outside the boundary r_L, and scaled to have the same integrated overpressure as the cubic function (Figure 7). We repeat the (constant ice thickness) experiments shown in Figure 2 with the new pressure distributions. There is some dependence of r_U/r_L on the pressure distribution (Figure 8), with r_U/r_L varying by about 0.005 across the different profiles for H = 1, 000 m and 100 kPa maximum overpressure. The quintic profile, which has the highest overpressure near the lake boundary (Figure 7), results in the highest r_U/r_L, while the cubic profile, the most concentrated toward the lake center, results in the lowest r_U/r_L. Conversely, the maximum nondimensional uplift rate ∂_tw_max is greatest for the cubic profile and smallest for the quintic profile, although the differences are small (Δ∂twmax<2.5×10-3 for H = 1, 000 and 100 kPa maximum overpressure, and an order of magnitude smaller for thicker ice).

Had we instead taken the elastic limit (η → ∞) in (5)–(7) and repeated the derivation of the plate bending model, we would have arrived at an equation of the same form as (15), except solving for the uplift w instead of the uplift rate ∂_tw and with D = Eh³/9 instead of ηh³/6. Solving the analogous obstacle problem shows that the uplift w depends inversely on the elastic modulus E, and that the solution for r_U/r_L depends on h and p in exactly the same way as in the viscous case. It follows that the full viscoelastic obstacle problem needs to be solved for only one radius.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.