A three-dimensional full Stokes model of the grounding line dynamics: effect of a pinning point beneath the ice shelf

. The West Antarctic ice sheet is conﬁned by a large area of ice shelves, fed by inland ice through fast ﬂowing ice streams. The dynamics of the grounding line, which is the line-boundary between grounded ice and the downstream ice shelf, has a major inﬂuence on the dynamics of the whole ice sheet. However, most ice sheet models use simpliﬁca-tions of the ﬂow equations, as they do not include all the stress components, and are known to fail in their representation of the grounding line dynamics. Here, we present a 3-D full Stokes model of a marine ice sheet, in which the ﬂow problem is coupled with the evolution of the upper and lower free surfaces, and the position of the grounding line is determined by solving a contact problem between the shelf/sheet lower surface and the bedrock. Simulations are performed using the open-source ﬁnite-element code Elmer/Ice within a parallel environment. The model’s ability to cope with a curved grounding line and the effect of a pinning point beneath the ice shelf are investigated through prognostic simulations. Starting from a steady state, the sea level is slightly decreased to create a contact point between a seamount and the ice shelf. The model predicts a dramatic decrease of the shelf velocities, leading to an advance of the grounding line until both grounded zones merge together, during which an ice rumple forms above the contact area at the pinning point. Finally, we show that once the contact is created, increasing the sea level to its initial value does not release the pinning point and has no effect on the ice dynamics, indicating a sta-bilising effect of pinning points.


Introduction
Most of the land ice which could contribute in the future to sea level rise belongs to Antarctica.Consequently, it is of significant relevance to investigate the amount of ice which could be discharged from Antarctica and further contribute to sea level rise.Lately, geological observations seemed to show that the last deglaciation ended with a rapid (around a century) meter-scale sea level rise, caused by ice-sheet instability (Blanchon et al., 2009).Some laser altimeter surveys carried out in West Antarctica indicated an acceleration of the ice discharge (Thomas et al., 2004), which has been confirmed recently with observations based on interferometry (Rignot et al., 2008), and also through two independent measurement techniques applied to the last 8 years (Rignot et al., 2011).Such rapid changes in Antarctica may arise from two causes which lead to instability when they are combined.The first is that a large part of the bedrock rests well below sea level, essentially due to isostatic adjustments of the crust below the Antarctic Ice Sheet.The second is that the bedrock is mainly upsloping seaward near the estimated location of the grounding line.These two statements are made by Lythe and Vaughan (2001) and by Le Brocq et al. (2010).As first mentioned by Weertman (1974), since confirmed by Schoof (2007), and verified lately by numerical computations (Durand et al., 2009a), the combination of these two characteristics is a source of ice sheet instability, and is called Marine Ice Sheet Instability (MISI).
In East Antarctica, it seems that the ice mass loss is near zero overall.However, while East Antarctica is less covered by data than West Antarctica, the existing sets of data Published by Copernicus Publications on behalf of the European Geosciences Union.
suggest a MISI for both eastern and western bedrocks (Lythe and Vaughan, 2001), and thus point out potential instabilities in East Antarctica as well.
The dynamics of marine ice sheets are mostly controlled by the dynamics of the transition zone between the grounded ice sheet and the floating ice shelf around the grounding line.Mechanically, grounded ice is mainly both vertically and sideways sheared, while floating ice is mainly longitudinally stretched and sideways sheared.These different stress components are superimposed within the vicinity of the grounding line, resulting in a complex stress pattern involving all components.This complexity is difficult to represent numerically.Until 2005, as demonstrated by Vieli and Payne (2005), none of the numerical models available was able to describe properly the grounding line migration.The IPCC Fourth Assessment Report (Solomon et al., 2007) also pointed out that dynamical effects in marine ice sheets were badly understood, and poorly represented in models.
The simplest of the current models are based on the socalled Shallow Ice Approximation (SIA), initially proposed by Hutter (1983), in which only vertical shear components are considered.The SIA is valid for a small aspect ratio between domain thickness and length, but not where the basal sliding contribution to the ice dynamics increases, or where the basal surface slope is too steep (Le Meur et al., 2004), like in the vicinity of the grounding line.More complex models were also developed, such as those including higher order stress gradients (e.g., Blatter, 1995;Pattyn, 2003;Saito et al., 2003), but these models neglect the resistive stress, whereas at the transition zone between grounded and floating ice, this term should be accounted for.Some models also adopt an SIA approach for the grounded ice where shearing dominates the flow, and a Shallow Shelf Approximation (SSA, e.g., Morland, 1984;MacAyeal, 1989) for the floating ice where stretching dominates.The most complete models resolve the full Stokes set of equations without any approximations in the stress tensor.So far, this approach has not been applied to the grounding line migration issues in 3-D, even if some attempts have been made in 2-D by Lestringant (1994), and more recently by Nowicki and Wingham (2008), Durand et al. (2009a) and Gagliardini et al. (2010).
A likely reason for modellers to simplify the full Stokes set of equations lies in the many degrees of freedom that need to be considered.Most of them are related to the high resolution required (in the case of a fixed grid) at the transition between grounded and floating ice, where it is crucial to capture correctly the spatial evolution of the stress pattern.However, these high numerical costs are increasingly adressed by parallelising the codes and running computations on computer clusters.
Some models take into account a curved grounding line in a 3-D geometry, as in Pattyn (2003) where the shelf is modelled, or in Katz and Worster (2010) where the shelf is replaced by a flux at the grounding line, computed from asymptotic theories from Schoof (2007), but not represented physically.However, little attention has been paid numerically on the effect of a pinning point touching the ice basal surface and creating a local grounded zone surrounded by floating ice, except for Goldberg et al. (2009) in which an SIA-SSA coupled model is used.These pinning points are known to stabilize marine ice sheets (Fricker et al., 2009), and have been detected in Antarctica beneath the Ross and the Amery ice shelves.
In this study, we present a 3-D numerical model based on a full-Stokes approach of a marine ice sheet.We are mainly interested in the mechanical effects, thus temperature is assumed to be uniform and constant.The model is derived from the work done by Durand et al. (2009a), and extended to the 3-D case.The purpose of this paper is mostly to present the model and its ability to represent correctly the behaviour of curved grounding lines and distinct grounded areas in 3-D.
To do so, we applied the model to a pinning point experiment, from which we obtained these grounding lines.
A general description of the model is provided is Sect.2, followed by a description of some numerical details in Sect.3, in particular how the position of the grounding line is solved numerically.Section 4 presents all the numerical experiments, including the verification of the model by testing its reversibility.We also investigate the dynamical effect of a pinning point in contact with the ice-shelf basal surface, and discuss its stabilising effect on a marine ice sheet.

Detailed description of the model
We solve the full-Stokes equations for a 3-D gravity-driven problem flow of isothermal, incompressible and nonlinear viscous ice, which slides over a rigid bedrock or floats over sea water.The grounding line, which separates the floating and the grounded part of the ice, falls into fixed grid points and migrates depending on the solution of a contact problem.

Geometry and main hypothesis
The geometry of the marine ice sheet is three-dimensional.It includes a grounded part that slides over a rigid bedrock and a floating part which undergoes the sea pressure.The main direction of the flow is aligned with the x axis, the z axis is the vertical pointing upwards, and the transversal y axis is perpendicular to the (x,z) plane (see Fig. 1 for geometry details).The domain is bounded transversally by two lateral boundaries, both parallel to the (x,z) plane.The ice divide (x = 0) is considered to be a symmetry plane and the length of the ice sheet remains similar over time, which means that the calving front has a fixed position throughout the simulation.
The constitutive law for the ice behaviour follows a Norton-Hoff type law (also named Glen's flow law in glaciology): The where S is the deviatoric stress tensor, and D the strain rate tensor, whose components are defined by: where u is the velocity vector.The effective viscosity η is expressed as: where εe is the second invariant of the strain rate defined as: In Eq. (3), A is the fluidity parameter, and Glen's flow law exponent is n = 3.

Governing equations
To determine the ice flow velocities and pressure, we solve the full Stokes equations over the ice volume, consisting of the momentum conservation equations without inertial effects: and the mass conservation equation in case of incompressibility: In Eq. ( 5), σ = S − pI is the Cauchy stress tensor with p = −trσ /3 the isotropic pressure, ρ i the ice density and g the gravity vector.The three-dimensional ice body is vertically limited by two free surfaces, namely the upper ice/atmosphere interface z = z s (x,y,t), and the lower interface z = z b (x,y,t) between ice and bedrock or sea (see Fig. 1).An advection equation is solved to determine the evolution of the two free interfaces; its general form is: where u x , u y , u z are functions depending on the triplet (x,y,z j ), and z j and a j depend on the triplet (x,y,t).The subscript j can refer either to the upper interface (j = s) or to the lower interface (j = b), u x (x,y,z j ), u y (x,y,z j ) and u z (x,y,z j ) are the velocity components for the considered interfaces, and a j (x,y,t) is the accumulation/ablation function.Here, the accumulation a s over the upper interface is a downward vertical flux uniform and constant, while the melting/accretion a b below the ice shelf is neglected (see Table 1).

Ice divide and calving front
At x = 0, the ice divide is assumed to be a symmetry plane, so that the prescribed velocity is u x (0,y,z) = 0.The end of the domain corresponding to the end of the ice shelf, called the calving front, is kept at a constant position.If it is sufficiently large, the distance between the calving front and the grounded ice has no influence on the upstream flow dynamics, because no friction is accounted for on the lateral boundaries (see Section 2.3.4).The front boundary undergoes normal stress due to the sea pressure p w (z,t), that evolves vertically as follows: where ρ w is the sea water density and l w the sea level.The stress imposed is thus σ nn = p w and σ nt = 0.

Upper interface
The atmospheric pressure exerted against the upper interface z = z s (x,y,t) is neglected, which implies that , where n is the unit normal vector of the interface pointing outward, t 1 and t 2 the unit tangent vectors at the interface.

Lower interface
The lower interface z = z b (x,y,t) may be in contact with either the sea or the bedrock, so two kinds of boundary conditions coexist on a single surface.Where the bottom surface is in contact with the sea (i.e.above the bedrock), the sea pressure is applied with no shearing, which means that: where the bottom surface is in contact with the bedrock, a nonlinear friction law is prescribed, which is summarised as: where the values of the parameters C and m entering the friction law are given in Table 1.The nodes belonging to the grounding line are defined as grounded nodes with at least one floating neighbour.

Lateral walls
The lateral boundaries of the domain are parallel planes.The first plane (y = 0) is an actual border of the domain, and the second one (y = W ) is a plane of symmetry, which makes it possible to model half of the geometry in the ydirection (see Section 3.2).In both cases, no flux is considered through the surfaces and the boundary condition prescribed is u y (x,0,z) = u y (x,W,z) = 0.The model gives us the possibility to prescribe a nonlinear friction law on this boundary, like the one discussed above in Equation ( 11).This possibility will not be explored in this paper and perfect sliding is assumed.

Numerics
This problem is solved using Elmer/Ice, the glaciological extension of the open-source finite-element code Elmer (downloadable from www.csc.fi/elmer).

General considerations for the mesh
One of the conclusions of the fourth assessment of IPCC (Solomon et al., 2007) was that dynamical effects in ice sheets were badly understood, and badly represented in models.Among these, the grounding line migration is the most important process that models do not represent well if applied to the Antarctic marine ice sheet.Durand et al. (2009a) came to the conclusion that along a flowline, a mesh refinement within tens of meters around the grounding line was necessary to yield consistent results and to provide reversible simulations, even if a minimum offset could not be removed completely for any of their experiments.
Here, the grid size in the y-direction is constant and equal to 2.5 km.The maximum refinement adopted in the vicinity of the grounding line is 50 m in the x-direction.As a fixed grid is used, this refinement is applied wherever the grounding line is expected to migrate, which means over tens of kilometers for experiments carried out in this study (see Section 4).Consequently, we end up with a large amount of mesh nodes: a little under 800 nodes in the x-direction, 21 nodes in the y-direction and 11 nodes in the z-direction, which corresponds to about 180 000 nodes.To optimise CPU time, this mesh is broken down into 48 partitions, each consisting of about 4000 nodes.

Construction of the initial geometry
The initial geometry of the ice sheet is first constructed in 2-D with the model of Durand et al. (2009a) which uses an adaptive mesh refinement, so the finest mesh extends over a small domain around the grounding line.Starting from a 200-m layer of ice extending as far as to the end of the domain, the geometry is evolved until a steady state is obtained.The grounding line always advances during the spin up procedure.The 2-D steady state geometry is then laterally extruded in the y-direction.The 3-D problem having a plane of symmetry (y = 50 km), this is only done over half of the domain in order to save computational resources.The resulting 3-D geometry is then meshed with a different x-distribution to be suitable for the planned 3-D experiments (see experiments detailed in Section 4).Consequently, the relaxation that is done after the extrusion may induce a slightly further migration of the grounding line, which however never exceeds the size of one element.The whole procedure is much faster than growing a marine ice sheet in 3-D from an initial slab meshed by a fixed grid of uniform thickness.

Model algorithm
A simplified algorithm of the model is presented in Fig. 2. The first step, carried out once at the beginning of the simulation, is to initialise the model with a geometry and a related mesh, plus initial pressure and velocity fields.
The non-Newtonian stress-strain relationship and the nonlinear sliding law introduce nonlinearities into the problem, and, after discretisation of the full Stokes equations by the finite element method, the system to solve adopts the nonlinear form A(X) • X = b(X).This system is linearised using the fixed point algorithm A(X i−1 ) • X i = b(X i−1 ) solved using direct methods.The numerical solution is stabilised using The Cryosphere, 6, 101-112, 2012 www.the-cryosphere.net/6/101/2012/ the residual free bubbles method (Baiocchi et al., 1993).The convergence criterion of the nonlinear system is based on the change in the Euclidian norm of the solution between the (i − 1)th and the ith iterations.Convergence is first reached when this criterion falls below 10 −3 , assuming a fixed grounding line at the position obtained during the previous time step.
After the initial convergence, the nodal force F w exerted by the sea on the ice lower surface is computed by integrating the sea pressure (see Eq. 8) over the element boundary faces according to the shape function.Then the contact between the grounded nodes and the bedrock is tested.If the contact force R, computed from the residual of the Stokes system (see Durand et al. (2009a) for details), exerted by the ice on the bedrock, falls below the force exerted by the sea F w , the node concerned is no longer subject to a Dirichlet condition but undergoes the sea pressure (see Fig. 3).This change in the basal conditions leads to a retreat of the grounding line.The finite element problem is thus altered by the change in the boundary condition.The nonlinear system after imposing the altered boundary condition is iterated with a convergence criterion of 10 −5 .
As mentioned before, Eq. ( 7) is solved to adjust the lower and upper surfaces.Technical and numerical details on the method are explained in Gagliardini and Zwinger (2008).The following update of vertical surface positions requires to check whether the difference between the bedrock and the ice lower surface evolves.At this stage, some floating nodes may be reattached to the bedrock following the condition already mentioned above.

Geometry and initial steady state
The bedrock is composed of a linear surface and a Gaussian surface, and defined in metres by: with which corresponds to a linear downsloping bedrock with a superimposed 3-D Gaussian curve acting as a pinning point (see Fig. 4).
The initial steady state is obtained by applying the parameters given in Table 1 (Sect.3.2 describes the spin up procedure).Note that the fluidity A is uniform and kept constant since the ice is assumed to be isothermal.The summit of the Gaussian bump is located about 0.8 m below the ice basal surface.At this stage, only one grounding line exists, it is perpendicular to the (x,z) plane and located at x = 536.8km.At the beginning of the simulation, the sea level is instantaneously lowered by 1 m, which creates rapidly (within one time step) a contact between the Gaussian bump and the ice basal surface, at the position (x,y) = (563,50) in km, about 27 km downstream from the initial grounding line.This is the starting point of all the three experiments presented below, which only differ by the way the sea level evolves and by whether the pinning point is maintained or not (see

Verification of the full Stokes model
The purpose of the following experiment is to verify the model reversibility, as proven by Schoof (2007).After the initial drop in the sea level, we let the geometry evolve over 350 a.The sea level is then increased back to its initial level, and the Gaussian bump is merely removed so that only one boundary between grounded ice and floating ice still remains.Of course, this abrupt change of the bedrock geometry is not geologically justified, but it allows us to test if the model is reversible, i.e. if the grounding line recovers its initial characteristics.The simulation is finally continued until the grounding line reaches a nearly stabilised position after 2000 years.Here, we are only interested in the model reversibility so we only focus on these aspects.More detailed explanations about the effect of the pinning point on the dynamics are given in Sect.4.3.Under the effect of the pinning point, the grounding line curves and advances (see Fig. 5a).Before the artificial removal of the bump, this advance ranges between 9 and 12 km, the largest advance being located directly upstream the top of the bump.Afterwards, the grounding line rapidly recovers its straightness within a century.
At this stage, given the uniform bedrock and the uniform grounding line in the y-direction, the 3-D problem is similar to a 2-D problem.In this case, according to Schoof (2007), the steady state has to be unique, so the grounding line should retreat back to its initial position.Nevertheless, we expect that initial and final steady positions should be similar within a residual offset, depending on the grid size at the grounding line.
After 2000 years of computation, velocities along the two flowlines y = 0 and y = 50 km nearly recover their initial pattern (see Fig. 6), and the grounding line is nearly stabilised (see Fig. 5b)).The offset in both positions equals 3 km.Experiments not reported here indicate that this grid-dependent offset is similar to the one presented in the paper even for larger displacements of the grounding line.This was also shown by Durand et al. (2009a) using a similar (but 2-D) model.Consequently, we are confident that if we had simulated a greater initial advance of the grounding line, the offset would nonetheless have been similar.This is also convenient because simulating these migrations in 3-D is merely not feasible within achievable times.

Experiment A: effect of the pinning point
For experiment A, we start from the same initial steady state as for the verification experiment, the pinning point is never removed, the sea level is kept constant and equal to -1 m.This new contact created between the ice and the bump has dramatic effects on the ice-sheet dynamics, much stronger than if only the sea level was decreased without any bump pinning the ice lower surface.Due to the supplementary amount of friction introduced by the pinning point, a rapid decrease of velocities is observed between both grounded areas (see Figs. 8 and 9), which is as high as the lateral coordinate is close to that of the bump summit.In the vicinity of the grounding line, the longitudinal velocity at y = 50 km is approximately twice its value at y = 0 km, where the influence of this new contact is less significant.Under the effect of this non uniform decrease of velocities, the grounding line advances seaward (see Fig. 8) 6. Vertically averaged velocity in x-direction ūx plotted over x for selected time-levels in the verification experiment, for the two flowlines y = 0 and y = 50 km.Initial steady state (without the pinning point) in black, last point in time before removing the pinning point in blue, and final state after 2000 years in green.at different rates in the order of tens of metres per year, the greatest ones being located directly upstream the pinning point.The ice volume within the domain is also increased (see Fig. 7).As a direct consequence of the non lateral uniformity of velocities, the grounding line bends during its advance.After about 800 a of simulation, both grounded areas start to merge into a single one, which lasts 200 a before the complete merge is achieved after 1000 a (see Fig. 10).
A similar experiment was carried out by Goldberg et al. (2009) with an SSA model, for a much larger gap (around 300 km) between the main grounding line and the pinning point.The authors also observed a merge between both grounded areas with advance rates of the grounding line www.the-cryosphere.net/6/101/2012/The Cryosphere, 6, 101-112, 2012 situated around 70 m a −1 , whereas in our case the rates are lower and between 15 and 30 m a −1 .This discrepancy is not surprising since both experiments were performed under different configurations and are hence not directly comparable.For example, the SSA experiments were conducted with a null velocity prescribed along the lateral walls, whereas we imposed free slip.When velocities drop under the fect of the pinning point, the ice volume is increased and consequently the lateral friction is greater, which acts as a positive feedback accelerating the grounding line migration, which could explain the difference with the SSA experiments of Goldberg et al. (2009).
The merge between both grounded areas is followed by oscillations of the last grounded point (see Fig. 10), which could be related to the fact that the grounding line lies on an upsloping part of the bedrock, which is a well known source of instability in marine ice sheets (Schoof, 2007).This instability stops after 1000 a when only one grounded zone remains.
At the beginning of the simulation, the last grounded point (blue curve in Fig. 9), which belongs to the grounding line part of the pinning point, retreats from its position at the top of the bump, towards the main grounding line.The pinning point actually grows and migrates upstream under the thickening of the ice located between both grounded zones.This migration adopts a linear profile during the first 100 a, to stop at a nearly stable position, around 5 km upward from its initial position, under the effect of the sea pressure around the pinning point.by the friction gap from positive to zero when moving seaward around the grounding line.This is immediately followed seaward by a stress oscillation, which is also observed inside the ice.A similar behaviour is observed at the secondary grounding line which belongs to the pinning point.This tensile behaviour is counterbalanced by a compressive behaviour between both grounded zones before they merge.These variations of the stress pattern concentrate in the central flowline, directly upstream the pinning point.It is also due to the jump in basal friction properties.At the end of the simulation, when both grounded areas merge into one, it seems that this compressive behaviour disappears.This corresponds to the fact that the intermediate floating zone becomes completely grounded, which removes the intermediate gap in the basal friction.The downstream shift between compressive and tensile behaviours forms an ice rumple on the upper surface.Its geometry is more or less stable between 100 and 900 a of simulation and its existence is a consequence of the altered stress field induced by the pinning point underneath.The ice rumple is 50 m high, which represents one tenth of the iceshelf thickness.It is vertically aligned with the position of the basal pinning point, and extends over a similar horizontal distance of around 3 km.

Experiment B: mechanical irreversibility
The purpose of experiment B is to check if the pinning point introduces irreversibility to the ice flow dynamics.This is tested by reverting the sea level back to 0 m, its initial value, a short time after the beginning of experiment A, at t = 1.5 a (this is also done for t = 100 a).In this case, the bedrock is not modified so the bumped part is still in contact with the ice lower surface.Both volume curves for experiment A and experiment B (see Fig. 7) show the same pattern, which indicates that the contact between the pinning point remains.This happens at the very beginning of the simulation so we expect it to continue throughout the following time steps.
To explain this mechanical irreversibility, we plotted in Fig. 12 the minimum critical value of sea level rise sl eq that would be required to release this contact.This value is computed for the grounded lower surface of ice as follows: This parameter corresponds to the difference between the normal stress undergone by the ice surface and the sea pressure, normalised by ρ w g, to interpret it in terms of sea level equivalent.Note that a few of the sl eq values are negative for the more advanced part of both grounding lines.There, the sea load against the ice is higher than the contact load exerted by the ice on the bedrock, indicating that these nodes are about to detach from the bedrock.It turns out that, even after only 1.5 a of experiment A, the value of sl eq is already high.To remove the contact between the seamount and the ice lower surface, it would be necessary to rise the sea level by more than 20 m.This is clearly an irreversibility induced by the presence of the new contact area with the seamount in the middle of the ice shelf.After 1.5 a, the contact with the seamount below the shelf extends over no more than 300 m.Afterwards, during a transition period of 100 a, already described and shown in Fig. 10, this contact area migrates upward and grows to finally extend over 3 km.Until the merge between both grounded areas, this zone keeps the same position and size.During this time period, the value of sl eq equals nearly instantaneously one half of the value taken after 100 years.This means that a newly formed pinning point below the shelf leads very rapidly to a configuration where the ice sheet grows in volume whatever the further physically acceptable increase that is applied to sea level.

Conclusions
Motivated by the lack of predictive ability of dynamical effects within marine ice sheets, particularly when the transversal direction is considered, we modelled the grounding line dynamics in 3-D using the finite element code Elmer/Ice.The 3-D full Stokes velocity field was computed and coupled with the evolution of both the ice-air and the ice-bedrock/sea interfaces.The algorithm of the model was presented in order to highlight the fact that grounding line retreat and advance are treated as different processes.The retreat is the possible result of a contact problem, evaluated by the deviation of the ice overburden from integrated sea water pressure.The advance is based on the geometrical comparison between altitudes, and is a result of the lower free surface evolution.The model's reversibility was verified.In the corresponding experiments, the grounding line was initially completely aligned with the lateral direction before the perturbation.During the perturbation, the grounding line curved and ad-vanced at maximum over 12 km, and recovered its initial alignment and position in the main flow direction with a 3 km offset when the pertubation was removed.This result is satisfactory since our full-Stokes model performs much better comparatively to the traditional fixed grid models which exhibit at least a numerical irreversibility of tens of kilometers (Gladstone et al., 2010).
The pinning point experiment was run for a physical time period of more than 1000 a.A merge between both grounded areas occurred between 800 a and 1000 a.This indicates that pinning points may have a dramatic influence on icesheet stability by, in that case, slowing down the discharge rates observed over the recent years in West Antarctica (as in Tinto and Bell. (2011)).On the opposite, losing the contact with a pinning point, by increasing basal melting for example, might conduct to a massively enhanced discharge rate.
The complex stress pattern occurring in the vicinity of the grounding line and pinning point was discussed.The shift between longitudinal compressive and tensile stress is particularly sharp and fast around grounding lines, and has a great impact on the ice-sheet geometry, by creating an ice rumple at the vertical position above the pinning point.This shows the relevance of applying a full-Stokes model solving properly the contact between ice and bedrock, to resolve the complex problem of grounding line dynamics.
Once the contact was created between the pinning point and the ice shelf, reverting the sea level back to its initial level The Cryosphere, 6, 101-112, 2012 www.the-cryosphere.net/6/101/2012/without removing the seamount is not sufficient to release the contact.Further investigations have shown that the necessary increase of the sea level to remove the pinning point would be at least 20 times the initial decrease.This indicates that a pinning point plays a stabilising role on marine ice sheets, even if the contact lasts only a very short time.While Schoof's flux boundary condition (Schoof, 2007) is beginning to be used by the glaciological community (Pollard and DeConto, 2009), it is important to test its validity for more complex ice-sheet / ice-shelf problems, particularly when the ice mass evolves in a transient state within a non fully stabilised geometry.The model presented here is ready to be used to carry out such studies.Melting and buttressing can be considered physically, and the bedrock topography can be modified to simulate real geometries, so a longerterm perspective is to apply the model to real Antarctic outlet glaciers.

Fig. 4 .
Fig. 4. view of the bedrock geometry (zoom around the initial grounding line and the bump, the 60 km width shown does not represent the full domain), consisting of a linearly inclined area and a Gaussian bump.Altitude isovalues are displayed every 25 m.The symmetry plane at y = 50 km is shown by the dashed line.

Fig. 5 .
Fig. 5. Migration of the grounding line over time (the corresponding averaged velocities are shown in Fig. 6) for the verification experiment (see Table 2).Advance and retreat processes are shown in blue and green, respectively.(a) Top view in the (x,y) plane.(b) Positions of the main grounding line for y = 0 and y = 50 km.The new grounded area created by the pinning point is not shown in this figure.

Fig. 7 .
Fig. 7. Evolution of the ice volume over time for experiment A (blue curve), and experiment B (green curve).The inset is a zoom into the first 100 years.

Fig. 8 .
Fig. 8. Effect of the pinning point on the grounding line migration, for selected time-levels between 0 and 900 a. Isovalues of the basal velocity norm are also represented every 20 m a −1 .

Fig. 9 .
Fig. 9.Effect of the pinning point on the vertically integrated velocity ūx .Depicted profiles are along the symmetry plane y = 50 km (dashed lines) and along the lateral wall y = 0 km (solid lines), for t =0, 100, 500, 900 years.The position of the most landward grounding line for similar y-coordinate is also shown.A separate marker for the pinning point is not contained in this figure.

Figure 11 Fig. 10 .
Figure11shows 3-D views (t = 0, 100, 500, 900 a) of the modelled marine ice sheet around grounding line areas.Longitudinal deviatoric stresses are superimposed to the illustrated geometries.At the main grounding line, the ice undergoes a longitudinal tensile stress, which can be explained

Fig. 11 .
Fig.11.The top series show 3-D snapshots for experiment A (after t = 0, 100, 500, 900 a), superimposed with the longitudinal deviatoric stresses.Each snapshot is associated in the bottom series with its corresponding (x,z) plane view of the flowline at y = 50 km.An ice rumple forms on the upper surface under the effect of the contact between the pinning point and the ice lower surface.Both are vertically aligned over all the simulation until both grounded parts merge together.

Fig. 12 .
Fig. 12. Difference between normal stress and sea pressure in terms of sea level rise equivalent (slre) at the beginning of experiment B (starting similarly to experiment A), at t = 1.5 and 100 a. (top) 2-D view around the grounding line areas.(bottom) Along the flowline y = 50 km, only the grounded points are represented.Horizontal lines represent the sea level rise to impose on the model to release the pinning point (Experiment B).

Table 1 .
Numerical values of the parameters adopted for the simulations.

Table 2 .
Flowchart of the experiments carried out and their characteristics.
Fig.Main retreat process of the grounding line, which takes place when the convergence criterion of the Stokes nonlinear system falls below 10 −3 .For a grounded node, if the sea load is higher than the ice/bedrock contact force (a), then this node (red circle) becomes floating for the rest of the Stokes system convergence (b).Boundary conditions prescribed: N for Neumann and D for Dirichlet.