Enthalpy benchmark experiments for numerical ice sheet models

We present benchmark experiments to test the implementation of enthalpy and the corresponding boundary conditions in numerical ice sheet models. Since we impose several assumptions on the experiment design, analytical solutions can be formulated for the proposed numerical experiments. The first experiment tests the functionality of the boundary condition scheme and the basal melt rate calculation during transient simulations. The second experiment addresses the steady-state enthalpy profile and the resulting position of the cold–temperate transition surface (CTS). For both experiments we assume ice flow in a parallel-sided slab decoupled from the thermal regime. We compare simulation results achieved by three different ice flow-models with these analytical solutions. The models agree well to the analytical solutions, if the change in conductivity between cold and temperate ice is properly considered in the model. In particular, the enthalpy gradient on the cold side of the CTS goes to zero in the limit of vanishing temperate-ice conductivity, as required from the physical jump conditions at the CTS.


Introduction
Ice sheets and glaciers can be distinguished by their thermal structure into cold, temperate and polythermal ice masses.While in cold ice the temperature is below the pressure melting point, in temperate ice the pressure melting point is reached.In temperate ice the heat generated by viscous deformation cannot give rise to temperature changes, but will be used for melting (Fowler, 1984;Blatter and Hutter, 1991).Thus temperate ice may contain a liquid water content (moisture).Polythermal ice masses contain both cold ice and temperate ice, separated by the cold-temperate transition surface (CTS, Greve, 1997a, b).The large ice sheets in Greenland and Antarctica show the polythermal structure of a Canadiantype glacier, which are mostly cold except for a temperate layer at the base (Aschwanden et al., 2012, and references therein).The liquid water inclusion in temperate ice makes this ice considerably softer than cold ice, resulting in a strong relationship between viscosity and moisture content (Duval, 1977;Lliboutry and Duval, 1985).The importance of this feature for the ice dynamics is obvious especially for temperate ice at the base where stresses are highest.
The enthalpy scheme presented in Aschwanden and Blatter (2009) and Aschwanden et al. (2012) describes temperature and water content in a consistent and energy-conserving formulation.Changes in the enthalpy are caused by changes of temperature in the cold ice part and by changes of the water content in the temperate ice part.The CTS position is implicitly given as the level-set of the pressure melting point and can be derived from the enthalpy field.Therefore, no restriction to the topology and shape of the CTS exists and there is no need to track it as in front-tracking models (e.g.Hutter et al., 1988;Blatter and Hutter, 1991;Greve, 1997a, b).Compared to the front-tracking models neither jump conditions nor kinematic conditions are required at the CTS.
The enthalpy scheme has already been used in model studies for the Greenland Ice Sheet.In the "referenceimplementation" of Aschwanden et al. (2012) the enthalpy scheme was compared to a cold ice scheme.A simplified version of the enthalpy scheme (regarding basal boundary conditions and ice rheology) was used to assess the effect of the initial thermal regime on century-scale simulations (Seroussi et al., 2013).Thus far we are lacking analytical solutions for thermo-mechanically coupled polythermal ice flow to test the enthalpy implementations in ice sheet models.
Published by Copernicus Publications on behalf of the European Geosciences Union.

T. Kleiner et al.: Enthalpy benchmark experiments
Here, two numerical experiments for the enthalpy field are presented for which analytical solutions exist.Similar to other studies on ice sheet modelling (Huybrechts et al., 1996;Bueler et al., 2005;Pattyn et al., 2012) we aim to verify the enthalpy method by comparing numerical solutions to analytical solutions under simplified boundary conditions.While artificially constructed exact solutions require additional compensatory terms to be incorporated in the numerical model (e.g.Bueler et al., 2005Bueler et al., , 2007)), the proposed experiments are chosen in a way that numerical models should be able to perform them with no or only minor modifications of their source codes.Therefore it is ensured that the models run through the same model components and execute the same code for the proposed numerical experiments and for real-world simulations.

Governing equations
Compared to thermodynamics usage, the enthalpy described in Aschwanden et al. (2012) is the specific internal energy.The work associated with changing the volume is not considered, since ice is assumed to be incompressible.The use of the name "enthalpy" is made to match other cryospheric applications (e.g.Notz and Worster, 2006).In the enthalpy approach temperature T and moisture ω are diagnostically computed from the modelled enthalpy field E (units: J kg −1 ).The following transfer rules are used: where p is the pressure, T ref is a reference temperature (to have positive values for the enthalpy for typical temperatures in glaciers), and L the latent heat of fusion.The enthalpy of the solid ice at the pressure melting point is defined as where T pmp (p) = T 0 − βp is the pressure melting point temperature, β is the Clausius-Clapeyron constant and T 0 is the melting point at standard pressure (see Table A1 for parameter values).
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 the ice density ρ i , the ice velocity vector v = (v x , v y , v z ), the conductive flux q i , and the heat source by internal deformation .The conductive flux in cold ice is represented by Fourier's law in enthalpy form with the conductivity K c = k i /c i .In temperate ice the conductive flux is composed of the sensible heat flux (caused by variations in the pressure melting point) and latent heat flux, thus (3) At the present state, K 0 is poorly constrained.To test the sensitivity of the models on this parameter, different values have been used according to Table A1.

Boundary conditions
At the upper ice surface the enthalpy is prescribed from the surface temperature with zero moisture content corresponding to a Canadian-type polythermal glacier (see Blatter and Hutter, 1991).In the following description of the basal conditions -Cold base (dry): if the glacier is cold at the base and without a basal water layer (i.e.E < E pmp and H w = 0), then The geothermal flux is the only source of heat as basal sliding and therefore frictional heating is forbidden for ice with temperatures below the pressure melting point.
The geothermal flux is assumed to be constant, thus changes of the heat storage in the underlying bedrock cannot affect the basal heat budget of the ice.
-Temperate base: if the glacier is temperate at the base without an overlying temperate ice layer, but with melting conditions at the base (i.e.E ≥ E pmp , H w > 0 and This condition applies to basal ice maintained at the enthalpy of the pressure melting point, when the geothermal flux and frictional heating (caused by sliding) exceed the heat flux away from the base into the ice.In this case the remaining heat is used for melting.
-Temperate ice at base: if the glacier is temperate at the base with an overlying temperate ice layer (i.e.E ≥ E pmp , H w > 0 and ∇T • n b = β/K c ), we let The proposed insulating Neumann boundary condition suppresses the diffusive enthalpy flux into the temperate ice layer even in the case of K 0 = 0.In this case the energy at the base is balanced by the basal melt rate calculation.
-Cold base (wet): if the glacier is cold, but has a liquid water layer at the base that is refreezing (i.e.E < E pmp and H w > 0), then It is assumed here that the subglacial water is at the pressure melting point and the heat stored in the water layer does not allow the basal enthalpy to be below the pressure melting point (continuity of temperature).As a consequence, the refrozen ice has a zero water content.
Note that, in addition to the temperate base condition, E ≥ E pmp , it is necessary to check if there is a temperate layer of ice above, ∇T • n b = β/K c .
Since we are dealing with polythermal glaciers, melting of ice or refreezing of liquid water at the base plays a role.The calculated melting/refreezing rate, a b (units: m s −1 ice equivalent), obeys with the frictional heating F b due to basal sliding, the heat flux in the ice q i , and the geothermal flux q geo entering the ice at the base.
Although explicit boundary conditions for CTS are not required in the enthalpy scheme, they are used to evaluate the numerical results and to derive analytical solutions later in the text.According to Greve (1997b) melting, freezing and parallel flow conditions must be distinguished depending on the CTS velocity.The enthalpy method allows for all three conditions in general.However the basal boundary conditions used in Aschwanden et al. (2012) only permit melting conditions, as Eq. ( 6) inhibits the increase of enthalpy towards the CTS.Further, the numerical models applied here do not allow a discontinuous enthalpy solution in the case of freezing or parallel flow conditions (ω > 0 at the CTS).
In the case of melting conditions at the CTS, the total enthalpy flux (advective and diffusive) at both sides of the CTS must be equal: where the superscripts "+" and "−" denote the cold and the temperate side of the interface, respectively, and n + CTS and n − CTS are the normal vectors pointing toward the CTS.This is based on the general assumption that the total heat flux leaving a representative volume through a particular face must be identical to the flux entering the next representative volume through the same face.
At the CTS, ice at its pressure melting point and without any moisture flows into the temperate layer.Hence, the enthalpy is continuous at the CTS: Assuming a horizontal CTS and according to Eqs. ( 9) and ( 10), the enthalpy derivative at the CTS is discontinuous in the given case of K c = K 0 : The condition further implies that for K 0 → 0 the enthalpy gradient on the cold side of the CTS (+) vanishes.

Numerical models
The numerical models used here are all three-dimensional flow models including a thermal component for ice.They all allow the evolution of the ice thickness, although this is not applied here.

TIM-FD 3 (finite differences)
In the Thermocoupled Ice-flow Model (TIM-FD 3 , Kleiner and Humbert, 2014) the relevant equations are discretized using finite differences in terrain-following (sigma) coordinates.For the advective terms in Eq. ( 2) the hybrid difference scheme of Spalding (1972) is used.This scheme switches between the second-order central-difference scheme and the first-order upwind-difference scheme according to the local cell Péclet number.It allows stable numerical solutions for the advection-dominated transport in the temperate ice layer.The conductive terms in Eq. ( 2) are discretized using the second-order central-difference scheme for the second derivative, where the conductivities are evaluated midway between the grid nodes (e.g.Greve and Blatter, 2009, chap. 5.7.3).The transport due to sensible heat flux in the temperate layer is assumed to be small and considered as a source term in the model.The time stepping is performed using a semi-implicit Crank-Nicolson scheme with a constant time step.
Special attention is required for the diffusion term, since the conductivity is discontinuous at the CTS.The most straightforward procedure for obtaining the interface conductivity would be to assume a linear variation of the conductivity between nodes (arithmetic mean).However, this approach cannot handle the abrupt changes of conductivity at the CTS.We use the harmonic mean of the conductivities, as suggested by Patankar (1980, chap. 4.2.3)not only at the CTS but for all conductivities evaluated between grid nodes.

ISSM (finite element)
The open source Ice Sheet System Model (ISSM; https: //issm.jpl.nasa.gov/) is applied here.A detailed model description can be found e.g. in Larour et al. (2012).It now implements the entire set of field equations and boundary conditions of the enthalpy formulation presented by Aschwanden et al. (2012).Since Seroussi et al. (2013) et al. (2012, Fig. 5).
The enthalpy field equation is discretized using a finiteelement method with linear elements.The steady-state equation and implicit time stepping scheme, respectively, give rise to a nonlinear system.It is solved using a parallelized solver.The numerical scheme can be stabilized using artificial diffusion or streamline upwind diffusion.For best comparison to the respective analytical solutions no numerical stabilization has been used here.Jumps in heat conductivity at the CTS are being accounted for by taking a volume-weighted harmonic mean of the heat conductivities over the element, see Patankar (1980).

COMice (finite element)
Numerical solutions are obtained using the COMice model (Rückamp et al., 2010) that is based on the commercial finiteelement software COMSOL Multiphysics © (www.comsol.com).The domain is approximated by a structured triangular mesh with vertical equidistant layers.Enthalpy (Eq.2) is solved with first-order Lagrange elements stabilized with streamline diffusion.The time derivatives are discretized using the implicit backward Euler scheme.An adaptive time stepping method according to Hindmarsh et al. (2005) controls the chosen time step with respect to a given tolerance.We apply Newton's method to solve the resulting system of nonlinear algebraic equations at each time step.
The step of the conductivity from K c (E) to K 0 at the CTS is implemented using Comsol's built-in operator circumcenter(expr): The operator interpolates the enthalpy solution to the circumcentre of the mesh element to which the point belongs.In doing so, circumcenter(E) is constant on each triangle and discontinuous along the edges.Therefore the conductivity jump is located on a mesh edge.This implementation shows better and faster convergence compared to other tested methods like a Heaviside function or a smoothed Heaviside function as used in the COMSOL implementation of Aschwanden and Blatter (2009) to compute the conductivity jump at the CTS.For post-processing, the CTS position is linearly interpolated between nodes.

Experiment A: parallel-sided slab (transient)
The simulation set-up is designed to test the implementation of the basal decision chart for boundary conditions and melting rates (Aschwanden et al., 2012, Fig. 5).Depending on the different thermal situations that occur at the base, the numerical code may have to switch between Neumann and Dirichlet boundary conditions for the enthalpy and the corresponding basal melt rate calculation.The main idea of this set-up is to test the reversibility during transient simulations.The conservation of water volume is also addressed here.An initially cold ice body that runs through a warmer period with an associated build-up of a liquid water layer at the base must be able to return to its initial steady state.This requires refreezing of the liquid water at the base.To test this behaviour we assume a simple heat conducting block of ice.
A parallel-sided slab of ice of constant thickness H is considered.The surface is parallel to the bed and has a constant inclination γ = 0 • to guarantee |v| = 0 and = 0. To make the set-up basically vertical 1-D, in order to be able to consider only vertical heat transport, we impose periodic boundary conditions at the sides of the block.Hence the horizontal extension does not play a role.The geothermal flux q geo at the base is constant.All parameters and their values are listed in Table A1.
The model run is as follows: Initial phase (I): starting under cold conditions with an imposed surface temperature of T s = T s,c = −30 • C and an isothermal initial temperature field T (0, z) = T s,c the simulation is run for 100 ka.
Warming phase (II): the surface temperature is switched to T s = T s,w = −10 • C and the simulation is continued for another 50 ka.
Cooling phase (III): the surface temperature is switched back to the initial value of T s = T s,c and the simulation is continued for a further 150 ka.
Since is zero, a temperate layer of ice at the base will not form and cold ice conditions hold everywhere inside the ice.The ice thickness and vertical alignment of the block is held constant over time although a significant water layer can be built up during the warming phase.Further, the water is stored at the base and no restriction of the maximum water layer thickness is applied.

Experiment B: polythermal parallel-sided slab (steady state)
To test the numerical solution for enthalpy in a vertical ice column with ice advection, we apply the "Slab with Melting Conditions at the CTS" set-up with a known analytical solution for K 0 = 0 kg m −1 s −1 (Greve and Blatter, 2009, chap. 9.3.6).However, the knowledge about latent heat flux in temperate ice is poorly constrained as laboratory experiments and field observations are scarce.We vary the values of K 0 to highlight the effect on the resulting polythermal structure.Similar to Experiment A, a parallel-sided slab of constant ice thickness H and a constant surface and bed inclination γ in the x-direction is considered (Table A1).Ice flow is decoupled from the thermal quantities by using a constant flow rate factor A. The velocity throughout the ice column is prescribed as Note that this set-up is not mass conservative, as there is no process considered that balances the accumulation rate required for a constant ice thickness.For simplicity we do not account for ice thickness evolution.The geothermal flux q geo is set to zero and basal sliding is neglected (F b = 0).Strain heating = 4 µε 2 eff is the only source of heat, where µ and εeff are the viscosity and the effective strain rate.The Glen-Steinemann power-law rheology (Glen, 1955;Steinemann, 1954) for the deformation of ice is used, thus The strain heating is largest at the base and reaches According to the assumptions in Greve and Blatter (2009, p. 246) the enthalpy conductivity K 0 in the temperate ice is zero, and the enthalpy flux at the cold site of the CTS (Eq.11) must vanish.The CTS in this experiment is uniquely determined because the vertical velocity is downward.At the ice surface (z = H ) the enthalpy is prescribed corresponding to the surface temperature T s = −3 • C and zero water content.At the ice base (z = 0) one of the boundary conditions given in Eqs. ( 4)-( 7) holds depending on the basal thermal conditions.All simulations start from a constant enthalpy corresponding to a temperature of −1.5 • C and zero water content.An analytical solution for the steady-state enthalpy profile based on the solution of Greve and Blatter (2009) is given in Appendix A2.The solution leads to a CTS position of approx.19 m above the bed.The conductivity ratio CR = K 0 /K c varies from CR = 10 −1 to 10 −5 for TIM-FD 3 and COMice and to 0 for ISSM, respectively for this set-up.The simulations are performed on vertically equidistant layers using different vertical resolutions z = (10.0,5.0, 2.0, 0.5) m.
Note that in both experiments outlined above no frictional heating at the base occurs.Drainage of moisture that exceeds a certain limit to the base needs to be considered, when a coupling of moisture to the ice viscosity is used, but is also ignored in this study.The implementation of a basal hydrology model is beyond the scope of this study, hence basal water is accumulated at the place of origin with no restriction to the water layer thickness.

Experiment A
The set-up does not allow for a temperate ice layer and therefore enthalpy variations are given only by temperature variations.The simulated basal temperatures, basal melt rates and the basal water layer thicknesses over time are shown in Fig. 1.As heat conduction is the only process of heat transfer, the vertical enthalpy profiles are linear in the steady states, which are reached at the end of each phase.At the steady states of the initial (I) and the cooling (III) phase the total vertical temperature gradient is given by the geothermal flux at the base and Eq. ( 4).This leads to the basal temperature of T (I,III) b = T s,c + H q geo /k i = −10 • C and zero melting at the base, revealed by all three models (| T | < 5 × 10 −2 • C).
In the warming phase (II) the basal temperature reaches the pressure melting point after a few thousand years and a basal water layer develops based on the basal melt rates.At the end of this phase temperatures reach the steady state (| T | < 5 × 10 −2 • C) and the basal melt rates can be calculated based on the steady-state temperature gradient between the surface and the base according to Eq. ( 8 For this setting the basal melt rate is a (II) b = 3.1 × 10 −3 m a −1 water equivalent (w.e.).The models agree well with Eq. ( 18) as shown in Fig. 1 (| a (II) b | < 10 −5 m a −1 w.e.).Phase III can be separated into two different parts: phase IIIa where the base is temperate because of the remaining basal water layer from phase II, and phase IIIb, where all subglacial water is refrozen and the base returns to cold conditions.As long as a basal water layer exists, the basal temperature is kept at pressure melting point independent of the applied surface temperature and temperature profile according to Eq. ( 7).At the end of phase IIIa, the basal melt rates can therefore be found by replacing T s,w with T s,c in Eq. ( 18).Due to the low surface temperature, refreezing conditions arise and reach steady-state values of a | < 10 −5 m a −1 w.e.).Since we have not included either a hydrology model or a reasonable upper limit for the subglacial water layer thickness, it is free to reach arbitrary thicknesses.That, in turn, is an advantage of the set-up, as we want to observe the system behaviour over longer time periods.The simulations lead to a maximum water layer thickness of ∼ 130 m that occurs a few thousand years after the end of the warming phase (II).A realistic liquid water layer thickness of about 2 m would vanish in a few time steps and would not allow for steadystate considerations at the end of IIIa.
We have chosen phase IIIa to compare not only the quasisteady-state solutions of the models at the end of each phase, but also the transient behaviour of the models compared to the analytical solution.For the comparison we use the basal melt rate instead of the temperature profile, since the correct melt rate requires a correct temperature profile and is easier to compare.In Fig. 2 the simulated basal melt rates for the first 20 ka of phase IIIa are compared to the analytical solution given in Appendix A1.After ∼ 1000 years the cold signal from the surface reaches the base and melting starts to decrease until the temperature gradient in the overlying ice does not allow for further melting and refreezing sets in.All models agree well with the analytical solution.The COMice solution is sometimes slightly below the analytical solution because of the very large time steps.The transition between melting and freezing occurs after ∼ 4684.7 years in the analytical solution.Model simulations show this transition at a comparable modelled time.
All model results clearly reveal reversibility: after the whole simulation period of 300 ka, the models return to the initial steady state at the end of phase I.

Experiment B
Here, model results of the steady-state simulations of experiment B are compared to the analytical solution given in Appendix A2.For TIM-FD 3 and COMice the steady state is assumed after 1000 model years, while in ISSM a thermal steady-state solver is applied.The final steady-state CTS positions for all simulations are shown in Fig. 3.
For the maximum value of temperate ice conductivity (CR = 10 −1 ) and the highest vertical resolution ( z = 0.5 m) the models result in CTS position slightly below 36 m.In these simulations the thickness of the temperate ice layer is almost doubled compared to the results achieved by using the smallest value of temperate ice conductivity (CR = 10 −5 ) with the same vertical resolution.The CTS positions decrease with decreasing CR and converge to the analytical so-lution.The models have approximately the same spread for the different vertical resolutions.The spread of the CTS position is smallest for CR = 10 −3 independent of the applied model.Compared to ISSM, TIM-FD 3 and COMice implementations do not allow for solving the case K 0 = 0 as in the analytical solution.
The steady-state enthalpy profiles and the corresponding temperature and moisture profiles are shown in Fig. 4 together with the analytical solution given in Appendix A2.The profiles are shown for the lowest (10 m) and highest (0.5 m) vertical resolution and the lowest conductivity ratio CR = 10 −5 used by all models.The results of all models agree well with the analytical solution for high resolutions.At coarser resolutions the simulated enthalpy profiles differ noticeably from the analytical solution.In the following we compare enthalpy differences as E = E analytic − E simulated .
In the ISSM simulation with the coarsest resolution ( z = 10 m), the enthalpy differs from the analytical solution by ∼ 1720 J kg −1 close to the CTS.This results in a temperature difference of ∼ 0.9 • C in the cold ice part.TIM-FD 3 and COMice reveal also a lower enthalpy at the cold side of the CTS compared to the analytical solution, but only to a minor extent (TIM-FD 3 : ∼ 0.2 • C, COMice: ∼ 0.1 • C).Note that the analytical solution only holds for K 0 = 0, thus small differences are expected here.
As the method chosen for interpolating heat conductivities in ISSM strongly favours the lower value, a quasi-isolating layer thicker than in the analytical solution is artificially created.Thus, heat flux into the upper cold ice column decreases, and that column cools.Vice versa, excess heat is accumulated in the lower temperate ice column, such that this part of the ice column heats up.The result is a negative temperature and positive water fraction offset.It scales with vertical mesh resolution, but stays detectable even on the highest mesh resolution tested here.
In the TIM-FD 3 simulation with the coarsest resolution ( z = 10 m), the enthalpy difference E is largest at the base (∼ 2530 J kg −1 ).As the base is temperate, this difference in the enthalpy corresponds to a difference in the basal water content of ∼ 0.8 %.With this resolution the temperate ice layer needs to be resolved within the lowermost three grid points.The slope in the profile is caused by second-order one-sided discretization (e.g.Payne and Dongelmans, 1997) of the basal boundary condition (Eq.6) in TIM-FD 3 .Compared to the FE models neither strain heating nor transport of heat is considered for basal grid nodes.
With increasing vertical resolution the maximum deviation from the analytical solution decreases for all models.For the highest resolution ( z = 0.5 m) and CR = 10 −5 the maximum differences are ∼ 150 J kg −1 , ∼ 100 J kg −1 and ∼ 10 J kg −1 for TIM-FD 3 , ISSM and COMice, respectively.The differences remain positive, thus the enthalpy is slightly underestimated.Only ISSM is able to perform the experiment with K 0 = 0 as in the analytical solution, but the maximum enthalpy difference does not further decrease.As ex- pected from Eq. (11) all models show small enthalpy gradients at the cold side of the CTS.
The observed order of accuracy measured as the rootmean-square deviation (RMSD) to the analytical solution has been obtained in the series of vertical mesh refinements from z = 10 to 0.5 m and is shown in Fig. 5.The models TIM-FD 3 and ISSM show approximately first-order convergence as z → 0, while in COMice the RMSD drops only for z below 2 m.The finite difference discretization scheme in TIM-FD 3 is formally second-order accurate in space (and time) and the finite element models ISSM and COMice use linear basis functions, thus one would expect second-order convergence as z → 0 for smooth problems.However, this is not the case here, since the observed order of accuracy depends on the strength of discontinuities (conductivity ratio between cold and temperate ice) and on the CTS implementation details.

Discussion
All three models are able to run the time-dependent experiment A and agree with the analytical solutions in terms of absolute values, timing and reversibility.However, not all types of basal boundary conditions have been tested here.Since the absence of strain heating suppresses the formation of a tem- perate ice layer at the base, the insulating boundary condition (Eq.6) could not be tested.
Beside the test of the implementation of the boundary conditions, this experiment addresses the importance of a basal water layer.Although the surface temperature changes, the basal temperature is kept at pressure melting point as long as a basal water layer exists.The amount of water at the base is crucial for the temperatures in the ice, because it acts as an energy buffer.It slows down the response of basal temperatures to surface cooling.The water layer thicknesses simulated here are unrealistically high compared to conditions under real ice masses.More realistic simulations would require a subglacial hydrology model, but this is beyond the scope of this paper.
Experiment B addresses whether the models are able to reproduce the steady-state analytical solution for certain polythermal conditions including advection, diffusion and strain heating.The models agree well with the analytical solution for K 0 → 0, if the vertical resolution is high.All models meet the transition conditions for the melting CTS, although no explicit boundary conditions are implemented.An adequate treatment of the abrupt change of conductivity at the CTS in the numerical discretization scheme is required to achieve this behaviour.The usage of an arithmetic mean (TIM-FD 3 ) or a Heaviside as well as a smoothed Heaviside function (COMice) for the conductivity jump leads to oscillations in the enthalpy solution that are visible e.g. in a timevarying CTS position.Consequently, no steady-state solution is reached under these conditions.The harmonic mean approach of Patankar (1980, chap. 4.2.3) for the conductivity (TIM-FD 3 and ISSM) leads to a continuous heat flux at the CTS and violates the condition of Eq. ( 11) (non-continuous).Nevertheless, the harmonic mean strongly favours the lower conductivity K 0 for small ratios CR = K 0 /K c and this leads to the apparent jump in ∂E/∂z.
TIM-FD 3 tends to underestimate the water content at the base of a temperate ice layer.This would result in stiffer ice at the base.In typical applications of the model the vertical layers are not equidistant as in this study, but refined towards the base.We therefore expect only a minor influence on the velocity field.ISSM simulations underestimate the temperature in the cold part accompanied by an overestimation of the water content in the basal temperate layer at coarse resolution.Implications for the overall stiffness are hard to obtain.Ice would deform more in the temperate part at the base, but less in the cold part above.
The understanding of moisture transport in the temperate ice is poor.If the latent heat flux can be represented as in Aschwanden et al. (2012), then it is crucial to consider the assumption made on the chosen value of K 0 .Simulations with a relatively high value of K 0 would lead to a much thicker temperate ice layer in contrast to simulations where K 0 ≈ 0. Stable numerical solutions could be obtained for temperate ice diffusivities in the chosen range of K 0 ≈ 10 −4 to 10 −8 kg m −1 s −1 and 0 for ISSM.The lower bound is therefore several magnitudes lower than K 0 = 10 −4 kg m −1 s −1 as the lowest value possible for a stable solution in Aschwanden and Blatter (2009).If one assumes a vanishing latent heat flux in the temperate part of a glacier, we would recommend using a value of K 0 ≈ 10 −6 kg m −1 s −1 (CR = 10 −3 ).For this value the CTS positions of all models are close to the analytical solution and show the smallest spread with varying vertical resolutions (Fig. 3).
The evolution equation for the enthalpy field is similar to the temperature evolution equation already implemented in thermomechanically coupled ice sheet models.Therefore an enthalpy scheme allows one to convert so-called "cold ice" models into polythermal ice models with only minor modifications, but with the restriction of melting conditions at the CTS.The question of whether exclusive melting conditions at the CTS are valid in an ice sheet is not conclusive.At least simulations of the Greenland Ice Sheet based on a two-layer front-tracking scheme performed with the polythermal ice model SICOPOLIS indicate that freezing conditions are relatively rare (Greve, 1997a, b).In the most recent version of SICOPOLIS the enthalpy scheme of Aschwanden et al. (2012) and a modified version of this scheme have been implemented as conventional one-layer enthalpy scheme and one-layer melting CTS scheme, respectively (Blatter and Greve, 2014).Thus comparisons to the two-layer front-tracking scheme can be performed for continental-scale ice sheets in the future.
The dynamics of glaciers, ice caps and ice sheets are strongly linked to the description of the rheology of temperate ice and its uncertainties.Besides the limited knowl-edge on the rheology of temperate ice, the current experimentally based relationship for the flow rate factor is only valid for water contents up to 1 % (Duval, 1977;Lliboutry and Duval, 1985).However, actual water contents found in temperate and polythermal glaciers are sometimes substantially larger (up to 5 %, Bradford and Harper, 2005).The advantage of deriving the water content by solving numerically for the enthalpy is limited by the use of a flow rate factor with a restricted validity range.Consequently, deformation experiments with temperate ice are urgently needed.

Conclusions
The proposed numerical experiments provide tests for the enthalpy implementation in numerical ice sheet models.All models applied here (TIM-FD 3 , ISSM, COMice) are able to perform these experiments successfully and agree to the analytical solutions.The enthalpy scheme determines the coldtemperate transition surface (CTS) and the vertical enthalpy profile in a polythermal glacier correctly without the need of tracking the CTS explicitly and applying additional conditions at this internal boundary.This is in particular the case for high vertical resolution for all three models.TIM-FD 3 and COMice also perform well for low vertical resolution, while the ISSM solution shows a significant enthalpy difference to the analytical solution although the analytical CTS position is met.There is a clear need for an empirical determination of the temperate ice conductivity K 0 and an improved description of the temperate ice rheology.
βp is the temperature relative to the melting point, H w is the basal water layer thickness and n b is the outward pointing normal vector.The type of basal boundary condition (Neumann or Dirichlet) is time dependent.The decision chart for local conditions given inAschwanden et al. (2012, Fig. 5) need to be evaluated at every time step.The chart encompasses four different situations:

Figure 1 . 29 Figure 1 .
Figure 1.Results for Experiment A simulated with TIM-FD 3 (blue), ISSM (red) and COMice (black) overlay each other.Phases I to III are described in the main text.The warming phase II is shaded in grey.

Figure 2 . 30 Figure 2 .
Figure 2. Simulation results compared to the analytical solution (thick solid grey line) for phase IIIa in Experiment A. TIM-FD 3 as blue solid line, ISSM as red dashed line, and COMice as black filled circles.
84 × 10 −3 m a −1 w.e. at the end of this phase as shown by the model solutions (| a (IIIa) b

Figure 3 .Figure 3 .
Figure 3.Comparison of simulated steady state CTS positions for different values of the temperate ice conductivity in Experiment B. The different models are shown as: TIM-FD 3 (blue), ISSM (red) and COMice (black).Results of different models are slightly shifted on the x axis to not overlay each other.The dashed black line indicates the CTS position of the analytical solution derived for K0 = 0. 31 Figure 3.Comparison of simulated steady-state CTS positions for different values of the temperate ice conductivity in Experiment B. The different models are shown as: TIM-FD 3 (blue), ISSM (red) and COMice (black).Results of different models are slightly shifted on the x-axis to not overlay each other.The dashed black line indicates the CTS position of the analytical solution derived for K 0 = 0.

Figure 4 .Figure 4 .
Figure 4. Simulated steady state profiles of the enthalpy, E, the temperature, T , and the water content, ω for TIM-FD 3 (blue), ISSM (red) and COMice (black) compared to the analytical solution (gray).ζ = z/H is the normalised vertical coordinate.The vertical resolution is ∆z = 10 m (upper row) and ∆z = 0.5 m (lower row), CR = K 0 /K c = 10 −5 .In the lower row the model results overlay the analytical solution.

Figure 5 .Figure 5 .
Figure 5.The root-mean-square deviation (RMSD) of the model results to the analytical solution (Experiment B) for different vertical grid resolutions ∆z.Model results for TIM-FD 3 (blue crosses), ISSM (red triangles) and COMice (black circles) are obtained for the lowest conductivity ratio CR = K0/Kc = 10 −5 applied to all models.