{'fraction_non_alphanumeric': 0.03484372345846781, 'fraction_numerical': 0.01847290640394089, 'mean_word_length': 4.504558204768584, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical approximation of viscous contact problems applied to glacial sliding

Viscous contact problems describe the time evolution of fluid flows in contact with a surface from which they can detach and reattach. These problems are of particular importance in glaciology, where they arise in the study of grounding lines and subglacial cavities. In this work, we propose a novel numerical method for solving viscous contact problems based on a mixed formulation with Lagrange multipliers of a variational inequality involving the Stokes equation. The advection equation for evolving the geometry of the domain occupied by the fluid is then solved via a specially-built upwinding scheme, leading to a robust and accurate algorithm for viscous contact problems. We first verify the method by comparing the numerical results to analytical results obtained by a linearised method. Then, we use this numerical scheme to reconstruct friction laws for glacial sliding with cavitation. Finally, we compute the evolution of cavities from a steady state under oscillating water pressures. The results depend strongly on the location of the initial steady state along the friction law. In particular, we find that if the steady state is located on the downsloping or rate-weakening part of the friction law, the cavity evolves towards the upsloping section, indicating that the downsloping part is unstable.


Introduction
Viscous contact problems are time-dependent fluid flow problems in which the fluid is in contact with a solid surface from which it can detach and reattach. Contact problems of this type arise when modelling glacial ice flow, which is typically treated as a viscous fluid flow (Schoof & Hewitt 2013). On a large scale, they are relevant to marine ice sheets with a grounding line (Schoof 2007(Schoof , 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 2 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 enables us to evolve the surface of the viscous flow in a 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 two-dimensional 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.

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. 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 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 non-zero tangential stress (a friction law), and a free-surface at the ice-atmosphere interface.
Subglacial cavitation occurs at the ice-bedrock 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 ( , ) > ( ).
Glacial ice is generally modelled with the Stokes equations. If we denote the velocity and pressure by = ( , ) and respectively, these can be written as (2.1b) Here, = (| ( )|) is the effective viscosity of the ice, which is usually modelled with Glen's law (Glen 1958), where ≈ 3 for ice and is a potentially temperature-dependent 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, For a given velocity and pressure field, we define the stress tensor = ( , ) by where is the identity tensor field. If denotes the outwards-pointing normal vector to the boundary of Ω, the normal and tangential stresses at the boundary are defined by = ( ) · and = − . (2.5) The Stokes equations (2.1) 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 For the subglacial cavity problem, we assume this water pressure to be uniform along the length of the bedrock because gravity is unimportant on the spatial scales under consideration. 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 = 0. The possibility of detachment is realised by enforcing the contact boundary conditions · 0, − and ( · ) ( + ) = 0 on Γ . (2.7) Notice that here we are enforcing the ice to either remain attached (i.e. · = 0) whenever the compressive normal stress is larger than the water pressure or to have the possibility of detaching ( · < 0) 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 = on Γ , (2.8) or the Neumann boundary condition 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). However, the value of in (2.9) coincides with the basal shear stress along the base of the domain (see (3.1) below) due to the overall balance of forces in the horizontal direction.
The Stokes equation must be coupled to a local advection equation which describes the evolution in time of the cavity roof: . (2.10b) Since is time-dependent, 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 (2.1) together with the contact boundary conditions (2.7) is equivalent to a variational inequality. This mathematical structure can be exploited by a finite element discretisation.

A Stokes variational inequality
The finite element method can be used to compute approximate solutions of the Stokes equation with the contact boundary conditions (2.7). As a first step, we must find a weak formulation by multiplying (2.1) by suitable test functions and integrating by parts. Weak formulations of partial differential equations are formulated in terms of weak derivatives and the solutions are generally sought in Sobolev spaces (Adams & Fournier 2003). In order to keep the notation simpler, we assume in this section and in Section 2.2 that we enforce the Neumann boundary condition (2.9) instead of (2.8). In this case, a suitable space of test functions is the vector-valued Sobolev space = 1, (Ω) for the velocity and the Lebesgue space = (Ω) for the pressure. Here, = 1 + 1/ , where is the coefficient in Glen's law, and = /( − 1) is the Hölder conjugate of .
Due to the contact boundary conditions (2.7), the associated weak formulation of the Stokes equation is a variational inequality. Given the convex set of admissible velocity fields the weak formulation of (2.1) is: find ( , ) ∈ × such that Here, the operators : × → R and : × → R are defined by (2.13) and the linear functional : → R by (2.14) As we assume that the water pressure is constant along the bedrock, the incompressibility of the flow velocity can be used to rewrite the right hand side as follows: where = − is the effective pressure. Consequently, the solution only depends on , not on the particular values of and . The variational inequality (2.12) is equivalent to the Stokes equation (2.1) 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 (2.1) with the boundary conditions introduced in the previous section. Of the three conditions in (2.7), only the kinematic condition · 0 is enforced explicitly in the definition of ; the remaining conditions are enforced implicitly in (2.12).
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 · 0 in the definition of is enforced via the addition of a nonlinear penalty term in (2.12a) and the variational inequality is transformed into a variational equality over 1, (Ω). This is the approach used in Stubblefield et al. (2021). We find that what works best in combination with the advection equation (2.10) is a mixed method with Lagrange multipliers that enforces an average-wise contact condition. We remark that such a condition could also be enforced via a penalty term, but in this case the contact boundary conditions would not be enforced exactly.
The contact boundary conditions (2.7) can be explicitly enforced in (2.12) 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 (2.12), we need to define a suitable set of admissible Lagrange multipliers. We define the range of the normal trace operator onto Γ by (2.16) 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 + 0 in (2.7): The mixed formulation of the variational inequality (2.12) may be written as: find ( , , ) ∈ × × Σ such that (2.19) Formulation (2.18) offers several advantages for finite element approximation in comparison to (2.12). 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 (2.18c), which are a weak analogue of the contact boundary conditions (2.7), can be rewritten as a nonlinear equation at the discrete level, see (2.22c) below. This means that the contact boundary conditions can be enforced exactly by solving the finite element system analogous to (2.18). 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 (2.12) 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 (2.18).

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 (2.1) to the time-dependent advection equation (2.10). 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/Firedrake-20211103.0 (2021). The code for the viscous contact solver presented in Algorithm 1 is openly available, as described in the data availability statement below.

Finite element approximation of the Stokes variational inequality
On each domain Ω , we seek a finite element approximation to (2.18). 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 velocity-pressure pair is known to be inf-sup stable in two dimensions (Boffi et al. 2013). Furthermore, as we explain in Section 2.2.2, working with a piecewise constant Lagrange multiplier enables the use of a very simple upwinding scheme for solving the advection equation (2.10) in a manner consistent with the discrete contact boundary conditions. We we denote the respective degrees of freedom (DoFs) in R , R and R by u, p and . In order to write an algebraic counterpart of (2.18c), we need the operator that returns the average normal components of a vector ℎ ∈ ℎ along the edges on Γ . The discrete counterpart to (2.18) is (2.22c) Here, we have introduced the matrices B ∈ R × and D ∈ R × , the vector f ∈ R and the nonlinear operators A : R → R and C : R × R → R . The matrices are given by the elements B = ( , ) and D = ∫ Γ ( · ) d and the vector f by The nonlinear operator A is given by where is a regularised form of (2.2) that prevents infinite viscosity at zero strain rates; it is defined by for a regularisation parameter > 0 which we set to = 10 −2 . The nonlinear function C is for an arbitrary > 0 which we set to = 1. We find the computations carried out in this paper to be insensitive to the value of > 0 when ∈ [10 −6 , 10 6 ]. In (2.26), the max operation is understood to be carried out for each of the elements in the vector − + ( u) ∈ R . The use of (2.26) is a common way of expressing contact conditions; a particular advantage is that the nonlinear system (2.22) can be solved with a semi-smooth Newton method (Hintermüller et al. 2002). Moreover, (2.22c) is equivalent to enforcing u 0, 0 and ( u) · = 0, (2.27) a discrete equivalent of (2.7). By solving (2.27), we see that on each edge on Γ we either have that the edge wishes to detach ( u < 0) or to remain attached to the bed ( 0 and u = 0). At the continuous level, one cannot write (2.18c) as in (2.22c) because in general the Lagrange multiplier is only an element of the dual space Σ and may not even be a function on Γ (Stadler 2007).

The discrete advection equation and a contact criterion
In order to write a discretisation of (2.10), 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 ( −1 , −1 ) and ( , ) is denoted by . We write , for the value of u 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 (2.10a) should satisfy two conditions. First, it must be stable, as described in standard textbooks (LeVeque 2007) (in Appendix A.2 we present results with and without stabilisation of (2.10a)). Additionally, the scheme should be compatible with the discrete contact boundary conditions (2.27), in the sense that if for two edges and +1 we have that , = 0 and , +1 = 0, then we should have that +1 = ; that is, the node in between the edges remains unchanged. This means that the discrete counterpart of · ( / , −1) should be defined in terms of u. If not, we generally find that all edges will detach within a few timesteps due to approximation errors. Here, we propose the following explicit scheme: 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 non-positive, see (2.27). However, a node on a previously detached edge can penetrate the bedrock. To avoid this and enforce (2.10b), we simply set +1 = ( ) whenever +1 < ( ). 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 ( −1 , −1 ) 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: belongs to Γ if and only if − ( ) > tol, (2.29) where the tolerance tol is set to 10 −9 . Consequently, given a cavity roof ℎ at time , the attached and detached regions Γ and Γ are fully determined.
Once we have computed a new cavity roof +1 ℎ , the Stokes problem is then solved over a new domain Ω +1 . This means that the mesh is deformed according to +1 ℎ . Since we use vertically extruded meshes for the computations in this paper, a simple algorithm is used which moves each node vertically in a linear manner with respect to the original position of the nodes at = 0. If we denote the mesh nodes over by and let = be the upper for a given 0 , 0 and +1 .

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 ice-bedrock interface in large-scale glacier models which do not resolve the smaller-scale 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 (2.2).
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 (after subtracting from the stress, the boundary conditions become = − on Γ , = 0 on Γ and 0 on Γ ). We non-dimensionalise the calculations by scaling lengths with the wavelength of the bed , velocities with a characteristic velocity scale , and stresses with the characteristic stress scale ( /(2 )) 1/ . The resulting problem depends only on the non-dimensional amplitude of the sinusoidal bed, the non-dimensional velocity prescribed along the upper boundary, and the non-dimensional effective pressure .

Steady subglacial cavities
The steady state of equations (2.1) and (2.10) 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 ( +1 − )/Δ is below a prescribed threshold of 10 −4 . In this section, we find the steady cavity for a Newtonian flow ( = 1 in Glen's law (2.2)) over a sinusoidal bed of small amplitude = 0.01 given by ( ) = cos (2 / ). The non-dimensional effective pressure is set to = 0.3 and the non-dimensional horizontal velocity to = 1. We first 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). A brief description of this method is included in Appendix B.
We use five different meshes with vertices uniformly distributed along the lower boundary. In Table 1 we present the non-dimensional basal shear stress and sliding speed along the cavities. These values are calculated with where = ( , ). The expression for can be derived from the expression for the force exerted by the ice on the bed as in Schoof (2005). Recall that the Lagrange multiplier represents the normal stress + along Γ , see (2.19); therefore can be calculated from via As we show in the next section, the values of and can be used to construct a sliding law.
The formula for presented in (3.1) 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 (2.2). 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 (3.1) 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 (2.7) are satisfied exactly at the discrete level for all of the meshes, because + 0. 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 (2.29), 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 non-zero values of + left of the reattachment node in Figure 3. This is particularly visible for the coarsest mesh with = 16.

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 power-law 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 1981Fowler , 1986Gudmundsson 1997;Schoof 2005;Gagliardini et al. 2007). The law proposed in Gudmundsson (1997) for the uncavitated case can be written as 2) is set to = 1 (green), 3 (blue) and 5 (red); for each , we compute the sliding law for beds of amplitude = 0.01 (circles, light), 0.04 (squares) and 0.08 (triangles, dark). For these computations, we set the scaled velocity at the top of the domain to = 1 and the scaled viscosity parameter to = 0.5. We use a mesh with 192 cells along the lower boundary. The parameter ( ) is computed from the slope of the curve near the origin for the lowest .
pressures, no cavitation occurs and a linear sliding law as in (3.3) with 0 = 1 is found. For effective pressures lower than a critical value 8 2 / , cavitation occurs and the sliding law becomes non-linear, varying with as well as .
We compute the sliding law over a sinusoidal bed of different amplitudes and for = 1, 3 and 5, and we plot the results in Figure 4 using the scaling suggested by (3.3). 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 0 can then be calculated from (3.4); these values can be found in Table 2 together with those obtained in Gudmundsson (1997) and Gagliardini et al. (2007). We see that the results obtained in these works are broadly similar to ours and that the value 0 = 1 obtained in the linearised theory is approached in all cases when = 1.
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  0) 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 (3.3) along the linear segment of the curves (where either little or no cavitation has occurred) decreases with increasing values of and . For example, when = 5, one can observe that the linear segment of the sliding law for = 0.08 clearly does not collapse onto the corresponding linear segment for = 0.01. 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 (3.3) 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 = 1, the numerical results computed with the finite element solver highly resemble those obtained with the linearised solution. For = 0.08 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.

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 & Bindschadler 1986;Sugiyama & 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 & Bindschadler 1986;Mair et al. 2002;Sugiyama & 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 speed 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 solving the Stokes problem 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 0 , a basal sliding speed ,0 , and a basal shear stress ,0 . Instead of prescribing the Dirichlet boundary condition = on Γ , we enforce the Neumann boundary condition = ,0 (1, 0) 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 boundary-layer treatment of the problem. We set = 3 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 = 0.08.
The effects of unsteady water pressures differ depending on the initial steady state from which we evolve the cavity. To illustrate this, we first evolve two different points along the upsloping component of the steady sliding law by oscillating the effective pressure with an amplitude of 0.1 0 and a fixed non-dimensional frequency of 0.4. As a reference, note that one non-dimensional 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. These results indicate that, with increasing cavitation, the amplitude of the sliding speed increases. For the case of small cavitation (Figure 5(a)), the sliding speed is slightly out of phase with the effective pressure. However, this phase difference disappears with larger cavitation, as observed in Figure 5(d) and also in Figures 6 and 7 below. This implies that the maximum sliding speed is most often reached when the effective pressure is lowest. Field measurements have also found maximum surface speeds to take place at moments of maximum water pressures (Iken & Bindschadler 1986;Sugiyama & Gudmundsson 2004). On the other hand, the phase difference between the sliding speeds and the cavity volume appears to change in each numerical test: in Figure 5(a), one can observe that the maximum sliding speed is reached when the cavity is still growing, while in Figure 5(d) it is reached at the time of maximum cavitation. 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.
As mentioned above, the non-dimensional time it takes a fluid particle to traverse the domain is of order ≈ 1. Therefore, the scaled frequency = 0.4 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 = 0.4 and = 2 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. In contrast, 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. More specifically, Figure 6(d) indicates that the maximum of the sliding speed and of the rate of change of the cavity volume occur simultaneously.
We can therefore expect a phase difference between the sliding speed and the cavity volume to arise for large frequencies. For very low frequencies, at each time instant , the cavity shape is approximately that of the corresponding steady state for the values of ( ), ( ) and ( ). Therefore, since increases with the cavity volume under steady conditions, this phase difference will disappear. When comparing the panels (a) and (d) of Figure 5, we see that this phase difference is larger when the cavity is smallest. This could indicate that the characteristic time scale with respect to which we measure the frequency increases for smaller cavity volumes.
Following Sugiyama & Gudmundsson (2004), in Figure 6 we also plot the sliding speed against the effective pressure throughout one cycle. For the lower frequency = 0.4, a clear loop arises in which the sliding speed is greater during cavity growth. For the higher frequency = 2, the loop nearly collapses into a single line. These plots can be compared with those obtained from field measurements in Sugiyama & 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). These plots also show the extent to which unsteady sliding can differ from its steady counterpart.
The downsloping section of the friction curve produces a so-called rate-weakening 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. Rate-weakening sliding has been observed in a laboratory setting for ice sliding over a sinusoidal bed (Zoet & Iverson 2015), although several authors have questioned whether such a sliding regime can arise for more realistic bed geometries (Fowler 1987;Schoof 2005;Helanow et al. 2021). An implication of rateweakening sliding is that the sliding law becomes double-valued, 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 & Hewitt 2013). In Figure 7 we perturb a steady state along the downsloping section of the curve with an oscillating effective pressure of nondimensional frequency = 0.4 with the amplitude set to 0.1 0 (left panels) and 0.01 0 (right panels). As shown in Figures 7(b) and 7(d), we observe that, for perturbations with both large and small amplitudes, the cavity quickly evolves towards the steady state along the upsloping section for a similar value ,0 / 0 . In fact, we find that this phenomenon continues to occur for different frequencies in the oscillations of the effective pressure and different steady states along the downsloping section of the sliding law. This finding could offer an additional reason to not use a sliding law with a rate-weakening regime: since such a regime is unstable in the sense described above, we could expect it to be unachievable under natural conditions.

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. Although we have only applied our method to the problem of subglacial cavitation, it can be extended to the ice-ocean grounding line problem. 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 (2.2). Finally, in Section 4, we have explored the temporal evolution of cavities under unsteady effective pressures and its effect on glacial sliding. One of the features of unsteady sliding studied in this work is the phase difference between the sliding speed, the cavity volume and the effective pressures. Our results show that with increasing frequencies, the phase difference between the sliding speed and cavity volume increases. They also seem to indicate that the maximum sliding speed occurs at the point of minimum effective pressure, at least for sufficiently cavitated states. Similar phase differences have been found in data obtained from field measurements (Iken & Bindschadler 1986;Sugiyama & Gudmundsson 2004;Andrews et al. 2014). Although our results could offer an explanation in terms of an idealised model, it should be noted that changes in measurements of surface elevation of an ice sheet can be the result of many cavities in different states. Finally, we also find that when we fix the value of as a Neumann boundary condition, the downsloping section of the sliding law, also known as the rate-weakening regime, is unstable under finite perturbations. In particular, if we perturb a steady state along the downsloping section, the cavity quickly evolves towards the corresponding point with a similar value of / along the upsloping part.
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 in Elmer when solving the advection equation, see (11) in Gagliardini & 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 Gagliardini et al. (2007, Figure 1)).
We also attempted to solve the subglacial cavity problem with the numerical method from Stubblefield et al. (2021). Although this method is suitable for solving the grounding line and subglacial lake problems presented in that reference, we found that it was not able to evolve the cavity correctly. This method solves the variational inequality via a penalty method, in which the constraint · 0 in (2.11) is enforced pointwise at the discrete level on the attached region via a penalisation term. However, due to the discretisation of the advection equation used in Stubblefield et al. (2021), enforcing this condition does not guarantee that the term · ( / , −1) is exactly zero where · = 0 at the discrete level. Moreover, according to our computations, when applying this penalty method to the subglacial cavity problem, the computed values of ℎ · along the attached region appear to be much less accurate than when using the method presented here, as shown in Appendix A.1. We also remark that the method in Stubblefield et al. (2021) does not include any stabilisation of the advection equation. As demonstrated in Appendix A.2, without stabilisation numerical oscillations appear along the cavity roof after some time steps. For these reasons, we expect that the numerical method from Stubblefield et al. (2021) should be modified for the subglacial cavity problem.
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 (2.18). 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 (2.10). We believe that suitable combinations of these methods could yield further improvements in terms of the speed and robustness of the algorithm. We are also exploring extensions of this method to three-dimensional problems; in particular, generalising the scheme for the advection equation is not obvious and requires further consideration. simpler, because we are enforcing · 0 on Γ in a pointwise manner at the discrete level. However, our numerical tests on the subglacial cavity problem indicate that this is not the case. This appears to be because the discretisation of the advection equation in (A 2) is not compatible with the discrete contact boundary conditions in the sense that ℎ · ( ) = 0 for all ∈ Γ does not necessary imply that − ℎ + = 0 (A 3) on a vertex ( , ) ∈ Γ . Moreover, the values of ℎ · computed with the penalty method appear to be less accurate then the corresponding values computed with the method presented in this work. Due to this lack of accuracy, we no longer have a clear distinction between regions that detach ( · < 0) and those that stay attached ( · = 0). This latter issue is illustrated in Figure 8. In this test, we evolve the cavity by solving the variational inequality equipped with the Dirichlet boundary conditions (2.8) with the penalty method described above with = 0.08 over a mesh with 16 elements along the lower boundary. We set = 1, = 1, Δ = 0.01, = 1, = 0.5 and = 10 −6 in (2.1a). Additionally, in order to avoid detachment, we increase the tolerance in (2.29) to tol = 10 −3 . In Figures 8(a)-(d) we present the evolution of the cavity from a fully attached state when computed our implementation of the method from Stubblefield et al. (2021). Our numerical results show that by = 0.5 the cavity shape is deforming excessively. In Figures 8(e) and (f) we also compare the computed values of ℎ · in the first time step with the penalty method from Stubblefield et al. (2021) and with our method, respectively. As seen in Figure 8(e), we have points where ℎ · > 0 in between points where ℎ · = 0.
We have carried out further computations with more refined meshes, and these indicate that the lack of accuracy in ℎ · persists when using the penalty method. We would like to emphasise that this numerical issue does not arise when solving the problems presented in Stubblefield et al. (2021). It must therefore be due to some particularity of the subglacial cavity problem.
A.2. On the use of upwinding when solving the advection equation In Section 2.2.2 we state that the advection equation (2.10a) should be stabilised with e.g. an upwinding method. Since no stabilisation is required for the numerical tests solved in Stubblefield et al. (2021), here we demonstrate that without upwinding, the subglacial cavity problem quickly becomes unstable.
In order to show this, we solve the same problem as in Appendix A.1 with a new method that is identical to the one presented in this paper except for the fact that the advection equation is discretised as instead of (2.28). That is, instead of taking the average values of ℎ · along the edge immediately upstream each node, the average of the edges upstream and downstream is taken. In this way, we still have compatibility with the discrete boundary conditions but no upwinding. In Figure 9 we can see the evolution of the cavity with and without upwinding. Oscillations clearly arise in the cavity roof at = 1 when no upwinding is used.