Response to reviewer # 1 of “ A new sub-grid surface mass balance and flux model for continental-scale ice sheet modelling : validation and last glacial cycle ”

first sentence I would not write like this in the abstract. Although I agree that typical grid resolution at the moment is around 10 to 50km for long-term computation, this is not always a necessary condition. Rather, I would state simply that this resolut ion is a current typical configuration (instead of ‘need to be run. . . .’). "need to be run" was modified to "are typically run". p3038, L26. better to delete ‘coarse’ (I feel it a bit subjective ) as the same reason above. I would just state the fact simply, at this stage. The following sentences n aturally drive us this resolution as ‘coarse’ one. "coarse" has been removed here. But we have added it in quotations in th e abstract to help define what we mean by coarse. p3039, L6 ‘the mean surface elevation’ of a coarse grid? "the mean surface elevation" has been changed to "the mean surface elev ation of a coarse grid cell". p3039, L8, citing Abe-Ouchi et al.: The first part is somewhat misleading and confusing. Van den Berg et al explicitly discuss the sensitivity of ice-sheet evolution to the g rid resolution, while Abe-Ouchi et


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 coarse resolutions of about 10 to 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 metres to a few kilometres, and therefore have to be parametrized.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 is below the equilibrium line altitude (ELA), ice ablation is overestimated (e.g.Tarasov and Peltier, 1997).Thus, lower grid resolution can lead to temporal and spatial errors in ice sheet inception (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 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 parametrized form (Seth et al., 1994;Leung and Ghan, 1995;Marshall and Clarke, 1999;Giorgi et al., 2003;Ke et al., 2013).For example, 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 parametrized with different hypsometric levels, representing a discrete number of elevations and their associated area.In addition to ablation and accumulation at each SG level, 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 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 Figures

Back Close
Full SG and coarse grid (CG) cells 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 Marshall and Clarke (1999) and Marshall et al.'s (2011) approach.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 accounts for SG slope statistics.An effective width is added for the representation of the ice fluxes between SG levels.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 level.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 levels, 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 the SG model on the Glacial Systems Model (GSM, formerly the MUNGSM) 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
This model expands the approach of Marshall and Clarke (1999) and Marshall et al. (2011).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.Introduction

Conclusions References
Tables Figures

Back Close
Full The elevation of a 3-D region is parametrized using a hypsometric curve.Differences between the new SG model and the Marshall et al. (2011) approach are summarized in Table 1.

Hypsometric curves
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 1 km resolution digital elevation model (DEM) (BODC, 2010).The curves are obtained in a two-step process.First, the region is divided into 10 levels (or bins) of equal altitude range and the number of grid cells present in each of these bins is computed.Then, the size of these bins is updated to avoid empty levels.The alternative of using equal areas in each level has been discarded as it smooths the results in regions of low peak density.10 levels 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 levels in their hypsometric curves.
Instead of setting the hypsometric slopes with an effective length proportional to the horizontal extent of the CG cell, the surface slope length accounts for the cubic dependence of ice flow on surface slope.Specifically, for each hypsometric level, we compute the cube root of the mean of the cube of the magnitude of the slopes.The effective length, L, used to update these slopes when ice starts to build up, is computed for each level as: where k represents the different hypsometric levels, from 1 to N, and 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 level is used at the lowest hypsometric level.A small effective length can generate unrealistically high velocities Introduction

Conclusions References
Tables Figures

Back Close
Full at that level.To avoid this, the lowest level effective length is set to the mean effective length of all the hypsometric levels when the altitude difference between the two lowest levels is less than 50 m.At any time step, t, the surface slopes are updated using: To compute the slope at the lowest level, the surface elevation h d ,N+1 is set to the basal elevation of the lowest hypsometric level h b,N .
We include a new width parameterization to improve the representation of flux between hypsometric levels.The effective width of each hypsometric level is set to the number of grid cells, multiplied by the spatial resolution, that are in contact with adjacent lower hypsometric levels grid cells.

Ice velocity
The vertically integrated ice velocity of the SG model, u, is computed using the shallow ice approximation (SIA): where 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 project (Payne et al., 2000).

Surface mass balance
Positive degree day methods have been widely used in surface mass balance models (Johannesson et al., 1995;Tarasov and Peltier, 1999;Hock, 2003;Shea et al., 2009;Barrand et al., 2013).Here, we use the PDD method described in Tarasov and Peltier (1999) to compute the ice ablation and accumulation from the temperature and precipitation fields.Ablation rates are derived from monthly mean temperature (T m ).To increase the accuracy, hourly temperatures are considered normally distributed, with a standard deviation (σ PDD ) of 5.5 • C, around the monthly mean.A lapse rate is also used to adjust the temperature forcing to the ice surface elevation.The number of days where the temperature is above 0 • C in a year is computed as: The amount of snow and ice are assumed to melt proportionally to the number of positive degree days.Snow is melted first and the remaining positive degree days are used to melt ice.The ablation rate factors for snow (γ snow ) and ice (γ ice ) have a mean June/July/August temperature (T jja ) dependence extracted from energy balance modelling (Braithwaite, 1995;Tarasov andPeltier, 2002, 2003):

Conclusions References
Tables Figures

Back Close
Full and In addition, the amount of superimposed ice for a year is computed as per Janssens and Huybrechts (2000): where P r is the rainfall in a year, P s is the snow fall in a year, M is the snow melt in a year, 2.2 is the capillarity factor, d is the active thermodynamic layer (set to 1 m), c i is the ice specific heat capacity (152.5 + 7.122T ) in J kg −1 K −1 , L is the latent heat fusion (3.35 × 10 5 ) in J kg −1 , and T surf is the surface temperature.
A normal distribution of the hourly temperature is also used to compute the amount of snow accumulation from the precipitation.A lower standard deviation σ RS = σ PDD − 0.5 is assumed in that case to account for the smaller temperature variability during cloudy days.Precipitation is assumed to fall as snow when the temperature is below 2 A parameterization of the elevation-desertification effect (Budd and Smith, 1981) reduces the precipitation by a factor of two for every kilometre increase in elevation.This exponential reduction is a function of the surface height difference to that of presentday with an ensemble parameter threshold for activation (Tarasov and Peltier, 2006;Tarasov et al., 2012).Figures

Back Close
Full

Ice thickness evolution
The ice thickness is computed from the vertically integrated continuity equation as: Ṁs represents the surface mass balance rate (basal melt is ignored).d represents the effective diffusivity given by: This prognostic equation 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.
At the highest level, we assume that no ice flows into the region.At the lowest level 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 Introduction

Conclusions References
Tables Figures

Back Close
Full flows from one SG level to another using an average slope.Our model configuration does not allow for ice at high elevations to flow in an adjacent coarse grid cell, or 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 as basal sliding is 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 level, slopes are computed assuming ice cliffs 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 seasonally changes in the surface temperature lapse rates over mountain regions and in the glaciers boundary layer with a mean annual value of about 3.7 to 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 .

Ice Sheet System Model (ISSM)
As 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 Sec.2.1.3,has been incorporated in ISSM.

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 three-dimensional 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 years interval.This module is based on the complete linear visco-elastic field theory for a Maxwell model of the Earth (Tarasov andPeltier, 2002, 2007).Introduction

Conclusions References
Tables Figures

Back Close
Full At the surface, the parametrized 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 I or II general circulation models results with weighting dependent on the maximum elevation of the Keewatin ice dome.The interpolation follows a glacial index derived from the GRIP δ 18 O record at the summit of the Greenland ice sheet (Dansgaard 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.3.A surface drainage solver is fully coupled asynchronously at 100 yr time step.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 annual surface temperature dependence.The inhibition of calving due to the presence of landfast sea-ice is also parametrized.To reduce misfits between the model results and geological evidences of the ice configuration, mass-balance forcings are nudged to promote compliance with geologically inferred deglacial margin chronologies (Tarasov and Peltier, 2004).
The GSM has been subject to a Bayesian calibration against a large set of paleo constraints for the deglaciation of North America, as detailed in Tarasov et al. (2012).
We use a high-scoring sub-ensemble of 600 runs from this calibration.

GSM and sub-grid model coupling
In this section, we describe the conditions applied to activate or deactivate the SG model in each synoptic cell and how the two models exchange information about ice thickness, surface mass balance, and surface temperature.Introduction

Conclusions References
Tables Figures

Back Close
Full 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 Alaskan Peninsula where synoptic 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 turned off when the lowest hypsometric level surface elevation reaches the bedrock elevation of the highest level.
This criterion keeps the SG model activated for a longer period of time than in Marshall and Clarke (1999) where the SG model is turned off when ice reaches the lowest level.
During deglaciation, 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 at the CG level is lower than half of the difference between the basal elevation of the highest hypsometric level and the basal elevation of the CG cell.This differs from Marshall and Clarke (1999) who uses only SG information to set the threshold to a fraction of the variation in SG basal elevation.

Interaction between the sub-grid model and the GSM
There is two way communication between the CG and SG models.CG ice is added to the SG level when the SG model is turned on, and 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.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 levels) from ice accumulation and ice flowing from surrounding mountain peaks.A SG level is classified as filled once its surface elevation reaches the basal elevation of the adjacent higher level.The surface mass balance rate and surface temperature of the synoptic cells are updated to the new elevations.When the SG model is turned off, the total SG ice volume is transferred to the CG cell.
Once the SG model is re-activated in a synoptic grid cell during deglaciation, the ice volume present at the CG level is distributed over the different hypsometric levels.To account for the higher volume of ice in valleys, represented by the lowest hypsometric levels, the average of the following two mass-conserving distributions is used for SG initialization.The first is even distribution across every level.The second keeps equal surface elevation for the lowest levels, starting from the lowest level and using as many levels 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 below the lowest unfilled level flows out of the coarse grid region; therefore, only a fraction of the CG flux is permitted.This ratio is computed as the area of the SG levels below the lowest unfilled level over the total area.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 level.Once that level reaches the elevation of the next higher level, the remaining ice is used to fill up the two levels at the same elevation.This process is repeated using as many levels as necessary to redistribute all the ice.For CG ice flux of the cell, the same amount of ice is removed from all the SG levels lower than the lowest unfilled level.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 levels one after another.Figures

Back Close
Full The SG model flux module is coupled asynchronously and runs at half the SG mass balance time step.Glacial isostatic adjustment from the CG model is imposed on the SG basal topography.

Sub-grid surface mass balance model performance
We compare the results of the SG model and ISSM before testing alternative parameterizations.

Comparison with ISSM
We use 2 kyr simulations, applying constant sea level temperature and precipitation over an inclined bed and 21 different regions in the Canadian Rockies.These regions have a dimension of 30 km by 60 km and the topographic data are available at a resolution of 1 km.To improve correspondence between ISSM and SG model, the minimum ice thickness allowed in the SG model is set to 10 m, and no basal sliding is allowed in ISSM.The boundary conditions at the ice margin in ISSM are computed as an ice-air interface.

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 750 m.For this case, the accuracy of the SG model correlates with the number of hypsometric levels as shown in Fig. 1 (ice and velocities profiles shown in the Supplement, Fig. S1).Reducing the number of SG levels increases the surface gradient between two hypsometric levels and thereby the computed ice velocities.With 10 hypsometric levels, the ice volume simulated by the SG model can be as low as 40 % of ISSM prediction.The misfits are not significant in simulations where no ablation is present (e.g. for a temperature set to −5 • C).Introduction

Conclusions References
Tables Figures

Back Close
Full

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. 2a), altitude (e.g.Fig. 2b) and slopes (e.g.Fig. 2c).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. 2 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 forcings as used in the inclined plane experiments, are not shown as they present similar misfits with ISSM.
In contradiction with the simplified inclined plane configuration, increasing the number of hypsometric levels does not reduce the misfits with ISSM simulations (Fig. 3).
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.Another complication for the "real" topography scenario comes from topographic "jumps" not addressed in the SG model.Some high resolution adjacent grid cells belong to nonadjacent hypsometric levels.The ice flow between these two locations is not accurately captured.The number of "jumps" increases with the number of levels used (Fig. S2 in the Supplement).10 hypsometric levels 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 Sec.3.2.

Test of alternative parameterizations
We examine the impact of including more topographic characteristics such as 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 ra-Introduction

Conclusions References
Tables Figures

Back Close
Full The ISSM and the sub-grid model were run until steady state (2 kyr) for simulations with a constant precipitation rate of 1 mm 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 slope squared do not reduce the misfits.The slope variance does not improve the results when combined with the other two parameters.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. 4) 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 std , 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 used in the stepwise regression fit section.The parameter P 1 , P 2 and P 3 (respectively 4.87, 0.016 and 2.8) are obtained using a least square approach that minimizes the differences between the velocities computed by ISSM and the SG model after one iteration (0.01 yr).Introduction

Conclusions References
Tables Figures

Back Close
Full The lowest hypsometric level have 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 level 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. 5).Ice thickness, velocities and slopes over the six regions analyzed 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 thousand 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. ( 3), were run.An extra parameter was added in the velocity equation to account for neglected stresses.Turning off the internal SG model flux term increased significantly the misfits with ISSM simulations as shown in Fig. 6.The basal elevation downstream of the terminus has been computed using a linear extrapolation of two or three upstream levels.The lowest hypsometric level effective length generated with these basal elevations did not reduce the misfits with ISSM results.

Conclusions References
Tables Figures

Back Close
Full We present results from a 600 member ensemble of simulations of the last glacial cycle.
We also examine model sensitivity to different coupling and flux parameters.For ease of interpretation, the ice volumes are presented as eustatic sea level (ESL) equivalent 2 .

Last glacial cycle simulations over North America
The SG model can significantly alter the pattern of ice accumulation and loss.Figure 7 shows an example of SG ice accumulating while melting in the CG model (Fig. 7a), and an example where CG ice is about 60 % greater than the SG ice (Fig. 7b).
An ensemble of simulations of the last glacial cycle over North America with the SG model activated generate, on average, between 0 and 1 mESL more ice than when the SG model is turned off (Fig. 8).
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 activated, ice accumulating in higher regions, flows downhill and accumulates in regions close to the ELA and in valleys (Fig. 9).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 run mean and standard deviation of the differences between runs with SG on and off at 110 ka, are respectively 0.4 and 1 mESL.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 can reach 5 mESL.Figure 10 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  Supplement).In extreme cases, differences can reach tens of mESL (Fig. S7 in the Supplement).

Sensitivity of the model to different flux and coupling parameters
Five of the best fit to calibration constraints parameter vectors obtained from the glacial cycle calibration are used in this section.We focus on results from one of the ensemble parameter vectors as all 5 runs, unless otherwise stated, display similar behaviour.
Turning off the SG fluxes has varying impacts over a glacial cycle simulation (Fig. 11).At 50 ka, for example, the flux module allows ice, accumulating in higher elevation regions, to flow into the ablation zone and reduces the total ice volume by 50 %.The same process is observed during deglaciation (Fig. S8b in the Supplement).During inception, the flux module transports ice to below the ELA at a greater rate than it can be melted, thereby increasing the total amount of CG ice (Fig. S8a in the Supplement).With Marshall et al.'s (2011) flux equation, differences between runs with SG fluxes turned on vs. off are smaller (Fig. S8 in the Supplement).
As expected from the comparison between ISSM and the SG model results, not allowing ice fluxes out of coarse grid cells with SG active generates ice volumes up to 60 % higher at 50 ka (Fig. 11).
As described in Sec.2.4.2, the CG ice thickness used by the GSM, conserves the ice volume of the SG levels under the lowest unfilled level (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 surface elevation to the maximum value between the surface elevation of the coarse grid cell and the lowest hypsometric level (Surface Conservation, SC, method).We also implemented a method using the maximum of the two former methods (Maximum Conservation, MC, method).During inception, the VC method generates up to 20 % more ice than the two other methods (Fig. 12).This method was used for the ensemble runs as it generates more ice over Alaska, thereby reducing misfits Introduction

Conclusions References
Tables Figures

Back Close
Full against geological inferences.After inception, the flux redistribution methods have different impacts (Fig. S9 in the Supplement) as they start from different ice configuration.The impact of the minimum altitude variation SG activation threshold is more complex.Figure 13 shows a non-linear dependence on the threshold.Simulations using different parameter vectors 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: -Accounting for accumulation and ablation at different SG levels alters the ice volume evolution in a GSM cell.
-Depending on the regional topographic characteristics, the new SG model simulates ice volumes 45 % lower to 15 % higher than ISSM (using 10 hypsometric levels).Increasing the number of hypsometric levels to more than 10 did not reduce misfits for simulation over rough topographic regions extracted from the Canadian Rockies.
-Turning off the SG internal fluxes significantly increase the ice volume misfits with ISSM simulations.
-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.The topographic characteristics tested in the alternative parameterization were: the flow direction, the terrain ruggedness, the sum of the squared slopes, the variance in the slopes, the number of local

Conclusions References
Tables Figures

Back Close
Full 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 mESL (with a standard deviation up to 5 mESL at 50 ka) of ice generated with inclusion of the SG model.At the end of inception, at 110 ka, the SG model increases the ice volume on average by 0.4 mESL (with a standard deviation up to 1 mESL) but still does not generate sufficient ice in the Alaskan peninsula.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 levels generates up to 50 % more ice 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 the total amount of ice by up to 60 % at 50 ka.
-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.Introduction

Conclusions References
Tables Figures

Back Close
Full 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.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.One issue we have not examined is the downscaling of the climatic forcing.Temperature and especially precipitation can exhibit strong vertical gradients in mountainous region.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).
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | H is the ice thickness and n represents the creep exponent parameter of Glen's flow law and 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 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | a detailed description of ISSM is given in Larour et al. (2012), only a brief description of the model components used in this study are presented here.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 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Marshall et al. (2011)  export SG ice to the CG level only when the lowest SG level is filled and the SG model for the given CG cell is deactivated in the timestep.In our model, SG/CG ice transfer is as follows.While the SG model is turned on, CG ice volume is set to that of the SG levels below the lowest unfilled level.The rationale for Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | dius sizes of 2, 6 and 10 grid cells) and the standard deviation of the surface elevation topography in the velocity parameterization.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 4 Behaviour of the sub-grid model in the GSM Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | maxima and the standard deviation of the surface elevation topography.The latter did improve the fit in the single time step tests.
The flux term used inMarshall et al. (2011)  study, with the ice rheology parameter representing ice at 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.
Discussion Paper | Discussion Paper | Discussion Paper | Author contributions.Kevin Le Morzadec and Lev Tarasov designed the experiments.Kevin Le Morzadec developed the SG model code and performed the simulations.Kevin Le Morzadec and Lev Tarasov coupled the SG model into the GSM.Mathieu Morlighem and Helene Seroussi supported the installation of ISSM and helped including the new module in ISSM.Kevin Le Morzadec prepared the manuscript with contributions from Lev Tarasov and the other co-authors.Lev Tarasov heavily edited the manuscript.Discussion Paper | Discussion Paper | Discussion Paper | Leung, L. R. and Ghan, S.: A subgrid parameterization of orographic precipitation, Theor.Appl.Climatol., 52, 95-118, 1995.3039 Marshall, S.: Modelled nucleation centres of the Pleistocene ice sheets from an ice sheet model with subgrid topographic and glaciologic parameterizations, Quaternary Int., 95, 125Discussion Paper | Discussion Paper | Discussion Paper |

FluxesFigure 3 .Figure 5 .
Figure 3. Ratio of the SG model over ISSM total ice volume for six different regions in the Rockies as a function of hypsometric levels.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 ISSM and the SG model (using 10 hypsometric levels) are presented in Fig.S3of the Supplement.

Figure 8 .
Figure 8.Average (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.Total ice volume evolution for a simulation using a better fitting (relative to calibration constraints) parameter vector.Different curves represent simulation where the SG model and the flux code are used ("flux on") or not ("flux off"), and if the flux between coarse grid cells is turned off ("NofluxOut").

Figure 12 .Figure 13 .
Figure 12.Ice volume evolution for a simulation over North America during inception with the SG model turned on.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.

Table 1 .
Differences between our new SG model and Marshall and Clarke (1999)/Marshall et al. (2011) models.Done as soon as the SG model is activated Done only when the lowest H CG = Volume of lowest SG levels SG level is filled (total volume when SG is turned off) CG to SG: H SG = Average between equal ice thickness and lowest levels filled