1 Introduction
Viscous contact problems are timedependent fluid flow problems in which the fluid is in contact with a solid surface from which it can detach. Contact problems of this type arise when modelling glacial ice flow, which is typically treated as a viscous fluid flow (Schoof and Hewitt, 2013). On a large scale, they are relevant to marine ice sheets with a grounding line (Schoof, 2007, 2011) and, on a smaller scale, to the formation of subglacial cavities when the ice slides over bedrock undulations (Fowler, 1986; Schoof, 2005; Gagliardini et al., 2007). These problems share a very similar mathematical structure and are of great importance for understanding ice sheet dynamics and predicting future sea level rise.
In this paper, we present a novel numerical approach for solving viscous contact problems in an accurate and robust way. This method relies on the formulation of the contact problem as a variational inequality. There exist many different approaches to solving variational inequalities with finite element methods. Of these, we find that solving the contact problem with a piecewise constant Lagrange multiplier is particularly well suited since it allows us to satisfy a discrete version of the contact boundary conditions exactly in the discretisation. We find that this latter property of the numerical scheme allows us to evolve the surface of the viscous flow in a particularly robust way. To the best of our knowledge, viscous contact problems have only previously been solved as variational inequalities in Stubblefield et al. (2021), where a penalty method is used.
Although the numerical method presented in this paper is suitable for any viscous contact problem, here we focus on subglacial cavitation and its application to glacial sliding. Glacial sliding represents a fundamental component of glacier dynamics and is one of the major uncertainties in ice sheet modelling (Ritz et al., 2015). It involves a variety of physical mechanisms among which is cavitation. Sliding with cavitation has been studied theoretically by means of linearised approaches (Fowler, 1986; Schoof, 2005) and these investigations are limited to relatively simple scenarios with twodimensional geometries, smooth bedrocks with small variations and steady conditions. Numerical methods, on the other hand, allow the extension to more complicated situations in a more straightforward manner (Gagliardini et al., 2007; Helanow et al., 2020). In this work, we attempt to more fully exploit the mathematical structure of the viscous contact problem, to improve the accuracy and robustness of the numerical methods employed.
In Section 2 of this paper, we introduce our algorithm for solving viscous contact problems, using the setup of the glacial cavitation problem to provide a concrete context. We then present two applications of the algorithm. First, in Section 3, we compute the steady sliding law for flow over a sinusoidal bed for linear and nonlinear rheologies, expanding on the results in Gagliardini et al. (2007). Furthermore, in this section we compare our results with those obtained from the linearised theory to validate the algorithm. Then, in Section 4, we explore the effects of unsteady water pressures on glacial sliding by calculating the basal sliding velocities and cavity shapes under oscillating water pressures.
2 Formulation and approximation of the problem of subglacial cavitation
We focus here on formulating the viscous contact problem of subglacial cavitation and describing a finite element scheme to approximate it. Subglacial cavitation is mathematically very similar to the problem of a marine ice sheet with a grounding line. Therefore, the extension of the finite element scheme presented here to grounding line problems requires minor modifications, such as the inclusion of nonzero tangential stress (a friction law). Note that the problem described here is the same, though with some different notation, as the subglacial cavitation problem studied in Gagliardini et al. (2007). A linearised version, assuming small amplitude of the bed bump, is equivalent to the problem studied in Fowler (1986) and Schoof (2005).
Subglacial cavitation occurs at the icebedrock interface, over length scales corresponding to the size of the bedrock obstacles. These length scales are generally several orders of magnitude smaller than those of the glacier. For this reason, the computational domain in which we model the formation of cavities is a thin layer of ice of finite height located under a larger mass of ice, see Figure 1. We assume the bedrock, and therefore also , to be periodic in the horizontal direction. The upper boundary represents a fictional boundary separating from the ice above. The height of the bedrock is given by the function and the height of the lower boundary by . The lower boundary is divided into an attached region where , and a detached (i.e. cavitated) region where .
The ice is generally modelled with the Stokes equations. If we denote the velocity and pressure by and respectively, these can be written as
 ∇⋅(2ηε(u) ) + ∇p &= 0,
∇⋅u&= 0.
Here, is the effective viscosity of the ice, which is usually modelled with Glen’s law (Glen, 1958),
(1) 
where for ice and
is a potentially temperaturedependent parameter, but which we treat as constant. The tensor
represents the strain rate of the ice and is given by the symmetric component of the velocity gradient,(2) 
For a given velocity and pressure field, we define the stress tensor by
(3) 
where is the identity tensor field. If
denotes the outwardspointing normal vector to the boundary of
, the normal and tangential stresses at the boundary are defined by(4) 
The Stokes equations eq:stokeseq must be complemented with a set of boundary conditions. Along the cavitated region of the lower boundary, we assume the ice to be in contact with water at a pressure , and we prescribe
(5) 
On the attached region , we assume the ice to be lubricated by a thin layer of water connected to the subglacial drainage system. For this reason, we allow the ice to slide freely and set . The possibility of detachment is realised by enforcing the contact boundary conditions
(6) 
Notice that here we are enforcing the ice to either remain attached (i.e. ) whenever the compressive normal stress is larger than the water pressure or to have the possibility of detaching () if the stress equals the water pressure.
On the top boundary , we enforce the boundary condition , where is the overburden ice pressure. Finally, we close the system with either the Dirichlet boundary condition
(7) 
or the Neumann boundary condition
(8) 
where is the basal shear stress. As explained further on in Section 3, we remark that the horizontal velocity is not the same as the basal sliding speed (though they become approximately equal in the linearised theory).
The Stokes equation must be coupled to a local advection equation which describes the evolution in time of the cavity roof:
∂θ∂t + u⋅( ∂θ∂x, 1 )^⊤&= 0,
θ&≥b.
Since is timedependent, so is the domain and the attached and detached regions and along the lower boundary.
Our goal in the rest of this section is to make use of the fact that the instantaneous Stokes problem eq:stokeseq together with the contact boundary conditions eq:contactbc is equivalent to a variational inequality. This mathematical structure can be exploited by a finite element discretisation.
2.1 A Stokes variational inequality
The finite element method can be used to compute approximate solutions of the Stokes equation with the contact boundary conditions eq:contactbc. As a first step, we must find a weak formulation by multiplying eq:stokeseq by suitable test functions and integrating by parts. In order to keep the notation simpler, we assume in this section and in Section 2.2 that we enforce the Neumann boundary condition eq:bcNeumann instead of eq:bcDirichlet. In this case, a suitable space of test functions is the vector Sobolev space for the velocity and the Lebesgue space for the pressure. Here, , where is the coefficient in Glen’s law, and is the Hölder conjugate of .
Due to the contact boundary conditions eq:contactbc, the associated weak formulation of the Stokes equation is a variational inequality. Given the convex set of admissible velocity fields
(9) 
the weak formulation of eq:stokeseq is: find such that
a(u,v u)  b(p, v u) &≥f(v u) && ∀v∈K,
b(q,u) &= 0 && ∀q∈Q.
Here, the operators and are defined by
(10) 
and the linear functional by
(11) 
If the water pressure is constant along the bedrock, as we shall assume, the incompressibility of the flow velocity can be used to rewrite the right hand side as follows:
(12) 
where is the effective pressure. Consequently, the solution only depends on , not on the particular values of and .
The variational inequality eq:vi is equivalent to the Stokes equation eq:stokeseq in the sense that if the velocity is at least twice differentiable and the pressure differentiable, then integration by parts can be performed in order to arrive at eq:stokeseq with the boundary conditions introduced in the previous section. Of the three conditions in eq:contactbc, only the kinematic condition is enforced explicitly in the definition of ; the remaining conditions are enforced implicitly in eq:vi.
There exist many different approaches for solving a variational inequality with the finite element method (Glowinski et al., 1981). For example, in the penalty method, the constraint in the definition of is enforced via the addition of a nonlinear penalty term in eq:vi1 and the variational inequality is transformed into a variational equality over . This is the approach used in Stubblefield et al. (2021). Of these approaches, we find that what works best in combination with the advection equation eq:advection is a mixed method with Lagrange multipliers.
The contact boundary conditions eq:contactbc can be explicitly enforced in eq:vi with a Lagrange multiplier which essentially represents the sum of the normal stress and the water pressure along the attached region . In order to introduce this reformulation of eq:vi, we need to define a suitable set of admissible Lagrange multipliers. We define the range of the normal trace operator onto by
(13) 
The space of admissible Lagrange multipliers is then the convex cone of elements in the dual space satisfying a weak equivalent of the dynamic contact boundary condition in eq:contactbc:
(14) 
The mixed formulation of the variational inequality eq:vi may be written as: find such that
a(u,v)  b(p, v)  λ(v⋅n—_Γ_a) = f(v) & && ∀v∈V,
b(q,u) = 0 & && ∀q∈Q,
μ(u⋅n—_Γ_a) ≥0 ∀μ∈Λ, λ∈Λ and λ(u⋅n—_Γ_a) = 0 &.
Equation eq:mixed1 can be obtained by multiplying eq:stokeseq1 by a test function , integrating by parts, and setting
(15) 
Formulation eq:mixed offers several advantages for finite element approximation in comparison to eq:vi. First, there is no longer any need to build an approximation of the convex set because we seek the velocity in the space . Second, the conditions eq:mixed3, which are a weak analogue of the contact boundary conditions eq:contactbc, can be rewritten as a nonlinear equation at the discrete level, see eq:algebraic3 below. This means that the contact boundary conditions can be enforced exactly by solving the finite element system analogous to eq:mixed. As we explain shortly, this leads to a robust and stable numerical method for solving the complete subglacial cavitation system.
Under certain conditions on the geometry of the domain and the linear functional , the variational inequality eq:vi has a unique solution. A proof can be found in de Diego et al. (2021), together with an analysis of the finite element approximation of the mixed system eq:mixed.
2.2 A numerical method for subglacial cavitation
In this section we present a numerical algorithm for solving the complete subglacial cavitation problem which results from coupling the Stokes equation eq:stokeseq to the timedependent advection equation eq:advection. We write the discrete counterparts to these equations and explain how they are coupled by deforming the domain according to a contact criterion. The resulting algorithm is summarised in Algorithm 1. The solver is implemented in Firedrake (Rathgeber et al., 2016), using the version available at zenodo/Firedrake20211103.0 (2021). The code for the viscous contact solver presented in Algorithm 1 is openly available in the archived repository de Diego (2021) and in https://bitbucket.org/gonzalogddiego/subglacialcavitysolver_2021/src/master/.
Subglacial cavitation is a timedependent problem, so its discretisation requires a partition of a given time interval into intervals of length . At each time step for , we must consider a discrete cavity roof and domain given by
(16) 
Additionally, the lower boundary is the union of the attached and detached regions and respectively. These are also timedependent and are determined at each time step by the contact criterion eq:criterion_contact below.
2.2.1 Finite element approximation of the Stokes variational inequality
On each domain , we seek a finite element approximation to eq:mixed. To this end, we consider a triangulation of and the finite element spaces , and given by piecewise continuous quadratic vector fields, piecewise constant scalar fields and piecewise constant scalar functions on respectively. This velocitypressure pair is known to be infsup stable in two dimensions (Boffi et al., 2013). Furthermore, as we explain in Section 2.2.2, the use of piecewise constant elements for the Lagrange multiplier enables the use of a very simple upwinding scheme for solving the advection equation eq:advection in a manner consistent with the discrete contact boundary conditions.
We write , and , where , and . For the functions
, we denote the respective degrees of freedom (DoFs) in
, and by , and . In order to write an algebraic counterpart of eq:mixed3, we need the operator(17) 
that returns the average normal components of a vector along the edges on . The discrete counterpart to eq:mixed is
A(u)  B p  Dλ &= f,
B^⊤u &= 0,
λ+ C(λ,u) &= 0 .
Here, we have introduced the matrices and , the vector and the nonlinear operators and . The matrices are given by the elements and and the vector by
(18) 
The nonlinear operators are given by
(19) 
and
(20) 
for an arbitrary . Here, the operation is understood to be carried out for each of the elements in the vector . The use of eq:complementarity is a common way of expressing contact conditions; a particular advantage is that the nonlinear system eq:algebraic can be solved with a semismooth Newton method (Hintermüller et al., 2002). Moreover, eq:algebraic3 is equivalent to enforcing
(21) 
a discrete equivalent of eq:contactbc. By solving eq:discretecontact, we see that on each edge on we either have that the edge wishes to detach () or to remain attached to the bed ( and ).
2.2.2 The discrete advection equation and a contact criterion
In order to write a discretisation of eq:advection, we introduce some further notation. We denote the points along the lower boundary by with the index increasing from left to right as in Figure 2. The edge between and is denoted by . We write for the value of associated to at time (by the definition of , this value corresponds to the average value of along ).
A numerical scheme that solves the advection equation eq:advection1 should satisfy two conditions. First, it must be stable, as described in standard textbooks (LeVeque, 2007). Additionally, the scheme should be compatible with the discrete contact boundary conditions eq:discretecontact, in the sense that if for an edge we have that , then we should have that both and remain unchanged. A naive discretisation of eq:advection will detach all edges within a few timesteps due to approximation errors. Here, we propose the following explicit scheme:
(22) 
where the velocity at each node is taken from the edge located immediately upstream. This results in an upwinding scheme with the property that penetration cannot occur along , because the average values of along the edges are nonpositive, see eq:discretecontact. However, a node on a detached edge can penetrate the bedrock. To avoid this and enforce eq:advection2, we simply set whenever .
Finally, given a cavity roof at , we must set the attached and detached regions and to solve the Stokes variational inequality at . This is carried out by assigning each edge with endpoints and to either or . We do this by looking at the node downstream of the edge, in accordance with the scheme for the advection equation:
(23) 
Consequently, given a cavity roof at time , the attached and detached regions and are fully determined.
3 Steady sliding with cavitation
The sliding of a glacier over its bedrock has been widely studied since Weertman’s seminal work in 1957 (Weertman, 1957). In general, these studies attempt to build a function known as the sliding law that captures the steady relationship between the basal sliding speed , the basal shear stress and other variables such as the water pressure . This sliding law can then be used to prescribe a boundary condition at the icebedrock interface in largescale glacier models which do not resolve the smallerscale shape of that interface.
As a first application of the algorithm described in the previous section, we build the sliding law for ice flowing over a sinusoidal bed. We first present detailed results for a particular steady state in Section 3.1 to evaluate the accuracy of the solver and the effect of mesh refinement. Then, in Section 3.2 we compute sliding laws for different values of the parameter in Glen’s law eq:glenslaw.
The computations presented in this section and in Section 4 have been carried out under the assumption that the water pressure is uniform along the lower boundary. As explained in the previous section, the problem then only depends on the effective pressure , rather than individually on or (the boundary conditions become on , on and on ). We nondimensionalise the calculations by scaling lengths with the wavelength of the bed , velocities with a characteristic velocity scale , and stresses with the characteristic stress scale . The resulting problem depends only on the nondimensional amplitude of the sinusoidal bed, the nondimensional velocity prescribed along the upper boundary, and the nondimensional effective pressure .
3.1 Steady subglacial cavities
mesh cells  detachment  reattachment  

16  96  0.014772  0.98667  0.7500  1.0000 
32  192  0.015143  0.98633  0.7188  1.0000 
64  768  0.015484  0.98598  0.7188  1.0000 
128  3072  0.015679  0.98577  0.7109  1.0000 
192  7296  0.015741  0.98570  0.7135  0.9948 
The steady state of equations eq:stokeseq and eq:advection with the boundary conditions and on the top boundary can be found with Algorithm 1 by evolving the cavity from an initial state until the norm of the discrete derivative is below a prescribed threshold of . In this section, we find the steady cavity for a Newtonian flow ( in Glen’s law eq:glenslaw) over a sinusoidal bed of small amplitude given by . The nondimensional effective pressure is set to and the nondimensional horizontal velocity to . We consider a linear rheology and a bedrock of small amplitude in order to compare our results with the analytical solution of the linearised cavitation problem considered in Fowler (1986) and Schoof (2005).
We use five different meshes with
vertices uniformly distributed along the lower boundary. In Table
1 we present the nondimensional basal shear stress and sliding speed along the cavities. These values are calculated with(24) 
where . Recall that the Lagrange multiplier represents the normal stress along , see eq:lambdameaning; therefore can be calculated from via
(25) 
As we show in the next section, the values of and can be used to construct a sliding law.
The formula for presented in eq:basal might seem strange if the subglacial cavity domain is interpreted as a boundary layer between an ice sheet and the bedrock. This would suggest we take the sliding speed to be the average value of the horizontal velocity along the top boundary . However, our computations indicate that, if the height of the domain is sufficiently large, we can expect the shear stress to approach a constant value and the horizontal velocity to vary with as approaches , where here is the exponent in Glen’s law eq:glenslaw. Therefore, the horizontal velocity along the top boundary depends strongly on the height of the domain. For this reason, and following Gagliardini et al. (2007), we use the equation in eq:basal to calculate . In this case, we find that is independent of for sufficiently large values of . In particular, throughout this paper we set . In agreement with Gagliardini et al. (2007), we find this value of to be sufficiently large.
In Figure 3 we present the steady cavity shape and normal stresses along the attached region for the different meshes. We can see from these figures that the cavity shape is accurately computed even with the coarsest mesh. Additionally, we also present the stress distribution obtained from the linearised theory, which is uniquely determined for an effective pressure and a sliding speed . The result from the linearised theory plotted in Figure 3 is computed with the value of calculated with the most refined mesh.
The plot for the normal stress distribution demonstrates that the contact conditions eq:contactbc are satisfied exactly at the discrete level for all of the meshes, because . This plot also exhibits the singularity of the normal stresses at the reattachment point. This singularity complicates the approximation of the normal stresses along the attached region and can lead to very inaccurate computations of the sliding law in largely cavitated states. However, Figure 3 also indicates that, with increasing mesh refinement, the solver appears to converge towards the linearised solution.
At the discrete level, the Lagrange multiplier , which we use to calculate , is piecewise constant on each edge along the lower boundary. In Figure 3, these values, the DoFs of , are plotted at the midpoints of each edge. Due to the contact criterion eq:criterion_contact, the edge immediately upstream of the first reattached node is treated as part of (this is required in order to allow for the possibility of subsequent detachment there), which explains why there are nonzero values of left of the reattachment node in Figure 3. This is particularly visible for the coarsest mesh with .
3.2 Computation of the linear and nonlinear steady sliding law
We next perform similar calculations to those of Section 3.1 but for varying effective pressure and powerlaw exponent . This allows us to map out a steady sliding law for ice sliding over a hard bed with cavitation as in Gagliardini et al. (2007). The steady sliding law can be mapped out by only varying and keeping fixed, since dimensional analysis of the steady problem shows that the scaled basal shear stress depends only on the ratio , and not independently on or (Fowler, 1986) (the same will not be true of the unsteady problem in Section 4).
Several previous studies have suggested what form the sliding law should take, both with and without cavitation (Kamb, 1970; Fowler, 1981, 1986; Gudmundsson, 1997; Schoof, 2005; Gagliardini et al., 2007). The law proposed in Gudmundsson (1997) for the uncavitated case can be written as
(26) 
where and is a function depending on . The function is related to the parameter (which also depends on ) considered in Gudmundsson (1997) by
(27) 
For a Newtonian flow, the complex analysis method presented in Fowler (1986) and Schoof (2005) yields an exact solution to the linearised problem. In particular, for high effective pressures, no cavitation occurs and a linear sliding law as in eq:nonlinearslidinggudmundsson with is found. For effective pressures larger than a critical value , cavitation occurs and the sliding law becomes nonlinear, varying with as well as .
We compute the sliding law over a sinusoidal bed of different amplitudes and for , 3 and 5, and we plot the results in Figure 4 using the scaling suggested by eq:nonlinearslidinggudmundsson. The mesh has 192 vertices along the lower boundary and 7296 cells. The location of the cavity endpoints is also plotted in Figure 4, along with the solution to the linearised problem calculated with the method from Fowler (1986) and Schoof (2005). For each , the parameter is computed by calculating the slope of the curve near the origin (where there is no cavitation) for the lowest value of . The corresponding values of the parameter can then be calculated from eq:c0; these values can be found in Table 2 together with those obtained in Gudmundsson (1997) and Gagliardini et al. (2007).
The computed sliding laws with cavitation in Figure 4 are multivalued for as expected (Fowler, 1986; Schoof, 2005). This aspect of the law justifies the use of the Dirichlet boundary condition instead of the alternative Neumann boundary condition. We find that if we use the Neumann boundary condition and initiate the cavity from a fully attached state, the solver always evolves to the steady state associated to the upsloping region of the curve (see also Section 4). From Figure 4, we also deduce that the validity of the sliding law eq:nonlinearslidinggudmundsson along the linear segment of the curves (where either little or no cavitation has occurred) decreases with increasing values of and . For example, when , one can observe that the linear segment of the sliding law for clearly does not collapse onto the corresponding linear segment for . As soon as the cavity size increases and the sliding laws cease to be linear, the aspect of these curves largely differ for different values of .
In Figure 4, we use a different scaling to the one used in Gagliardini et al. (2007). In Gagliardini et al. (2007), the computed maximum value reached by is included in the scaling for the sliding law. In this way, the maximum value reached by the scaled sliding law equals 1 by design. However, we preferred the scaling based on eq:nonlinearslidinggudmundsson because it contains fewer terms that are unknown a priori. It is also worth mentioning that, for different values of , the curves in Figure 4 do not collapse into a single curve when plotted with the scaling from Gagliardini et al. (2007).
For the linear case with , the numerical results computed with the finite element solver highly resemble those obtained with the linearised solution. For
a slight difference with the linearised solution can be seen near the peak of the sliding law. This difference is probably a consequence of nonlinear effects that are accentuated with increasing amplitudes of bedrock roughness.
4 Unsteady sliding
In the previous section, the sliding law was constructed by computing steady cavity states. However, field measurements from alpine glaciers and from the Greenland Ice Sheet have found short term variations in the water pressure, on timescales down to hours (Iken, 1981; Iken and Bindschadler, 1986; Sugiyama and Gudmundsson, 2004; Andrews et al., 2014; Hoffman et al., 2016). In these studies, variations of water pressure have been correlated with variations in surface speeds, vertical strain and uplift. Subglacial cavitation has been considered a possible mechanism causing these correlations (Iken and Bindschadler, 1986; Mair et al., 2002; Sugiyama and Gudmundsson, 2004). These observations motivate an investigation of glacier sliding under unsteady conditions. In this section, we therefore compute the evolution in time of subglacial cavities under oscillating water pressures and calculate the corresponding unsteady basal sliding speeds and shear stresses. The study published by Iken of the transient stages between steady cavity shapes (Iken, 1981) is the only numerical investigation of unsteady cavitation known to the authors.
We initialise the computations from a steady state corresponding to a point in the sliding law determined by an effective pressure , a basal sliding speed , and a basal shear stress . Instead of prescribing the Dirichlet boundary condition on , we enforce the Neumann boundary condition on . We consider it more physically realistic to have the basal shear stress fixed rather than the sliding speed because we can expect the basal stresses to balance the gravitational driving stresses, which are essentially fixed on these timescales (in practice, if water pressure variations are spatially localised, the driving stress can be transferred to neighbouring regions of the bed, but it is not easy to account for this within the current boundarylayer treatment of the problem). We set in Glen’s law to model the nonlinear rheology of ice. Algorithm 1 is used with a mesh with 192 elements along the lower boundary over a sinusoidal bed of amplitude .
The effects of unsteady water pressures differ depending on the initial steady state from which we evolve the cavity. To illustrate this, three different points along the steady sliding law are evolved by oscillating the effective pressure with an amplitude of and a fixed nondimensional frequency of 0.4 (note that one nondimensional time unit is approximately the time taken for ice at the top of the domain to traverse one wavelength of the bed). The results are plotted in Figure 5
. In the top and middle panels, the steady states are points on the upsloping section of the curve. These results indicate that, with increasing cavitation, the amplitude of the sliding speed increases. Moreover, the phase difference between the oscillating sliding speeds and the cavity volumes also changes: in the top figure, one can observe that the maximum sliding speed is reached when the cavity is still growing, while in middle one it is reached at the time of maximum cavitation. On the other hand, maximum sliding speeds occur at the times of lowest effective pressure in all cases. Field measurements have also found maximum surface speeds to take place at moments of maximum water pressures
(Iken and Bindschadler, 1986; Sugiyama and Gudmundsson, 2004). There are slight oscillations in the computed sliding speed when the cavity volume is at its largest. These are numerical artefacts due to the stress singularity at the reattachment point of the cavity having an increasing effect on the overall solution of the problem as the cavity volume grows. In these situations, a small displacement of the reattachment point has a large effect on the stress distribution along the bed and therefore also on the computed sliding speed.The downsloping section of the curve produces a socalled rateweakening sliding regime in which, for a supposed fixed effective pressure, an increase in the sliding speed is accompanied by a decrease in the basal drag. Rateweakening sliding has been observed in a laboratory setting for ice sliding over a sinusoidal bed (Zoet and Iverson, 2015), although it has been questioned whether such a sliding regime can arise for more realistic bed geometries (Fowler, 1987; Schoof, 2005). Under rateweakening sliding, increases in the sliding velocity lead to decreases in the basal drag; therefore, a glacier whose weight is balanced only by basal friction could accelerate unboundedly. For this reason, rateweakening sliding has been considered a potential cause of glacier surges (Lliboutry, 1979; Fowler, 1987). Another implication of rateweakening sliding is that the sliding law becomes doublevalued, as seen in Figure 4. This invalidates the commonly used shallow ice approximation of the Stokes equations which requires the friction law to be invertible (Schoof and Hewitt, 2013).
We find that, when a steady state along the downsloping section of the curve is perturbed with a timedependent water pressure, the cavity quickly evolves towards the steady state along the upsloping section for a similar value . This is depicted in the lower panel of Figure 5, where we can see that the volume of the cavity and the sliding speed quickly decrease and settle around the other steady state associated to . We observed that this phenomenon continues to occur for lower amplitudes and different frequencies in the oscillations of the effective pressure. Therefore, our numerical results indicate that under finite perturbations of the effective pressure, the ice evolves away from the steady states along the rateweakening section of the curve. This finding could offer an additional reason to not use a sliding law with a rateweakening regime: since such a regime is unstable in the sense described above, we could expect it to be unachievable under natural conditions.
The nondimensional time it takes a fluid particle to travel across the domain is of order . Therefore, the scaled frequency can be considered a relatively slow frequency that allows the cavity to approximately follow the steady shapes associated to the effective pressure at each instant in time as calculated in Section 3. In Figure 6, we compare results obtained with frequencies and to examine the effect of faster oscillations in the water pressure. With a higher frequency, the magnitude of the change of cavity volume is significantly lower. Contrastingly, the amplitude of the sliding speed increases slightly with a higher frequency. The phase difference between the velocity and the cavity volume also changes when the frequency is increased. For a high frequency, the maximum velocity is no longer attained when the cavity reaches its largest extent, but before, when the cavity is still growing.
Following Sugiyama and Gudmundsson (2004), in Figure 6 we also plot the sliding speed against the effective pressure throughout one cycle. For the lower frequency , a clear loop arises in which the sliding speed is greater during cavity growth. For the higher frequency of , the loop nearly collapses into a single line. These plots can be compared with those obtained from field measurements in Sugiyama and Gudmundsson (2004, Figure 5); a qualitative similarity between both is that higher speeds are reached when decreases. Additionally, similar plots are presented in Andrews et al. (2014, Extended Data Figure 4).
5 Discussion and conclusions
In this work, we have presented a novel numerical method for solving viscous contact problems by formulating the Stokes variational inequality as a mixed problem with Lagrange multipliers. We have also introduced a numerical scheme for the advection equation specifically designed to ensure both stability and coherence with the contact conditions of the Stokes problem. This combination leads to a robust and simple method for solving viscous contact problems. In Section 3 we have used this method to compute steady cavity configurations. In particular, we have validated the numerical method by comparing our results over bedrocks with small amplitudes to the linearised approach in Fowler (1986) and Schoof (2005) and we have reconstructed steady sliding laws for different values of the parameter in Glen’s law eq:glenslaw. Finally, in Section 4, we have explored the temporal evolution of cavities under unsteady effective pressures and its effect on glacial sliding.
Although not included in this paper, we have compared our method with those from Gagliardini et al. (2007) and Stubblefield et al. (2021). When computing the points along the steady sliding law in Figure 4 with our method, the number of time steps required to converge to a steady state can become very large (of order 1000) in the highly cavitated stages along the downsloping region of the curve. This is due to very small scale oscillations that travel across the cavity but have a significant effect on the calculated values of due to the stress singularity at the reattachment point. Contrastingly, when using the method from Gagliardini et al. (2007), these oscillations seem to dampen and the method can converge in about 100 iterations for highly cavitated states. We speculate that this is due to the use of a numerical stabilisation used in Elmer when solving the advection equation, see (11) in Gagliardini and Zwinger (2008). Despite this difference in computational times, we find that the basal stress computations carried out with our method appear to be more accurate due to the exact enforcement of discrete contact conditions (see Figure 3 and compare with Figure 1 from Gagliardini et al. (2007)).
We also attempted to solve the subglacial cavity problem with the method from Stubblefield et al. (2021). Although this penalisation method appears to be suitable for solving the grounding line and subglacial lake problems presented in that reference, we found that this method was not able to evolve the cavity correctly and total detachment would occur after a few time steps. In this method, via a penalisation term, the constraint in eq:convex_set is enforced pointwise at the discrete level on the attached region. However, due to the discretisation of the advection equation used in this method, enforcing this condition does not guarantee that the term is exactly zero where at the discrete level. It is this lack of “compatibility” between the discretisations of the Stokes contact problem and the advection equation that causes problems when solving the subglacial cavity problem, where the geometry of the cavity can undergo rapid changes under strong variations of the effective pressure.
In the numerical method introduced in this paper, we are essentially solving a Stokes variational inequality by enforcing the average values of to be less than or equal to 0 along the attached region. When formulated as a mixed problem, we arrive at system eq:mixed. However, many other possibilities exist at the discrete level; for example, we could enforce that the projection of onto a space of piecewise linear functions along the attached region of the boundary be less than or equal to zero. Alternatively, we could also enforce this constraint to hold for the values of at the midpoints of the edges. We are currently exploring these different discretisations in combination with different time stepping schemes for the advection equation eq:advection. We believe that suitable combinations of these methods could yield further improvements in terms of the speed and robustness of the algorithm.
[Acknowledgements]We acknowledge helpful discussions with O. Gagliardini.
[Funding] This work was supported by the Engineering and Physical Sciences Research Council (P. E. F., grant numbers EP/V001493/1 and EP/R029423/1). GGdD was supported by the University of Oxford Mathematical Institute Graduate Scholarship.
[Declaration of interests]The authors report no conflict of interest.
[Data availability statement]The code for the viscous contact solver presented in Algorithm 1 is openly available in the archived repository de Diego (2021) and in https://bitbucket.org/gonzalogddiego/subglacialcavitysolver_2021/src/master/. The solver is implemented in Firedrake Rathgeber et al. (2016), using the version available at zenodo/Firedrake20211103.0 (2021).
[Author ORCID]G. G. de Diego https://orcid.org/000000029896024X; P. E. Farrell, https://orcid.org/0000000212417060; I. J. Hewitt https://orcid.org/0000000291676481
References
 Direct observations of evolving subglacial drainage beneath the Greenland Ice Sheet. Nature 514 (7520), pp. 80–83. Cited by: §4, §4.
 Mixed finite element methods and applications. Springer. Cited by: §2.2.1.
 On the finite element approximation of a semicoercive Stokes variational inequality arising in glaciology. External Links: 2108.00046 Cited by: §2.1.
 Finite element solver for subglacial cavities, https://doi.org/10.5281/zenodo.5647206. Zenodo. External Links: Document, Link Cited by: §2.2, §5.
 A theoretical treatment of the sliding of glaciers in the absence of cavitation. Philos. T. R. Soc. A 298 (1445), pp. 637–681. Cited by: §3.2.
 A sliding law for glaciers of constant viscosity in the presence of subglacial cavitation. P. Roy. Soc. AMath. Phy. 407 (1832), pp. 147–170. External Links: Document Cited by: §1, §1, §2, §3.1, §3.2, §3.2, §3.2, §3.2, §5.
 Sliding with cavity formation. J. Glaciol. 33 (115), pp. 255–267. Cited by: §4.
 Finiteelement modeling of subglacial cavities and related friction law. J. Geophys. Res.Earth 112 (F2). Cited by: §1, §1, §1, §2, §3.1, §3.2, §3.2, §3.2, §3.2, Table 2, §5.
 The ISMIPHOM benchmark experiments performed using the finiteelement code Elmer. The Cryosphere 2 (1), pp. 67–76. Cited by: §5.
 The flow law of ice: a discussion of the assumptions made in glacier theory, their experimental foundation and consequences. In Physics of the Movement of Ice: Symposium at Chamonix 1958, pp. 171–183. Cited by: §2.
 Numerical analysis of variational inequalities. North Holland. Cited by: §2.1.
 Basalflow characteristics of a nonlinear flow sliding frictionless over strongly undulating bedrock. J. Glaciol. 43 (143), pp. 80–89. Cited by: §3.2, §3.2, Table 2.
 Sliding relations for glacier slip with cavities over threedimensional beds. Geophys. Res. Lett. 47 (3). Cited by: §1.
 The primaldual active set strategy as a semismooth Newton method. SIAM J. Optimiz. 13 (3), pp. 865–888. Cited by: §2.2.1.
 Greenland subglacial drainage evolution regulated by weakly connected regions of the bed. Nat. Commun. 7 (1), pp. 13903. Cited by: §4.
 Combined measurements of subglacial water pressure and surface velocity of Findelengletscher, Switzerland: conclusions about drainage system and sliding mechanism. J. Glaciol. 32 (110), pp. 101–119. Cited by: §4, §4.
 The effect of the subglacial water pressure on the sliding velocity of a glacier in an idealized numerical model. J. Glaciol. 27 (97), pp. 407–421. Cited by: §4.
 Sliding motion of glaciers: theory and observation. Rev. Geophys. 8 (4), pp. 673–728. Cited by: §3.2.

Finite difference methods for ordinary and partial differential equations: Steadystate and timedependent problems
. SIAM. Cited by: §2.2.2.  Local friction laws for glaciers: a critical review and new openings. J. Glaciol. 23 (89), pp. 67–95. Cited by: §4.
 Evidence for basal cavity opening from analysis of surface uplift during a highvelocity event: Haut Glacier d’Arolla, Switzerland. J. Glaciol. 48 (161), pp. 208–216. Cited by: §4.
 Firedrake: automating the finite element method by composing abstractions. ACM Trans. Math. Softw. 43 (3), pp. 1–27. Cited by: §2.2, §5.
 Potential sealevel rise from Antarctic icesheet instability constrained by observations. Nature 528 (7580), pp. 115–118. Cited by: §1.
 Icesheet dynamics. Annu. Rev. Fluid Mech. 45 (1), pp. 217–239. Cited by: §1, §4.
 The effect of cavitation on glacier sliding. P. Roy. Soc. AMath. Phy. 461 (2055), pp. 609–627. External Links: Document Cited by: §1, §1, §2, §3.1, §3.2, §3.2, §3.2, §4, §5.
 Ice sheet grounding line dynamics: Steady states, stability, and hysteresis. J Geophys. Res.Earth 112 (F3). Cited by: §1.
 Marine ice sheet dynamics. Part 2. a Stokes flow contact problem. J. Fluid Mech. 679, pp. 122–155. Cited by: §1.
 Variational formulation of marine icesheet and subglaciallake groundingline dynamics. J. Fluid Mech. 919, pp. A23. External Links: Document Cited by: §1, §2.1, §5, §5.
 Shortterm variations in glacier flow controlled by subglacial water pressure at Lauteraargletscher, Bernese Alps, Switzerland. J. Glaciol. 50 (170), pp. 353–362. Cited by: §4, §4, §4.
 On the sliding of glaciers. J. Glaciol. 3 (21), pp. 33–38. Cited by: §3.
 Software used in ‘Numerical approximation of viscous contact problems applied to glacial sliding’, https://doi.org/10.5281/zenodo.5643646. Zenodo. External Links: Document, Link Cited by: §2.2, §5.
 Experimental determination of a doublevalued drag relationship for glacier sliding. J. Glaciol. 61 (225), pp. 1–7. Cited by: §4.
Comments
There are no comments yet.