Comparison of hybrid schemes for the combination of shallow approximations in numerical simulations of the Antarctic Ice Sheet

. The shallow ice approximation (SIA) is commonly used in ice-sheet models to simplify the force balance equations within the ice. However, the SIA cannot adequately re-produce the dynamics of the fast ﬂowing ice streams usually found at the margins of ice sheets. To overcome this limitation, recent studies have introduced heuristic hybrid combinations of the SIA and the shelfy stream approximation. Here, we implement four different hybrid schemes into a model of the Antarctic Ice Sheet in order to compare their performance under present-day conditions. For each scheme, the model is calibrated using an iterative technique to infer the spatial variability in basal sliding parameters. Model results are validated against topographic and velocity data. Our analysis shows that the iterative technique compensates for the differences between the schemes, producing similar ice-sheet conﬁgurations through quantitatively different re-sults of the sliding coefﬁcient calibration. Despite this we observe a robust agreement in the reconstructed patterns of basal sliding parameters. We exchange the calibrated sliding parameter distributions between the schemes to demonstrate that the results of the model calibration cannot be straightfor-wardly

Abstract. The shallow ice approximation (SIA) is commonly used in ice-sheet models to simplify the force balance equations within the ice. However, the SIA cannot adequately reproduce the dynamics of the fast flowing ice streams usually found at the margins of ice sheets. To overcome this limitation, recent studies have introduced heuristic hybrid combinations of the SIA and the shelfy stream approximation. Here, we implement four different hybrid schemes into a model of the Antarctic Ice Sheet in order to compare their performance under present-day conditions. For each scheme, the model is calibrated using an iterative technique to infer the spatial variability in basal sliding parameters. Model results are validated against topographic and velocity data. Our analysis shows that the iterative technique compensates for the differences between the schemes, producing similar ice-sheet configurations through quantitatively different results of the sliding coefficient calibration. Despite this we observe a robust agreement in the reconstructed patterns of basal sliding parameters. We exchange the calibrated sliding parameter distributions between the schemes to demonstrate that the results of the model calibration cannot be straightforwardly transferred to models based on different approximations of ice dynamics. However, easily adaptable calibration techniques for the potential distribution of basal sliding coefficients can be implemented into ice models to overcome such incompatibility, as shown in this study.

Introduction
Accurate projections of ice-sheet-driven sea level changes require the use of numerical models that are capable of capturing the dynamics of rapidly flowing regions and groundingline zones (Pattyn et al., 2013). This requirement can be best accommodated using the most complete models currently available for modelling the ice dynamics, referred to as full Stokes (FS) models (e.g. Gagliardini et al., 2013). However, the timescales over which an ice sheet builds up and disintegrates in response to variations in the climatic forcing typically involve many thousands of years. Numerical experiments over such time spans are necessary to separate the long-term transient component from relatively fast fluctuations in the ice volume during the observational record. These long-term, continental-scale palaeo-simulations are currently infeasible using FS models due to the computational expenses triggered by the non-linearity of the model equations and the complex interdependence of the involved quantities.
To overcome the contemporary spatio-temporal limitations of FS models, a hierarchy of approximations has been developed over the last decades (e.g. Hindmarsh, 2004). The shallow ice approximation (SIA; Hutter, 1983) is a zerothorder approximation of the momentum balance equations that keeps only the gravity-driven vertical shear stress, predicting reasonably well the behaviour of grounded ice masses which are characterised by a thickness much smaller than their horizontal length scales. Ice floating in the sea wa-ter experiences almost no friction at the base, and its behaviour is often described by the shallow shelf approximation (SSA; Morland, 1987), which omits the vertical shear stress from the FS equations and neglects the basal drag. The transition between grounded and floating ice regions exhibits areas where ice flow is often enhanced by basal conditions favourable for sliding, generating rapid ice flow features known as ice streams. In these ice-sheet sectors membrane stresses become increasingly important, sharing many similarities with the floating ice shelves, and thus the SIA is no longer appropriate to describe the ice dynamics. It is also important to note that the absence of a membrane stress transfer in the SIA renders this approximation invalid for modelling the grounding-line migration, i.e. the migration of an interface between grounded and floating ice sectors (Pattyn et al., 2012).
More sophisticated methods have been designed to overcome the limitations of SIA models when reproducing the dynamics of ice streams, which still aim at low computational costs. The approaches used by this new generation of continental-scale ice-sheet models include depth-integrated Blatter-Pattyn models (Blatter, 1995;Pattyn, 2003) based on the asymptotic analysis by Schoof and Hindmarsh (2010), algorithms that detect and use the SIA where it is applicable and the FS elsewhere (Ahlkrona et al., 2016), and so-called hybrid models utilising heuristic combinations of the SIA and the shelfy stream approximation (SStA, which is the SSA including basal drag; Bueler and Brown, 2009;Winkelmann et al., 2011;Pollard and DeConto, 2012a). These hybrid models utilise a number of algorithms to identify zones of potential fast flow where ice streams may operate and then combine contributions of each approximation based on predefined criteria. The use of hybrid models enables simulations over hundreds of thousands of years on continental scales, yet showing a reasonable performance in idealised scenarios and intercomparison tests compared to higher-order models (e.g. Pattyn et al., 2013;Feldmann et al., 2014). Since the combination of the SIA and the SStA is based on heuristics, the approaches used to combine the two approximations vary from model to model, ranging from weighted averages of both velocity solutions to a simple summation over the entire domain.
Despite the above differences among existing models, all of them are subject to common limitations when applied to the present-day Antarctic Ice Sheet (AIS). These limitations include the scarcity of observational data needed to reduce the errors introduced by poorly constrained model parameters and boundary conditions, e.g. the flow enhancement factors introduced to account for anisotropy of ice flow, geothermal heat flux, glacial isostatic adjustment, and distribution of water-saturated sediments at the ice-sheet base. The latter has the potential to enhance basal sliding and is currently considered to be a major source of large, widespread misfits between the observed and modelled elevations of the AIS (e.g. de Boer et al., 2015). Recent studies have attempted to quantify potential distributions of these intrinsic bed properties using sophisticated inverse methods (e.g. Joughin et al., 2009;Morlighem et al., 2010;Arthern and Gudmundsson, 2010;Pralong and Gudmundsson, 2011;Arthern et al., 2015). These diagnostic methods focus mainly on the fit between the modelled and observed ice velocities. Pollard and DeConto (2012b) presented a much simpler algorithm, aiming to fit the observed surface elevations instead of velocities. The prognostic model is run forward in time, and the local elevation error is used to periodically adjust the basal sliding parameters until the best fit between the observed and modelled elevations is attained. This procedure has the ability to drastically reduce large elevation errors during the calibration and initialisation of ice-sheet models, which is an important requirement for simulations that would otherwise be undermined by poor parameter choices.
In this paper, we evaluate the performance of four different hybrid schemes implemented as part of the same continentalscale ice-sheet model, applied to the entire AIS. To this end, a calibration procedure based on the aforementioned iterative technique of Pollard and DeConto (2012b) is applied to each hybrid scheme. For comparison purposes, the same procedure is carried out using a SIA-only model. The results of these experiments are validated against an independent, observational data set of surface ice velocities (Rignot et al., 2011). Additionally, we provide insights into the relative contributions of the shallow approximations in different hybrid schemes. By exchanging the inferred distributions of basal sliding parameters between the applied hybrid schemes, we test the applicability of the model calibration results in different types of ice flow models. For hybrid approaches involving adjustable parameters, we also explore the sensitivity of the results to parameter variations. First, the ice-sheet model and the hybrid schemes are described in Sect. 2, where we also detail the iterative technique for the calibration of the basal sliding parameters. The observational and model-based data sets used in our simulations are described in Sect. 3. The set-up of the numerical experiments can be found in Sect. 4. The results are presented and discussed in Sect. 5, followed by the summary and conclusions provided in Sect. 6.

Model overview
In this study, the simulations of the AIS are carried out using the open source, three-dimensional, thermomechanical ice sheet-shelf model SICOPOLIS (SImulation COde for POLythermal Ice Sheets) version 3.2-dev, revision 619 (Greve, 1997;Greve and Blatter, 2009;Sato and Greve, 2012). It uses finite differences to solve the numerically approximated SIA and SSA equations for grounded and floating ice, respectively. Relevant modifications introduced in the model specifically for this study are presented in a greater detail in Sect. 2.2 and 2.3. SICOPOLIS is applicable for the modelling of polythermal ice sheets; i.e. it explicitly identifies potential temperate regions in which the modelled ice temperature is at the pressure-melting point (Greve, 1997). Within these regions, ice and small amounts of liquid water can coexist, and the water content is used as an additional input for a regularised Glen's flow law (Glen, 1955) utilised in our experiments, following Greve and Blatter (2009), in the form where η denotes ice shear viscosity, T m is the temperature difference relative to the pressure-melting point, σ e is the effective shear stress, σ 0 = 10 kPa is a small constant used to prevent singularities when σ e is very small, n = 3 is the flow law exponent, and A is a temperature-and pressuredependent rate factor (Cuffey and Patterson, 2010). In temperate ice regions, A is modified to account for the aforementioned water content, following Lliboutry and Duval (1985). The empirical coefficient E is a flow enhancement factor, which is introduced to account for the effects of anisotropic ice fabric. The value of E depends on the deformation regime of ice, which is different for grounded ice, where horizontal shear prevails, and floating ice, dominated by longitudinal stretching (Ma et al., 2010). In general, large-scale marine ice-sheet models use a homogeneous, constant value ranging between 1 and 8 for the grounded ice and between 0.2 and 1 for ice shelves (e.g. de Boer et al., 2015). Within these ranges, the computed age of ice is often used to assign different values for glacial and interglacial ice. Here, we use E = 1 and E = 0.5 for grounded and floating ice, respectively. These values are smaller than those chosen in previous studies using SICOPOLIS (e.g. Sato and Greve, 2012) and are based on our initial tests and the sensitivity analysis by Pollard and DeConto (2012b). At the base of the grounded ice sectors, stress conditions at the bedrock and the associated potential for sliding are linked to the basal velocity, u b , used as a boundary condition for the computation of the SIA velocities, through an empirical Weertman-type sliding law (Weertman, 1964;Dunse et al., 2011), in the form where τ b is the basal shear stress, and p = 3 and q = 2 are the sliding law exponents. N b is the effective basal pressure, computed as where ρ ice and ρ sw are the density of ice and sea water, respectively, g is the gravitational acceleration, H is the modelled ice thickness, and H sw is the difference between the mean sea level and the ice base topography. The parameter C b depends on the basal temperature and pressure conditions: where the exponential function controls the amount of sliding, depending on the temperature below the pressuremelting point T m and the sub-melt-sliding parameter γ , ensuring a smooth transition across different basal thermal regimes and, thus, prohibiting discontinuities in the velocity field (Bueler and Brown, 2009). A spatially varying factor C 0 is introduced to account for differences in the bedrock material properties affecting sliding (e.g. hard bedrock vs. soft sediments). Potential distributions of C 0 have been explored using different iterative and inverse approaches and a variety of sliding laws, aiming to find spatially varying values that minimise the discrepancy between the modelled and observed quantities such as ice thickness, ice surface velocity, and elevation change (e.g. Joughin et al., 2009;Morlighem et al., 2010;Arthern and Gudmundsson, 2010;Pralong and Gudmundsson, 2011;Pollard and DeConto, 2012b;Arthern et al., 2015). A particular iterative technique implemented in our model is described in Sect. 2.3. Other model components include evolution equations for ice temperature and ice thickness, with the latter forced by independent modules for the computation of the surface and basal mass balances (Greve and Blatter, 2009;Sato and Greve, 2012). A summary of the model parameters used in this study is provided in Table 1.

Hybrid schemes
For this study, four hybrid approaches have been implemented into the model. Each of them offers a different way to identify the fast flowing zones and combine the horizontal SIA velocity, u, and the horizontal SStA velocity, v. For each scheme, individual velocity solutions from the shallow approximations are calculated independently. It is important to note that SStA velocities in grounded ice regions include basal drag. This implies a major difference from SSA velocities computed in the floating ice shelf sectors, for which the friction at the ice-ocean interface is negligible. For consistency, the basal drag term that enters the SStA equations is computed using the same sliding law as described in Sect. 2.1 (Eq. 2), in the form where v b is the basal SStA velocity and the drag coefficient, β drag , is computed as Here, v b x and v b y are the horizontal components of the SStA velocity at the ice base, and v 0 = 0.01 m yr −1 is a small regularisation quantity introduced to prevent singularities at the www.the-cryosphere.net/11/247/2017/ The Cryosphere, 11, 247-265, 2017 locations where there is no basal sliding (Bueler and Brown, 2009). The first hybrid scheme (henceforth HS-1) is the original implementation in SICOPOLIS v3.2-dev (revision 619) based on the slip ratio of grounded ice, computed as where u b is the Weertman sliding velocity (Eq. 2) and u s is the surface SIA velocity. At each iteration, and for each velocity component, the local slip ratio r is compared to a prescribed threshold r thr ranging from 0 to 1 (Sect. 4). If r is larger than the threshold, the grid point is flagged as streaming ice where the SStA velocities should be computed. Once SStA velocities are computed, the individual contributions of the SIA and SStA velocities at each streaming grid point are determined using the weight Then, for each streaming grid point, the hybrid horizontal velocity U is computed as recalling that u and v are the horizontal SIA and SStA velocities, respectively. The second approach (henceforth HS-2) is based on the idea by Bueler and Brown (2009), in which SStA velocities are calculated over the entire ice sheet and used as a sliding velocity complementing a non-sliding SIA model. SIA and SStA velocities are combined as in Eq. (9), and the weighting function is adopted from Bueler and Brown (2009, Eq. 22): where v ref is a reference ice velocity (Sect. 4). The velocity v ref marks the point for which the SIA and SStA contributions are equally weighted; i.e. the resulting hybrid velocity is a standard mean of both solutions. The weighting function is smooth and monotone, and its value converges towards 0 for small velocities and towards 1 when v is large compared to the reference velocity v ref . As in the HS-1, w is used to compute respective contributions from the SIA and SStA, with the difference that Eq. (10) uses v ref as the only criterion to determine the SStA contribution. The SStA velocities are calculated over the entire ice sheet, and an a priori identification of fast flow zones is not required. As described in Sect. 2.1, the SIA solution in SICOPO-LIS is computed using the Weertman sliding law (Eq. 2) as a boundary condition. To assess the influence of a SIA solution including the Weertman sliding, we have split this hybrid scheme into two: a sub-scheme (HS-2a) that replicates the idea of Bueler and Brown (2009) with no basal velocity prescribed in the computation of the SIA and a sub-scheme (HS-2b) that has a Weertman sliding component and uses it to compute a slightly modified weight.
where u b is the basal sliding velocity as in Eq.
(2). Thus, in the HS-2b the SStA solution does not serve as a replacement of a sliding law. It is rather used to determine to what extent the computed SStA velocity should replace the basal velocity used to compute the SIA solution.
The third approach (henceforth HS-3) simply adds up the non-sliding SIA and SStA solutions: This superposition of approximations has been employed in recent studies using SIA models in combination with a SStA solution as a sliding law (e.g. Winkelmann et al., 2011). It bypasses the need for additional free parameters, such as r thr and v ref in the HS-1 and HS-2, respectively. This approach is based on the assumption that on ice shelves the SIA contribution is negligible due to low surface gradients, and therefore the modelled ice flow is dominated by the SStA solution, whereas in the continental interior the modelled ice flow is dominated by the SIA solution (Winkelmann et al., 2011).
Since the SIA and SStA solutions are computed over the entire domain, their superposition enables a smooth transition across different flow regimes, ranging from slow ice motion in the interior to a characteristic fast flow of ice shelves, thereby allowing for stress transmission across the grounding line. As in the HS-2, an identification of fast flowing zones is purely diagnostic and not required during the computation of U . Table 2 presents a summary of the hybrid schemes applied in this study.

Calibration of basal sliding coefficients
We have implemented an iterative method following Pollard and DeConto (2012b) in order to determine the distribution of sliding coefficients C 0 that minimises the difference between the modelled and observed ice thickness. The method starts from a spatially uniform guess value for the distribution of C 0 and runs the model forward in time, as described in Sect. 4. At a given time step, t inv , the method uses the basal temperature below the pressure-melting point, T m , to identify grounded grid points where basal sliding may occur. If the absolute value of T m is smaller than the parameter γ from Eq. (4), i.e. close to the pressure-melting point, the method computes the difference between the modelled and observed ice thickness, which is then used to locally adjust C 0 at this grid point according to where C * 0 is an updated sliding coefficient and H = (H − H obs )/H inv . Here, H is the difference between the modelled and observed ice thickness, scaled by a factor, H inv , in order to prevent overshoots. For the same reason and following the implementation by Pollard and DeConto (2012b), variations in the value of the multiplicative factor 10 H are further limited by a range of ∼ 0.03 to 30. In contrast to previous studies using SICOPOLIS where γ = 1 K, here we set this parameter to a value of 3 K, allowing for a more frequent calibration of the sliding coefficients (Pollard and DeConto, 2012b).
Studies using this iterative technique and inversion methods have shown that potential distributions of sliding coefficients C 0 are highly heterogeneous, with values spanning several orders of magnitude (e.g. Pollard and De-Conto, 2012b;Arthern et al., 2015). To ensure numerical stability, we limit our inferred values to a range of 1 to 10 5 m yr −1 Pa −1 during the calibration procedure. Additionally, we have implemented the following condition: when the computed surface ice velocity reaches an ancillary speed limit at a certain grid point, the adjustment of C 0 for that point is halted. This prevents the method from overadjusting the sliding coefficients when the velocity limit has been reached and no noticeable changes occur in response to further adjustments of C 0 . This additional constraint is applied in order to ensure numerical stability and keep the modelled ice velocities within the range of observations. For the experiments presented in this study, the lower speed limit is defined as 0.1 m yr −1 , whereas the upper limit is set to 4000 m yr −1 . These values are based on the observed surface velocities of Rignot et al. (2011).
www.the-cryosphere.net/11/247/2017/ The Cryosphere, 11, 247-265, 2017 The iterative technique involves an additional limiting condition that prevents overadjustments of C 0 . For each individual grid point, if the difference between the modelled and observed ice thickness reduced at the previous time step, the adjustment at the current time step is deactivated. This allows previous adjustments to fully develop their effects over the following time steps and prevents the technique from adding unnecessary extra adjustments that often result in overshoots. The calibration is activated again as soon as the time derivative of the modelled ice thickness drops to zero (i.e. the difference between the modelled and observed ice thickness is not reduced anymore) or the misfit starts increasing (e.g. due to increased influx from surrounding areas). Our experiments have shown that this additional feature enables the use of a smaller t inv (50 years used here compared to 500-10 000 years in Pollard and DeConto, 2012b) because further adjustments will only be applied when and where strictly necessary. A further benefit is that it indirectly allows nonlocal adjustments of C 0 to influence the local ice dynamics: if an adjustment applied in the vicinity of a grid point reduces the misfit, further adjustments at this grid point will still be halted.

Data sets
The calibration procedure takes advantage of the improved quality of the modern, continental-scale Antarctic data sets, such as climatic forcing (Van Wessem et al., 2014) and topography (Fretwell et al., 2013). The forcing data serve as time-invariant boundary conditions for our equilibrium (steady-state) model simulations. It should be noted that the modern AIS is unlikely in a steady state and a transient simulation, e.g. of the entire last glacial cycle, would provide a more realistic scenario for the calibration procedure. However, existing reconstructions of the Antarctic palaeoclimate and past ice-sheet configurations still contain large uncertainties, with in situ data being sparse in space and time. Keeping this in mind, we think that using modern data sets and assuming equilibrium conditions is a valuable first-order approximation to a more complex model calibration. Furthermore, such equilibrium set-up can serve as an initial guess for transient deglaciation simulations, which include time-dependent processes not considered here (Fyke et al., 2014).  (Fretwell et al., 2013), including the location of the sites mentioned in the text. Initial modern conditions for surface topography, ice shelf thickness, and bedrock elevations relative to the present-day sea level are derived from the BEDMAP2 data set (Fretwell et al., 2013). BEDMAP2 is a compilation of 24.8 million ice thickness data points obtained from a variety of sources including airborne and over-snow radar surveys, satellite altimetry, seismic sounding data, and satellite gravimetry (Fretwell et al., 2013). This compilation is complemented by surface elevation data from several digital elevation models to derive previously unknown bedrock features and allows for a detailed modelling of the AIS. Figure 1 shows the bedrock topography data from BEDMAP2, together with the locations mentioned in the text. The main sources of uncertainty in the ice thickness and bedrock elevation maps are the errors in surface digital elevation models and ice thickness measurements, as well as the applied regridding, which produce overall uncertainties ranging from 59 m across areas with smooth landscapes to 1000 m in regions where only gravimetric data are available (e.g. south of Coats Land).
At the base of the thermal bedrock, geothermal heat flux is prescribed according to the map of Shapiro and Ritzwoller (2004). This map is derived from a global seismic model of the upper mantle and the crust assuming a relation between seismic velocities and mantle temperatures and uses observations from regions with similar structures to infer heatflow probability distributions where such observations are scarce or non-existent. In the inferred heat-flux map, West Antarctica is characterised by an average heat flow that is nearly three times higher than in East Antarctica. Although the resulting map depends on the accuracy of available observations as well as on the choice of the seismic model and similarity functional, the inferred distributions are robust to internal parameter changes, especially for continental areas such as Antarctica (Shapiro and Ritzwoller, 2004). However, the method cannot reproduce small-scale patterns caused by local variations in crustal heat production, which may exert an important control on the dynamics of the Antarctic ice streams.
Boundary conditions at the surface include Antarctic precipitation rates and near-surface air temperature model output from the regional climate model RACMO2.3 (Van Wessem et al., 2014), averaged over the period of 1979 to 2010. RACMO is forced at its boundaries by the reanalysis data from ERA-Interim (Dee et al., 2011). In the interior of the domain the Antarctic climate conditions are modelled with a horizontal resolution of 27 km and 40 levels in the vertical direction. RACMO contains modules that are specifically developed for glaciated regions, including a multilayer snow model. The model output from RACMO compares well with 3234 in situ observations of surface mass balance used for its validation, displaying a correlation of r 2 = 0.77 and a particularly good fit for the dry East Antarctic plateau (Van Wessem et al., 2014).
A simple lapse-rate correction of 0.008 • C m −1 is used to account for the discrepancies between the modelled and observed surface elevations. Surface melt is computed with a positive degree-day (PDD) scheme (Reeh, 1991;Calov and Greve, 2005) using the factors α ice = 8 mm d −1 • C −1 and α snow = 3 mm d −1 • C −1 (ice equivalent) for ice and snow, respectively (Ritz et al., 2001), and a standard deviation of α std = 5 • C for the statistical fluctuations in air temperature.
All input fields are projected onto a regular, rectangular, polar stereographic grid covering the entire Antarctic continent and the surrounding Southern Ocean, with a nominal horizontal resolution of 20 km, corresponding to 301 × 301 grid points. Our choice of resolution has been motivated by the large number of experiments presented in Sect. 5, each spanning 400 000 years (400 kyr), and the fact that Pollard and DeConto (2012b) have shown that the results remain essentially unchanged when the horizontal resolution is increased to 10 km (in a nested simulation), even in rapidly flowing sectors. However, our initial tests with a model resolution of 40, 20, and 10 km have identified areas where the modelled ice flow is more sensitive to the resolution used (see the Supplement). These mainly occur close to the icesheet margins where the grounding zone of a glacier is often represented by only one grid cell at lower resolution. In such regions the use of a finer grid allows for a more detailed treatment of the topographically constrained glacial flow. In the context of this comparison study, however, these limitations equally affect the performance of all hybrid schemes and do not impact our conclusions. In the vertical direction, ice columns consist of 91 layers (11 equidistant grid points for temperate ice and 81 grid points for "cold" ice densifying towards the base, sharing the grid point at their interface), mapped to a [0, 1] interval using a σ transformation (Greve and Blatter, 2009).

Experimental set-up
The modified version of SICOPOLIS described in Sect. 2 is applied to the entire AIS and the fringing ice shelves. The experiments performed during the calibration of the ice sheet-shelf system aim at quantifying the differences and similarities between the hybrid approaches. Default values for the parameters of the hybrid schemes are r thr = 0.5 (HS-1) and v ref = 100 m yr −1 (HS-2a and HS-2b). As mentioned in Sect. 2.2, the HS-3 does not include any free parameters. Additionally, and for comparison purposes, the same experiments are performed using a SIA-only scheme for grounded ice (henceforth SoS).
The bulk of the calibration consists of an iterative technique used to derive the distribution of sliding coefficients C 0 (Sect. 2.3) that exerts a dominant control on the resulting ice distribution and its fit to observations. The calibration procedure starts from a homogeneous guess value of C 0 = 1 m yr −1 Pa −1 , which is the lower limit of the considered range (see Sect. 2.3). The time steps between adjustments and the scaling factor in Eq. (13) are set to t inv = 50 years and H inv = 5000 m, respectively. We follow the method by Pollard and DeConto (2012b) and allow for a free evolution of both the ice sheet and ice shelf thickness, but their interface (the grounding line) is kept at its present-day observed position. Free evolution is needed because the calibration of the sliding coefficients requires an evolving ice thickness that will be routinely compared to observations. The reasons for a constrained grounding line are twofold: this approach (1) enables a comparison with observations in coastal regions during the calibration and (2) prevents artificial transitions between grounded and floating areas caused by equally artificial effects of unrealistic initial thermal regimes and C 0 distributions that evolve from initial guess values. Our tests show that such artefacts can produce feedbacks that are difficult or impossible to reverse. For the same reason, glacial isostatic adjustment is not included in the simulations, and ice shelf fronts are constrained to the observed locations. As mentioned above, the ice shelf thickness is allowed to evolve, but basal melt rates are adjusted at each time step in order to www.the-cryosphere.net/11/247/2017/ The Cryosphere, 11, 247-265, 2017 keep the modelled ice shelves as close as possible to observations. This ensures a consistent computation of mass fluxes across all flow regimes and does not overlap with the calibration of sliding coefficients, since these are not applied in floating ice sectors. In this study, the calibration procedure is divided into four steps. In the first three steps, a relaxation scheme is applied to the evolution of the modelled ice thickness, H . Here, the difference between the current solution of the ice thickness equations, H new , and the solution from the previous time step, H old , is scaled at every time step by a factor h rlx ranging between 0 and 1, as follows: For a time-invariant forcing, our tests have shown that different values of h rlx will result in very similar equilibrium states. The relaxation simply delays the time point at which this state is reached. However, such relaxation procedure allows for bigger time steps at which the topography evolution is computed, without affecting the internal temperature evolution. In this way, an equilibrium with the boundary conditions will be reached faster by the temperature field than by the topography. This effectively minimises transient effects in the closely associated ice thickness and velocity fields, especially at the beginning of the calibration runs when model parameters follow the initial guess values. More importantly, it allows us to simultaneously apply the iterative calibration technique described in Sect. 2.3, in contrast to approaches in which the topography is fixed (e.g. Sato and Greve, 2012). This ensures that the modelled ice thickness distribution has the closest possible match to the initial observed distribution throughout the initial calibration stages. The value of the relaxation factor is set to h rlx = 0.001 in the first stage of the calibration driven over a simulation time of 100 kyr, using a time step of 5 years. The following second and third stages replicate the first one, but with larger scaling factors of h rlx = 0.01 and h rlx = 0.1, respectively, for an additional simulation time of 100 kyr in each stage. In the final stage, the relaxation is deactivated and a smaller time step of 1 year is applied over an additional simulation time of 100 kyr, which provides a long enough time span to attain a dynamic equilibrium.

Results and discussion
In this chapter we present an ensemble of simulations of the AIS aiming to comprehensively evaluate different hybrid schemes combining the SIA and SStA. Our evaluation uses the degree of agreement between the modelled and observed ice-sheet geometries and surface velocities as a measure of model performance, allowing for a point-by-point quantification of the model errors.
Keeping in mind that each hybrid scheme builds upon the SIA solution, partially or entirely replacing it at variable lo-cations with the SStA solution, we have also included the results from the SIA-only scheme in our comparison. In particular, this enables a qualitative separation of relative contributions of the SIA and SStA, providing a new insight into the internal differences between the schemes and their applicability to ice-sheet areas with diverse dynamical characteristics.
By applying an automated model calibration against the observed ice thickness to each of the hybrid schemes we infer spatial distributions of poorly constrained sliding coefficients as a proxy for mechanical conditions at the ice-bedrock interface and assess their sensitivity to the choice of a particular hybrid scheme. In addition, the influence of variations in the parameters controlling relative contributions of SIA and SStA velocities in each scheme is assessed for a wide range of parameter values.

Comparison of equilibrium states
As described in Sect. 3, our experiments use the BEDMAP2 observational data set (Fretwell et al., 2013) as an initial icesheet configuration, running the model forward in time under a relaxation scheme (Eq. 14) until the temperature distribution within the ice sheet reaches an equilibrium state. After this initialisation, the relaxation is deactivated and the model runs until full thermal and dynamic equilibria are attained. The equilibrium is defined as the point in time in which the variations of the total grounded ice volume over a prolonged time (> 10 kyr) are smaller than 0.01 %.
Starting from the initialised states after a 300 kyr relaxation procedure, the time required to reach an equilibrium varies from scheme to scheme (Fig. 2). The HS-1, HS-2a, and HS-3 attain a virtually invariable state after only 20 kyr, as opposed to the SoS and HS-2b that display significant oscillations in the total volume around a mean equilibrium value even after a period of 100 kyr, suggesting unstable equilibrium states. Compared to the SoS, the computation of SStA velocities in the grounded ice sectors requires an extra computational effort for each iteration in the numerical solvers, with the computing time increasing by a factor of ∼ 4 for the applied hybrid schemes. The characteristic computing time of the HS-1 is somewhat shorter than that of the other hybrid schemes (but still longer than for the SoS) due to the prognostic identification of ice streams that obviates the need for the computation of SStA velocities over the entire ice sheet. However, we have observed that the iterative solvers in the model require a substantially smaller number of iterations when the hybrid schemes are used, making them numerically more stable compared to the SoS.
The calibration procedure applied to all schemes yields total grounded ice volumes which are in a close agreement with the reference value of ∼ 2.56 × 10 7 km 3 from the observational data (Fretwell et al., 2013), with the maximum deviation being below 1.5 %. Individual values of the modelled ice volumes and their respective deviations for most of our  Figure 2. Overview of the calibration procedure. Main: the evolution of the total grounded ice volume during the calibration procedure, in ×10 7 km 3 . At the end of each 100 kyr stage, the relaxation coefficient h rlx increases by 1 order of magnitude, from 0.001 in the first stage to 1 (free evolution), until a dynamic equilibrium is attained. Inset: mean absolute differences between the modelled and observed ice thickness at the end of the simulations, in metres. For each scheme, a mean absolute error is calculated for the entire ice sheet (left bars) and separately over the areas where basal sliding is identified (right bars). Fractions of the mean absolute error arising from under-and overestimations are shown in blue and pink, respectively. experiments are summarised in Table 3. The best fit to observations is obtained using the HS-2b, which corresponds to an overestimation of the total grounded ice volume by 0.1 %. The other schemes produce relatively larger misfits, with the HS-2a simulation yielding the greatest deviation of 1.46 % arising from an overestimation of the total ice volume. The smallest ice sheet is produced by the SoS that underestimates the total grounded ice volume by 0.72 %. The inset of Fig. 2 shows the total ice thickness errors for all schemes. The errors are computed as an average of the absolute misfit values in all grounded grid points and remain below 50 m for all hybrid schemes (Table 3). Among these, the largest average error of 46.6 m is produced by the HS-2a, while the smallest error of 40.0 m is obtained using the HS-2b. For the former, 83 % of the misfit is due to overestimations of ice thickness, while the latter shows nearly even contributions of under-and overestimations. This is in accordance with the results shown in Fig. 2, where the misfits obtained from the two schemes using the SStA as a sliding law are dominated by an excessive ice thickness. The hybrid schemes that include a Weertman basal sliding component during the computation of SIA velocities are not necessarily dominated by underestimations of the ice thickness. For example, less than a quarter of the misfit produced by the HS-1 can be attributed to areas with an ice thickness deficit. For comparison, the averaged absolute error produced by the SoS is 52.4 m, with 65 % of the misfit coming from an underestimation of the ice thickness.
The results depicted in Fig. 2 reflect only a timedependent, continental-scale information about the fit between the modelled and observed ice-sheet geometries, providing no insights into local model performance. In order to have a more detailed overview of the results, Fig. 3 (left column) shows the corresponding spatial maps of ice thickness errors. All schemes provide a reasonably good fit to the observational data in the continental interior, with larger discrepancies mainly occurring at the ice-sheet margin. It can be readily observed that the modelled ice sheet is too thick over mountainous regions and across the region between the Shackleton Range and the Pensacola Mountains (south of Coats Land), which has some of the largest uncertainties in the topographic data (Fretwell et al., 2013). These are common features for all schemes, regardless of the approach chosen for the sliding component, and could originate from insufficient geothermal flux, too high precipitation rates, and/or an unrealistic smoothness of the bedrock in this area that prevents the formation of topographically driven ice streams. The distribution of zones where the modelled basal temperatures are far below the pressure-melting point are depicted as white-coloured areas in Fig. 3 (right column). In general, these areas coincide with the locations where the largest ice thickness errors occur. Areas where the ice thickness is underestimated are mainly located at and around the ice margins. In many of these areas, however, sliding is identified and the calibration of C 0 reaches its lower limit, implying that the surface mass input is insufficient to minimise the misfit. The SoS produces too thin ice along most of the Table 3. Summary of the results at the end of the calibration procedure for each of the applied hybrid schemes (HS, see Sect. 2.2) and SIA-only scheme (SoS), including the total grounded ice volume, V grd (km 3 ); a deviation from the total grounded ice volume, V grd (%); a mean absolute error in the ice thickness, H (m); a fraction of the total area where basal sliding occurs, A sld (%); a mean ice thickness error only where basal sliding operates, sld H (m); a mean surface velocity error, v s (m yr −1 ); a fraction of the grounded ice area dominated by the SIA (only where w < 0.25), A SIA 0.25 (%); and a fraction of the grounded ice area dominated by the SStA (only where w > 0.75), A SStA 0.75 (%). ice-sheet margin, which explains the high percentage of the ice thickness errors arising from underestimations (inset of Fig. 2). In contrast, the hybrid schemes produce error patterns that differ only slightly between each other. The only exception is the HS-2b, which is closest to the SoS, albeit exhibiting smaller underestimations at the ice-sheet margins. Furthermore, we evaluate the performance of the calibration procedure over the areas where basal sliding is identified and the calibration is applied. For this purpose we have calculated averaged absolute errors in the ice thickness across regions where basal sliding operates (inset of Fig. 2). It should be kept in mind that the calibration procedure also affects the ice masses located in the immediate proximity to sliding areas through surface elevation gradients and/or stress transmission. Mean ice thickness errors over the sliding areas are smaller than those estimated for the whole ice sheet, although the degree of relative improvement varies from scheme to scheme (Table 3). For example, the misfit produced by the SoS decreases by only ∼ 20 % if calculated over the sliding ice-sheet sectors, while more than 60 % of the errors resulting from the HS-2a occur over the areas where no sliding is identified and C 0 is not calibrated. For all hybrid schemes, the percentage of the errors associated with an underestimation of the ice thickness increases substantially when their performance is assessed across areas where basal sliding coefficients are calibrated. This supports the observation that the modelled ice is excessively thick in regions where C 0 is not calibrated.
The inferred distributions of C 0 are shown in Fig. 3 (right column). In general, the areas where the calibration is performed are similar for all schemes (Table 3), although there is a significant spread in the retrieved values. The HS-1 scheme predicts a minimum corresponding to 45 % of the total grounded ice area, contrasted by higher percentages for the other schemes ranging between 54 and 62 %. The upper limit of 10 5 m yr −1 Pa −1 for the sliding coefficient is reached by all schemes across up to 18 % of the ice-covered land, with the highest concentration occurring in the Siple Coast region, where ice streams flow rapidly over a smooth and deformable bed due to lubrication by water saturated subglacial sediments (e.g. Blankenship et al., 1986;Alley et al., 1987;Kamb, 2001). The upper limit is also reached in Coats, Mac-Robertson, and Ellsworth lands (Fig. 1). Similarly, all calibration runs estimate the C 0 value to reach its lower limit of 1 m yr −1 Pa −1 over ∼ 10 to ∼ 32 % of the ice-sheet-covered area. The SoS and HS-2b infer this value throughout most of the West AIS and over vast parts of the East AIS, particularly at the ice margins. Its incidence is relatively lower for the other hybrid schemes, especially for the HS-1 and HS-2a. Quantitatively, a good agreement between the inferred coefficient distributions is mainly restricted to the areas where C 0 reaches an upper or lower limit. In other areas the values generally vary across orders of magnitude.
Direct comparison of our results to those from Pollard and DeConto (2012b) is hindered by the differences in the sliding laws and hybrid schemes. Nevertheless, we have found a good qualitative agreement in the inferred distributions of C 0 , with similar patterns of high vs. low values. For instance, both studies identify the Siple Coast and Coats Land regions as areas where the ice streaming flow is driven by basal conditions favourable for sliding. A similar agreement is found for the Thwaites and Pine Island Glacier areas, as well as in the MacRobertson Land (Fig. 1). Low values of C 0 are predicted in the continental interior of West Antarctica and over most of the East AIS.

Analysis of the SIA and SStA contributions
As described in Sect. 2.2, different hybrid approaches are not expected to produce exactly the same equilibrated ice velocity fields. We demonstrate this in Fig. 4, where the modelled steady-state surface velocities derived from the SoS and the hybrid schemes are compared to a continental-scale observational data set (Rignot et al., 2011). This high-resolution (900 m) data set contains many small-scale features that are unresolved by the model due to its lower resolution. Nevertheless, all schemes are able to reproduce the observed range   Rignot et al. (2011) (bottom), regridded to the model resolution of 20 km. Right column: ratios of the modelled to observed surface velocities, plotted on a logarithmic scale. Velocities smaller than 2 m yr −1 are excluded. Colour-code saturates at ratios larger than 100 or smaller than 0.01. of ice flow velocities, distinguishing between the ice-sheet areas with low velocities near the ice divides and fast flowing ice streams reaching the ice-sheet margins. In the transition zones between the continental interior and the ice-sheet margins, all schemes reproduce to some extent the fast flowing ice streams identified by observations. However, in contrast to the results from the hybrid schemes, surface ice velocities derived from the SoS are characterised with noise-like patterns, which are especially visible in the areas of rapid ice flow. We attribute these artefacts to a combination of lacking stress transmission in the SIA, which allows for steep gradients in the modelled velocities, and the calibration procedure, which can potentially amplify these gradients through local adjustments of C 0 .
Although the overall character of the observed surface ice velocities is qualitatively well reproduced by all the hybrid schemes, modelled ice flow is clearly too fast at and around several ice stream locations, such as at the Siple Coast. Furthermore, modelled surface velocities are generally overestimated close to the grounding lines of most outlet glaciers, in particular due to the resolution-related limitations discussed in Sect. 3. These overestimations are particularly large in the SoS simulation, even though the discrepancies between the modelled and observed ice thickness are small (Fig. 3, left  column). Overestimations of the ice flow velocity by the hybrid schemes at the grounding zone are smaller than those derived from the SoS simulation, but they cover a wider area and reach further upstream. Similar observations have been made by Pollard and DeConto (2012b) using a different hybrid ice sheet-shelf model and a different set of topographic ice observations and external forcing. They proposed that the overestimation of surface velocities in these areas may be caused by a coarse model resolution, exaggerated snowfall rates, or an excessive internal deformation compared to sliding near the ice margins.
In order to quantitatively evaluate the model fit to observations, we have calculated point-by-point ratios between the modelled and observed surface velocities (Fig. 4, right column). As mentioned above, the modelled velocities are in a good agreement with observations in the continental interior characterised by slow ice motion. In contrast, at the margins the ice velocity simulated by the hybrid schemes sometimes reaches values that are several hundred times higher than in the observational data set. However, this mostly happens across areas where our model generates a non-existent fast flow. One of the best examples of such model artefacts is the former Ice Stream C in the Siple Coast, which is predicted by the model but has been stagnant for ∼ 150 years in reality (Hulbe and Fahnestock, 2007;Engelhardt and Kamb, 2013). In some cases, the locations of the modelled ice streams are shifted relative to the observed ones, thereby generating adjacent zones of under-and overestimations of the surface velocity. In other cases, the modelled rapid ice flow follows a different route compared to observations, sometimes merging with neighbouring ice streams. These shifts may origi-nate from local errors in the bedrock topography data accentuated by their projection onto the coarse horizontal grid we use here.
The mean errors in the absolute surface ice velocity fall within the range of ∼ 16 to ∼55 m yr −1 (Table 3), with the HS-2b and SoS producing the minimum and maximum misfits, respectively, analogously to the results for the mean ice thickness errors discussed in Sect. 5.1. In the SoS simulation a general underestimation of the ice thickness near the icesheet margins coincides with areas where sliding coefficients tend to reach the lower limit prescribed for the calibration. In turn, the use of such low values triggers a slowdown of the ice flow in the transition zone.
Although the results of the HS-2b simulation presented in Sect. 5.1 are in many aspects similar to those from the SoS, their abilities to reproduce the observed ice flow patterns are very different. The HS-2b uses the basal velocity from Eq. (2) (utilised by the SoS) to compute the weights of relative contributions of the shallow approximations, thereby adding the SStA contribution where sliding velocities from Eq. (2) are high. In these rapidly flowing sectors, we attribute a better performance of the HS-2b compared to the SoS to the inclusion of the stress transmission by the SStA. It also fosters a SStA-dominated modelled ice flow in the surroundings of ice streams, particularly upstream, thereby improving the overall fit to observations. Compared to the HS-2b, the HS-2a has a similar performance with the second smallest mean velocity error of 17.6 m yr −1 , followed by higher mean errors of 23.2 and 29.7 m yr −1 produced by the HS-3 and HS-1, respectively. Thus, the use of hybrid schemes allows for a twoto threefold reduction in the surface velocity misfit.
As described in Sect. 2.2, the hybrid schemes used in this study mainly differ in how the SStA weight, w, is computed. This controls not only relative contributions from the shallow approximations but also the locations where such combinations are implemented. In order to provide a deeper insight into how each hybrid scheme combines the SIA and SStA velocities, Fig. 5 illustrates the equilibrated distributions of w. It can be immediately observed that the distributions of w produced by the hybrid schemes are very different, both in spatial coverage and inferred values. For example, the percentage of the grounded ice area where the modelled ice flow is dominated by the SStA (w > 0.75) is ∼ 30 % for the HS-1, but only 2.4 and 8.5 % for the HS-2a and HS-2b, respectively (Table 3). Furthermore, the transition between the SIA-and SStA-dominated ice-sheet sectors appears patchy in the HS-1, whereas it is smooth and collocated with the present-day ice streams in the HS-2a, thereby resembling the observed surface velocities. Such transitions are sharp in the HS-2b, implying a simple differentiation between fast and slow ice flow areas. This particular scheme exhibits a SIA-dominated ice flow regime (w < 0.25) over 82 % of the grounded ice area that may explain a high degree of similarity with the SoS, especially keeping in mind that it also includes basal sliding in the SIA. In contrast, the HS-2a sets basal ice ve-locities to zero in the computation of the SIA velocities, still producing a similar percentage of SIA-dominated area. In these sectors, differences in the ice thickness derived from the HS-2a and HS-2b (Fig. 3, left column) can be attributed to a presence or absence of basal sliding in the computation of the SIA solution that may prevent overestimations e.g. in interior East Antarctica and cause underestimations in other sectors, such as the surroundings of Dronning Maud and MacRobertson Lands.
It is clear that the degree of agreement between the modelled and observed ice flow is sensitive to the choice of a particular hybrid scheme and the way it measures relative SIA and SStA contributions. This can be best visualised using scatter plots of the modelled vs. observed surface velocities for each scheme, colour-coded for the values of the corresponding w distributions (Fig. 5). In general underestimations of the modelled velocities occur across slowly flowing (< 10 m yr −1 ) areas, which are dominated by the SIA, while the areas of fast flow dominated by the SStA are responsible for most of the overestimation. At first glance it may seem that overestimations are caused by an excessive contribution of the SStA, but a comparison with the SoS scatter plot shows that this is not necessarily the case. The largest overestimations occur in the SoS simulation, with surface velocities reaching the upper permitted limit clustered in the upper part of the scatter plot (Fig. 5). As mentioned above, the SoS and HS-2b share many similarities, and a comparison between their respective scatter plots shows that the use of the hybrid scheme reduces both underestimations in the lower velocity range and overestimations in the fast flowing areas.

Intercomparison of the inferred basal sliding coefficients
Although the performance of the hybrid schemes is quantitatively similar when evaluated against observations, the inferred values of basal sliding coefficients (C 0 ) vary by orders of magnitude in many regions of Antarctica (Fig. 3, right column). It is important to keep in mind that these large differences between the inferred values of C 0 arise solely from the differences in the hybrid schemes, namely in their ways to combine the shallow approximations, since all other model components are exactly the same for all experiments.
To demonstrate the effects of prescribing a distribution of C 0 derived from one ice model (using a specific hybrid scheme) into a different ice model (using a different hybrid scheme or other level of approximation), we performed an additional set of experiments, which utilises the equilibrium states described in Sect. 5.1 as initial conditions. In these experiments, we exchange the distributions of C 0 inferred from the HS-2a and HS-2b, and then run the model over a period of 100 kyr. As described in Sect. 2.2, these schemes slightly differ in how they compute the SStA contribution, mainly due to different techniques used to account for basal sliding. Although both schemes identify similar locations of rapidly flowing ice (Sect. 5.2), their inferred distributions of basal sliding coefficients contain the highest variability among all hybrid schemes implemented in this study. At the end of these additional 100 kyr runs, the mean misfit between the modelled and observed ice thickness is above 200 m for the HS-2a and almost 350 m for the HS-2b (not shown). This represents increments of ∼ 330 and ∼ 775 % in the deviations from observations, respectively.
To further exemplify the significance of the associated uncertainties in the retrieved basal sliding parameters, we prescribe the median of the inferred distributions of C 0 , computed from all hybrid schemes, at the base of the modelled ice sheet. The choice of a median over an average is motivated by our initial tests in which generally larger values of C 0 inferred from the HS-2a tend to produce an average biased towards this particular distribution. Figure 6 shows differences between the modelled and observed ice thickness at the end of these additional runs. Comparison with Fig. 3 reveals a general degradation of the fit between the model and observations. The HS-2a and HS-2b exhibit the largest sensitivity to the change in the basal sliding parameters, with a significant amplification of over-and underestimations of the ice thickness occurring across most of the ice sheet, respectively. For the HS-2a, an absolute ice thickness misfit increases by ∼ 180 % to a mean value of 131 m (Table 4), of which 95 % is due to overestimations. The misfit increases by ∼ 355 % to a mean value of 182 m in the HS-2b simulation, with underestimations of the ice thickness accounting for 88 % of the total difference. The degree of degradation is less pronounced when the HS-1 and HS-3 are used. The misfit for the HS-1 increases by only ∼ 50 %, displaying a mixture of areas where the modelled ice thickness is either too large or too small relative to observations. The HS-3 experiment shows an intermediate degree of sensitivity to a change in basal parameters and exhibits similar misfit patterns as the HS-2b, albeit the magnitude of the underestimation is smaller and the mean absolute error in the ice thickness remains around 100 m. These experiments show that care is needed when using quantifications of the basal conditions obtained from external sources (e.g. Joughin et al., 2009;Morlighem et al., 2010;Arthern and Gudmundsson, 2010;Pralong and Gudmundsson, 2011;Pollard and DeConto, 2012b;Arthern et al., 2015) as input in ice flow models. This concerns not only differences between the basal sliding approaches implemented in each model, but also models using the same sliding law as part of different hybrid schemes, as shown in this study.

Exploration of the hybrid parameter space
The disparity in the results from the hybrid schemes presented in the previous sections showcases the impacts of slight changes in the model representation of fast flowing zones of ice sheets, even though it corresponds to a small fraction of the parameter space in ice-sheet models. In order Figure 6. Difference between the modelled and observed ice thickness for all hybrid schemes, in metres, after additional 100 kyr runs in which a median of the inferred distributions of C 0 is prescribed at the ice-sheet base. The values of C 0 derived from the SoS are not included in the median distribution. Table 4. Summary of the results for each hybrid scheme at the end of additional 100 kyr runs where the median of the inferred basal sliding coefficients from all implemented hybrid schemes is prescribed as an external forcing, including the total grounded ice volume, V grd (km 3 ); a deviation from total grounded ice volume, V grd (%); a mean absolute error in the ice thickness, H (m); and a mean ice thickness error only where basal sliding operates, sld H (m). to explore the sensitivity of the results to parameter variations within this parameter space, we perform an additional series of experiments where we vary the somewhat arbitrary threshold and reference quantities used by some of the hybrid schemes.
As described in Sect. 2.2, the hybrid schemes HS-1, HS-2a, and HS-2b use the weight w to calculate relative contributions of the SIA and the SStA. For the HS-1, the computation of w involves a prescribed threshold for the slip ratio of The effects of variations in these parameters on the evolution of the total grounded ice volume during the calibration procedure are demonstrated in Fig. 7. Here we test parameter values within a range that contains almost every possible scenario because values outside the range are non-physical (HS-1) or exhibit no noticeable differences compared to the range limits (HS-2a and HS-2b).
For the HS-1, the upper limit for the slip ratio threshold, r thr = 1, implies the use of a SIA-only solution, whereas the lower limit, r thr = 0, implies that a combination of SIA and SStA solutions is applied everywhere, determined by the weight w computed using the slip ratio. Within the range of tested values, maximum deviations from the observed grounded ice volume are below 1.2 %. Only the use of r thr = 1 produces an ice sheet that is smaller than observed. An analysis of the mean absolute differences between the modelled and observed ice thickness (Fig. 8) reveals that the use of larger values of r thr leads to a larger misfit with observations, with the maximum error corresponding to r thr = 1. Figure 8. Mean absolute differences between the modelled and observed ice thickness at the end of the equilibrium simulations, in ×10 7 km 3 , for different values of free parameters included in the computation of the SStA weight. As in Fig. 2, a mean error is calculated for the entire ice sheet (left bars) and separately over the areas where basal sliding is identified (right bars). Each bar is divided into relative errors arising from under-(blue) and overestimations (red).
In the HS-2a, all tested values produce an ice sheet that is larger than observed. This is caused by a general overestimation of the ice thickness in the continental interior. Here, the modelled ice flow is slow and dominated by the SIA. Thus, the weights w are small, implying a negligible contribution of the SStA, independently of the value of v ref . Values smaller than v ref = 100 m yr −1 are found to produce total grounded ice volumes and mean ice thickness misfits that are close to the results of the reference run (v ref = 100 m yr −1 ), although with a slight improvement in the model fit to observations. Values of v ref > 100 m yr −1 lead to a larger misfit, with the highest tested value of v ref = 1000 m yr −1 pro-ducing a grounded ice volume that displays a deviation of 2.3 % from observations. The use of even larger values of v ref asymptotically decreases the contribution of the SStA. It is important to note, however, that since the HS-2a does not include sliding in the SIA component, its solutions will never approach that of the SoS.
In contrast, the HS-2b activates the sliding law even in the regions where ice streams have not been identified by the scheme, as long as the conditions for sliding described in Sect. 2.1 are fulfilled. This may explain the similarity between the grounded ice volumes produced by the SoS and the HS-2b solutions, using a reference velocity set to the upper limit of v ref = 1000 m yr −1 . This parameter value produces a deviation of −0.9% from the observed ice volume, which is similar (but with an opposite sign) to the misfit obtained from the simulation using the lower limit of v ref = 5 m yr −1 tested in this study. For values higher than the default value of v ref = 100 m yr −1 , the simulations produce an ice sheet that is smaller than observed. The use of v ref = 60 m yr −1 leads to the best fit between the modelled and observed ice thickness among all schemes and parameter values, reaching ∼ 39 m, with a small deviation from the observed grounded ice volume of 0.27 %.
However, it is important to keep in mind that the calibration procedure has a limited power to reduce the misfits between the model and observations in areas where no sliding occurs. This is why some deviations are expected from our numerical experiments, in particular due to overestimations of the ice volume over mountainous regions where the calibration procedure does not operate. Therefore, a perfect fit to the observed ice volume obtained by any of the schemes would likely involve underestimations of ice thickness in other regions.

Summary and conclusions
We have implemented and compared the performance of four hybrid schemes for the combination of the shallow ice and shelfy stream approximations as part of the ice-sheet model SICOPOLIS. The use of shallow approximations enables continental-scale, long-term palaeo-simulations of ice sheets, preventing the restrictive computational expenses of more complex models. Moreover, hybrid schemes overcome the limitations of simpler SIA-only models in regions characterised by rapid ice flow driven by basal sliding. The hybrid schemes in this study differ in the ways (1) relative contributions of the shallow approximations are computed, (2) areas where a combination of shallow approximations is applied are identified, and (3) basal sliding is accounted for.
By adapting a simple iterative technique to infer the distribution of basal sliding parameters from observational topographic data sets (Pollard and DeConto, 2012b), we show that all the applied hybrid schemes produce dynamic equilibrium states that are in a good agreement with observations. For all the schemes, mean ice thickness misfits range from 40 to 52.4 m, and the total grounded ice volume deviations do not exceed 2.5 % for a wide range of the associated parameter uncertainties. For optimal parameter choices, mean errors in the ice thickness remain below 40 m. For comparison, present-day simulations in continental-scale ice-sheet models typically produce widespread errors in the ice thickness reaching hundreds of metres, with deviations in total grounded ice volume that can exceed 10 % (e.g. de Boer et al., 2015). A clear limitation of the calibration method presented here is the possibility of an error cancellation during the adjustment of C 0 . At present there is no easy way to differentiate between the effects of an artificial compensation of, for example, the lacking model physics or errors in the input data sets through the model calibration and the model parameters that are representative of the actual conditions at the base of the ice sheet. This shortcoming will hopefully be overcome in the future once the necessary observational data sets become available and the use of more sophisticated ice models is more feasible. The degree of applicability of the calibrated C 0 distributions for long-term palaeoclimate simulations is currently unclear, and its quantification requires further modelling studies of the past ice-sheet geometries. However, we believe that the retrieved variability in basal sliding coefficients is a good first-order approximation and is certainly a better guess than commonly used spatially uniform values of C 0 over the entire domain. We believe that the method used to derive C 0 could be potentially extended to the use for other time periods or climate regimes (or even adapted for the transitions between largely dissimilar regimes), providing that the necessary topographic and climate data are available. Considering the recent efforts to reconstruct the past geometries of the AIS Bentley et al., 2014), this may well become possible in the upcoming decades.
We also computed a non-hybrid, SIA-only solution to allow for a qualitative separation of relative contributions of the SIA and the SStA. Such direct comparison is possible because all the hybrid schemes analysed in this study share identical model components. Although the calibration procedure applied to the SIA-only scheme also produces a reasonable fit to observations, the misfits are overall larger than those from the hybrid schemes. Moreover, the modelled surface velocities exhibit noise-like patterns at and around locations where rapid ice flow exists. Since these patterns are not present in the hybrid solutions, we attribute their occurrence to the lack of stress transmission in the SIA flow model.
We find that individual weights assigned to the SIA and SStA solutions vary significantly from scheme to scheme. For the schemes in which the SIA and SStA weights are computed, the modelled ice flow is dominated by the SIA over ∼ 50 to ∼ 80 % of the AIS, with the predominant contribution of the SStA generally limited to the locations and surroundings of the observed ice streams.
Comparison of our results with an independent, observational data set of surface ice velocities reveals a reasonable agreement in the continental interior of Antarctica. However, the modelled ice flow appears to be too fast at the ice-sheet margins in all hybrid schemes, especially close to the observed ice stream locations. These misfits between the modelled and observed surface velocities are contrasted by relatively small differences between the modelled and observed ice thickness in these areas. Such misfits may originate from other factors such as uncertainties in the observational ice thickness and surface velocity data sets (Fretwell et al., 2013;Rignot et al., 2011, see Sect. 3), limited model resolution, the assumption of an ice-sheet equilibrium with present-day climate conditions, and/or errors in the model-based geothermal heat flux and/or climatic forcing (Pollard and DeConto, 2012b).
Although all hybrid schemes produce a comparably good fit to the observational topographic data set, the inferred values of basal sliding coefficients largely differ between different schemes. Therefore, the discrepancies in the representation of the ice dynamics by different hybrid schemes are compensated through adjustments of this poorly constrained parameter. Nevertheless, the schemes mostly agree on the areas where the ice base is at or close to the local pressuremelting point. Furthermore, there is a qualitative agreement in the patterns of low vs. high values of C 0 obtained from each calibration run. Although any attempt to quantify a local potential for sliding would require new and improved observational data sets describing basal conditions, the robustness of the inferred patterns provides an initial guess for their real distribution.
Furthermore, we have evaluated the performance of each hybrid scheme based on both its fit to observations and its numerical stability for a range of model set-ups and have found that at grid resolutions higher than 20 km (e.g. 15 and 10 km) the HS-1 and HS-2b become numerically unstable and stop before the experiments are completed, owing to large gradients in basal sliding coefficients arising from the use of basal velocities as boundary conditions for the SIA solutions in conjunction with the calibration of sliding coefficients. In addition, variations in the free parameter r thr of the HS-1 generate artefacts in the resulting basal temperature field. The HS-2a and HS-3, which utilise the SStA as a sliding law, are numerically more stable to variations in model parameters and changes in grid resolution. The drawback of the HS-2a is in its limited ability to improve the fit between the modelled and observed ice thickness in ice-sheet sectors where the SStA velocities are low (< 100 m yr −1 ), which is the case for large parts of the Antarctic interior. This limitation is independent of the choice of the parameter v ref (Fig. 7). The HS-3 overcomes this limitation by accounting for the SStA contribution everywhere. However, the HS-3 also utilises the SIA velocities over the entire ice sheet, which in certain areas, such as the steep ice-sheet margins (Fig. 4), are excessively high, leading to an increased root-mean-square error compared to HS-2a (Fig. 5). To improve the performance of both schemes, our future work will reconcile the drawbacks of HS-2a with the advantages of HS-3, providing a very stable and flexible hybrid scheme.
Finally, we assessed the effects of differences between the calibrated sliding parameters derived from the four hybrid schemes by performing additional experiments in which the inferred distributions of C 0 were either averaged or exchanged between the schemes. Our results show that the parameter distributions are not exchangeable between the schemes, since this leads to a strong degradation of the model fit to observations. This suggests that results of a model calibration and/or initialisation cannot be straightforwardly transferred to a model that uses a different level of approximation of the Stokes equations. Given an increasing number of studies attempting to quantify the basal conditions under ice sheets through a variety of methods including ice flow models, our experiments show that one needs to be careful when using the inferred basal sliding parameters as external boundary conditions in glaciological models.

Data availability
The ice-sheet model SICOPOLIS is available as free and open-source software (under the GNU General Public License) via (Greve, 2017). The data produced by SICOPOLIS for this study can be obtained by contacting the corresponding author.
The Supplement related to this article is available online at doi:10.5194/tc-11-247-2017-supplement.