Extended enthalpy formulations in the ice flow model ISSM version 4.17: discontinuous conductivity and anisotropic SUPG

The thermal state of an ice sheet is an important control on its past and future evolution. Some parts of the ice sheet may be polythermal, leading to discontinuous properties at the cold–temperate transition surface (CTS). These discontinuities require a careful treatment in ice sheet models (ISMs). Additionally, the highly anisotropic geometry of the 3D elements in ice sheet modelling poses a problem for stabilization approaches in advection dominated problems. Here, we present extended enthalpy formulations within the finite-element Ice Sheet System Model (ISSM) that show a better performance to earlier 5 implementations. In a first polythermal-slab experiment, we found that the treatment of the discontinuous conductivities at the CTS with a geometric mean produce more accurate results compared to the arithmetic or harmonic mean. This improvement is particularly efficient when applied to coarse vertical resolutions. In a second ice dome experiment, we find that the numerical solution is sensitive to the choice of stabilization parameters in the well-established Streamline Upwind Petrov–Galerkin (SUPG) method. As standard literature values for the SUPG stabilization parameter do not account for the highly anisotropic 10 geometry of the 3D elements in ice sheet modelling, we propose a novel Anisotropic SUPG (ASUPG) formulation. This formulation circumvents the problem of high aspect-ratio by treating the horizontal and vertical directions separately in the stabilization coefficients. The ASUPG method provides accurate results for the thermodynamic equation on geometries with very small aspect ratios like ice sheets. Copyright statement. 15

Ice sheets and glaciers can exhibit a polythermal state that includes both cold (below the pressure melting point) and temperate (at the pressure melting point) domains, separated by the cold-temperate transition surface (CTS) (Blatter and Hutter, 1991). In temperate ice, the heat generated by viscous deformation leads to a change of phase (Fowler, 1984;Blatter and Hutter, 1991), hence temperate ice contains liquid water. The decrease of the ice viscosity with increasing content of liquid water 25 in temperate ice in turn enhances ice flow (Duval, 1977), especially if the temperate ice is present in basal layers, where shear deformation is largest.
Modern state-of-the-art ice sheet models (ISMs) simulate the thermal state according to the enthalpy method originally formulated in Aschwanden et al. (2012) and further developed and verified in Kleiner et al. (2015), Blatter and Greve (2015), Greve and Blatter (2016) and Hewitt and Schoof (2017). The main advantage of this formulation is the elimination of tracking 30 the CTS, as both cold and temperate ice domains are handled within one equation for the enthalpy E; temperature T and liquid water fraction ω are diagnostically computed from enthalpy. An increasing number of ice flow models are adopting an enthalpy scheme (e.g. Aschwanden et al., 2012;Brinkerhoff and Johnson, 2013;Seroussi et al., 2013;Kleiner et al., 2015;Greve and Blatter, 2016;Hoffman et al., 2018).
In ISMs, the governing thermodynamic equation are discretized, e.g. using the finite element method (FEM). Special care 35 has to be taken to the parabolic thermodynamic equation as numerical instabilities inherent to the advection component of this equation tend to occur without stabilization. When employing the FEM the standard Galerkin finite element method is often stabilized with the popular Streamline Upwind Petrov-Galerkin (SUPG) method (Brooks and Hughes, 1982). Although the SUPG method is well-established for advection-dominated problems, the optimal parameter choices are still subject of extensive research (e.g. Tezduyar and Osawa, 2000;John and Knobloch, 2007). Low aspect ratio mesh elements in the FEM 40 are particularly problematic, and error analysis is often restricted to two dimensions (e.g. John et al., 2018). Moreover, current mathematical and numerical analyses are not always general enough to apply to real-world applications (John et al., 2018).
ISMs are dealing with low aspect ratio, since the ice vertical extent (up to ∼4 km) is much smaller than its lateral extent (up to several thousands of kilometres). As a consequence, 3D elements are frequently highly anisotropic and pose a challenging problem in order to maintain stabilization. A non-optimal choice of stabilization parameters could result either in under-or 45 over-stabilization of the numerical solution. As a consequence of increasing computer power and modern models frequently relying on the FEM, Helanow and Ahlkrona (2018) investigated the accuracy and robustness of linear equal order finite elements discretization with Galerkin least-squares (GLS) stabilization on the Stokes equation system with anisotropic meshes.
They found that common literature values for this stabilization scheme perform well on simple domains. However, on more complex geometries, in particular, at the ice margin of outlet glaciers, the choice of standard parameters results in significant 50 oscillations in the vertical component of the surface velocity.
Beside the need for efficient stabilization in FEM, the phase change in the enthalpy formation leads to discontinuous thermal properties. This feature needs to be handled with care when seeking a numerical solution. Of particular concern are discontinuities of the thermal conductivity (Patankar, 1980;Voller and Swaminathan, 1993;Voller, 2001;Nield and Bejan, 2013). Kleiner et al. (2015) mentioned, that treating the discontinuous conductivity at the CTS as an arithmetic mean causes non-plausible 55 oscillations in the enthalpy solution that are visible, e.g. in a time-varying CTS elevation. Our work addresses the current lack of accuracy of the simulated vertical enthalpy profile to the analytical solution obtained with the ice flow model ISSM with a coarse vertical resolution (∆z=10 m, Kleiner et al., 2015, see Fig. 4 (upper row) therein).
We describe and analyze here recent developments designed to obtain an enthalpy formulation within the finite-element model ISSM (Ice Sheet System Model, Larour et al., 2012) that performs well over a wide range of grid aspect ratios in 60 advection dominated problems. The focus of this work is twofold: on the one hand, we focus on treatments of discontinuous conductivities at the CTS. Here, we test three formulations for the discontinuous conductivity proposed in Nield and Bejan (2013) for porous medium. On the other hand, we test SUPG formulations on thin geometries like ice sheets. Therefore, we run sensitivity experiments to test distinct parameter choices. One component of this study is the presentation of a novel anisotropic SUPG (ASUPG) method in ice sheet modelling that decouples the vertical from the horizontal direction to account for their 65 different scales. The formulations presented are extensions of the current implementations within the ice flow model ISSM (version 4.17) compared to Seroussi et al. (2013) and Kleiner et al. (2015).

Mathematical model
Let Ω(t) ⊆ R 3 be a three-dimensional domain with t ∈ [0, t max ]. The equations are given in Cartesian coordinates, in which x 70 and y are in the horizontal plane, and z is positive upward. The enthalpy balance equation reads with the specific enthalpy (internal energy) E, the ice velocity vector v = (v x , v y , v z ), the ice density i , the non-advective enthalpy flux q, and the heat source Ψ. The enthalpy field equation of the ice-water mixture depends on whether the mixture is cold (E < E pmp ) or temperate (E ≥ E pmp ), with E pmp the enthalpy at the pressure melting point. The non-advective enthalpy 75 flux in cold ice is represented by Fourier's law but replacing temperature T by E. In the temperate domain, the non-advective enthalpy flux is the sum of sensible and latent heat fluxes (e.g. Greve and Blatter, 2009, p. 239 (Hutter, 1982) with K 0 = ν/L, the latent heat of fusion L, the liquid water fraction ω and liquid water diffusivity ν. The diffusivity is assumed to be constant although it could depend on ω (Hutter, 1982). However, other approaches for the water mass flux, e.g. transport according to Darcy's law, are equally feasible (e.g. Fowler, 1984;Hewitt and Schoof, 2017). Sensible heat flux within the 85 temperate ice is assumed to be small compared to heat production due to deformation and considered as a source term in Eq. 1. Thus, with K c = k i /c i , where k i is the temperature conductivity and c i the specific heat capacity of ice and

90
where Φ is the heat production term due to deformation. The temperature dependence of the heat conductivity and specific heat capacity is neglected as well as the contribution of the liquid water conductivity to the ice/water mixture (Eq. 71 in Aschwanden et al., 2012).
In most cases, the liquid water fraction is increasing but temperature is decreasing towards the base, because of the Clausius-Clapeyron relation. Therefore, the transport of latent heat down the liquid water fraction gradient (Eq. 2) occurs against the 95 temperature gradient. However, the temperate ice conductivity K 0 remains poorly constrained as laboratory experiments and field observations are scarce. In the polythermal sided slab experiment proposed in Greve and Blatter (2009, sec. 9.3.6) the liquid water diffusivity is neglected and thus K 0 = 0. Nevertheless, numerical implementations will automatically generate some numerical diffusion (Greve, 1997). Sometimes a small diffusivity is used for numerical stability rather than physical reasons, e.g. K 0 = 10 −6 kg m −1 s −1 (Greve and Blatter, 2016). In ISMs typical ratios for K 0 /K c are between 10 −1 (Aschwanden 100 et al., 2012) and 10 −3 (Greve and Blatter, 2016). In this study, K 0 is simply varied to test its sensitivity on the polythermal structure.
Dirichlet boundary conditions are imposed at the upper surface in all setups. The type of basal boundary condition (Neumann or Dirichlet) is time dependent and follows the decision chart for local basal conditions given in Aschwanden et al. (2012).
However, the boundary conditions for the conducted experiments in this study are specified below.

Finite element formulation
In ISSM (Larour et al., 2012;Seroussi et al., 2013), the enthalpy equation (Eq. 1) is discretized with piecewise linear P1×P1 elements and stabilized using the SUPG method according to Franca et al. (2006). The stabilized finite element methods for Eq. 1 can be written as: where (·, ·) is the inner product of the Hilbert space H 1 0 (Ω) of square integrable functions and derivatives, and are zero on the domain boundary. The term S(E, w) is added to the standard variational formulation such that consistency is preserved and 115 numerical stability enhanced. There are different stabilization schemes that are usually considered (Franca et al., 2006); here we rely on the SUPG method: where K denotes an arbitrary element of the triangulation T h , τ K is a stability coefficient and (·, ·) K denotes integration over K. Please note, that for bilinear elements ∆E = 0.

120
The stabilization parameter, τ K is formulated as follows (Brooks and Hughes, 1982;Franca et al., 2006) h K is a characteristic dimension of element K (referred to as local mesh parameter), ξ is an upwind function and Pe K is the 125 local Peclet number. The usual Peclet definition is modified by including m k , which takes into account the effect of the degree of interpolation, k. For linear interpolations, m k=1 is 1/3 (Franca et al., 1992). For the velocity norm |v| we use the euclidean norm.

Anisotropic SUPG
The standard stabilization techniques were initially developed for isotropic meshes, which essentially require that the elements 130 have a similar size in all spatial directions. Once the elements become anisotropic, the local mesh parameter plays an important role in the calculation of stabilizing coefficients. Various definitions have been utilized based on e.g. the maximum edge length, minimum edge length, circumradius of an element, and the element length aligned with the upwind direction (e.g., Mittal, 2000;Knobloch, 2008;Brinkerhoff and Johnson, 2015). Apart from that, Becker and Rannacher (1995) and Blasco (2008) introduced stabilization coefficients for GLS diffusion that cover geometrical information from different spatial directions.

135
These definitions do not cover the element characteristic that stems from thin 3D elements. In ice sheet modelling, 3D meshes are generally formed by extruding vertically triangular meshes, leading to prismatic elements that are highly anisotropic since the vertical extent is typically one or two orders of magnitude smaller than the horizontal extent. Typically, 15 to 20 horizontal layers are used, with thinner layers close to the base. Considering a one-kilometer thick ice sheet, that is discretized in the horizontal direction between 0.5 km and 20 km, aspect ratios could exceed 100. Taking the maximum edge length as the local mesh parameter h K , which is a default choice for isotropic elements, would lead to over-stabilization, while taking the minimum edge length as h K would result in under-stabilization.
In order to develop a new SUPG stabilized method for anisotropic meshes, which accounts for geometrical information from the mesh, we consider a Cartesian three-dimensional mesh with prismatic elements. In doing so, we split the traditional SUPG formulation into a horizontal and vertical direction with the stabilization parameters τ horizontal K and τ vertical K , respectively.

145
Relying on the ideas for stabilization parameters in different spatial direction by Becker and Rannacher (1995) and Blasco (2008), the anisotropic SUPG (ASUPG) stabilization term S(E, w) is written as The stabilization parameters τ horizontal K and τ vertical K are similar to those calculated in Eqs. 9, 10, and 11, but the ASUPG 150 approach replaces the local mesh parameter h K with the characteristic horizontal and vertical dimension of the element K.
That means h k is replaced by h horizontal K and h vertical K in the two spatial directions. Here, both are calculated as the maximum extent of the element K in the respective directions.

Treatment of discontinuous conductivity
Since the conductivity is discontinuous at the CTS, special attention must be paid to the treatment of the effective conductivity 155 K eff in Eq. 3. The effective thermal conductivity of the solid-fluid system is related to the conductivity of the solid (ice), K c , and to the conductivity of the fluid (water), K 0 , and depends in a complex way on the geometry of the medium. In Nield and Bejan (2013), three models are proposed: 1. The effective thermal conductivity is the weighted arithmetic mean: 160 2. The effective thermal conductivity is the weighted harmonic mean: 3. The effective thermal conductivity is given by the weighted geometric mean: The weighting term θ ∈ [0, 1] indicates the volume fraction occupied by liquid water in a grid cell K. The volume fraction 165 of K is defined as the sum of the enthalpy in the temperate phase, E t = i E if E ≥ E pmp , divided by sum of the enthalpy E m = i E, with i the number of nodes of K. It follows that (a) θ is zero if E t = 0 and (b) θ is one if the whole grid cell is temperate (i.e. E m = E t ). The discontinuous conductivity model is only evaluated for elements that contain a CTS.
The applicability of the three models is controversial in the literature and depends strongly on the problem (e.g. Midttømme and Roaldset, 1999;Wang et al., 2006;Reddy and Karthikeyan, 2010;Jorand et al., 2011;Nield and Bejan, 2013;Ghanbarian 170 and Daigle, 2016). However, Nield and Bejan (2013) recommend the arithmetic mean if the heat conduction in the solid and fluid phases occurs in parallel. On the other hand, the harmonic mean is appropriate if the structure and orientation of the porous medium is such that the heat conduction takes place in series, with all of the heat flux passing through both solid and fluid.
Since heat conduction through porous media is likely a combination of both structures, a geometric mean can be interpreted as accounting for both processes as it always results in a value in between an arithmetic and harmonic mean (assuming K c = K 0 ).

175
Instead of employing a geometric mean a combination of the arithmetic and harmonic mean models may reveal comparable results for the effective conductivity (e.g. combinatory rules are used by Wang et al., 2006;Reddy and Karthikeyan, 2010).
When K c and K 0 are equal, the three models give the same effective thermal conductivity. For the limit case, where K 0 → 0, the harmonic and geometric mean imply insulating properties as K eff → 0 and no heat flux occurs across the interface; the arithmetic mean retains a non-zero flux.

Experiments
We ran several experiments with the emphasis to test our modifications in ISSM on accuracy and on stability. The discontinuous conductivity treatments are verified against an analytical solution within a polythermal slab experiment. As this experiment results effectively in a one-dimensional vertical experiment, it is not suitable to test the SUPG parameter choices. Therefore, we setup a synthetic second ice dome experiment with variations in the topography. Constants and model parameters used in 185 the experiments are summarized in Tab. 1.

Polythermal slab
We repeat the well-established polythermal sided slab experiment proposed in Greve and Blatter (2009)  in Eq. 6 is ignored. Please note that the analytical solution considers K 0 = 0 kg m −1 s −1 . In this experiment, we apply a thermal steady-state solver (i.e. ∂E/∂t = 0 in Eq. 1). Comparisons of results when applying a transient solver or a steady-state solver revealed no difference in the steady-state enthalpy profile.

Ice dome
In this experiment, a more realistic set-up than the polythermal slab experiment is considered with a three-dimensional ice dome based on the Vialov profile (Vialov, 1958). Other settings and parameters are borrowed from the EISMINT Phase 2 benchmark (Payne et al., 2000). The surface z s and bedrock z b of the entire ice sheet are defined as: with the ice thickness h(x, y), the maximum ice thickness h max , the radius r = x 2 + y 2 , the maximum extent r max , and the Glen exponent n. The summit of the ice dome is located at (x, y) = (0, 0).
In this experiment, a thermo-mechanical coupling is considered. The Glen-Steinemann power-law rheology (Steinemann, 210 1954;Glen, 1955) is used for the deformation of ice. The ice viscosity reads where A is the flow rate factor andε eff the effective strain rate (considered as the second invariant of the strain-rate tensor). A is assumed to be dependent on the temperature T * (temperature relative to the pressure melting point T pmp ) and liquid water fraction ω: where A 0 and A t 0 are constants, Q a is the activation energy for creep, and R is the gas constant. The constant A t 0 is equal to A(T * = T pmp , ω = 0). The upper bound of the water fraction ω is 0.01 to ensure validity of the flow rate factor parameterization in the temperate part with the experimental dataset (Duval, 1977;Lliboutry and Duval, 1985). For the dynamical model, we employ the higher-order Blatter-Pattyn approximation (Blatter, 1995;Pattyn, 2003). Basal 220 sliding is allowed everywhere and the basal drag, τ b , is written as: where v b,i is the basal velocity component in the horizontal plane and i = x, y and k 2 the friction coefficient. The effective pressure is defined as N = i g h. At the ice front a zero pressure boundary condition is applied as all the ice is above sea level.
A traction-free boundary condition is imposed at the ice/air interface.

225
For the thermal model, we impose a Dirichlet condition at the surface: T (x, y) = 238.15 K + 1.67 × 10 −5 K m −1 r.
The ice sheet base is subject to the decision chart presented in Aschwanden et al. (2012). In this implementation, the basal boundary condition is allowed to switch between Neumann and Dirichlet type dependending on the thermal basal conditions.
The geothermal flux, q geo , is considered spatially constant.

230
To investigate the sensitivity of over-and under-stabilization, we perform experiments with three different stabilization formulations (Tab. 2). The setup SUPG maxK is the standard SUPG setup based on the maximum edge length of an element K for the local mesh parameter h K . In contrast, the SUPG minK uses the minimum edge length as, however, recommended for anisotropic 2D meshes (Harari and Hughes, 1992;Mittal, 2000). Finally, the ASUPG is employed.
To study whether the stabilization is dependent on different mesh resolutions and the amount of advection, we vary the

Polythermal slab
The final steady-state CTS elevations for all simulations are shown in Fig. 1. For the maximum value of temperate ice conductivity (K 0 /K c = 10 −1 ), the models result in a CTS elevation around 36-39 m. With decreasing K 0 /K c , the temperate ice layer 245 thickness consistently decreases for the harmonic and geometric mean models and is almost halved for the lowest conductivity ratio K 0 /K c = 10 −5 ; the solution converges to the analytical CTS elevation for the high mesh resolution. However, for the harmonic mean, we detect a larger spread over the grid-resolutions at low K 0 /K c compared to the geometric mean. The simulations with the arithmetic mean yield a completely different picture. The range in the CTS elevation increases considerably with decreasing K 0 /K c and the analytical CTS elevation is met for the highest mesh resolution, below 2 m.

250
The steady-state results of the three conductivity models are verified with the analytical solution of the vertical enthalpy profile. Figure 2 shows the simulated vertical enthalpy profiles for ∆z = 10 and 0.5 m and the lowest conductivity ratio The results of all models agree well with the analytical solution for high resolutions. At coarser resolutions however, the simulated enthalpy profiles differ noticeably from the analytical solution for the arithmetic and the harmonic mean, while the geometric mean coincides well with the analytical solution. Please note, that the results for the harmonic mean 255 are similar to those presented in Kleiner et al. (2015) for ISSM.
The accuracy of the simulations with the lowest conductivity ratio is measured with the root-mean-square error (RMSE) to the analytical solution. The RMSE as a function of vertical resolution is shown in Fig. 3. All three models exhibit different behaviors. The arithmetic mean reveals a somewhat inconsistent behavior, while the harmonic mean shows approximately   first-order convergence as ∆z → 0. Overall, the geometric mean shows low errors, and the error remains on a similarly low 260 level even for coarse resolutions.
The different behaviors highlight the dependency of the solution on the CTS implementation details. As already identified by Kleiner et al. (2015) the usage of an arithmetic mean leads to oscillations in the enthalpy solution that are visible e.g. in a time-varying CTS elevation. Consequently, no steady-state solution is reached under these conditions. Here, when applying a steady-state solver, the solver does not converge and the CTS elevation flips between the non-linear iterations.

Ice dome
In this experiment, we explore the impact of the parameter choices in the SUPG formulation on the reliability and accuracy of the results. In Fig. 4 the simulated basal enthalpy field is shown for the lowest resolution l min = 10 km and high sliding case k 2 = 50 a m −1 for the three employed stabilization formulations. Due to symmetry, only the upper-right part of the domain is shown. As expected, the SUPG minK produces unphysical oscillations in the simulated enthalpy field. SUPG maxK and

270
ASUPG reveal a smooth result with merely minor oscillations a the ice front, where the surface slopes becomes singular. The same picture is observed along a cross-section of the ice sheet interior (Fig. 5). For the SUPG minK, the numerical oscillations in the enthalpy field are visible in the whole ice profile. The same qualitative behavior among the SUPG formulations is detected  for all employed grid resolutions and sliding cases (Fig. 6). Increasing the mesh resolution leads to a significant reduction in upstream oscillations. However, oscillations still occur close to the ice margin. This is in line with the theory that τ k must 275 vanish as grid refinement increases, and no stabilization may be necessary for sufficiently fine meshes. The amount of basal sliding, which controls the amount of advection, plays a secondary role.
Surprisingly, SUPG maxK and ASUPG are visually indistinguishable and result in qualitatively similar results. However, when re-running the polythermal slab experiment with the three SUPG formulations, distinct differences in the simulated enthalpy are obtained (Fig. 7). The simulations with ASUPG and SUPG minK both match the analytical solution 280 with RMSE=0.01 and 0.01 kJ kg −1 , respectively. The simulation with SUPG maxK deviates considerably from the analytical solution with RMSE=0.48 kJ kg −1 . Overall, we find that (1) using SUPG maxK as the local mesh parameter results in an oscillation-free enthalpy field but tends to produce far more diffusion than the other choices, (2) using SUPG minK as the local mesh parameter results in unphysically large oscillations for more complex geometries, and (3) ASUPG provided realistic solutions in all conducted experiments.

285
Our results demonstrate that choosing the stabilization parameter in a heuristic or ad-hoc manner, without knowledge of the possible effects, can impact the solution significantly. Choosing a sub-optimal value for the stabilization parameter can affect  e.g., by coupling to the evolution of the ice thickness.
Since the above-presented solutions for the ASUPG method are excellent, the parameter choices for the local mesh param- , and the velocity norm |v| are not further investigated. The velocity norm is here treated equally in both directions (Eq. 9), and no differentiation is made between the horizontal and vertical direction. Some test runs (not shown here) applying direction-dependent Euclidean norms of the velocity revealed no discernible differences to the above-presented 295 results. Additionally, in the current implementation, the local mesh parameter in the horizontal direction, h horizontal K , does not cover anisotropy of elements in the horizontal plane. However, these simplifications have so far not led to numerical problems, but might be subject to future work.

Conclusions
We presented extended enthalpy formulations within the ice flow model ISSM compared to Seroussi et al. (2013) and Kleiner 300 et al. (2015). Treating the discontinuous conductivities at the CTS as a geometric mean results in a good solution for coarse resolutions compared to the analytical solution. This treatment is an improvement compared to earlier ISSM results presented in Kleiner et al. (2015) and based on a harmonic mean.
Additionally, we tested various SUPG stabilization formulations on their ability to deal with the high aspect ratio of 3D elements in glaciological applications. We found that the traditional parameters in the SUPG stabilization coefficients are 305 susceptible to stabilization parameter choices, here the local mesh parameter which is easily adjustable. We propose a novel anisotropic SUPG (ASUPG) method that circumvents the high aspect-ratio problem in ice sheet modelling by treating the horizontal and vertical direction separately in the stabilization coefficients. The ASUPG method provides accurate results for the thermodynamic equation on geometries with very small aspect ratios like ice sheets. Competing interests. The authors declare that they have no conflict of interest.
Acknowledgements. Martin Rückamp acknowledges support of the Helmholtz Climate Initiative REKLIM (Regional Climate Change). We 315 thank Stephen Cornford and Alexander Robinson for the reviews that helped to improve the manuscript. For discussions and suggestions we thank Vadym Aizinger (University Bayreuth), Yonqi Wang (University Darmstadt) and Luca Wester (University of Erlangen).