Parameterising the grounding line in flow-line ice sheet models

Realistic predictions of the behaviour of marine ice sheets require that models are able to robustly simulate grounding line migration. Fixed-grid ice sheet models have been shown to exhibit inconsistent and hence unreliable grounding line migration, except at very high resolution not yet achievable in whole ice sheet simulations. In this study we present several different approaches to parameterising the grounding line. These are distinguished by choices regarding the ice thickness profile from the last grounded to the first floating grid point, and how this profile impacts the gravitational driving stress and basal drag. We demonstrate that the most obvious choice of thickness parameterisation, linear interpolation from the last grounded to the first floating grid point, is not the most effective. We show that use of a grounding line parameterisation greatly improves performance, and that choice of a better grounding line parameterisation over a simpler one can bring further improvements, in terms of both accuracy and self consistent behaviour, comparable to increasing the grid resolution by factor two (i.e. doubling the number of grid points). The approach presented here to parameterising the grounding line does not completely solve the grounding line problem, however it reduces the resolution required. The parameterisations are presented in the context of a one dimensional “shelfy-stream” flow-line model, but could be extended to cope with more than one dimension and other model formulations.


Introduction
The potential for marine ice sheets such as the West Antarctic Ice Sheet (WAIS) to undergo rapid collapse (sometimes referred to as "marine ice sheet instability"), and the possible Correspondence to: R. M. Gladstone (r.gladstone@bristol.ac.uk) resulting sea level rise, has been discussed since the 1970s (e.g.Mercer, 1978).An introduction to marine ice sheet instability and recent developments in the area is given by, for example, Schoof (2007a) or Katz and Worster (2010).
In order to make model based predictions of the behaviour of marine ice sheets, ice sheet models must include a realistic representation of the motion of the grounding line.Vieli and Payne (2005) demonstrated that the grounding line in models where computations are carried out at fixed horizontal locations exhibits strongly resolution-dependent behaviour.More recently, Gladstone et al. (2010) showed that this problem can be overcome at very high resolution when using a simple parameterisation for the grounding line.However, the resolution required -a grid cell size of O(100 m) or finer -makes full ice sheet simulations prohibitive in terms of computational resource.High resolution can be achieved through adaptive mesh refinement near the grounding line (Gladstone et al., 2010;Goldberg et al., 2009;Durand et al., 2009).However, the computational cost is still significant, as is the programming time required to implement adaptivity in an existing non-adaptive model, especially for a full, three-dimensional ice sheet model.
In the current study we investigate whether adaptivity can be avoided, or at least its computational cost reduced, through implementation of a parameterisation to determine the grounding line position at sub-grid scale precision.The Grounding Line Parameteriastions (GLPs) presented here build on those of Pattyn et al. (2006) and Gladstone et al. (2010), adding not only further variations to the approach taken in those studies but also further corrections to both the gravitational driving stress and basal drag.The GLPs are intended to be usable in existing full ice sheet models, whether adaptivity is present or not.
The GLP design rationale is given below.A brief summary of the model is given in Sect.2, followed by a detailed description of the different GLPs (Sect.3).Results from a series of grounding line migration experiments (described in Published by Copernicus Publications on behalf of the European Geosciences Union.606 R. M. Gladstone et al.: Parameterising the grounding line Sect.4) using these GLPs are presented in Sect. 5 and discussed in Sect.6.

GLP design rationale
The motivation behind each of the steps involved in implementing the GLPs is discussed here, followed by a detailed description in Sect.3. The GLPs all use the flotation condition to define the grounding line position, where H is the ice thickness, b is the bedrock depth (positively downwards from sea level), ρ is the ice density, and ρ w is the density of sea water.Ice with thickness greater than flotation (ρH > ρ w b) is considered to be grounded, and ice with thickness below flotation (ρH < ρ w b) is considered to be floating.For simplicity variations in the ice density (e.g.low density firn layer) are ignored, but such variations could easily be incorporated.
In fixed grid models without a GLP, the grounding line is typically assumed to lie at the last grounded grid point (e.g.Vieli and Payne, 2005).In the current study the grounding line is allowed to lie exactly at the point of transition from grounded to floating ice, irrespective of whether this point lies at a model grid point.This sub-grid scale grounding line position is used to apply a correction to the force balance in the grid cell containing the grounding line.This is achieved using (a subset of) the following steps: 1. Thickness and bedrock profiles (i.e.values defined as a function of position) are constructed across the grid cell containing the grounding line (Sect.3.1).
2. These profiles are used along with the flotation condition to calculate the grounding line position with subgrid scale precision (Sect.3).
3. A velocity profile across the grid cell containing the grounding line is constructed by using the thickness profile from step 1 and the assumption of a linear flux profile (Sect.3.2).
4. The above profiles and grounding line position are used to apply a correction to the basal drag and the gravitational driving stress in the grid cell containing the grounding line position (Sect.3.2).
Assessment of this approach is based on convergence with grid resolution of grounding line behaviour, and on comparison of steady state grounding line position against analytic solutions (Schoof, 2007a).This provides an overall performance assessment.The individual steps outlined above are not assessed directly, however we would expect errors to be largely attributable to choice of thickness profile across the grid cell containing the grounding line (explained below).The study can be viewed as a test of the validity of the thickness profiles, and indeed a test of whether the approach of choosing a single function to calculate thickness profiles at all timesteps can be justified in practice.
The determination of the grounding line position and the modification to the forcing terms are fully consistent with the equations governing the model (Sect.2).Hence any inaccuracy in these calculations must be attributed mainly to inaccuracy in the thickness, bedrock and velocity profiles defined as part of the parameterisation in the grid cell containing the grounding line.We argue that choice of thickness profile (step 1) is most important.The default assumption that the grounding line lies at the last grounded grid point is clearly wrong, however we cannot give a robust mathematical justification for the particular thickness profiles presented here.Instead several different profiles are tested, without advance confidence that they are accurate.
The assumptions of linearity for bedrock and flux profiles are expected to cause less error than the choice of thickness profile.The linear bedrock assumption is correct in the current idealised study given that a linear bedrock profile is prescribed for the whole domain (Sect.3) and so does not contribute to error in the current study.The linear flux assumption is certainly true at steady state for the current idealised study due to the surface mass balance (SMB) being constant in time and space, hence any error due to the flux assumption is only applicable during spin up.The velocity profile is a function of thickness and flux profiles and incorporates no other source of error.

Model description
All the simulations presented here are carried out using the fixed grid ice stream ice shelf (FGSTSF) model of Gladstone et al. (2010).This is identical to the FGSTSF model of Vieli and Payne (2005) except that the higher order piecewise parabolic method (PPM) is used for thickness evolution (see Gladstone et al. (2010) for a description of the PPM method).It is a vertically integrated (vertical shear is not represented) flow-line model.The governing equations are presented below.
Conservation of mass for ice sheets, streams and shelves in the case of a single dimension, x, is given by where u is the horizontal velocity, a is the SMB and H is the ice thickness.Conservation of momentum for ice stream and shelf flow in the current study is given by here s is the height of the upper ice surface above sea level, g is acceleration due to gravity, β 2 is a basal drag coefficient, m is a constant determining the power law for basal drag, and v is the vertically averaged effective viscosity.Except where stated otherwise a linear drag law is used (m = 1).v is given by The force balance terms modified by the GLPs (Sect.3.2) are basal drag (second term on the left side of Eq. ( 3) and gravitational driving stress (right side of Eq. 3).
For the ice shelf, basal drag is removed by setting β 2 = 0.The left hand boundary of the domain represents the ice divide and has zero velocity and zero surface slope boundary conditions.The right hand boundary represents the calving front of the floating ice shelf, and a force balance boundary condition is used.Eq. ( 2) is solved explicitly using finite differences.See Gladstone et al. (2010) for a full description of the boundary condition implementation and how the above equations are solved.

Parameterising the grounding line
The GLPs (outlined in Sect.1.1) are applied at every timestep to each grid cell containing a grounding line (i.e. each grid cell that is grounded on one side and floating on the other).
Note that these GLPs are used to modify terms in the force balance equation.They only indirectly impact on thickness evolution, hence they do not affect mass conservation.The distinct GLPs are described below, but first some notation is introduced.
The subscript i is used to denote the grid point at the landward (i.e.grounded) side of a grid cell containing a grounding line, and i +1 for the grid point at it's seaward side.In the experiments presented here there will always be exactly one grid cell containing the grounding line but the GLPs all generalise without modification to the case of multiple grounding lines.
The GLPs are named (Table 1) according to the choice of one of six different thickness interpolation functions (Sect.3.1) and one of effectively four forcing corrections (Sect.3.2), giving 24 different GLPs.
To prevent the GLP equations from becoming unwieldy, a scaled dimensionless variable, λ (∈ R[0,1]), is used to express distance from the ith grid point (i.e. the last grounded grid point) in the x-increasing (i.e.seaward) direction.λ is given by where x is distance in km from the inland edge of the domain (i.e.ice divide), x i is the distance in km of the ith grid point from the inland edge of the domain, and x is the grid cell size in km.Hence λ = 0 at the last grounded grid point and λ = 1 at the first floating grid point.Using this notation, the dimensionless grounding line position is given by where x g is the grounding line position in km from the landward edge of the model domain.
The bedrock profile b(λ) is assumed to be linear across the grid cell containing the grounding line, though higher resolution bedrock data could be used if available.

Parameterising the thickness profile
Six methods for constructing a thickness profile across the the grid cell containing the grounding line are presented below.These methods are summarised in Table 1 (which also summarises forcing parameterisations), and an example illustration of them is shown in Fig. 1. Figure 1 demonstrates how closely these six thickness profiles match a very high resolution thickness profile given the following assumptions: the coarse resolution grid points have thicknesses at the same accuracy as the higher resolution simulation; the grounding line lies at the midpoint of a grid cell.Given that neither of these assumptions are true in general the performance of the different profiles cannot be predicted from Fig. 1.Instead Fig. 1 serves to illustrate the approach, and to emphasize the inaccuracy of the default assumption that the grounding line lies at the last grounded grid point.
The bedrock profile Eq. ( 7), the thickness profile Eq. (see below), and the flotation condition (Eq. 1) are solved simultaneously at the grounding line to find grounding line ice thickness, H g , grounding line bedrock depth, b g , and grounding line position, λ g .

Linear interpolation
The simplest reasonable assumption that can be made about the thickness profile across the grid cell containing the grounding line is that of linearity between the known values at grid points i and i + 1, Substituting Eqs. 7 and 8 into 1 at λ = λ g gives the grounding line position This parameterisation is abbreviated as LI, see Table 1.

Pattyn's parameterisation
Instead of making explicit assumptions about both thickness and bedrock profiles across the grid cell containing the grounding line, Pattyn et al. (2006) constructed a function of both thickness and bedrock depth, www.the-cryosphere.net/4/605/2010/The Cryosphere, 4, 605-619, 2010 and used interpolation of this function to calculate a grounding line position.With the assumption of linear bedrock, this implies a thickness profile of The grounding line equation, equivalent to Eq. ( 8) in Pattyn et al. (2006), is then This parameterisation is abbreviated as PA, see Table 1.

Linear extrapolation
From visual inspection of the thickness profile across the grounding line in very high resolution simulations (e.g.see Fig. 1) the thickness gradient changes abruptly in the vicinity of the grounding line.This and the next (cubic interpolation, Sect.3.1.4)choice of thickness profile make use of the gradients landward and seaward of the grounding line in addition to the thicknesses.
Here, linearly extrapolated thickness is used from both the grid points to the landward (i.e.upstream in simulations presented here) of the grounding line, H [up] , and to the seaward (downstream), H [do] : Substituting Eqs. 7 and 13 into 1, and Eqs.7 and 14 into 1, at λ = λ g gives two expressions for grounding line position where λ g[up] and λ g [do] are potential grounding line positions predicted by landward and seaward extrapolation respectively.Assuming that H [up] and H [do] intersect in the grid cell containing the grounding line, the landward and seaward thickness equations are combined to give the thickness profile where λ × is the point of intersection of the two extrapolation functions The grounding line position is then given by In the case that H [up] and H [do] do not intersect in the grid cell containing the grounding line, no sensible thickness profile can be constructed from H [up] and H [do] , and so linear interpolation is used instead (LI, Sect.3.1.1).LI is also used in the case that two potentially viable grounding line positions are given (i.e. ).This linear extrapolation parameterisation is abbreviated as LE, see Table 1.

Cubic interpolation
This thickness profile implements higher order interpolation using thicknesses from two grid points landward and two grid points seaward of the grounding line position instead of just one (i.e.grid points i − 1 and i + 2 are used in addition to i and i + 1).
A cubic equation for thickness is fitted across the grid cell containing the grounding line, where four constraints are required to determine the four coefficients of the cubic, a, b, c and d.Two of these are provided by setting the thickness to H i and H i+1 at grid points i and i + 1 respectively, as in the other parameterisations.
The other two are provided by setting the thickness gradients at i and i + 1 to those of the neighbouring grid cells, Substituting Eqs. 7 and 20 (with the above expressions for the coefficients) into 1 at λ = λ g gives an expression for the grounding line position where Note that upper case letters are used for the coefficients simply to emphasize that these are not the same coefficients as in equation 20 above.The cubic equation is solved as in Tuma and Walsh (1998), p7.If no real roots are found, or if more than one root is found within the grid cell containing the grounding line, linear interpolation (LI) is used instead.This parameterisation is abbreviated as CI, see Table 1.

Harmonic mean based parameterisation
The harmonic mean of two numbers a and b is given by 2ab/(a + b).In numerical heat transfer problems a function based on the harmonic mean is used to represent the effect of step changes in conductivity on heat flux at sub-grid scale precision (Patankar, 1980).Here we adopt the approach of Patankar (1980) to construct a thickness profile Substituting Eqs. 7 and 22 into 1 at λ = λ g gives an expression for the grounding line position where which is solved using the quadratic reduction formula.If no real roots are found, or if more than one root is found within the grid cell containing the grounding line, linear interpolation (LI) is used instead.This parameterisation is abbreviated as HM, see Table 1.

Second order harmonic mean based parameterisation
Replacing H with H 2 in Eqs.22 also gives a tractable thickness profile Substituting Eqs. 7 and 24 into 1 at λ = λ g gives an expression for the grounding line position where This cubic equation is solved as in Tuma and Walsh (1998), p7.If no real roots are found, or if more than one root is found within the grid cell containing the grounding line, linear interpolation (LI) is used instead.This parameterisation is abbreviated as H2, see Table 1.

Parameterising the forcing terms
In order to allow the thickness parameterisation to affect evolution, it must be allowed to influence the way in which the forcing terms are implemented in the grid cell containing the grounding line.In previous studies this impact has been implemented via the basal drag.Pattyn et al. (2006) imposed a transition zone in their model by setting the drag coefficient to be a function of distance from the grounding line.Gladstone et al. (2010) scaled the drag coefficient linearly in the grid box containing the grounding line according to the proportion of the grid box that was grounded.Here we introduce new approaches to modifying both basal drag and gravitational driving stress.
The FGSTSF model used in the current study employs a staggered grid for calculation of velocity.The forcing terms are defined on the staggered grid.This means that the forcing terms for the grid cell containing the grounding line are defined mid way between the ith and (i + 1)th grid points, which we will denote by (i + 1 2 ).

Gravitational driving stress
The gravitational driving stress, G, is given by the right hand side of Eq. 3.For the typical case that both H and s are linear across the grid box, G at grid point i + 1 2 is given by For the more general case that H and s are not linear across the grid box (and note that this is the case even for the linear thickness profile LI due to the discontinuity in s across the grounding line), G at grid point i + 1 2 can be calculated more accurately by This calculation is carried out numerically by dividing the grid cell containing the grounding line into 1000 equally sized segments and using the approximation of Eq. ( 26) for each segment.This number was chosen through experimentation.Below 100 segments errors due to numerical integration start to become measurable and above 10 4 the computation starts to impact on model run time.
It should be noted that while cumbersome (unwieldy 6th order polynomials are required in places), all the thickness profiles presented in Sect.3.1 are tractable to analytical solutions of the above integral.In practice, a computational implementation of the analytical solution was in some cases found to be highly inaccurate, due to catastrophic cancellation.
Note that this modification to the gravitational driving stress forcing term need be carried out only in the grid cell containing the grounding line (so it doesn't have a measurable impact on run time).
This profile scaling parameterisation for gravitational driving stress is abbreviated as "G" in Table 1.For example "LI GB1" uses linear interpolation to calculate a thickness profile across the grid cell containing the grounding line, uses the method described above to modify gravitational driving stress in this grid cell, and the linear basal drag correction described below.

Basal drag
All the GLPs in the current study involve modification of the basal drag term in the grid cell containing the grounding line, and assume that the basal drag is zero for the floating part of the grid cell.
The simplest parameterisation for basal drag is to scale the basal drag coefficient β 2 linearly with the fraction of grounded ice in the grid cell containing the grounding line, This linear scaling is referred to as B1, see Table 1.B1 gives a basal drag force in the grid cell containing the grounding line of −β 2 u i+ 1 2 (1 − λ g ).Given that the true velocity profile in the vicinity of the grounding line is not expected to be linear this scaling is not in general correct.If the velocity profile u(λ) across the grid cell containing the grounding line were known then a more appropriate scaling could be used, Although u(λ) is not known, given that the GLPs presented here all involve prescribing a thickness profile, and that the assumption of linear flux across the grid cell is a safer assumption than that of linear velocity a profile for u(λ) can be calculated, where q is the flux given by This profile scaling parameterisation for basal drag (Eqs.29, 30, and 31) is abbreviated to "B2", see Table 1.This approach can be taken with both linear (m = 1) and non linear (m = 1/3) drag laws used in the current study.

Experiments
The impact of the different grounding line parameterisations (GLPs) is investigated in idealised simulations.
where x is the distance from the ice divide, all distances are in metres, and b is measured positively upwards from sea level.
Determination of approach to steady state is by visual inspection of grounding line evolution plots.The simulation lengths are 35 kyr and 80 kyr for advance and retreat experiments (described below) respectively, and this is sufficient for the final grounding line position to be close to steady state in all cases.
As discussed by Gladstone et al. (2010), fixed grid grounding line models can exhibit a region containing multiple locally stable grounding line positions, and the limits of this region can be determined by 'advance' simulations (in which the grounding line must advance by more than x as steady state is approached) and "retreat" simulations.This region is a numerical artifact and converges towards zero as resolution increases (Schoof, 2007a;Gladstone et al., 2010).Both advance and retreat simulations are used in the current study, and their implementation is described in detail in Appendix A. Gladstone et al. (2010) demonstrated that when using the linear thickness GLP (FGSTSF GI in Gladstone et al. (2010), identical to LI B1 in the current study) the steady state grounding line position approaches the analytical solution (Schoof, 2007a) as resolution increases, at least to within a few kilometres, for both advance and retreat simulations.

Assessing performance
In the current study convergence of the final (close to steady state) grounding line position with resolution is quantified and plotted for performance assessment.Also, two metrics are defined that give a measure of error.The values of these error metrics with increasing resolution are assessed for all GLPs.
The convergence of final grounding line position, x gs , is assessed by plotting the change in final grounding line position with increasing resolution, x gs .For a given resolution, x, this is given by where the superscript denotes resolution.
x gs is plotted against resolution.This can be done independently for both advance and retreat simulations.
The first of the two error metrics is a quantification of the size of the region of locally stable grounding line positions (Gladstone et al., 2010).This is referred to as "Retreat minus advance" (RMA) and is defined as RMA = x gr − x ga (34) where x gr is the final grounding line position from a retreat experiment and x ga is the final grounding line position from the corresponding advance experiment.It is worth noting that x gr ≥ x ga for all simulations in the current study.
The second metric, ACC, is an attempt to measure model accuracy.Accuracy is the discrepancy between simulated and theoretical steady state grounding line positions, but the fact that there are multiple viable modelled steady state grounding line positions (the advance and retreat simulations give different predictions) makes this problematic to quantify.Here we have made the choice that our "best" prediction for a given model setup is the mid point between the predictions from advance and retreat simulations.Thus ACC is defined by where x gs is the analytic steady state grounding line position given by Schoof (2007a).This assumption is not "correct" as a measure of accuracy, but it does give a quantifiable metric that converges to a correct measure of accuracy as the RMA metric approaches zero.These metrics should not be confused with the "convergence" and "accuracy" errors defined by Gladstone et al. (2010).
Since only one steady state solution can exist for the grounding line position in a shelfy-stream model with a linear downsloping bed (Schoof, 2007a), an ideal model solution would have RMA = 0 and ACC = 0.
For each GLP, an advance and retreat simulation is carried out at each resolution, where resolution varies from x = 4.8 km and t = 0.4 years, to x = 0.3 km and t = 0.025 years.x and t decrease by a factor of 2 each time giving a total of 5 different resolutions.The GLPs are assessed by comparison of final grounding line position with the analytic solution (Schoof, 2007a), convergence of final grounding line position with resolution, and behaviour of the metrics RMA and ACC with increasing resolution.

Results
The time evolution of simulated grounding line position is analysed in Sect.5.1.A comparison is presented of the simplest GLP (LI B1, see Table 1) against the default assumption that the grounding line lies at the last grounded grid point (i.e.no parameterisation is used, henceforth referred to as "no-GLP").Aspects of the time evolution of the grounding line are then compared across GLPs.In Sect.5.2 final The Cryosphere, 4, 605-619, 2010  grounding line positions and convergence with resolution are compared across GLPs using metrics RMA and ACC.Several of the thickness parameterisations, specifically LE, CI, HM and H2, are designed to resort to LI in the case that no valid solution can be found.It is worth noting that in practice this happens extremely rarely, and we do not consider it significant, except in the case of LE.LE frequently fails to find a valid solution within the grid cell containing the grounding line, and hence frequently reverts to LI.In practice that LE gives results virtually identical to LI.

Time evolution
The time evolution of the grounding line both for the no-GLP case and for the simplest GLP (LI B1, see Table 1) is shown in Fig. 2. The orange lines indicate the analytical steady state grounding line positions (Schoof, 2007a).The grounding line position in a good advance simulation would be expected to approach the lower orange line at steady state, whereas the grounding line in a good retreat simulation should approach the upper orange line towards the end of the first phase, and the lower line towards the end of the second phase.
In all the advance simulations initially rapid advance gradually slows towards steady state (except for the no-GLP x = 4.8 km simulation, which becomes unstable and fails to complete).In both no-GLP and LI B1 cases the higher resolution final grounding line positions are closer to the ana-lytic solution than the lower resolution simulations.As found by Gladstone et al. (2010), the no-GLP simulations show errors O(100 km) whereas the LI B1 simulations show errors O(10 km) or less (errors are defined here as the difference between final grounding line position and analytic solution).The first phase of the retreat simulations shows behaviour similar to the advance simulations.In the second phase of the retreat simulations, initially rapid retreat gradually slows towards steady state, but the onset of retreat is delayed at lower resolutions.This delay can be better understood after considering the finer details of simulated grounding line evolution (see below).Most of the no-GLP simulations become unstable in retreat, with only the x = 0.3 km simulation completing successfully.This numerical instability relates to the interaction between basal drag, gravitational driving stress and the grounding line position (Gladstone et al., 2010).The errors seen in the LI B1 simulations reduce from O(100 km) to O(10 km) as resolution increases from x = 4.8 km to x = 0.3 km.
Of the ten no-GLP simulations (both advance and retreat for five different resolutions), the x = 0.3 km retreat simulation is the only one to run to completion with a smaller error (only by O(10 km)) than the equivalent LI B1 simulation.Given that most no-GLP simulations either become unstable and fail to complete or show much greater errors than the equivalent LI B1 simulations, the no-GLP choice is not a viable option and will not be considered further in this study.
A close up of grounding line evolution in an advance simulation using the LI B1 GLP is shown in Fig. 3.Although the mean rate of advance is very similar across different resolutions, the advance appears to occur in steps of size x (Fig. 3 upper panel).This behaviour would be expected of simulations without a GLP where the grounding line must always lie at a grid point.A closer inspection (Fig. 3 lower panel) shows that this behaviour is due to sudden accelerations of the grounding line as the grounding line passes a grid point, followed by gradual deceleration as the next grid point is approached.This suggests that the LI B1 GLP, whilst allowing for grounding line positions anywhere within the grid cell, does not allow for a continuous, smooth response of the grounding line position to the changing state of the system.The 'state' of the system is essentially the thickness profile of the whole simulated ice sheet, which determines the gravitational driving stress and basal drag.In other words, the grounding line resists advance (i.e.advances very slowly) until a threshold (corresponding to the grounding line passing a grid point) is passed in the evolving thickness profile, after which very rapid advance occurs.The lower frequency, higher amplitude accelerations seen in the lower resolution simulations indicate that a larger change is needed in the thickness profile to trigger grounding line accelerations.
These accelerations are also seen in both the first (not shown) and second (Fig. 4 the grid point is approached, but in both retreat and advance simulations the steepest part of the curve occurs immediately after a grid point has been passed.We suggest that the delayed onset of retreat seen in the second phase of the lower resolution retreat simulations is due to the greater change in thickness profile needed to reach the threshold for the first grounding line retreat acceleration.This is a numerical artifact closely related to the existence of a region of locally stable grounding line positions (Gladstone et al., 2010).
The grounding line evolution over the range of GLPs with resolution fixed at x = 2.4 km is shown in Fig. 5.
Use of the different GLPs does induce a spread in the results, but this spread is smaller than that induced by resolution for the LI B1 GLP.The time evolution and final positions from the advance simulations vary little (within O(10 km) of the analytic solution in all cases), but the retreat varies considerably, by O(10 2 km).The simplest GLP con- ceptually, LI B1, is one of the worst in terms of final grounding line position from the retreat simulation.The best GLP by this measure is H2 GB2.A close up of the retreat behaviour for these two GLPs is shown in Fig. 6.The sudden accelerations in grounding line motion can be seen in both LI B1 and H2 GB2 (and indeed in all the GLPs, not shown).The better performing GLP, H2 GB2, shows slightly smoother grounding line motion than the poorer LI B1.Although the time averaged retreat speed of the H2 GB2 grounding line is greater than that of the LI simulation over the time interval shown in Fig. 6, the maximum magnitude of the retreat velocity is greater in the LI B1 simulation (Fig. 6, lower panel).So the H2 GB2 simulation has smaller peak speeds but a higher mean.However, none of the GLPs completely flatten out these velocity spikes, just as none of the GLPs give accurate, matching final grounding line positions from both advance and retreat experiments.
The Cryosphere, 4, 605-619, 2010  Results are shown for all GLPs (see Table 1).Example first order and second order convergence (grey bars) are shown for comparison (note that the starting point of the grey bars is arbitrary, it is the gradient that defines the order of convergence).

Steady state grounding line position
The error metrics are plotted against grid resolution for all GLPs in Fig. 7.Both metrics (ACC and RMA, described in Sect.4.1) decay as resolution increases, typically linearly or slightly slower (by comparison to grey bars in Fig. 7).Convergence of the final grounding line position approaches first order as resolution increases (Fig. 8).
ACC appears to be converging faster at higher resolutions (Fig. 7).However, this may be due to the definition of the metric rather than being indicative of faster convergence with resolution.There are a number of possible explanations for this.The ACC metric is based on comparison to an analytic solution, which may itself contain minor errors due to assumptions made obtaining the solution (Schoof, 2007a).It is also possible that the model is converging to a location close to the analytic solution but not an exact match.It is possible that the final modelled grounding line positions are not at exact steady state (though they are close to steady state).RMA is a more reliable measure of convergence than ACC as it is purely a measure of internal consistency.
Figure 9 shows how the different forcing term corrections (B1, B2 and G, Sect.3.2) impact on performance for a specific thickness interpolation (in this case H2, Sect.3.1.6).Although the more sophisticated handling (i.e.H2 GB2) does show smaller errors according to both metrics, the impact is small, and RMA and ACC appear to converge at similar rates for the different forcing term corrections.This result is similar for other thickness interpolations (not shown), with the GB2 corrections generally giving the smallest errors and the B1 correction giving the largest errors.The differences are not large and convergence of RMA and ACC does not vary greatly.
Figure 10 shows convergence of RMA and ACC for the different thickness interpolations (Sect.3.1) when the simplest basal drag correction (B1, Sect.3.2.2) is used.The linear interpolation, LI, shows the greatest error (except at the lower resolutions where CI is worse) and the second order harmonic mean based interpolation, H2, shows the lowest error.The cubic interpolation GLP, CI, appears to converge slightly faster than the others.
The "best" GLP is H2 GB2.This gives the lowest errors at all resolutions and for both error metrics.The previously published LI B1 (Gladstone et al., 2010) gives poor performance.PA B1, based on the parameterisation of Pattyn et al. (2006), gives mid range performance.,Sect. 3.1.6).Results are shown for all forcing term corrections (see Table 1).

Non-linear drag law
The results presented so far use a linear drag law (m = 1 in Eq. 3).We now consider the impact of choice of drag law on GLP performance.Figure 11 shows RMA and ACC against resolution for all the GLPs presented in the current study, but with a non-linear basal drag law given by m = 1 3 in Eq. 3. The drag coefficient is given by β 2 = 7.624×10 6 Pa m − 1 3 s 1 3 .The results are broadly similar to using the linear drag law, though the errors are smaller by approximately a factor of two.As with the linear drag law, the metrics appear to approach zero approximately linearly with resolution (by comparison against grey bars in Fig. 11), indeed convergence may be slightly faster with the non-linear drag law.The ranking of GLPs (not shown) is generally similar to the linear drag law, with the more sophisticated forcing parameterisations giving smaller errors.This indicates that the basal drag formulation can impact on performance but not greatly on the choice of suitable GLP.

Discussion
The aim of the current study is to provide an easily implementable and computationally efficient approach to parameterising the grounding line that can reduce grounding line errors in full ice sheet models, and to justify this approach through experimentation.The GLPs presented in the current study could all be extended to two horizontal dimensions, though this might not be trivial in the case of the more sophisticated parameterisations.
It is clear that the difference between not using a GLP and using the simplest GLP (namely LI B1) is large (Sect.5.1, see also Gladstone et al. (2010)).Given the large errors and the unstable nature of grounding line retreat in a fixed grid shelfy-stream model without a GLP, use of a GLP is necessary, though which of the present GLPs to use is less clear.
In Sect.5.2 the more sophisticated GLPs were shown to give better performance than the simpler ones, but this performance difference is not as large as the difference between no GLP and the simplest GLP.The best GLP in the current Error metrics ACC and RMA (Sect.4.1) against resolution when using the non-linear drag law (Sect.5.3).Results are shown for all GLPs (see Table 1).Example first order and second order convergence (grey bars) are shown for comparison.study, H2 GB2 (Sect.3.1.6),gives errors comparable to the worst GLP, LI B1, run at twice as fine a resolution (i.e.double the number of grid points).This result holds for both the linear and non-linear drag laws.When implemented in an ice sheet model with two horizontal dimensions, use of H2 GB2 instead of LI B1 would represent a significant (at least factor 8) saving in computational resource.However, LI B1 would be easier to implement than H2 GB2 in two horizontal dimensions.Although errors at a given resolution are reduced in more sophisticated GLPs, the rate of convergence does not vary significantly across GLPs.None of the GLPs presented here can fully overcome the grounding line problem inherent to fixed grid models (Vieli and Payne, 2005): very high resolution is still needed.
The inability of the current approach to fully solve the problem suggests that either the correct interpolation function has not been found, or that the approach itself is limited.We suspect the latter.Given the excellent fit of the cubic interpolation, CI, to the high resolution profile in Fig. 1, the CI GLPs might be expected to perform better than the other GLPs.However, this is not the case, due to the quality of fit of the CI interpolation varying during model evolution.This suggests that the approach of choosing one interpolation function for thickness over the the grid cell containing the grounding line is fundamentally limited, and that such a function would itself need to evolve as the model evolves.
Another way of viewing this problem is in terms of the step like behaviour in grounding line evolution (Sect.5.1).The GLPs are intended to solve the grounding line problem by allowing the grounding line to move smoothly across the grid cells.But grounding line movement still exhibits rapid accelerations as grid points are passed, demonstrating that the grounding line problem is only partially solved using the approaches in the current study.
This behaviour is not surprising -there is no a priori reason why the thickness profile over the grid cell containing the grounding line should match one particular interpolation function.However, the default assumption that the grounding line always lies at the last grounded grid point is clearly incorrect.The capacity of the GLPs presented here to allow the grounding line to lie at any point within a grid cell is not only a conceptual improvement, but gives demonstrably better results then the default assumption.
A more accurate method of parameterising the grounding line would therefore need to use a function that evolves as the model state evolves, possibly parameterised based on detailed studies of high resolution simulations.However, given that mesh adaptivity gives a true representation of the underlying equations and has been shown to address the grounding line problem well (Gladstone et al., 2010;Durand et al., 2009) adaptivity may provide a better solution than very complicated parameterisations.
An alternative approach to parameterising the grounding line was implemented by Pollard and DeConto (2009).Two separate models for grounded and floating ice were connected across the grounding line using an ice flux boundary condition.Cross grounding line ice flux was calculated as a function of ice thickness, rate factor, basal drag, and a scaling factor to represent buttressing (Eq.29 in Schoof (2007a), see also Schoof, 2007b).This specification of flux is valid in the special case of a flow-line model for plug flow where "ice is not too cold, sliding is slow, or the ice sheet is wide" (Schoof, 2007a).Errors associated with this flux prescription method in the case of actual ice streams and ice shelves have not yet been quantified, though the assumptions are more likely to be invalid away from steady state (Schoof, 2007b).A comparison against the GLPs described in the current study, and against very high resolution simulations (possibly using adaptivity) in a real world context would form a useful further study.
The flux prescription approach described above does not address the restriction imposed by fixed grid grounding line models that the grounding line must advance or retreat in steps of one grid cell at a time (which in turn causes step changes in the basal drag).The solution of Pollard and De-Conto (2009) to this limitation was to apply the prescribed flux either at the grounding line or downstream of it, depending on a flux criterion (details in supporting materal, Pollard and DeConto (2009)).The criterion overcomes the inconsistency between advance and retreat simulations but is without rigourous physical justification.

Conclusions
A general approach to parameterising the grounding line in fixed grid ice sheet models has been presented, expanding on previous work (Pattyn et al., 2006;Gladstone et al., 2010).The approach, centred on interpolating ice thickness over the grid cell containing the grounding line, shows greater reliability and an order of magnitude improvement in simulated grounding line position compared to the default assumption that the grounding line lies at the last grounded grid point.
Twenty four grounding line parameterisations (GLPs) have been presented, and tested in a fixed grid shelfy-stream model.The performance difference between the best and worst is comparable to a doubling of resolution.The GLPs are amenable to adaptation to two horizontal dimensions, where a doubling of resolution has a large (at least factor 8) impact on computational resource.
Two of the GLPs have been previously published.The simplest GLP, LI B1 (Gladstone et al., 2010), gives poor performance compared to the other GLPs.PA B1, based on the work of Pattyn et al. (2006), gives mid range performance.The new parameterisation H2 GB2, which includes a correction to the gravitational driving stress, gives the best performance.
None of these GLPs fully solve the grounding line problem, very high resolution is still needed.This is consistent with the conclusion of Schoof (2007a) that adaptivity (or high resolution) near the grounding line is essential.A combination of adaptive mesh refinement and a GLP would provide the most computationally efficient approach to minimising grounding line errors.

Appendix A Advance and retreat simulations
Simulations are carried out in pairs: an advance simulation in which the steady state grounding line position is approached from landward, and a retreat simulation in which the steady state grounding line position is approached from seaward.
It has been shown (Gladstone et al., 2010;Durand et al., 2009) that steady state grounding line position of retreat experiments is in general seaward of the steady state grounding line position for the corresponding advance experiment, except at very high resolution when the two steady states converge.The pair of simulations is needed to compute the metrics that are used to assess the performance of the GLPs (Sect.The retreat simulations have two phases.In the first phase advance occurs and in the second phase retreat occurs.The first phase of a retreat simulation has enhanced forcing under which the grounding line will advance much further than in the corresponding advance simulation.The second phase has forcing identical to that of the corresponding advance simulation as steady state is approached in retreat.As with the advance simulations, the final grounding line position is close to steady state.
The details of the forcing modification for retreat experiments are now described.During the first phase of 30 kyr the SMB and rate factor are both modified, and these are returned linearly to their standard values over the first 10 kyr years of the second phase.A further 40 kyr with forcing constant and identical to the corresponding advance simulation are then run in order to reach the final steady state, giving a total run length of 80 kyr for the retreat simulations.
It is not strictly necessary to approach steady state in the first phase of the retreat simulations, so long as significant retreat occurs in the second phase, to fulfil the requirement that the final grounding line position in a retreat simulation is approached from seaward.
The rate factor is given by A R (t * ) = where a R is the (time varying) SMB for the retreat simulation and a is the SMB used in the corresponding advance simulation, both measured in m yr −1 Note that the forcing in both advance and retreat experiments is identical (i.e.A R = A and a R = a) and constant as the final grounding line position is approached.

Fig. 1 .
Fig. 1.Example illustration of the different thickness interpolation functions used in the grounding line parameterisations.The solid grey lines show the ice sheet profile (bedrock, lower ice surface and upper ice surface from bottom upwards) from a snapshot during the evolution of a very high resolution simulation ( x = 0.15 km).The black lines show each of the different thickness profiles (Sect.3) at lower resolution ( x = 2.4 km) for the case where the high resolution simulated grounding line position lies near the centre of the lower resolution grid box.Low resolution grid point positions are shown with vertical grey dashed lines.The LE profile is not shown as it defaults to LI in this case.The default profile corresponds to no parameterisation -the grounding line is assumed to rest at the last grounded grid point.

Fig. 3 .
Fig. 3. Close up of the LI B1 advance simulations shown in Fig. 2 showing the step like nature of grounding line advance in detail.Resolutions are 0.3 km, 0.6 km, 1.2 km, 2.4 km and 4.8 km.The horizontal dashed lines in the lower plot indicate grid point locations at x =4.8 km resolution between 1010 km and 1030 km.

Fig. 4 .
Fig. 4. Close up of the LI B1 retreat simulations shown in Figure 2. The horizontal dashed lines indicate grid point locations at x =4.8 km resolution between 1645 km and 1660 km (upper plot) and at x =1.2 km resolution between 1398 km and 1404 km (lower plot).

Fig. 5 .Fig. 6 .
Fig. 5. Time evolution of grounding line position for all GLPs at a resolution of 2.4 km.The final grounding line positions from the advance simulations have been extended with dashed lines to facilitate comparison against the retreat simulations.The retreat simulations for the simplest (LI B1) and "best" (H2 GB2) GLPs are shown with dashed lines.

Fig. 7 .
Fig. 7. Error metrics ACC and RMA (Sect.4.1) against resolution.Results are shown for all GLPs (see Table1).Example first order and second order convergence (grey bars) are shown for comparison (note that the starting point of the grey bars is arbitrary, it is the gradient that defines the order of convergence).
Fig. 11.Error metrics ACC and RMA (Sect.4.1) against resolution when using the non-linear drag law (Sect.5.3).Results are shown for all GLPs (see Table1).Example first order and second order convergence (grey bars) are shown for comparison.
4.1).The Cryosphere, 4, 605-619, 2010 www.the-cryosphere.net/4/605/2010/R. M. Gladstone et al.: Parameterising the grounding line 619 Both retreat and advance simulations are initialised from a uniform slab of ice 200 m thick.Advance simulations are simply spun up for 35 kyr using a constant forcing.The final grounding line position is close to steady state.
t * /10 4 ) + At * /10 4 0 ≤ t * is the (time varying) rate factor for the retreat simulation, and t * is the time in years (measured positively forward in time) from the start of the second phase of the simulation (i.e.t * = t− 35 kyr where t is the time in years from the start of the simulation).The SMB is given by a R (t * ) =    a + 0.4 t * < 0 a + 0.4(1 − t * /10 4 ) 0 ≤ t *

Table 1 .
Summary of grounding line parameterisations (GLPs) used in this study.

Table 2 .
Model inputs and parameter values.
The experimental setup is similar (though not identical) to the Marine Ice Sheet Model Intercomparison Project (MISMIP) of www.the-cryosphere.net/4/605/2010/The Cryosphere, 4, 605-619, 2010 R. M. Gladstone et al.: Parameterising the grounding line Schoof et al. (2009) experiments 1b and 2b and to Gladstone et al. (2010).The domain size is 2112 km from ice divide (left boundary of domain) to ice front (right boundary of domain).The grid point spacing, x, and the timestep, t, vary as described below (Sect.4.1).SMB is prescribed and is spatially and temporally uniform over the domain (except for the first part of the retreat experiments, see below).The rate factor A, drag coefficient β 2 , and SMB are given in Table 2.The bedrock, b, is linear and downsloping with the same gradient as in the MISMIP experiments.
Convergence of steady state grounding line position (described in Sect.4.1).x gs is plotted against resolution for all GLPs (see Table1).Example first order and second order convergence (grey bars) are shown for comparison.