A new sub-grid surface mass balance and flux model for continental-scale icesheet modelling : testing and last glacial cycle

Abstract. To investigate ice sheet evolution over the timescale of a glacial cycle, 3-D ice sheet models (ISMs) are typically run at "coarse" grid resolutions (10–50 km) that do not resolve individual mountains. This will introduce to-date unquantified errors in sub-grid (SG) transport, accumulation and ablation for regions of rough topography. In the past, synthetic hypsometric curves, a statistical summary of the topography, have been used in ISMs to describe the variability of these processes. However, there has yet to be detailed uncertainty analysis of this approach. We develop a new flow line SG model for embedding in coarse resolution models. A 1 km resolution digital elevation model was used to compute the local hypsometric curve for each coarse grid (CG) cell and to determine local parameters to represent the hypsometric bins' slopes and widths. The 1-D mass transport for the SG model is computed with the shallow ice approximation. We test this model against simulations from the 3-D Ice Sheet System Model (ISSM) run at 1 km grid resolution. Results show that none of the alternative parameterizations explored were able to adequately capture SG surface mass balance and flux processes. Via glacial cycle ensemble results for North America, we quantify the impact of SG model coupling in an ISM. We show that SG process representation and associated parametric uncertainties, related to the exchange of ice between the SG and CG cells, can have significant (up to 35 m eustatic sea level equivalent for the North American ice complex) impact on modelled ice sheet evolution.


Introduction
The resolution used in any model of complex environmental systems (e.g.ice sheet models (ISMs), general circulation models or hydrological models) limits the processes that can be represented.For continental-scale glacial cycle contexts, ISMs are currently run at resolutions of about 10-50 km (Pollard and DeConto, 2012;Tarasov et al., 2012;Colleoni et al., 2014).Processes such as surface mass balance on mountain peaks, iceberg calving, and ice dynamics in fjords are sensitive to scales of about 100 m to a few kilometres and therefore have to be parameterized.For example, even at 10 km grid resolution, mountain peaks are smoothed to bumps in a plateau (Payne and Sugden, 1990), inducing errors in computed surface mass balance (Marshall and Clarke, 1999;Franco et al., 2012).If the mean surface elevation of a coarse grid cell is below the equilibrium line altitude (ELA), ice ablation is overestimated (e.g.Tarasov and Peltier, 1997).Thus, coarser grid resolution can lead to temporal and spatial errors in ice sheet inception (Abe-Ouchi and Blatter, 1993;Abe-Ouchi et al., 2013) and subsequent evolution (Van den Berg et al., 2006;Durand et al., 2011).
Any model of complex environmental systems will have sub-grid (SG) processes that are, by definition, not dynamically resolved.Accurate modelling of such systems must therefore determine whether the SG processes variability is relevant for the given context.If it is, some of the impact of this SG variability may be captured in a parameterized form (Seth et al., 1994;Leung and Ghan, 1995;Marshall and Clarke, 1999;Giorgi et al., 2003;Ke et al., 2013).For ex-Published by Copernicus Publications on behalf of the European Geosciences Union.
ample, to improve surface mass balance in continental-scale ice sheet models, Marshall and Clarke (1999) used hypsometric curves, which represent the cumulative distribution function of the surface elevation.In this method, each individual glacier is not explicitly represented.Instead, 2-D topographic regions are parameterized with different hypsometric bins, representing a discrete number of elevations and their associated area.In addition to ablation and accumulation at each SG bin, there is SG ice transport from high elevation regions to valleys where the average altitude is below the ELA.Starting with ice-free conditions, Marshall and Clarke (1999) found an increase in the total ice volume over North America after a 3 kyr1 simulation when this hypsometric parameterization is coupled to an ice sheet model.The impact and accuracy of this SG model have yet to be quantified.The model was only validated against observations of a glacier located in the region used for tuning the parameterization (Marshall et al., 2011).Moreover, the communication between the SG and coarse grid (CG) models was identified as a potentially important source of error (Marshall and Clarke, 1999), but its impact has not been documented.
In this paper, we develop a new SG model extending the approach of Marshall and Clarke (1999) and Marshall et al. (2011).We use hypsometric curves that account for a much larger set of topographic information than just the maximum, minimum and median elevation.We present a new slope parameterization to compute the velocities that account for SG slope statistics.An effective width is added for the representation of the ice fluxes between SG bins.In contrast to the one-way communication used in the past, another modification to the original model is a two-way exchange of ice between the SG and CG cells.The CG ice thickness updating accounts for SG ice thickness, and the SG model accounts for ice flux out of the CG cell.For the first time, we evaluate the accuracy of the SG model against high resolution simulations by a higher order ice sheet model (ISSM; Larour et al., 2012).Sensitivities to the SG model configuration, such as the number of hypsometric bins, are assessed.We examine the extent to which the inclusion of further topographic statistics (e.g. the peak density in a region or the variance of the slopes) can improve computed sub-grid fluxes.We also evaluate the impact of embedding the SG model in the Glacial Systems Model (GSM, formerly the MUNGSM), our coarse grid model, for last glacial cycle simulations of the North American ice complex, using an ensemble of parameter vectors from a past calibration of the GSM (Tarasov et al., 2012).Special attention is given to the impact of the coupling between the SG model and the GSM.

Sub-grid model
The sub-grid model is a finite difference flow line model composed of a diagnostic equation for the ice velocities and a prognostic equation for ice thickness evolution.The surface mass balance is calculated using a positive degree day (PDD) method.The elevation of a 3-D region is parameterized using a hypsometric curve.Differences between the new SG model and the Marshall et al. (2011) approach are summarized in Table 1.Marshall and Clarke (1999) built their hypsometric curves, representing the basal elevation of a region, synthetically from the minimum, maximum and median elevation of the topography.We generate the hypsometric curves from the GEBCO (General Bathymetric Chart of the Oceans) 1 km resolution digital elevation model (DEM) (BODC, 2010).To select a region fitting the coarse resolution grid cell of the GSM (degrees), the GEBCO Cartesian coordinates are converted in degrees assuming the Earth as a perfect sphere with radius of 6370 km.The curves are obtained in a two-step process.First, the region is divided into N bins of equal altitude range.Then, to avoid empty bins, the surface elevation range of each empty bin is expanded (consequently decreasing the elevation ranges of the higher and lower adjacent bins) using as many adjacent bins as necessary until these bins represent approximately the same surface area.This process is repeated from the highest bin to the lowest as many times as necessary.

Hypsometric curves
We use 1 km resolution gridded data, so that the area of each bin is proportional to the number of high resolution grid cells assigned to that bin.The alternative of using equal areas in each bin has been discarded as it smooths the results in regions of low peak density.A total of 10 bins have been selected in this study, based on the comparison against high resolution modelling (see Sect. 3.1).Marshall and Clarke (1999) and Marshall (2002) used, respectively, 10 and 16 hypsometric bins in their hypsometric curves.
At any time step, t, the surface slopes, S t k , for SG bin k, from 1 (highest) to N (lowest), are computed from the surface elevation h t d,k and an effective length L k : To compute the slope at the lowest bin we assume an ice cliff boundary condition.The surface elevation h d,N+1 is set to the basal elevation of the lowest hypsometric bin h b,N .Instead of setting the hypsometric slopes with an effective length proportional to the horizontal extent of the CG cell (Marshall and Clarke, 1999), we account for the cubic dependence of ice flow on surface slope (see Eqs. 3 and 4).Specifically, for each hypsometric bin we compute the slope, S 0 k , as the cube root of the mean of the cube of the magnitude of the slopes from the GEBCO data.The effective length, L k , for SG bin k is computed from the basal elevation h b,k : where h b is the basal elevation.As no information is extracted about the basal elevation downstream of the terminal SG cell, the effective length at the first upstream bin is used at the lowest hypsometric bin.A small effective length can generate unrealistically high velocities in that bin.To avoid this, the lowest bin effective length is set to the mean effec-tive length of all the hypsometric bins when the altitude difference between the two lowest bins is less than 50 m.The flow line model requires an effective width, W , for the representation of flux between hypsometric bins.W k of each hypsometric bin is set to the total contact length of the SG cells assigned to the bin with adjacent lower hypsometric bin grid cells as detailed in Fig. 1.

Surface mass balance
We use the positive degree day method described in Tarasov and Peltier (1999) to compute accumulation and ablation from monthly mean temperature and precipitation.A constant environmental lapse rate adjusts the temperature to the ice surface elevation.A parameterization of the elevationdesertification effect (Budd and Smith, 1981) reduces the precipitation by a factor of 2 for every kilometre increase in elevation.Snow is melted first and the remaining positive degree days are used to melt ice with allowance for the formation of superimposed ice.The Supplement includes a more detailed description of the surface mass balance module.
The GSM and ISSM compute the surface mass balance using the same PDD method.

Ice thickness evolution
The prognostic equation for the ice thickness (H ) is computed, at each hypsometric bin, from the vertically integrated continuity equation as S is the surface slope and Ṁs is the surface mass balance rate (basal melt is computed in the CG GSM but ignored in the SG model).u is the vertically integrated ice velocity of the SG model derived using the shallow ice approximation (SIA).The effective diffusivity D is given by The creep exponent n of Glen's flow law is set to 3. A 0 is the creep parameter in Pa −3 s −1 , ρ = 910 kg m −3 and g = 9.81 m s −2 .Ice flow is insignificant when the ice thickness is on the order of 10 m.To avoid potential numerical instabilities, velocity is set to 0 if ice thickness is less than 20 m.
In their most recent experiments, Marshall et al. (2011) tuned their revised model against the present day total ice volume (encompassing 27 % uncertainties) in the eastern slopes of the Canadian Rockies.This tuning sets the ice rheology parameter for an ice temperature equivalence of approximately −40 • C. As the SG model is used for regions that are either starting to accumulate ice or else deglaciating, basal ice temperature (where most deformation occurs) is likely close to freezing.The creep parameter is therefore fixed to a value corresponding to an ice temperature of 0 • C using the Arrhenius relation from the EISMINT (European Ice Sheet Modelling Initiative) project (Payne et al., 2000).
Equation ( 3) is solved semi-implicitly using a central difference discretization as The superscripts t and t +1 represent respectively the current and the subsequent time step.x is the effective length L and y is the effective width W defined in Sect.2.1.1.At the highest bin, we assume that no ice flows into the region.At the lowest bin ice is allowed to flow out of the region.

Model limitations
The shallow ice approximation, used to compute fluxes, is formally invalid for high surface slopes such as present in mountain ranges like the Rockies.Simulating ice evolution over a 3-D terrain using a flow line model limits the ice flow representation.Ice flows from one SG bin to another using an average slope.Our model configuration does not allow for ice at high elevations to flow into an adjacent coarse grid cell.Nor does it allow for ice present at low elevations, in isolated regions having a closed drainage basin, to stay in a coarse grid cell.Moreover, the Arrhenius coefficient is computed with a constant ice temperature of 0 • C. High velocities processes, such as periodical surges (Tangborn, 2013;Clarke, 1987), cannot be represented since basal sliding and basal hydrology are not present in the current study.
The hypsometric length parameterization inferred from the surface slopes are correct for ice free regions, but it is only an approximation once the ice starts building up.At the lowest hypsometric bin, slopes are computed assuming ice cliff boundary conditions.
For the comparison against ISSM results, the surface temperature is downscaled with a lapse rate of 6.5 • C km −1 .This typical value used in glacial modelling represents the average free-air lapse rate observed in the troposphere which need not match the impact of changing surface elevation.Studies over Iceland, Greenland, Ellesmere Island and the Canadian high Arctic report seasonal changes in the surface temperature lapse rates over mountain regions and glaciers, with a mean annual value of about 3.7-5.3• C km −1 (Marshall and Losic, 2011).Rates as low as 2 • C km −1 are measured in the summer (Gardner et al., 2009).These values are tested in the GSM ensemble simulations where the lapse rate ranges between 4 and 8 • C km −1 .

GSM
The core of the GSM is a 3-D thermomechanically coupled ice sheet model.The model incorporates sub-glacial temperatures, basal dynamics, a visco-elastic bedrock response, climate forcing, surface mass balance, a surface drainage solver, ice calving and margin forcing.The grid resolution used for this study is 1.0 • longitude by 0.5 • latitude.
The thermomechanically coupled ice sheet model, described in detail in Tarasov and Peltier (2002), uses the vertically integrated continuity equation and computes the 3-D ice temperature field from the conservation of energy, taking into account 3-D advection, vertical diffusion, deformation heating, and heating due to basal motion.Velocities are derived from the SIA equations.The sub-glacial temperature field is computed with a 1-D vertical heat diffusion bedrock thermal model that spans a depth of 3 km (Tarasov and Peltier, 2007).If the base of the ice is at the pressure melting point, basal motion is assumed to be proportional to a power of the driving stress.The exponent for this Weertman-type power law is set to 3 for basal sliding and 1 for till deformation (detailed description in Tarasov andPeltier, 2002, 2004).The geographic location of the sediment cover is determined from different data sets (Laske and Masters, 1997;Fulton, 1995;Josenhans and Zevenhuizen, 1990).Ice shelf flow is approximated with a linear function of the gravitational driving stress.At the base, ice melt is also computed from the energy balance.
The visco-elastic bedrock response is asynchronously coupled to the GSM with a 100-year interval.This module is based on the complete linear visco-elastic field theory for a Maxwell model of the Earth (Tarasov andPeltier, 2002, 2007).
At the surface, the parameterized climate forcing (Tarasov and Peltier, 2004, 2006, 2007) is based on a linear interpolation between the present day climatology, derived from a 14-year average (1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995) of the 2 m monthly mean reanalysis; Kalnay et al., 1996), and a last glacial maximum (LGM) climatology.The LGM climatology field is derived from a linear combination of PMIP (Paleoclimate Modelling Intercomparison Project)I and II general circulation model results with the linear combination dependent on the maximum elevation of the Keewatin ice dome (PMIP I boundary conditions lacked a major Keewatin ice dome, while PMIP II had a large dome).The interpolation follows a glacial index derived from the GRIP (Greenland Ice-core Project) δ 18 O record at the summit of the Greenland ice sheet (Dans-gaard et al., 1993;World Data Center-A for Paleoclimatology, 1997).The surface mass balance is derived from this climatology using the same methodology as described in Sect.2.1.2.A surface drainage solver is fully coupled asynchronously at 100-year time steps.It diagnostically computes downslope drainage, filling any depressions (lakes) if drainage permits (Tarasov andPeltier, 2005, 2006).
The calving module, described in detail in Tarasov and Peltier (2004), is based on a height above buoyancy criterion with added mean summer sea surface temperature dependence.The inhibition of calving due to the presence of landfast sea ice is also parameterized.To reduce misfits between the model results and geological evidences of the ice configuration, the mass-balance forcing is nudged to promote compliance with geologically inferred deglacial margin chronologies (Tarasov and Peltier, 2004).

GSM and sub-grid model coupling
In this section, we describe how the SG model is embedded in the GSM and the conditions applied to activate or deactivate the SG model in each CG cell.The GSM is run, at all times, over all the CG cells and the ice thickness is updated in cases where the SG model is activated.Figure 2 gives a summary diagram of the coupling between the GSM and the SG model.

Interaction between the sub-grid model and the GSM
There is two-way communication between the GSM and SG models to exchange information about ice thickness, surface mass balance, and surface temperature.Ice in a CG cell is added to the SG level2 when the SG model switches from deactivated to activated for a given CG cell.The information about the ice evolution at the SG level is used to update the ice thickness, surface mass balance rate and surface temperature at the CG level.Marshall et al. (2011) export SG ice to the CG level only when the lowest SG bin is filled and the SG model for the given CG cell is deactivated in the time step.In our model, SG-CG ice transfer is as follows.While the SG model is activated, CG ice volume is set to that of the filled SG bins.The rationale for this is the assumption that over a large mountainous region, such as the Rockies, an ice sheet grows by building up ice in major valleys (represented by the lowest hypsometric bins) from ice accumulation and ice flowing in from surrounding mountain peaks.A SG bin is classified as filled once its surface elevation reaches the basal elevation of the adjacent higher bin.The surface mass balance rate and surface temperature of the CG cells are updated to the new elevations.When the SG model switches from activated to -SG model inactive.
-Surface elevation of the CG cell drop below a thresold.
-CG cell ice distributed over the SG hypsometric levels.deactivated, the total SG ice volume is transferred to the CG cell.

OFF
Once the SG model is reactivated in a CG cell during deglaciation, the ice volume present at the CG level is distributed over the different hypsometric bins.To account for the higher volume of ice in valleys, represented by the lowest hypsometric bins, the average of the following two massconserving distributions is used for SG initialization.The first is even distribution across every bin.The second keeps equal surface elevation for the lowest bins, starting from the lowest bin and using as many bins as necessary.Marshall and Clarke (1999) have no ice flux to adjacent CG cells when the SG model is active.In our model, ice transport between CG cells, computed with the GSM, is modified using SG information.We assume that only the ice present in the filled bins flows out of the coarse grid region; therefore, only a fraction of the CG flux is permitted.This fraction is computed as the area of the filled SG bins divided by the total CG cell area.To avoid double counting of this inter CG flux, the SG model does not compute flux out of the lowest bin through Eq. (3) when coupled to the GSM.At every iteration, the SG model accounts for the CG ice flux.For CG ice flux into a cell with active SG, the ice fills the lowest hypsometric bin.Once that bin reaches the elevation of the next higher bin, the remaining ice is used to fill up the two bins at the same elevation.This process is repeated using as many bins as necessary to redistribute all the ice.For CG ice flux out of the cell, the same amount of ice is removed from all the filled SG bins.If the total volume of ice to be removed is not reached using that region of the SG cell, the excess remaining is used to empty higher bins one after another.
The SG model flux module is coupled asynchronously and runs at half the SG mass balance time step.Glacial isostatic adjustment from the CG level is imposed on the SG basal topography.

Sub-grid model activation/deactivation
Unlike Marshall and Clarke (1999), the SG model is activated only in cells above sea level with rough topography.A terrain is considered rough when the differences between the maximum and minimum basal elevation is higher than 500 m.To account for regions such as the Alaska Peninsula where CG cells represent regions including basal topography both above and below sea level, cells where at least half of the area is above sea level are treated at the SG level.During inception, ice accumulates and can flow into valleys, filling them and thereby reducing the surface elevation variation.The SG treatment becomes less critical and is deactivated when the lowest hypsometric bin surface elevation reaches the bedrock elevation of the highest bin.This criterion keeps the SG model activated for a longer period of time than in Marshall and Clarke (1999) where the SG model is deactivated when ice reaches the lowest bin.During deglacia-tion, mountain peaks become uncovered and surface elevation variations increase, reaching a point where both ablation and accumulation are present.The SG model is reactivated when the ice thickness in the CG cell is lower than half of the difference between the basal elevation of the highest hypsometric bin and the basal elevation of the CG cell.This differs from Marshall and Clarke (1999), who use only SG information to set the threshold to a fraction of the variation in SG basal elevation.

Ice Sheet System Model (ISSM)
As a detailed description of the ISSM is given in Larour et al. (2012), only a brief description of the model components used in this study are presented here.The ISSM is a finite element 3-D thermomechanically coupled ice flow model.The mass transport module is computed from the depth-integrated form of the continuity equation.Using the ice constitutive equation, the conservation of momentum provides the velocities.The model offers the option of computing the velocities using full Stokes, higher-order Blatter-Pattyn, shelfy-stream or shallow ice approximation equations.The higher-order Blatter-Pattyn approximation is used in this study.As the velocity equations depend on the temperature, this field is computed from conservation of energy, including 3-D advection and diffusion.For this study, a new surface mass balance module identical to the one present in the sub-grid model, and detailed in Sect.2.1.2,has been incorporated into the ISSM.

Sub-grid model performance and tests
The SG model computation time for a 3000-year simulation, using 10 hypsometric bins, is about 0.02 s.At a resolution of 1 km and using 10 cpus, ISSM run time is about 2-5 h (depending of the topographic region used).The sub-grid model adds 3-6 h (depending of the parameter vector used) to the glacial cycle runtime over North America.

Comparison with ISSM
We compare 2 kyr ISSM and SG simulations, applying constant sea level temperature and precipitation over an inclined bed and 21 different test regions in the Canadian Rockies.These regions, for both the ISSM and SG simulations, have a dimension of 30 km by 60 km and we use a DEM of 1 km resolution.To improve correspondence between the ISSM and the SG model, the minimum ice thickness allowed in the SG model is set to 10 m.The boundary conditions at the ice margin in the ISSM are computed as an ice-air interface.To isolate the impact of using the SIA to represent fluxes in a mountainous region containing steep slopes in the hypsometric parameterization, our current experiments have no basal sliding.As glaciers can experience sliding in this type of region, the next stage of this project will include sliding.

Inclined plane test
The bed topography for this test is an inclined plane topography with a constant slope of 0.014 and a maximum basal elevation of 800 m.For this case, the accuracy of the SG model correlates with the number of hypsometric bins as shown in Fig. 3 (ice and velocity profiles are shown in Fig. S1 in the Supplement).Reducing the number of SG bins increases the surface gradient between two hypsometric bins and thereby the computed ice velocities.With 10 hypsometric bins, the ice volume simulated by the SG model can be as low as 40 % of the ISSM prediction.The misfits are not significant in simulations where no ablation is present (e.g. for a temperature set to −5 • C).

Rocky Mountains test
The SG model is tested on 21 regions from the Canadian Rockies, representing a wide range of topographic complexity (e.g.Fig. 4a), altitude (e.g.Fig. 4b) and slopes (e.g.Fig. 4c).The slopes of these regions are higher than in the inclined plane case.We focus on the results for simulations over the six test regions in Fig. 4 forced with sea level temperature of 0 • C and a desertification effect factor of 0.5.The results of other simulations, using different regions and with similar forcing as used in the inclined plane experiments, are not shown as they present similar misfits against ISSM results.
In contradiction with the simplified inclined plane configuration, increasing the number of hypsometric bins does not reduce the misfits with ISSM simulations (Fig. 5).The SG model does not account for the build-up of ice in closed drainage basins where no flow is permitted out of the region before a threshold elevation is reached.tion for the "real" topography scenario comes from topographic "jumps" not addressed in the SG model.Some high resolution adjacent grid cells belong to non-adjacent hypsometric bins.The ice flow between these two locations is not accurately captured.The number of "jumps" increases with the number of bins used (Fig. S2 in the Supplement).A total of 10 hypsometric bins are then used to limit this effect.Even so, the SG model generates 45 % less to 15 % more ice than ISSM simulations (25 % less on average), depending on the regional topographic characteristics.No relation was found between the geographic complexity and the performance of the model, as explained in Sect.3.2.

Test of alternative parameterizations
We examine the impact of including more topographic characteristics in the velocity parameterization.Characteristics considered include the flow direction, the terrain ruggedness (measured as the variation in three-dimensional orientation using a radius of 5 grid cells around the grid cell of interest), the sum of the squared slopes, the variance in the slopes, the number of local maxima (tested with radius sizes of 2, 6 and 10 grid cells) and the standard deviation of the surface elevation topography.
The ISSM and the sub-grid model were run until steady state (2 kyr) for simulations with a constant precipitation rate of 1 m yr −1 and a sea level temperature forcing of 0 • C. The parameters minimizing ice volume differences were selected using a stepwise multilinear regression fit.The flow direction and the mean of the squared slope do not reduce the misfits.The slope variance does not improve the results when combined with the remaining two parameters (elevation standard deviation and terrain ruggedness).When used alone, it does reduce the errors but not as well as when the standard deviation of the topography is used.The terrain ruggedness and the peak density both represent the same physical characteristics and do not improve the results when used alone.Improvements are obtained when combined with the standard deviation of the topography.However, the improvement is not greater than with the standard deviation alone.The standard deviation of the topography is the parameter that correlates the most with the misfits.The average absolute value of the differences between the SG model and ISSM average ice thickness is 61 m.This difference is reduced to 21 m (see Fig. 6) when the regression model generated using the standard deviation of the topography is used.More details about the results of the stepwise regression fits are provided in the Supplement.
To explore potential improvement from accounting for the standard deviation of the high resolution topography, S SD , we test the following parameterization of the velocity, u 1 : This equation is used in a simulation initialized with the ice thickness, velocities and slopes of ISSM values at steady state.The parameters P 1 , P 2 and P 3 (respectively 4.87, 0.016 and 2.8) are obtained using a least-squares approach that minimizes the differences between the velocities computed by ISSM and the SG model after one iteration (0.01-year).
The lowest hypsometric bin has the most significant misfits (e.g.Fig. S4 in the Supplement).This is likely related to the margin ice cliff slope parameterization.To try to correct this, we test the following parameterization for the lowest hypsometric bin velocity: Using the same least-squares approach as above, the parameters P 4 and P 5 are respectively set to 5924.4 and −1.6383.These two parameterizations do not reduce the ice thickness differences with ISSM transient results (see Fig. 7).Ice thickness, velocities and slopes over the six regions analysed are presented for the different parameterizations in Fig. S5 of the Supplement.As the model is highly non-linear, the improvement generated by the least-squares fit method for an initialization with ISSM steady state conditions does not persist over 1000-year runs.
The following modifications of the current version of the SG model have been explored but did not improve the model.The central difference discretization of the ice thickness in the effective diffusivity coefficient was replaced by an upwind scheme.Simulations with different values of the Arrhenius coefficient, the power of the ice thickness and the slope, in Eq. ( 4), were analysed.An extra parameter was added in the velocity equation to account for neglected stresses.Turning off the internal SG model flux term increased the misfits with ISSM simulations by a minimum of 100 % (as shown in Fig. 8).The basal elevation downstream of the terminus has been computed using a linear extrapolation of two or three upstream bins.The lowest hypsometric bin effective length generated with these basal elevations did not reduce the misfits with ISSM results.

Behaviour of the sub-grid model in the GSM
We present results of simulations over the last glacial cycle.The 39 "ensemble parameters" of the GSM (attempting to capture the largest uncertainties in climate forcing, ice calving, and ice dynamics) have been subject to a Bayesian calibration against a large set of palaeoconstraints for the deglaciation of North America, as detailed in Tarasov et al. (2012).We use a high-scoring sub-ensemble of 600 parameter vectors from this calibration to compare the GSM behaviour when the SG model is turned on and off.The primary supplement of Tarasov et al. (2012) includes a tabular description of the 39 ensemble parameters as well as input data sets.For the purposes of clarity and computational cost, we examined model sensitivity to different coupling and flux parameters using five parameter vectors (of the 600 members  ensemble) that gave some of the best fits to the calibration constraints.As these five parameter vectors display similar behaviour, we present sensitivity results using the parameter vectors for the two runs described in detail in Tarasov et al. (2012) (identified in that paper as runs nn9894 and nn9927).
For ease of interpretation, the ice volumes are presented as eustatic sea level (e.s.l.) equivalent3 .

Last glacial cycle simulations over North America
The SG model can significantly alter the pattern of ice accumulation and loss.Figure 9 shows an example, for one of the parameter vectors of the ensemble of simulations, where SG ice accumulates while it melts at the CG level (Fig. 9a), and an example where CG ice is about 60 % greater than the SG ice (Fig. 9b).The ensemble of simulations of the last glacial cycle over North America with the SG model activated generates, on average, between 0 and 1 m e.s.l.more ice than when the SG model is turned off (Fig. 10).
The impact of the SG model depends, however, on the climate forcing and the ice sheet extent and elevations.During inception, when the SG model is turned on, ice accumulating in higher regions flows downhill and accumulates in regions close to the ELA and in valleys (Fig. 11).This allows, for example, ice to build up in the northern part of Alaska.For typical runs, the ice generated by the SG model in the Alaskan Peninsula is, however, insufficient as compared to geological inferences (Dyke, 2004).The ensemble mean and standard deviation of the differences between runs with SG on and off at 110 ka, are respectively 0.4 and 1 m e.s.l.However, at specific time slices, the differences can be much larger.Once the ice sheet has grown to a sizeable fraction of LGM extent, for example at 50 ka, the standard deviation of the ensemble-run differences (between SG on and off) reaches 5 m e.s.l. Figure 12 shows an example where ice in a region of low altitude in the centre of Canada is not allowed to grow when the SG model is used.On the other hand, a simulation using different ensemble parameters generates ice in this region only when the SG model is turned on (Fig. S6 in the Supplement).In extreme cases, differences can reach tens of m e.s.l.(Fig. S7 in the Supplement).We could not identify a reason for the strong sensitivity of ice volume around 50 ka other than the inherent non-linearity of the GSM.

Sensitivity of the model to different flux and coupling parameters
The accounting of SG fluxes has varying impacts over a glacial cycle simulation (Fig. 13).At 50 ka, for example, the total ice volume with parameter vector nn9894 is reduced by a.
b. 50 % when SG fluxes are included.During inception, on the other hand, inclusion of SG fluxes increases the total amount of CG ice (Fig. 14, again with nn9894).
To better understand the range of responses to CG ice flow between grid cells that have SG activated, three case scenarios can be considered.Case 1: ice flows out of the lowest SG bins located above the ELA into the lowest SG bins located above the ELA of another CG cell.There is limited impact of not allowing ice to flow out of the CG cell as in both cases ice accumulates.Case 2: ice flows out of the lowest SG bins located above the ELA into the lowest SG bins located below the ELA of another CG cell.In that case, turning off the fluxes between CG cells tends to reduce the total melt.Case 3: ice flows out of the lowest SG bins located below the ELA into the lowest SG bins located below the ELA of another CG cell.Ice flowing into lower SG bins generates higher melting rates, so permitting fluxes between CG cells will in this case tend to increase ice mass loss.In cases 2 and 3, the combination of ice flowing below the ELA from the adjacent CG cell and from the bins above the ELA can raise the surface elevation of lower bins above the ELA and reduce the melt.Depending on the proportion of each of these cases, not al-  lowing ice fluxes out of coarse grid cells with SG activated generates higher or lower ice volumes (Fig. 13).Moreover, 50 ka is an example of a 60 % increase of the total ice volume when the fluxes out of coarse grid cells (with SG activated) are not allowed.As a counter-example, 35 ka presents a case where turning off the fluxes out of (SG activated) coarse grid cells decreases the total ice volume.
With Marshall et al.'s (2011) flux equation, differences between runs with SG fluxes turned on versus off are negligible over the full glacial cycle (Fig. 14).
As described in Sect.2.3.1, the CG ice thickness used by the GSM conserves the ice volume of the filled SG bins (volume conservation, VC, method).As this ice is redistributed over the total area of the coarse grid cell, the surface elevation of the ice, and consequently the fluxes, are underestimated.The surface gradient between adjacent cells is then lower than the gradient at the SG level.We tested setting the CG surface elevation to the maximum value between the surface elevation of the coarse grid cell and the lowest hypsometric bin (surface conservation (SC) method).We also implemented a method using the maximum surface elevation generated by the two former methods (maximum conservation (MC) method).During inception (between 118 and 114 ka) the VC method generates between 10 and 20 % (which is equivalent to 0.5-1 m e.s.l.) more ice than the two other methods (Fig. 15).During the first 60 kyr of simulation, the difference in total ice volume stays under 1 m e.s.l.independently of the flux redistribution methods (Fig. S8 in the Supplement).Between 60 ka and the LGM, the SC method generates between 1 and 12 m e.s.l.less ice than the two other methods (Fig. S8).The VC method was used for the ensemble runs as it generates more ice over the Alaska Peninsula, northern and southern mountain ranges, thereby reducing misfits against geological inferences.Figure 16 shows the results of the glacial cycle simulation when the SG model is turned off and when the minimum altitude variation SG activation threshold is set to 50, 150, 300 and 500 m.A non-linear dependence on the threshold can be observed.At 50 ka, for example, setting the threshold to 50 m generates the lowest total ice volume while a threshold of 150 m leads to the highest ice volume.The difference between these two runs is 34.5 m e.s.l. at 50 ka.Thresholds of 300 and 500 m generate intermediate total ice volumes.Moreover, simulations using different parameter vectors (not shown) result in different behaviours.No conclusion could be drawn about the optimal threshold.

Conclusions
Our new sub-grid surface mass balance and flux model extends the initial work of Marshall and Clarke (1999) and Marshall et al. (2011).The evaluation of the model, done for the first time against results from a high resolution higher order model (ISSM), demonstrates that -Depending on the regional topographic characteristics, the new SG model simulates ice volumes 45 % lower to 15 % higher than simulated by the ISSM (using 10 hypsometric bins).Increasing the number of hypsometric bins to more than 10 did not reduce misfits for simulations over rough topographic regions extracted from the Canadian Rockies.
-Turning off the SG internal fluxes increases the ice volume misfits with ISSM simulations by a minimum of 100 %. -Increasing the number of topography characteristics used in the SG model, as suggested by Marshall and Clarke (1999), did not reduce the misfits with the high resolution model during transient runs.
An ensemble of simulations over the last glacial cycle of the North American ice complex shows, on average, an increase of ice generated with inclusion of the SG model.The ensemble mean for each time step is between 0 and 1 m e.s.l., with a standard deviation of a minimum of twice the mean and reaching 5 m e.s.l. at 50 ka).At the end of inception, at 110 ka, the increase of ice volume from SG model inclusion is still insufficient over the Alaska Peninsula when compared to geological inferences.Over the glacial cycle, the SG model generates different patterns of ice extent.In some instances, the SG model prevents ice growth, while in others it enables extra ice build up over thousands of square kilometres.
Simulated ice evolution is sensitive to the treatment of ice fluxes within the SG model and between the SG and CG levels.
-The flux term has an important impact on the SG model.
Not allowing ice to flow between hypsometric bins increases the total ice volume with a maximum increase of 50 % at 50 ka (in a glacial cycle run).During inception, however, the flux module can generate more ice.Different parameterizations of the flux term impact the results.A SG ice rheology parameter corresponding to ice at about −40 • C (as used in Marshall et al., 2011) generates the same amount of ice during inception as when the flux term is off.
- -Not allowing ice to flow out of a CG cell where SG is activated increases or decreases the total ice volume depending of the ice configuration.At 50 ka, the total increases by 60 %.
-The ice configuration from simulations over the last glacial cycle of North America is sensitive to the choice of SG to CG ice redistribution scheme.
We have identified the representation of SG fluxes between CG cells to be a challenging issue that can significantly impact modelling ice sheet evolution.
We have shown that the above geometric and ice dynamics factors can have significant impacts on modelled ice sheet evolution (with up to a 35 m e.s.l.difference in North American ice volume at 50 ka).Therefore, significant potential errors may arise if sub-grid mass-balance and fluxes are not accounted for in the coarse resolutions required for glacial cycle ice sheet models.Other alternatives to the hypsometric parameterization, such as running a high resolution SIA model in the region of rough topography, could be considered.One issue we have not examined is the downscaling of the climatic forcing.Temperature and especially precipitation can exhibit strong vertical gradients in mountainous regions.Whether this can have significant impact on CG scales is unclear.Improvements of the precipitation representation are possible using, for instance, a linear model of orographic precipitation for downscaling climatic inputs (Jarosch et al., 2012).

Figure 1 .
Figure 1.Schematic representation of the effective width of the 7th hypsometric bin for a region of 10 km by 10 km.Each square represents a high resolution (1 km) grid cell.The numbers define the hypsometric bin these SG grid cells belong to.The total length of all red lines (14 km) represents the effective width for the 7th bin.

Figure 2 .
Figure 2. Communication between the GSM and the SG model for one CG cell.

Figure 3 .
Figure 3. SG model vs. ISSM differences over an idealized inclined plane terrain.Average ice thickness differences (SG model − ISSM) are presented for simulations using different temperatures, desertification effect factors and number of hypsometric bins.

Figure 4 .
Figure 4. Topography characteristics for six regions over the Canadian Rockies.(a) summarizes surface elevations, (b) the hypsometric curves, and (c) the mean slope for each hypsometric bin.

Figure 5 .Figure 6 .
Figure5.Ratio of the SG model over ISSM total ice volume for six different regions in the Rockies as a function of hypsometric bins.The simulations were run until steady state with a constant sea level temperature of 0 • C and a desertification effect factor of 0.5.The steady state ice thicknesses, velocities and slopes from the ISSM and the SG model (using 10 hypsometric bins) are presented in Fig.S3in the Supplement.

Figure 7 .
Figure7.Average ice thickness root mean square error (RMSE) between the ISSM and the SG model for different topographic regions.Simulations are run over 2 kyr using a constant precipitation rate of 1 m yr −1 and a sea level temperature forcing of 0 • C. Different SG parameterizations are presented.Para 1 is the standard deviation of the topography parameterization (Eq.6) and Para 2 the lowest hypsometric slope parameterization (Eq.7).

Figure 8 .
Figure 8. Surface elevation generated by the ISSM (solid blue line), the SG model with no flux term, using 5 and 10 hypsometric bins, (dotted lines) and the SG model including the flux term (solid thin red line).These simulations use a constant sea level temperature of 0 • C and a desertification effect factor of 0.5.Results are shown at steady state after 2 kyr for six different regions with different topographic characteristics.

Figure 9 .
Figure 9. Elevation comparisons when the SG model is turned on (blue) or off (red) at different time steps using the parameter vector nn9894.h d 10 years is the CG surface elevation after 10 years.h dhyps 10 years is the SG surface elevation.h b is the basal elevation.(a) and (b) represent cases where the ELA is above and below the coarse grid surface elevation.

Figure 10 .
Figure10.Ensemble mean (solid red line) and standard deviation (dotted blue line) eustatic sea level equivalent of the total ice volume differences when the SG model is turned on and off, for an ensemble run over the last glacial cycle.

Figure 11 .
Figure 11.Ice field during inception at 115 ka for a simulation using one of the parameter vectors that generates best fits to the calibration constraints (nn9894).(a) Ice thickness with SG turned on.(b) Ice thickness differences between simulations with the SG model turned on and off.Zero differences are presented in the same colour as the continent.

Figure 12 .
Figure 12.Ice field at 50 ka for a simulation using parameter nn9894.(a) Ice thickness with SG turned on.(b) Ice thickness differences between simulations with the SG model turned on and off.

Figure 13 .Figure 14 .
Figure13.Total ice volume evolution for a simulation using parameter vector nn9894."flux on" and "flux off" both include the SG surface mass balance calculations but the latter has no SG ice fluxes."NofluxOut" has SG on, but no SG ice flux between coarse grid cells.The "SG OFF" line is most of the time hidden under the "flux off" line.

Figure 15 .
Figure15.Total ice volume evolution for a simulation over North America during inception with the SG model turned on (SG on) using the parameter vector of run nn9927.Different methods of ice redistribution at the CG level are compared.VC is for ice volume conservation, SC for surface elevation conservation and MC uses the maximum of the previous two methods."SG off" represents a run where the SG model has been turned off.

Figure 16 .
Figure16.Total ice volume evolution for a simulation using parameter vector nn9894.Different curves represent simulations with different minimum altitude variation thresholds used for the SG activation.
The flux term used in theMarshall et al. (2011)  study, with the ice rheology parameter representing ice at www.geosci-model-dev.net/8/3199/2015/Geosci.Model Dev., 8, 3199-3213, 2015 about −40 • C, generates an ice volume higher than when a flux parameterization with a rheology value representing ice at about 0 • C is used.

Table 1 .
Differences between our new SG model and the Marshall and Clarke (1999)/Marshall et al. (2011) models.