A moving-point approach to model shallow ice sheets: a study case with radially symmetrical ice sheets

. Predicting the evolution of ice sheets requires numerical models able to accurately track the migration of ice sheet continental margins or grounding lines. We introduce a physically based moving-point approach for the ﬂow of ice sheets based on the conservation of local masses. This allows the ice sheet margins to be tracked explicitly. Our approach is also well suited to capture waiting-time behaviour efﬁciently. A ﬁnite-difference moving-point scheme is derived and applied in a simpliﬁed context (continental radially symmetrical shallow ice approximation). The scheme, which is inexpensive, is veriﬁed by comparing the results with steady states obtained from an analytic solution and with exact moving-margin transient solutions. In both cases the scheme is able to track the position of the ice sheet margin with high accuracy.


Introduction
Ice sheets are an influential component of the climate system whose dynamics lead to changes in terms of ice thickness, ice velocity, or migration of ice sheet continental margins and grounding lines.Therefore numerical modelling of ice sheets needs accuracy not only of the physical variables but also in the position of their boundaries.However, simulating the migration of an ice sheet margin or a grounding line remains a complex task (Huybrechts et al., 1996;Vieli and Payne, 2005;Pattyn et al., 2012Pattyn et al., , 2013)).This paper introduces a moving-point method for the numerical simulation of ice sheets, especially the migration of their boundaries.In this paper we focus on the migration of continental ice sheet margins.
At the scale of an ice sheet or a glacier, ice is modelled as a flow which follows the Stokes equations of fluid flows (Stokes, 1845), even though the flow is non-Newtonian.Solving this problem at that scale is costly.A 3-D finite element model called Elmer/Ice has been developed for this purpose numerically (see Gagliardini et al., 2013, for a detailed description of Elmer/Ice).Other models take advantage of the very small aspect ratio of ice sheets and use a thin layer approximation differing only in the order of the approximation.The oldest and numerically least expensive model used for ice flow is the Shallow Ice Approximation, or SIA (Hutter, 1983).It gives an analytical formulation for horizontal velocities of ice in the sheet and for their vertically averaged counterpart.Although simple and fast, the SIA captures well the nonlinearity of the system due to shearing at large scales.However, the SIA is not designed to include basal sliding and is a poor approximation for small scales especially at the ice divide and the ice sheet margin.The SIA is, nevertheless, an excellent resource for testing numerical approaches, since moving-margin exact solutions exist in the literature (Halfar, 1981(Halfar, , 1983;;Bueler et al., 2005).
Significant efforts have been invested in ice sheet modelling.These have led ice sheet modellers to compare results obtained by various models for the same idealistic test problems.They first started by comparing results obtained with fixed-grid models for grounded ice sheets using the SIA (European Ice Sheet Modelling INiTiative (EISMINT): Huybrechts et al., 1996;Payne et al., 2000).Then the focus shifted to the simulation of grounding-line migration.The MISMIP (Marine Ice Sheet Model Intercomparison Project) and MISMIP3d projects (Pattyn et al., 2012(Pattyn et al., , 2013) ) have shown that fixed-grid methods perform poorly without high resolution around the grounding line or by enforcing the Published by Copernicus Publications on behalf of the European Geosciences Union.
flow at the grounding line using asymptotic results from a boundary layer theory (Schoof, 2007(Schoof, , 2011)).This has led ice sheet modellers to develop adaptive and moving techniques to overcome this issue.
One approach to gain high resolution is to use automated adaptive remeshing (Durand et al., 2009;Gudmundsson et al., 2012), forcing the resolution to stay high around the ice sheet margin.A related approach is to use adaptive mesh refinement (AMR) techniques, which allow improved resolution to be achieved in key spatially and temporally evolving subregions (Goldberg et al., 2009;Gladstone et al., 2010;Cornford et al., 2013;Jouvet and Gräser, 2013).However, even with AMR, the ice sheet margin still falls between grid points, although by adapting the grid to increase the resolution near the margin the accuracy is kept high.Adapting the grid is, nevertheless, an expensive procedure, as areas where refinement is needed have to be regularly re-identified.
Another possibility is to transform the moving domain.The number of grid points is kept constant in time, but the accuracy is kept high by the explicit tracking of the position of the ice sheet margin.This is done by transforming the ice domain to a fixed coordinate system via a geometric transformation.This approach has been successfully applied by Hindmarsh (1993) and Hindmarsh and Le Meur (2001) to an ice sheet along a flowline.However, it is not easily translated into two dimensions.
We consider here intrinsically moving-grid methods.As in the case of transformed grids, these methods allow explicit tracking of the ice sheet margin.There exist a number of techniques for generating the nodal movement in movinggrid methods.They can be classified into two subcategories: location-based methods and velocity-based methods (Cao et al., 2003).In location-based methods the positions of the nodes are redefined directly at each time step by a mapping from a reference grid (Budd et al., 2009).This is generally done by choosing a monitor function.This approach has been used by Goldberg et al. (2009), the main difficulty being the definition of the monitor function.In velocity-based methods, on the other hand, the movement of the nodes is defined in terms of a time-dependent velocity, which allows the nodes to be influenced by their previous position (Baines et al., 2005;Lee et al., 2015).Currently, this approach has not been applied to the dynamics of ice sheets.
In this paper, we apply a particular velocity-based movingpoint approach based on conservation of local mass fractions to continental ice sheets.The method is in the tradition of ALE (arbitrary Lagrangian-Eulerian) methods (Donea et al., 2004), with the difference that, instead of seeking a velocity intermediate between Lagrangian and Eulerian, the method uses both Eulerian and Lagrangian conservation to deduce the velocity and solution, respectively (see Baines et al., 2011, and references therein).We derive a finite-difference moving-point scheme in a simplified context and verify the approach with steady states obtained from an analytic solution and with exact moving-margin transient solutions in the case of radially symmetrical ice sheets.We show in particular that the scheme is able to track the position of the ice sheet margin accurately.The paper is organised as follows: in Sect. 2 we recall the SIA and detail the simplified context of our study, in Sect. 3 we describe our velocity-based movingpoint approach, and in Sect. 4 we verify our approach by comparison with exact solutions before concluding in Sect. 5.

Ice sheet geometry and Shallow Ice Approximation
We consider a single solid-phase ice sheet whose thickness at position (x, y) and time t is denoted by h(t, x, y).We assume that the ice sheet lies on a fixed bedrock and denote by b(x, y) the bed elevation.The surface elevation, s(t, x, y), is then obtained as (1) The evolution of ice sheet thickness is governed by the balance between the ice gained or lost on the surface, snow precipitation and surface melting, and ice flow draining ice accumulated in the interior towards the edges of the ice sheet.This is summarised in the mass balance equation where m(t, x, y) is the surface mass balance (positive for accumulation, negative for ablation), U (t, x, y) is the vector containing the vertically averaged horizontal components of the velocity of the ice, and (t) is the area where the ice sheet is located.Formally derived by Hutter (1983), the SIA is one of the most common approximations for large-scale ice sheet dynamics.Combined with Glen's flow law (Glen, 1955), the SIA provides (in the isothermal case) an analytical form for U as follows: Parameters involved in this formulation are summarised in Table 1.Regarding the exponent n > 1, its fixed value is classically set to 3 (see Cuffey and Paterson, 2010, for more details).

Radially symmetrical ice sheets
As a first step, we confine the study to limited-area ice sheets with radial symmetry, in other words, (t) = [0, r l (t)] × [0, 2π ].The ice sheet is centred on (0, 0), and r l (t) denotes the position of the ice sheet margin (edge of the ice sheet) at time t (see Fig. 1, which shows a section through the ice sheet).The radial symmetry implies that the geometry of the sheet depends only on r, so h(t, x, y) = h(t, r), s(t, x, y) = where r is the unit radial vector, and the mass balance Eq. ( 2) simplifies to (5) A symmetry condition is added at the ice divide (r = 0): and the ice sheet margin r l (t) is characterised by the Dirichlet boundary condition: We also assume that the flux of ice through the ice sheet margin is zero (no calving).Under hypotheses regarding the regularity of the ice thickness near the margin (see Calvo et al., 2002), we can differentiate Eq. ( 7) with respect to time.Using the mass balance equation Eq. ( 5) and h = 0 at the margin, ∂h ∂r and at the margin.Realistically, since h = 0 at the margin, ∂h/∂r will not be zero there, and hence This velocity (Eq.10) will be used in the moving-point approach described in the next section.

A moving-point approach
In the following paragraphs we describe the moving-point method that we use to simulate the dynamics of ice sheets in the context of Sect.2.2.This method is essentially a velocitybased (or Lagrangian) method relying on the construction of velocities for grid points at each time step.This allows the grid to move with the flow of ice.Moving points cover the domain only where the ice sheet exists, so that no grid point is wasted.Adjacent points move to preserve local mass fractions, and the movement is thus based on the physics (Blake, 2001;Baines et al., 2005Baines et al., , 2011;;Scherer and Baines, 2012;Lee et al., 2015).This conservation method has been applied to a variety of problems and is perfectly suitable for multidimensional problems (different examples are summarised in Baines et al. (2011) and references therein; see also Partridge (2013) for the special case of ice sheet dynamics).The key points of the method are given in the next paragraphs, and the numerical verification of the method is carried out in Sect. 4.

Conservation of mass fraction
Moving-point velocities are derived from the conservation of mass fractions (CMF).To apply this principle we first define the total mass of the ice sheet θ (t) as where ρ i is the constant density of ice.Since only mass fractions are considered in this paper and ρ i is assumed constant, we can omit ρ i without loss of generality.
Since the flux of ice through the ice sheet margin is assumed to be zero, any change in the total mass over the whole ice sheet is due solely to the surface mass balance m(t, r), We now introduce the principle of the conservation of mass fractions.Let r(t) be a moving point and define µ(r), the relative mass in the moving subinterval (0, r(t)), as The rate of change of r(t) is determined by keeping µ(r) independent of time for all moving subdomains of [0, r l (t)].

Trajectories of moving points
We obtain the velocity of a moving point by differentiating Eq. ( 13) with respect to time, giving Carrying out the time differentiation using Leibniz's integral rule and substituting for ∂h/∂t from the mass balance Eq. ( 5) gives with boundary conditions (Eq. 6) at r = 0. From Eqs. ( 14), (15), and ( 12), we can determine the velocity of every interior point as The point at r = 0 is located at the ice divide, which does not move during the simulation.The point at r l (t) is dedicated to the ice sheet margin, which moves with the velocity obtained in Eq. ( 10).We verify in Appendix A that the interior velocity calculated by Eq. ( 16) coincides with the boundary velocities calculated directly from the boundary conditions (see Eq. 10).

Determination of the ice thickness profile
Once the velocities dr/dt of the moving points r(t) have been found from Eq. ( 16), the points are moved in a Lagrangian manner.In addition, the total mass θ (t) is updated from Eq. ( 12).The ice thickness profile is then deduced from Eq. ( 13) as follows.Differentiating Eq. ( 13) with respect to r2 , we obtain which allows the ice thickness profile at time t to be constructed since dµ(r)/d(r 2 ) is constant in time and therefore known from the initial data.Note that the positivity of the ice thickness is preserved since µ is by definition a strictly increasing function (see Eq. 13).

Asymptotic behaviour at the ice sheet margin
As pointed out by Fowler (1992) and Calvo et al. (2002), singularities can appear with the SIA at the margin of grounded ice sheets.The singularity arises because of the vanishing of h at the margin and the steepening of the slope ∂h/∂r.Nevertheless the ice velocity U defined by Eq. ( 4) can remain finite even if the slope is infinite.We give more details on this subject in this subsection.We also detail the influence of the singularity on the movement of the ice sheet margin.At a fixed time and for points r sufficiently close to r l , we can write the ice thickness profile h(r) as the first term in a Frobenius expansion, to leading order, where φ l = O(1).If γ = 1, then h(r) is locally linear with slope φ l , but if γ < 1 the slope ∂h/∂r is unbounded.Hence in the asymptotic region near the margin, in the case where the bedrock topography b(r) is constant, from Eq. ( 4) which vanishes as r tends to r l if γ > n/(2n+1) and remains finite if γ = n/(2n + 1).Suppose that, in the evolution of the solution over time, γ (t) > n/(2n + 1) initially so that r l (t) is constant and the boundary is stationary (waiting).If r(t) follows a CMF trajectory, then, in the absence of accumulation/ablation, the velocity of the moving coordinate r(t) is given by Asymptotically, except at the boundary itself, this velocity is finite and positive, since U > 0 and its spatial derivative The Cryosphere, 10, 1-14, 2016 www.the-cryosphere.net/10/1/2016/∂U/∂r < 0 sufficiently close to the boundary, showing that the distance r l (t) − r(t) decreases with time.
It is a technical exercise to show that this property extends to cases with accumulation/ablation and with a general bedrock with a finite slope ∂b/∂r at the margin (see Partridge, 2013).The key point to notice is that the asymptotic behaviour depends on an infinite slope of h at the margin whereas b(r) always has a finite slope.

Numerics
We now implement a numerical scheme using a finitedifference method.The complete algorithm is detailed in Appendix B. In addition, we explain in Appendix B6 why our implementation respects the asymptotic behaviour of the ice sheet at its margin.

Numerical results
This section is dedicated to the verification of the numerical scheme derived from the moving-point method detailed in Sect. 3 and to the study of its behaviour.Every numerical experiment is performed with the parameter values given in Table 1.

Accurate estimation of steady ice thickness profile
We consider a surface mass balance m(r) independent of time.The steady state of an ice sheet occurs when the temporal change in ice thickness ∂h/∂t is zero.In that case, from Eq. ( 5), the following relationship is valid: with h ∞ (r) being the thickness of the steady ice sheet and U ∞ (r) its ice velocity.By integrating the previous equation and by including the boundary conditions (Eqs.6 and 7), the position of the margin r ∞ l can be obtained from If the bedrock is flat, the profile of the steady ice sheet, from Eqs. ( 21) and ( 4), is This approach was already in use in the EISMINT intercomparison project (Huybrechts et al., 1996) with the following constant-in-time surface mass balance: Eq. ( 22) has an analytical formulation with this surface mass balance.Therefore, r ∞ l is determined with machine precision by numerical root-finding algorithms (r ∞ l ≈ 579.81 km), and the profile of the steady state is accurately estimated from Eq. ( 23) by a single numerical integration ( r 0 m(s) s ds in Eq. ( 23) has an analytical form) using a composite trapezoidal rule (we take enough grid points to ensure that the error of the estimates is smaller than 0.01 m).

Runs with different initial profiles
We check the ability of the CMF method to track either advancing or retreating ice sheet margins by performing three different model runs.In each case, the numerical model has a grid with 21 points, uses the EISMINT surface mass balance, and is initialised using the following profile: For each of the three different runs, we take (a) a uniformly distributed initial grid with r l (0) = 450 km, h 0 = 1000 m, and p = 3/7; (b) an initial grid with r l (0) = 500 km and with higher resolution near the margin, h 0 = 1000 m, and p = 1; (c) a uniformly distributed initial grid with r l (0) = 600 km, h 0 = 4000 m, and p = 1/4.
The model is run for 25 000 yr with a constant time step t = 0.1 yr.
Figure 2 shows the evolution of the geometry and the overall motion of the grid points for each run.In run (a) the margin is staying at its initial position until the ice sheet is large enough and the sheet front steep enough to make it advance.Run (b) shows a retreating margin at an early stage before advancing, and run (c) captures the opposite behaviour.We also note that run (b) has no difficulty with a nonuniform initial grid and keeps the resolution high close to the margin.This stresses the flexibility of the CMF method to deal with various resolutions at the same time.
We then check the convergence of the three initial states to the same steady state.The calculated ice thickness at the ice divide and the position of the margin at the final time are compared with reference values in Table 2.In each case our numerical model has been able to approach the position of the margin with high accuracy (less than 400 m) at low resolution, as only 21 grid points have been employed.

EISMINT moving-margin experiment
We now perform the moving-margin experiment described in the EISMINT benchmark in order to both verify our numerical model in this case and compare our results with Table 2. Comparison between reference steady state described in Sect.4.1.1 and results obtained after a 25 000 yr run using 21 moving points with the EISMINT surface mass balance and initial profile described by Eq. ( 25).Exp 1(a): initial uniform grid with r l (0) = 450 km, h 0 = 1000 m and p = 3/7; Exp 1(b): initial grid with higher resolution near the margin with r l (0) = 500 km, h 0 = 1000 m, and p = 1; Exp 1(c): initial uniform grid with r l (0) = 600 km, h 0 = 4000 m, and p = 1/4.

Ice thickness at
Position of r = 0 (in m) the margin (in km)  those obtained by 2-D fixed-grid models used in Huybrechts et al. (1996).Compared with the experiments performed in Sect.4.1.2,the only differences are that we use an initial uniformly distributed grid with 28 nodes, an initial domain of length r l (0) = 450 km, and an initial ice thickness profile h(0, r) = t ×m(r), where t = 0.1 yr is the constant model time step and m(r) is given by Eq. ( 24).Then we run the model as in the EISMINT experiment for 25 000 yr to reach the steady state.
We first verify the result of our run with the steady state obtained in Sect.4.1.1.As shown in Fig. 3, absolute errors in the ice thickness profile mostly occur near the ice sheet margin, rising to 58.23 m at the last grid point (compared to an rms error of 15.71 m and an absolute error at the ice divide of 18.81 m).Regarding the ice sheet margin, its position is again well estimated (with an absolute error of only 138.5 m).
We next compare these results with results from fixed-grid models involved in the EISMINT intercomparison project.We confine our comparison to 2-D fixed-grid models as we only use radial symmetry (see Huybrechts et al., 1996).Regarding the ice thickness at the ice divide, our model result of 3005.76 m is within the range of estimation given by the intercomparison: 2982.3±26.4 m.These results are summarised in Table 3, showing that our moving-point method is able to achieve as good an equivalent estimation as classical fixed-grid methods with a small number of nodes while providing accurate tracking of the movement of the margin, in the context of a shallow grounded ice sheet.

Rates of convergence with EISMINT moving-margin experiment
We now study the rate of convergence of our method towards the reference solution in the EISMINT experiment.Rates of convergence are generally expressed in the form O(( r) γ ), with r being some mesh spacing.However, this approach is not appropriate in our case since the moving-point method has mesh spacings varying in time and space.Instead we present our estimated rate of convergence as a function of the number of grid points.The steady state from the EISMINT moving-margin experiment compared with our 25 000 yr model run with 28 nodes, uniformly distributed at the initial time.The reference profile is obtained by a numerical integration of Eq. ( 23) using a composite trapezoidal rule.The error in the ice thickness occurs mostly near the ice sheet margin, as in other experiments (rms error is 15.71 m and maximum error is 58.23 m).The position of the margin is well determined as the absolute error is only 138.5 m.
We calculate the absolute error for both the margin position and ice thickness at the ice divide from the results obtained in the EISMINT framework using an initial uniformly spaced grid with n r = 20, 30, 40, 60, and 80 grid points.From those results we estimate the rate of convergence for both errors.Results are summarised in Table 4.We observe that the error for the margin position decreases at an almost quadratic rate O(n −1.95 r ) and the error in the ice thickness at the ice divide at a linear rate O(n −1.16 r ).This confirms that our CMF method is well able to track the ice sheet margin without losing accuracy in the ice thickness profile.

Steady states with non-flat bedrock
The steady-state approach of Sect.4.1.1 is still valid for an ice sheet lying on a non-flat bedrock.However, the experiments in such cases are quite limited as we only have the position of the steady margin from Eq. ( 22).Nevertheless we carry out a few experiments in this context in order to demonstrate that the CMF moving-point approach is perfectly suitable for non-flat bedrock.
We consider the following fixed bedrock elevation: As in the previous section, experiments are performed with the EISMINT surface mass balance (Eq.24).At an initial time t = 0 we prescribe a uniformly distributed grid with a margin located at r l (0) = 450 km and an initial ice thickness h(0, r) = t × m(r) for the constant time step t = 0.1 yr.The model is run for 25 000 yr.The resulting evolution of the geometry and the overall motion of the grid points are shown for a grid of 20 points in Fig. 4. Regarding the position of the margin at steady state, our run has an absolute error of 127.7 m.This is even better than the previous result obtained for a flat bedrock.
www.the-cryosphere.net/10/1/2016/The Cryosphere, 10, 1-14, 2016  We also check the convergence of the estimated margin position at steady state towards its reference value by performing the same experiment with an initial uniformly spaced grid and n r grid points, n r = 20, 30, 40, 60, and 80. Absolute errors are summarised in Table 5.As in Sect.4.1.4we observe that the absolute error for the margin position decreases at a nearly quadratic rate O(n −1.83 r ).This corroborates the ability of the moving-point method to track the ice sheet margin even for non-flat bedrocks.

Verification with time-dependent solutions
In the previous paragraphs, steady states were used to verify our numerical CMF moving-point numerical method.However these experiments did not verify the transient behaviour of the ice sheet margin.To do so, we use exact timedependent solutions.

Similarity solutions
Few exact solutions for isothermal shallow ice sheets have been derived in the literature.Most are based on the similarity solutions established by Halfar (1981Halfar ( , 1983) ) for a zero surface mass balance.Bueler et al. (2005) extended this work to non-zero surface mass balance and established a new fam-ily of similarity solutions by adopting the following parameterised form for the surface mass balance: with ε being a real parameter in the interval −1 2n+1 , +∞ .Assuming that t > 0, this leads to the following family of similarity solutions: and The total mass of such ice sheets, as defined in Eq. ( 11), is where W 1 is a constant independent of ε: ds. (32)

Results
We study in this section the accuracy of transient model runs in comparison with the time-dependent exact solutions.The initialisation of every experiment is done by using the exact time-dependent solution (Eq.28), and, at each time step, The Cryosphere, 10, 1-14, 2016 www.the-cryosphere.net/10/1/2016/Rapid changes occur in the state of the sheet at the beginning of the simulation; then the dynamics dramatically slow.The ice thickness at the ice divide decreases at a rate t −1/9 , and the position of the margin increases at a rate t 1/18 .the surface mass balance is evaluated at each moving node by using the relationship m = ε t h from Eq. ( 27).When ε is non-zero, some feedback between the surface mass balance and the ice thickness is expected (Leysinger Vieli and Gudmundsson, 2004).Each model run in this section uses a fixed time step of t = 0.01 yr.
The first experiment is conducted with the constant mass similarity solution (ε = 0) between t = 100 yr and t = 20 000 yr for the reference period.Rapid changes occur in the state of the similarity solution between t = 100 yr and t = 1000 yr; then the dynamics dramatically slow (see Fig. 5 for the evolving ice thickness profile of the similarity solution).The ice thickness at the ice divide decreases at a rate t −1/9 , and the position of the margin increases at a rate t 1/18 .We begin by analysing the results obtained with a grid made up of 100 nodes, uniformly distributed at the initial time.In terms of thickness, errors mostly occur near the ice sheet margin as is the case with fixed-grid methods (see Bueler et al., 2005).For example, at the final time t = 20 000 yr, a maximum error of 134 m in the ice thickness is obtained at the computed margin while the interior of the sheet has errors less than 10 m (see Fig. 6).We also notice that errors in the ice thickness (both maximum and rms errors) decrease as the ice sheet slows down (see Fig. 7).Regarding the margin, even if the absolute error in its position increases in time, it is kept under one kilometre (880 m at the final time t = 20 000 yr).This confirms the combined ability of our method to model accurately the evolution of the ice thickness profile and to track precisely the movement of the ice sheet margin in transient behaviour.
We then study the convergence of our scheme at a final time t = 20 000 yr when the number of grid points is increased.We perform the same analysis for ε = −1/8, 1/4, and 3/4.Rates of convergence for different errors (rms error and maximum error for ice thickness profile, absolute error for the position of the margin and the volume of the ice sheet) are summarised in Table 6.These demonstrate the ability of the scheme to achieve accurate results for the position of the margin and the ice thickness profile for transient behaviour even with a small number of nodes.

Conclusions
In this paper, we have introduced a moving-point approach for ice sheet modelling using the SIA (including non-flat bedrock) based on the conservation of local mass.From this principle we derived an efficient finite-difference movingpoint scheme.The scheme was verified by comparing results with steady states from the EISMINT benchmark (Huybrechts et al., 1996) and time-dependent solutions from Bueler et al. (2005).Accurate results have been achieved with a small number of grid points in both cases.In particular our approach has been able to track the position of the ice sheet margin with high accuracy without compromising the www.the-cryosphere.net/10/1/2016/The Cryosphere, 10, 1-14, 2016 estimation of the ice thickness profile.Hence the comparison shows that the approach has considerable potential for future investigations.
Whilst this paper uses a vertically averaged horizontal ice velocity given by the shallow ice approximation, the movingmesh scheme is independent of the form of the ice velocity used here and could be used as a solver for mass balance alongside more complex vertically integrated approximations (see e.g.Schoof and Hindmarsh, 2010).
As mentioned earlier, the conservation approach is suitable not only for 1-D cases (flowline or radial) but also for 2-D scenarios.A first application has been demonstrated in Partridge (2013) and will be the subject of a new paper.The conservation approach can also be applied to marine ice sheets.In these cases, different kinds of boundaries have to be considered: e.g.grounding line, shelf front, and continental margin.Preliminary results with the moving-point method have been obtained in Dodd (2013) The limit is consistent with the velocity of the moving margin obtained in Eq. ( 10).The same approach can be used to show that dr/dt tends to 0 when r(t) tends to the ice divide r = 0.

Appendix B: A finite-difference algorithm
The moving-point method is discretised on a radial line using finite differences on the grid ri , i = 1, . .., n r , where 0 = r1 (t) < r2 (t) < . . .< rn r −1 (t) < rn r (t) = r l (t). (B1) The approximation of h(t, r) at ri (t k ) = rk i is written h k i , and that of the ice velocity U (t, r) as U k i .The velocity of the points is represented by v k i .The symbol θ k designates the numerical approximation of the total mass, and the constant mass fractions are represented by µ i for every µ(r k i ).Before giving the formula for every quantity calculated, we give the structure of the finite-difference algorithm in Algorithm 1.

B1 Initialisation
At the initial time the user needs to provide the location of each grid point {r 0 i } and the initial ice thickness {h 0 i } there.By definition, we assume that r0 1 = 0 and h 0 n r = 0. We estimate the total mass of the ice sheet at the initial time by using a composite trapezoidal rule approximating Eq. ( 11).This gives We derive the numerical approximation for the mass fractions µ i by discretising Eq. ( 13) following the same principle:

B2 Ice velocities
We confine the algorithm to n = 3 for the creep exponent in the Glen flow law.Then Eq. ( 4), giving the ice velocity, can be expanded by using the binomial theorem 27 343 ∂(h 7/3 ) ∂r

3
. (B4) We choose to rewrite the radial form of Eq. ( 4) in this way in order to ensure that the ice velocity at the ice sheet margin computed with a finite-difference scheme can be non-zero as noted in Sect.3.4.The bedrock elevation b and its derivative are known for every location of the domain.The sign of U k i (U k 1 = 0) is obtained by calculating the sign of s k i − s k i−1 (approximating the sign of the surface slope by an upwind scheme).We also approximate the derivatives of h p for any p > 0 by an upwind scheme:

B3 Approximate nodal velocities
The velocity of interior nodes is obtained by discretising Eq. ( 16) as where the integrals in Eq. (B6) are approximated by a composite trapezoidal rule.For the velocity of the ice sheet margin, Eq. ( 10) is discretised by using a first-order upwind scheme, namely,

Table 1 .
Parameters involved in the computation of the vertically averaged horizontal components of the velocity of the ice.r), and b(x, y) = b(r).The vector U can then be written in the radial coordinate system as

Figure 1 .
Figure 1.Section of a grounded radially symmetrical ice sheet.
et al.: Moving-point approach to model shallow ice sheets in radially symmetrical cases and hence the rate of change of the total mass, θ , is given by

Figure 2 .
Figure2.Evolution of the geometry (on the left) and overall motion of the grid points (on the right) for three experiments with the EISMINT surface mass balance and initial profile described by (25).Top: initial uniform grid with r l (0) = 450 km, h 0 = 1000 m, and p = 3/7; middle: initial grid with higher resolution near the margin with r l (0) = 500 km, h 0 = 1000 m, and p = 1; bottom: initial uniform grid with r l (0) = 600 km, h 0 = 4000 m, and p = 1/4.
Figure3.The steady state from the EISMINT moving-margin experiment compared with our 25 000 yr model run with 28 nodes, uniformly distributed at the initial time.The reference profile is obtained by a numerical integration of Eq. (23) using a composite trapezoidal rule.The error in the ice thickness occurs mostly near the ice sheet margin, as in other experiments (rms error is 15.71 m and maximum error is 58.23 m).The position of the margin is well determined as the absolute error is only 138.5 m.

Figure 4 .
Figure 4. Evolution of the geometry and overall motion of the grid points for the non-flat bedrock (topography given in Eq. 26) with the EISMINT surface mass balance.At steady state, the observed error for the position of the margin is 127.7 m.

Figure 5 .
Figure 5.The reference ice sheet profile (ε = 0) obtained from transient similarity solutions (see Sect. 4.3.1) is displayed for t = 100 years, for t = 1000 years, and at 1000-year intervals thereafter.Rapid changes occur in the state of the sheet at the beginning of the simulation; then the dynamics dramatically slow.The ice thickness at the ice divide decreases at a rate t −1/9 , and the position of the margin increases at a rate t 1/18 .

Figure 6 .
Figure 6.The result obtained at final time t = 20 000yr with 100 nodes equally distributed at initial time t = 100 yr, and a fixed time step t = 0.01 yr is compared to the reference transient similarity solution with ε = 0 (see Sect. 4.3.1).A maximum error of 134 m on the ice thickness is obtained at the computed margin, while the interior of the sheet has errors less than 10 m.The position of the margin is obtained with an error of 880 m.

Figure 7 .
Figure7.Evolution of the rms error and maximum absolute error in the ice thickness, and absolute error in the position of the margin between the run obtained with 100 nodes equally distributed at initial time t = 100 yr and a fixed time step t = 0.01 y and the reference transient similarity solution with ε = 0 (seeSect.4.3.1).Errors in the ice thickness decrease as the ice sheet slows down.The errors in the position of the margin increase in time, but their evolution is slower when the dynamics are slower.

Table 3 .
Huybrechts et al., 1996)comparison results for the EIS-MINT moving-margin experiment in steady state (see Table5inHuybrechts et al., 1996)and results obtained for the same experiment with the moving-point method with an identical number of grid points n r = 28.The reference values are obtained from the accurate evaluation described in Sect.4.1.1.

Table 4 .
Estimation of absolute errors from results obtained in the EISMINT framework using the moving-point method with an initial uniformly spaced grid with n r = 20, 30, 40, 60, and 80 grid points.Rates of convergence are estimated directly from the calculated absolute errors.

Table 5 .
Estimation of absolute error from results obtained with the non-flat bedrock described by Eq. (26) using the movingpoint method with an initial uniformly spaced grid with n r = 20, 30, 40, 60, and 80 grid points.Rates of convergence are estimated directly from the calculated absolute errors.

2016/ Appendix A: Consistency of the moving-point approach at boundaries
. However, the problem of initialisating such a model for use in real applications remains open.The incorporation of various data assimilation procedures is currently being investigated in this context.We now verify that dr/dt tends to the velocity obtained from Eq. (10) at the ice margin when r(t) tends to r l (t).Assuming the continuity of ∂h/∂r and m in the vicinity of the ice sheet The Cryosphere, 10, 1-14, 2016 www.the-cryosphere.net/10/1/