Performance and applicability of a 2.5-D ice-flow model in the vicinity of a dome

. Three-dimensional ice ﬂow modelling requires a large number of computing resources and observation data, such that 2-D simulations are often preferable. However, when there is signiﬁcant lateral divergence, this must be accounted for (2.5-D models), and a ﬂow tube is considered (volume between two horizontal ﬂowlines). In the absence of velocity observations, this ﬂow tube can be derived assuming that the ﬂowlines follow the steepest slope of the surface, under a few ﬂow assumptions. This method typically consists of scanning a digital elevation model (DEM) with a moving window and computing the curvature at the centre of this window. The ability of the 2.5-D models to account properly for a 3-D state of strain and stress has not clearly been established, nor their sensitivity to the size of the scanning window and to the geometry of the ice surface, for example in the cases of sharp ridges. Here, we study the applicability of a 2.5-D


Introduction
Computing performance has continued to increase through the last decade, and 3-D numerical simulations that could not be performed a few years ago are nowadays affordable. In particular, the Stokes equations can be solved on large 3-D data sets (e.g. Gillet-Chaulet et al., 2012), such that the complete state of stress is accounted for. Nevertheless, 2-D (x, z) models are still very useful; since they are easier to handle, the computing time is at least 2 orders of magnitude less, and they need fewer observations to describe the boundary conditions. Two-dimensional models are assumed to be unchanging in the transverse direction, such that they apply well when the topography is always the same in the missing y direction. This is a reasonable assumption for the case of large ice streams, or on the flanks of a symmetric ice cap (e.g. Martín et al., 2006). Under this assumption, the ice undergoes a strain in a vertical plane only, and all the components of strain in the transverse direction vanish (plane strain case).
In this study, we describe the volume between theoretical stream lines as a "flow tube" (Fig. 1, red lines). If the ice flow is locally converging or diverging, 2-D models are no longer valid, as they do not conserve the mass in this case, and cannot account for lateral mechanical stresses imposed by the width variations of the flow tube. For example,  Figure 1. The x axis is taken along the ridge stream line, and the y axis tangent to the surface elevation contour. The scanning window is used to determine the value of R in its centre. Lateral stream lines are represented to show the non-linear widening of the flow tube in this particular case, and the corresponding accumulation surface is coloured in grey. Sergienko (2012) showed that the 2-D (x, z) models could not account for the flow of ice on a bumpy bed, because they consider no lateral variability, and recommended using such models on a width-averaged topography. It is however possible to account for variations in the width of the flow tube in the model and still express the equations with (x, z) coordinates only (Reeh, 1988); we hereafter call these models 2.5-D models. Under a few flow assumptions (the bed elevation varies slowly in the transverse direction, the walls of the flow tube are vertical, and the lateral curvature of the flowline is small), such models represent an improvement on 2-D models because they account for lateral variability while maintaining the computational speed which full 3-D models lack.
A particularly simple case is that of an axisymmetric circular dome, where the width of the flow tube increases linearly along the flowline. In any other case, there are several possible methods for determining the width of the flow tube. At a large scale, the ice velocity can be determined via interferometric synthetic-aperture radar data (Rignot et al., 2011), but this is only reliable when ice velocity is large compared to the associated error. In the interior of ice sheets, where the ice velocity is low, velocity measurements can be made using stake monitoring. If correctly positioned, the direction and position of the stakes may give the width of the flow tube (Waddington et al., 2007). Without any available velocity observations, the only possibility is to use a digital elevation model (DEM) describing the shape of the surface. The velocity measurement method is preferable when possible, since the ice does not always flow perpendicularly to the surface contour lines. However, along a divide, or along the centre line of a drainage basin, the longitudinal surface velocity reaches its minimum or its maximum value compared to the lateral surroundings; moreover, if the flowline has a slight horizontal curvature, the horizontal shear strain rate can be neglected (Reeh, 1988). The ice velocity is then oriented along the steepest surface slope and the curvature of the surface contour lines linked to the width of the flow tube; for a dome or a 2-D ridge, determining the curvature of the surface contour line is sufficient to describe the widening of the flow tube. In this case, to determine the local curvature of a surface, the neighbourhood of each point of a DEM is approximated via interpolation. The size of this neighbourhood (scanning window) can be freely chosen, but no objective criteria are currently available to make this choice. This approach has been used to compute the surface curvature of an ice sheet (Rémy and Minster, 1997;Rémy et al., 1999), but the size of the scanning window was not discussed. Due to the noise affecting the DEM and due to the regional curvature, the local computed curvature may be ambiguous.
The 2-D models that account for the width of the flow tube vary in the complexity of their physics and their underlying assumptions (e.g. Martín et al., 2009;Koutnik and Waddington, 2012;Sergienko, 2012). The 2.5-D model proposed by Reeh (1988) assumes that the streamlines are perpendicular to the surface contour lines (Fig. 1), so that the mass conservation equation directly depends on the radius of curvature of the surface contour lines (here called R). The model imposes a vertical profile for the horizontal velocity, by the use of a shape function, and the momentum conservation equation was thus not given for any case. This vertical profile is assumed to be unchanging through the x direction (columnflow model). Later on, Hvidberg (1996) improved upon this approach by setting the stress equilibrium equations without any assumption about the shape of the velocity field. Additionally, when considering a width-varying flow tube, the heat equation should be modified as well, but only a few authors have considered this necessary for their purpose (Hvidberg, 1993(Hvidberg, , 1996Pattyn, 2002).
Different authors have used the modelling approach proposed by Reeh (1988), or similar ones, at different scales and for various geometries: for example, mountain glacier (Salamatin et al., 2000;Pattyn, 2002), ice-sheet domes (Reeh and Paterson, 1988;Hvidberg et al., 1997a), ice-sheet ridges (Hvidberg et al., 2002) or even for the whole of Greenland (Reeh et al., 2002). The shape and size of the ice bodies are quite diverse, but unfortunately the validity of the 2.5-D approach, and its assumptions, for these different cases has received little attention. In particular, the transverse deformation of the ice is assumed to be constant through depth (i.e. the walls of the tube are vertical), which may depend on the surface geometry. Furthermore, the error in the computation of R is only discussed by Hvidberg et al. (1997b), who estimated the error in the calculation of R to be 15 %, and Hvidberg et al. (2001) by about 50 %. The method used to measure R is not detailed either, and we doubt that its influence is negligible.
No assumptions of the 2.5-D model specifically prohibit its use for a highly diverging tube. However, Reeh (1989) determined that the model was incapable of handling such a flow regime, but this was established for a model as- , 9, 2301-2313, 2016 www.geosci-model-dev.net/9/2301/2016/ suming a constant horizontal velocity profile. Furthermore, Reeh (1989) accounted for the spatial evolution of R with a simple linear model, based on the surface observation. These considerations suggest that the domain of applicability of 2.5-D models with regards to flow divergence and surface geometry remains to be determined. As a consequence, the applicability of the 2.5-D model should particularly be examined on dome geometries, where simple 2-D models would be unable to account for flow divergence. The flow tubes may widen by several orders of magnitude on a few tens of kilometres, especially on the sharpest ridge of the dome. The goal of this study is to perform several comparisons between 2.5-D and 3-D models, for various dome geometries and temperature conditions, to answer the following questions.
1. What is the error associated with the computation of R from the DEM? This geometric error depends largely on the size of the scanning window.
2. How well can a 3-D state of stress be accounted for by the single parameter R(x) along the flowline? This is related to the inherent error of the 2.5-D model and its underlying assumptions.
To investigate the performance of the 2.5-D model, we compare the velocity fields resulting from the 2.5-D and 3-D models, the latter being taken as a reference. To our knowledge, no such systematic comparison between the results of 3-D and 2.5-D models has been carried out. As such, we hope the present study will guide future research using 2.5-D models in different scenarios. In the following, we present the equations of the 3-D and 2.5-D models. We then run the simulations on several domes of different shapes: a circular one (axisymmetric), a slightly elongated one, and a very elongated one. We first consider the isothermal case, before moving on to investigate the effects of temperature.

Description of the 3-D model
In future work, we hope to investigate a small dome on the East Antarctic plateau. As such, we presently consider a synthetic case with similar geometric and thermal conditions.

Geometry and mesh
In order to investigate the influence of flow divergence, we model a ridge of a dome, as this results in significant divergence. We perform the present model comparison on a synthetic geometry which consists of a 15 km-radius domain, whose shape is a quarter of cylinder only, for reasons of symmetry. The initial thickness of the ice is 3239 m at the summit, the mean surface slope is around 0.6/1000 and the underlying bed is flat. The space coordinate is a (x, y, z) Cartesian system. The 3-D mesh is horizontally unstructured and vertically extruded on 10 levels. The horizontal mean spacing between the nodes is 1 km (Fig. 2).

Conservation equations
We denote the velocity vector u, with components (u, v, w) t . The stress and strain rate tensors are denoted σ and˙ respectively, and their components, σ ij and˙ ij . The deviatoric part of σ is denoted τ , and its components, τ ij . The 3-D mechanical model consists of a Stokes problem for incompressible ice of density ρ, in which the mass and momentum conservations equations are written where g is the gravitational acceleration vector. The values of the different parameters are given in Table 1. The ice is assumed to deform following Glen's generalized flow law (Glen, 1958): where τ e is the second invariant of τ . We choose a value of n = 3. The rate factor A(T ) non-linearly depends on temperature, following an Arrhenius law (Cuffey and Paterson, 2010), reflecting the fact that warm ice is softer. In Antarctica, the temperature range between the surface and the bottom can be large, and the pressure melting point is frequently reached at the bedrock. As a consequence, the viscosity of the upper ice can be 2 orders of magnitude larger than near the bedrock. To study the influence of the temperature T (expressed in Kelvin) on the performance of our 2.5-D model, we first consider isothermal ice at 245 K, and then a nonisothermal ice, for which where s and b are the elevation of the surface and the bedrock respectively. This linear temperature profile is simple but realistic enough to show the effect of a warmer ice at the bottom. For convenience, and as the bed is flat, b = 0 in the following experiments.

Boundary conditions
Since the 3-D mesh is a quarter of a cylinder, the conditions have to be set on five different boundaries, numbered from BC1 to BC5 (Fig. 2). We consider frozen ice at the bed and no flow through the lateral boundaries (symmetry condition).
Considering no sliding at the bottom and neglecting the atmospheric pressure at the surface, the boundary conditions are written as follows: BC2 : u.n| x=0 = 0, BC4 : σ .n| z=s = 0, where n is the outward-facing normal vector on the surface. Since the surface (BC4) is let free, a kinematic boundary condition for the surface s(x, y, t) has to be solved as well: where a is the accumulation rate. The velocity boundary condition on BC5 is set to control the shape of the steady-state dome. The method to create an elongated dome is similar to that of Gillet-Chaulet and Hindmarsh (2011), in which the shallow ice approximation is used to prescribe a profile of horizontal velocity at the boundary. As such, BC5 is not a physical boundary, but it allows the domain to be reduced to a reasonable extent. The velocity vector is oriented along the surface slope, and pointing outwards. It is variable with depth to the power n + 1 and vanishes at the bed. Its norm is tuned depending on its orientation, so that the shape of the dome can be elongated in one preferred direction ( Fig. 3): where θ is the angle to the edge of the domain, and α a shape parameter controlling the elongation of the dome. The following results will correspond to a stabilized steady-state geometry, for a uniform and constant accumulation a in time and space. To do so, the output mean velocity ω is tuned to balance the surface accumulation, so that it can be expressed as where is the surface area of the accumulation zone on the top boundary (grey area in Fig. 1), and W L the width of the ice flow at the downstream position x = L. In this particular case, and W L are simply equal to 1/4 of the dome surface and dome perimeter, and L is the radius of the dome. Three different values were taken for the shape parameter α: 2 (axisymmetric case, circular geometry), 3, and 6 ( Fig. 3).
The axisymmetric domain will test the performance of the model for a perfectly known case. The latter two cases will test the ability of the model to account for divergence of varying magnitudes.
3 Description of the 2.5-D model The coordinate system used by Reeh (1988) is a curvilinear coordinate system with a right-handed oriented coordinate axis. The x axis is oriented along the flowline, the z axis is vertically oriented, and the y axis is transverse to flow and tangential to a surface contour line. As we only consider here straight flowlines (linear ridge of an ice divide), the coordinate system is locally Cartesian (Fig. 1). We refer to the 2.5-D model in the (x, z) coordinate system. We now recall the assumptions made by Hvidberg (1996), partly inherited from Reeh (1988).
1. The flowlines are perpendicular to the surface contour lines.
2. The directions of the horizontal velocity components are constant with depth, which implies that the walls of the flow tube are vertical.
3. There are no shear stresses on the vertical boundaries defined by the flow tube.
4. The ice deforms according to Glen's flow law.
These assumptions together mean that the surface horizontal strain is transferred to the bottom, so that the surface contour lines and the horizontal velocity in the flow direction Geosci. Model Dev., 9,[2301][2302][2303][2304][2305][2306][2307][2308][2309][2310][2311][2312][2313]2016 www.geosci-model-dev.net/9/2301/2016/ impose the transverse stresses. Such assumptions are reasonable in the centre of an ice sheet for a slowly varying bed (Reeh, 1988). If the bedrock spatial variations are too steep, they will warp the ice-free surface so that the velocity at the base may not be parallel to the velocity at the surface (Hvidberg, 1993;Sergienko, 2012).

Geometry
The 2-D domain is taken as a vertical slice of the 3-D domain, on one of its lateral boundaries (Fig. 2). The dome being elongated, we will run the 2.5-D model along the sharpest ridge of the 3-D dome (y = 0) or perpendicular to it (x = 0).

Mechanical model
We denote the width of the considered flow tube W (x). The radius of the surface contour lines R(x) is taken as positive for diverging flow and negative for converging flow. In this model, the assumption (2) implies that W has no dependence on z, so that geometrical considerations show that W (x) is directly linked to R(x) by An axisymmetric dome leads to the simple relations R(x) = x and W (x) ∝ x. If the flow tube is diverging more than the axisymmetric flow (on a ridge for example), the corresponding tube surface is narrowed for a given output width, and leads to lower output velocities.
The following sections present the equations of mass and momentum conservation, modified to account for the divergence of the tube. As the velocity field mainly depends on the input/output balance, some authors only conserve the mass (e.g. Parrenin et al., 2004;Todd and Christoffersen, 2014) but not the momentum; this approach can in particular be sufficient for a vertically integrated model (Hvidberg et al., 1997b). Other authors do not modify the Stokes equations at all, but instead add an extra-surface mass balance term (Cook et al., 2014;Gladstone et al., 2012) which depends on the divergence of the tube. This approach has the advantage of simplicity and results in a correct output flux, but neglects the true horizontal advection of the ice. However, this can be justified for ice-sheet margins, where the ice mainly undergoes sliding. For all these mass-only conservation models, the normal lateral stress of the surrounding ice is not accounted for, since the force equilibrium is not properly modified.

Mass conservation
At the ice divide, the velocity component v and its spatial derivatives vanish for reasons of symmetry, so that there is no dependence of the strain rates on the transverse coordinate. Under the above assumptions and considering a flow tube of width W (x) (and corresponding radius R(x)), the normal strain rates in the curvilinear system are then written as (Jaeger, 1969, p. 45 If the flow tube has a constant width, the value of R is infinite and the equation corresponds to the plane strain case. For the more complete form of these expressions, see the discussion of Reeh (1988). The mass conservation then follows:

Momentum conservation
For a straight flowline, the horizontal shear strain rate is written as (Reeh, 1988) Along a divide between drainage basins, v vanishes and u attains a local maximum so that ∂u/∂y = 0. In this particular www.geosci-model-dev.net/9/2301/2016/ Geosci. Model Dev., 9, 2301-2313, 2016 case, there is no horizontal shear strain, and the only nonzero components of strain are the normal ones and the vertical strain rate˙ xz . For dome geometries, Hvidberg (1993) shows that the curvilinear coordinate system is equivalent to a cylindrical coordinate system distorted in the y direction, for which the radius of curvature R(x) describes the local distortion. The force equilibrium equations, expressed in the (x, z) Cartesian coordinate system, are inherited from their formulation in cylindrical coordinates, and are written as (Hvidberg, 1996;Jaeger, 1969, p. 123) where σ yy is known in terms of u and R(x) (Eqs. 3 and 13).

Boundary conditions
The boundary conditions are inherited from the 3-D case: no sliding, free surface, vanishing velocity at the ice divide and an imposed horizontal velocity profile downstream. Note that the value of the mean output velocity ω directly depends on , which is now given by where L is now the length of the flowline. Equations (11) and (18) together mean that the errors in the calculation of W result in errors in the prescribed output velocity. Since we consider a straight flowline, there is no transverse flow across the considered plane. The free surface equation is thus derived from Eq. (8), and is equivalent to a simple 2-D case, here given for the (y = 0) plane:

Implementation in Elmer/Ice
The modified mechanical equations are implemented in the Elmer/Ice finite element software (Gagliardini et al., 2013). The correct implementation of the mass conservation was checked by comparing different 2.5-D simulations with the Vialov-type profiles (Vialov, 1958) computed for different diverging tubes. The expression of the Vialov profile in the case of a power-law varying flow tube is presented in Appendix A.

Determination of the contour radius
To determine the radius of curvature of the surface contour lines, we first export a DEM from the surface nodes of the 3-D model. These nodes are fitted using an inverse distance weighting, with a power of 4 in order to ensure a good smoothing of the computed surface, representative of a real ice sheet. For comparability with a real case, the spatial resolution of the DEM is taken equal to 400 m, which is the resolution of the DEM resulting from the ICESat mission on Antarctica (Schutz et al., 2005). The DEM surface curvature is computed for each pixel, using a scanning window for fitting (Fig. 1). The local surface inside the window is curvefitted by a bivariate quadratic polynomial function, and the analytical curvature of this function taken as the local curvature, which depends on the size of the scanning window used. For example, along an elongated ridge (Fig. 3, middle and right), the surface contours are close to ellipses, but only a sufficiently large scanning window will be able to account for this global shape. A small scanning window may compute an approximately circular curvature. By contrast, even though a large window may lead to a more accurate value of the curvature, geometrically speaking, it is not necessarily the case that it will also lead to a more accurate velocity field, since it is difficult to know which surrounding environment is mainly influencing the ice flow. As a consequence, three different sizes of scanning window are tested here, 2.8, 6 and 10 km, thus corresponding to a width of 7, 15 and 25 pixels. For each pixel, the value of R(x) is taken as the inverse of the curvature of the contours of the fitted surface in the centre of the window. This is done using GRASS GIS software.

Protocol of comparison
We first run the 3-D transient isothermal simulation, and stop when a steady state is reached (∂s/∂t < 10 −6 m a −1 ). We then use the resulting DEM to compute the profile of R to initialize the 2.5-D model. For each dome configuration (α = 2, 3, 6) we compare the different runs chosen: (a) 3-D, (b) 2.5-D for the three scanning windows, and a fixed geometry, (c) 2.5-D for the three scanning windows, with a free surface. For the axi-symmetric case we add a true 2.5-D axisymmetric run (imposing an analytical value R = x). Then we compare the 3-D and 2.5-D results for non-isothermal ice, to investigate the influence of the temperature, especially near the base of the ice sheet.

3-D/2.5-D axisymmetric comparison
The absolute error in the ice velocity for an axisymmetric 2.5-D model (R = x) is of the order of 10 −4 m a −1 (Fig. 4a, where the black and yellow curves are almost superimposed). The observed error is a result of discretization and should tend to zero as the element sizes decrease.

3-D/2.5-D comparison
The computation of the radius of the surface contour lines is strongly influenced by the size of the scanning window. For a circular geometry, the variation of W along x should be linear, which is almost the case for the two wider windows (Fig. 5, solid lines). With this regular geometry, the larger the window, the more precise the radius, since the fitted surface will be more accurate. By contrast, the width value computed with the smaller window is less regular and underestimated by about 30 %, meaning that we cannot evaluate a certain curvature from too small a sample. When choosing the intermediate or larger scanning window, the root mean square error (RMSE) in the velocity is 9.9 or 3.1 % respectively (Fig. 4a), and is a consequence of the error in the calculation of the radius.

3-D/2.5-D comparison, along the ridge
Along a ridge, the flow tube is non-linearly diverging. For a given output width, the accumulation area is smaller than in the axisymmetric case, thus leading to lower output velocities.
With fixed geometries, it clearly appears that the velocity is underestimated for elongated domes (Fig. 4b and c, dashed lines), meaning that the local surface slope cannot explain the ice motion by itself: the ice along the ridge is also, if not mainly, pulled by the surrounding lateral ice, which moves due to a steeper surface slope. The case of the small scanning window appears to be different (Fig. 4b, red dashed line), simply because of a really poor estimation of W , and thus of the output velocity.
The downstream velocities are always quite accurate (10 % error), since they mainly depend on the tube surface calculation, incorporated into the velocity boundary condition. When releasing the surface, the surface slope slightly increases to accommodate the velocity boundary condition, and the computed velocity field is then closer to the 3-D reference. The relative error made in the downstream part of the flow is comparatively higher near the divide since the velocities are very small.
In the case of a sharp ridge, the ice surface has the shape of a circus tent (Fig. 6, dotted line), i.e. the ice slope increases when going towards the divide. This shape is clearly not physical, and is the result of a numerical artefact. Since the vertical strain rate is always of the order of a/H , the conservation equation leads to a balance between ∂u/∂x and u/R close to the divide. For highly diverging flows, u/R is much smaller than for an axisymmetric case at the same x position (Appendix B). As a consequence, ∂u/∂x should have a com- paratively much higher value. The only way for u to increase over a short distance near the divide is for the surface elevation to decrease sharply. To handle this artefact, we increased the mesh resolution near the dome, but without any successful results. This artefact does not appear for α = 3 (slightly elongated dome).
For α = 6, the tube surface is better estimated with an intermediate window, whose size is closer to the local value of R. The RMSE in the velocity is 12.1 % for the intermediate window and 44.3 % for the large window. Too large a window would consider the whole shape of the dome and lead to an underestimation of R. The amplitude of the error between the different runs shows that for sharp ridges (or highly diverging tubes) the choice of the window size is not straightforward, as a wider window increases the regularity of the velocity field, but decreases the ability to capture the local curvature. Reeh (1989), using a 2.5-D model for dating purposes, explained the large errors for the diverging tube of Camp Century partly by the simplicity of his model, especially his linear model of R. The computed velocities show important discrepancy with the oberved ones, which leads to a bad estimation of the origin position of the ice, by 200 %. The error made in the velocity field with this more complete model seems now to be small enough for tracking ice particles correctly, for example with the intermediate window size. However, for a real case of highly curved surface, it is difficult to know a priori the best window size to use and whether the divergence can be properly accounted for, in the absence of a reference as is done here.

3-D/2.5-D comparison, perpendicular to the ridge
As the divergence is much smaller perpendicular to the sharpest ridge, the velocity field is much smoother, and its spatial evolution closer to the 3-D reference than the alongridge case. The RMSE is 11.9 and 7.5 % for the intermediate and large window size respectively, for α = 3, and 13.7 and 11.1 % for α = 6; this error is slightly higher than for the circular geometry (Fig. 4d). In the case of large radius values, towards the exterior of the dome, the wider window gives the more accurate results.

Non-isothermal ice
A supplementary comparison is carried out on the sharp ridge (y = 0) for temperature varying linearly through depth. The computed velocity field towards the divide (low x values) shows a reversed vertical profile, i.e. the basal ice goes faster than the upper ice (Fig. 7, bottom). This non-physical result in 2.5-D can be explained this way: as soon as a 3-D tube diverges more than for an axisymmetric flow, the warmer basal ice is more easily laterally strained than the colder surface ice. As a consequence, the walls of the 3-D flow tube are no longer vertical, and using the 2.5-D model in such a case would violate assumption (2). This effect is particularly pronounced close to the divide, where the tube is narrower, and can be seen in the 3-D simulations as follows (Fig. 8). The stream lines going through a flux gate at x = 1000 m are tracked on a few hundred metres (x = 1200 m for α = 6, x = 1500 m for α = 3). The divergence of the stream lines depends on their depth -the tube is larger near the bottom (Fig. 8, blue curves) -and this dependence is stronger for high diverging tubes. To accommodate the lateral strain in this case, the 2.5-D model computes high horizontal velocities in the bottom, whereas the real motion is in fact mainly laterally oriented. No mesh refinement has been able to correct this problem. Since it does not happen for a constant temperature, it certainly originates from the lower viscosity of the basal ice (Fig. 8, red curves). This result also suggests that, on sharp ridges and with non-isothermal ice, working with a fixed vertical profile of velocity will prevent such unintended behaviour. This artefact may affect the results of Hvidberg et al. (2002), who study a flowline between GRIP (Summit) and North GRIP, because their model accounts for temperature. Their flowline stretches along a ridge which can be quite sharp, and for which the flow divergence is probably higher than axisymmetry. However, it is much longer than ours, and the ridge is the sharpest far from the summit, whereas it is the contrary in our case. The artefact that appears in our simulations may not be as sensitive in their case. Nevertheless, care must be taken in such cases, since the basic assumptions may not be justifiable and the model is likely to be outwith its application domain.
For reasons of continuity, the walls of the flow tube cannot be vertical in the direction perpendicular to the ridge either, but the effect is too weak to impact the computed velocity field.

Conclusions
A systematic comparison between 2.5-D and 3-D models has been presented in order to evaluate the ability of the former to accurately compute the velocity field on a small dome of an ice sheet. The error made when estimating the value of the radius of the surface contour lines is of the order of 10 % if the computation window is well chosen, though it can be comparatively higher close to the divide. The radius of curvature of the surface elevation contour lines should be determined with a sufficiently large computation window, but choosing the optimum size is not completely straightforward; in any case we suggest it should not be less than one-third of the maximum measured radius, and several windows should be tested to ensure the robustness of the results.
The 2.5-D model can be used without any specific restriction for tubes diverging less than and up to an axisymmetric flow. For isothermal ice, the model can be used with tubes diverging more than an axisymmetric flow, if the divergence is not too high. For very high divergence, the ice is in our study mainly pulled by the output boundary condition, and the resulted velocity field may be somewhat irregular, and surface geometry unphysical close to the divide. This means that, in the case of a sharp ridge, accounting for a certain state of stress via a single parameter R(x) is clearly not possible. In practice, the examples of surface topographies given in this study may serve as a reference to evaluate whether using a 2.5-D model for a real case is justifiable.
For non-isothermal ice, the tube should not diverge more than axisymmetry, because the softer basal ice would be much more easily laterally strained in the case of an elongated dome. The walls of the flow tube are therefore not vertical, which violates the model assumptions, and the corresponding horizontal velocity profile may be not physical near the divide. This has significant consequences for dating purposes: as the computed velocity field shows too small values in the upper layers, the age of the ice is overestimated, and the modelled isochrone layers are too high. In the absence of a reference 3-D solution as produced here, the velocity field may not necessarily appear to be unphysical, but this does not mean that the numerical artefact does not significantly affect the age calculation.
This study shows that the use of a 2.5-D model, which is a trade-off between 2-D and 3-D, must avoid several pitfalls, and its applicability domain appears to be significantly narrower than initially thought. The geometric error, resulting from miscalculation of R, should have less impact on the quality of the results than the inherent errors related to the violation of the 2.5-D model's assumptions, especially when the temperature varies through the ice column. The 2.5-D model was designed for divides between drainage basins, but real flowlines along sharp divides should not be modelled this way. Finally, for a real case, it has not been established how flat the bedrock should be to ensure that the assumption of verticality of the walls of the flow tube is reasonably respected, even for a slightly divergent flow. The applicability domain of the 2.5-D model for slightly horizontally curved flowlines has not been investigated here, but further uncertainties are expected in this case, for which some secondorder terms are neglected in the mass conservation.

Code availability
The presented simulations were performed using the Elmer/Ice v.7.0 rev. 7016 finite element model. The source code of the 2.5-D model has been available in the distribution since v.8.0 rev. d9d4a2f, implemented in the AIFlow solver. The link to the source code is: https://github.com/ElmerCSC/ elmerfem.
Geosci. Model Dev., 9,[2301][2302][2303][2304][2305][2306][2307][2308][2309][2310][2311][2312][2313]2016 www.geosci-model-dev.net/9/2301/2016/ Appendix A: Vialov profile for a power-law diverging tube To check the correct implementation of the mass conservation in the 2.5-D model, we hereafter compute the height of a Vialov profile corresponding to a regularly diverging flow tube. Note that such a surface is only representative of a single flowline, and not of a whole surface, as is usually done for a Vialov profile in plane strain (Vialov, 1958) or axisymmetry (Ritz, 1992). Figure 5 shows that we may approximate the shape of the flow tube by a power law depending on the x coordinate. Let consider a flow tube of width W = W L x L β . For plane flow, β = 0, for axisymmetry, β = 1, and for sharp ridges, β > 1. The volume outflow q * for a certain coordinate x may be expressed in one of two ways (Cuffey and Paterson, 2010, p. 388): where H is the ice thickness, τ b the basal shear, a the accumulation rate, L the length of the glacier, and A the rate factor of Glen's flow law. Equating the two expressions yields a · x = 2A(β + 1) n + 2 τ n b H 2 .
The following reasoning is then similar to that of Cuffey and Paterson (2010), simply modified by a (β +1) multiplier. The final expression for the ice thickness is unchanged: except for the height of the ice sheet H * at the dome (x = 0), which is now H * = 2(n + 2) 1/n ρg n 2n+2 √ L a 2A(β + 1) 1 2n+2 . (A5) This expression is consistent with the one previously derived for axisymmetry. We then use this expression to control the 2.5-D model by comparing the value of H * computed by the model with its above theoretical value.

Appendix B: Radius and surface of a power-law diverging flow tube
We consider the same flow as in Appendix A. The value of the radius R(x) is then expressed as The surface area of the tube upstream of x can be expressed as As u is more or less proportional to the upstream surface area , u/R is expected to be proportional to W L x L β .
On the contrary, one can consider that the value of ∂u/∂x should be of the order of ω/L, i.e. simply proportional to (L)/L = W L /(β + 1). Near the divide, u/R is then comparatively much smaller than ∂u/∂x for sharp ridges than for axisymmetric flows, and imbalances the corresponding mass conservation.
Author contributions. The experiments were designed by Olivier Gagliardini, Frédéric Parrenin and Joe Todd. Olivier Passalacqua carried them out, helped by Catherine Ritz for analytical developments, and by Fabien Gillet-Chaulet for the Elmer/Ice implementation. Olivier Passalacqua prepared the manuscript with contributions from all co-authors.