A glacial systems model configured for large ensemble analysis of Antarctic deglaciation

This article describes the Memorial University of Newfoundland/Penn State University (MUN/PSU) glacial systems model (GSM) that has been developed specifically for large-ensemble data-constrained analysis of past Antarctic Ice Sheet evolution. Our approach emphasizes the introduction of a large set of model parameters to explicitly account for the uncertainties inherent in the modelling of such a complex system. At the core of the GSM is a 3-D thermo-mechanically coupled ice sheet model that solves both the shallow ice and shallow shelf approximations. This enables the different stress regimes of ice sheet, ice shelves, and ice streams to be represented. The grounding line is modelled through an analytical sub-grid flux parameterization. To this dynamical core the following have been added: a heavily parameterized basal drag component; a visco-elastic isostatic adjustment solver; a diverse set of climate forcings (to remove any reliance on any single method); tidewater and ice shelf calving functionality; and a new physically motivated, empirically-derived sub-iceshelf melt (SSM) component. To assess the accuracy of the latter, we compare predicted SSM values against a compilation of published observations. Within parametric and observational uncertainties, computed SSM for the present-day ice sheet is in accord with observations for all but the Filchner ice shelf. The GSM has 31 ensemble parameters that are varied to account (in part) for the uncertainty in the ice physics, the climate forcing, and the ice–ocean interaction. We document the parameters and parametric sensitivity of the model to motivate the choice of ensemble parameters in a quest to approximately bound reality (within the limits of 31 parameters).

The article describes the Memorial University of Newfoundland/Penn State University (MUN/PSU) glacial systems model (GSM) that has been developed specifically for large-ensemble data-constrained analysis of past Antarctic Ice Sheet evolution.We document the parameters and parametric sensitivity of the model to motivate the choice of ensemble parameters in a quest to approximately bound reality (within the limits of C1144 31 parameters).
One other key sentence in the abstract is the following: "Our approach emphasizes the introduction of a large set of model parameters to explicitly account for the uncertainties inherent in the modelling of such a complex system."This along with the novel subshelf melt(SSM) model was the reason that we chose to submit it to TC/TCD and not GMD.We believe this submission provides a good example of model configuration that really considers the underlying parametric uncertainties in paleo ice sheet models.
In other words, a key question we are trying to answer with this paper is: "does this model configuration (ie set of ensemble parameters) probe the relevant glacial system uncertainties to a major/significant extent" (and implicitly to a much greater extent than previous studies)?
Large ensemble analyses are the focus of our next submission.This was the best division we could come up with.Our next submission is at least 54 pages long and combining the two into one paper would result in a manuscript that few would have the patience to read (or review).
Other papers in TC/TCD that introduce models show a few to moderate number of results for model validation.The point to consider here is that this model largely includes/subsumes the already well documented and validated PSU model (in the sense that the appropriate choice of ensemble parameters can in large part bring the model within close configuration of the PSU model).So traditional verification of the whole model is not needed to the same extent.The most novel component is the SSM model, and that is explicitly validated.For our context, full validation entails the ensemble analysis in our next submission.
List of specific reviewer comments and our responses: The names of parameters make the reading difficult.I understand very well that C1145 these names may be used in the code but more readable or even pronounceable (try zcrhslid) names would help.Moreover, I doubt that many people are able to remember the names of 31 ensemble parameters, more than 50 other variables, and 6 metrics.
In the text, please recall a short definition as often as necessary to make the reading fluent.
This MUN/PSU GSM is derived from what is becoming a community model (the PSU).This version is intended to eventually become a community model as well.A published model description with variable names that do not correspond to those in the model has little value to those who will most need it.Furthermore, having reviewed the variable names, we find few that we would want to change in hind-sight.Rational variable name choice for one person will often tend to be idiosyncratic to another.For instance, most ice sheet modellers use "H" for ice thickness, and "h" for elevation.Only a French speaker would understand the origin of the the latter.Given this, the standard and we believe near-optimal solution in the literature is a clear tabulation of variables as is present in this submission.In the future, electronic publishing will hopefully enable clickable link descriptions for all variables/acronyms.
We have also gone through the text with appropriate edits to ensure that variable usage far from initial introduction/context has its definition appropriately recalled.We have also now added a table of all acronyms.
## It would be useful to have a typographic way (name, font, Bold ::: ) to make the difference between the 31 ensemble parameters and the characteristics that are derived from these parameters (Slk, Se, ::: ).In particular, when these characteristics are maps, please mention it (Se(X) for instance).It is done in the climate forcing part, but not in the basal drag one.Same comment in table 2 Distinguishing ensemble parameters in bold is a good idea and that is now implemented in the revised text.Incorporation of explicit spatial and temporal dependence would add a lot of superfluous clutter since only 2 out of 20 equations in our revised C1146 submission lack spatial dependence).The reason explicit temporal and spatial dependence are shown for precip and temperature was to clarify which fields provide spatial and which fields provided temporal dependence (this is much more obvious for other fields).
## To better illustrate the dependency of variables with parameters, an example would be a figure with fstd as a function of Se for various roughness.This is only an example, the same type of plot should be given for most parameterisations.
The relationships are clearly specified by the given equations.To graphically detail all or most of the parametric dependencies in the model would be a major undertaking, that few would have the patience to read (who would want to go through 50 or more extra plots?) and would also be problematic given various interdependencies.
We have provided the sensitivity table to summarize the response of the model to the ensemble parameters.Previous model descriptions by others (eg The Cryosphere, 5, 727-740, 2011) also do not provide such requested plots of introduced parametrizations.We have graphically detailed the SSM component as this we believe is quite novel and physically motivated.Perhaps the reviewer is also reading too much into the parametrizations.They are largely a synthesis of previous parametrizations with a focus on making explicit at least some of the underlying parametric sensitivities.Their validation stems from a combination of the documented source origin, physical plausibility/motivation, and resultant impact on characteristics/statistics of generated ice sheet chronologies.
## A map of basal drag at present time should be given because it would help comparison with other works.This can be done for the run that gives the smallest misfit with present topography, (or for the baseline run).For the same run, please give a map of velocity compared with observations.
A basal drag coefficient map is a good idea for demonstration purposes and we have added that to the supplement.As before, the velocity map is beyond the purpose of C1147 this submission and is included in our next submission.
## A table with the metrics would be useful.There are not that many types of observation over the ice sheet and I regret that the observed present velocities are not used as an additional metric.
The topic of metrics is fully addressed in a companion paper (Briggs and Tarasov, QSR 2013).To the extent the observed velocities match balance velocities, then velocities do not add constraint beyond the constraint of fit to present-day ice thickness (assuming that present-day model accumulation is reasonably accurate).As for the residual discrepancy (ie non-balance component), then there is little apriori basis to believe that a forward model should capture this residual without local tuning given the somewhat random nature of ice-stream switching.As such, we believe it's better to compare a posteriori velocities from a calibrated ensemble against observations than to calibrate or tune the model to match current velocities, especially given our uncertainty focus and glacial cycle context.
## The parametrisation of climate forcing is very clever.To illustrate the time evolution would it be possible to have the envelope of the climate record in ice core locations (for instance Dome C, Dome Fuji...) so it could be compared with ice records ?
This is addressed more appropriately in our next submission that examines ensemble results.As the temperature depends on surface elevation, an envelope can't be well defined for the present focus on model description and a small number of runs.Furthermore, as the parametrization spans previous approaches, there is no need to validate it.
## At the end a sensitivity study is presented, and here also it could be better illustrated.
For instance a plot of the envelope of the grounded volume evolution, with the baseline run.As a conclusion, the scientific content is interesting and useful but the presentation must be strongly improved.

C1148
Again, as detailed above, this is beyond the purpose of the present paper.The analysis of ensemble results is about to be submitted for publication.
# Detailed comments ## p. 1535-line 24.The phase space is larger but the dispersion between reconstructions might be reduced if the complexity comes from adding a physical process.This is usually not the case.Dispersion can be reduced with a more physically-based process representation of an already incorporated process that thereby reduces parametric uncertainty, but the addition of new physical processes tends to increase the degrees of freedom in the model.The possible exception is a process that added a clear negative feedback to the system (which dominated over any positive feedback added) with an appropriate level of parametric constraint.But none of us can come up with an example of such a process that would pertain to the current state of ice sheet (and associated systems) modelling.
## p. 1540.equation 1.The exponent of the flow law could have been a parameter as well .As Cuffey and Paterson mention, it could depend whether basal drag is governed by sticky spots and sediments properties.Can the authors justify their choice ?
In short no.We make no claim as to have fully captured all uncertainties, only to have considered uncertainties in this context to a much greater magnitude than any previous study.As we already stated in the submission, the present evidence leans towards a Coulomb plastic relation governing subglacial till deformation.Various relationships have also been proposed for hard-bedded sliding but as yet we lack the observations to decisively choose between them.The current lack of sub-glacial hydrology component in the model (soon to be rectified) adds further uncertainty in the parametric representation.Whether fully accounting for such uncertainties will have a significant effect would need further study and this point is now made explicitly clear in the text.In terms of the eventual Bayesian calibration of the model, a key challenge is to capture as many uncertainties as possible while minimizing redundancy.LT's calibration ap-

C1149
proach has much weaker convergence when multiple parameters are largely offsetting each other.So having both the coefficient and exponent of the basal drag law (and one would need separate ones for both soft and hard beds) as ensemble parameters could be a problem.
## p. 1542 line 12.Considering the non local aspect of both the SSA equation (even in its hybrid form) and mass conservation equation, the impact of crh on ice thickness is non local so the observation-model misfit might not be related to local crh.Can you comment on this point ?
We fully acknowledge that adjustments to basal sliding coefficients are local, whereas the effects on the dynamics and ice thickness are non-local (and non-linear).Nevertheless, Pollard and DeConto (2012, TC) show that good results are achieved with a similar "local" inverse method, i.e., the iterative procedure yields realistic ice elevations despite the non-local nature of the dynamics.
## p. 1542.equation 4. All these details are very useful if someone wants to reproduce the experiment or use this parameterization but it seems a bit "magical".Why for instance, zcrhsed/zcrhslid is at the power Se and not multiplicated by Se.
We have modified the text somewhat to better motivate this equation.and DeConto (2012, GMD).The focus of this submission is on detailing the differences from the PSU model.Other references provide more background but do not provide any more critical material that would be required by an expert modeller to reconstruct the full model.As such, there would be no purpose served in copying a previous model description (i.e., Pollard and DeConto, 2012, GMD) and doubling the size of this paper.
## -While I understand that it might be computationally too costly to generate large ensembles of 2000-3000 runs on a finer grid, I am concerned about the very coarse resolution of 40km.It is very probable that a change in resolution will have a much larger impact on the results than the choice of some of the 31 parameters discussed in the manuscript.The implications and limitations of using this low resolution need to be thoroughly discussed.
As the reviewer suggests, the 40 km grid was necessary to make the many long-term runs feasible for the project, and enabled a much more extensive large-ensemble suite of experiments than would otherwise have been possible.We feel that this is justified, and the 40-km results are robust, because previous tests with much the same model have shown very little sensitivity to resolution.Results at 40, 20 and 10 km resolutions are shown in Pollard and DeConto (2012, GMD, Fig. 6) and (2012, TC, App.C), with very little differences in modern and inverse ice distributions.The primary model feature that allows this insensitivity to resolution is the grounding-line flux parameterization of Schoof ( 2007), as discussed in the above 2 references.The last 2 sentences have been inserted into the text.
## -It seems as if some of the chosen metrics are more sensitive to the chosen parameters than others -could the authors comment on this observation and its implications?
This judgement depends on how you compare the values of the different metrics (eg rescaling the sensitivity plots could provide different visual impressions).The choice of scale we chose (for the plots) was the maximum range over the whole parameter ensemble.This choice negates any direct comparison of maximum range, but allows C1152 comparison of say number of parameters that cover say 20% or 50% of the max range.
For this context, the present-day grounded ice volumes have the least number of high impacting (say 30% range coverage) parameters.But that's not surprising given that all the climate parameters converge to provide the same present-day climate (aside from elevation effects).One could speculate about the impact of past evolution on presentday configurations, but it is much more direct and complete to carry out this analysis for a large ensemble (part of future work).
Good point.This is admittedly a qualitative statement, but we took it to mean at least one model output value changes by ∼20% or more relative to the maximum range of results over the whole sensitivity ensemble.We'll make the revised submission more explicit instead of the vague "significant".
## -I agree that the appropriate choice of the metrics depends on the scientific question (cf.page 1564, line 21).However, it remains unclear to me why the authors have chosen these 6 diagnostics?Why not include, e.g., present-day ice velocity as a metric?
The 6 diagnostics are basic and standard, and allow meaningful comparisons between model and observed fields.The difficulty with surface velocities are that they vary over several orders of magnitude, and have fine spatial structure around the margins (streaming valleys vs. slow ridges) where simple direct point-to-point comparisons would yield considerable scatter even if the pverall pattern is very good.This is seen in scatter plots of observed vs. model velocities (e.g., Pollard and DeConto, GMD, 2012, Fig 9d, and other publications).
## -How does the introduction of the grounding line flux condition affect the sensitivity of,e.g., the Ross and Filchner-Ronne shelf area?Does this affect the interpretation of the results with respect to metrics 4, 5 and 6?

C1153
It's hard to say how the grounding line flux condition affects the model, because it is a fundamental part of the model physics and is present in all model runs.Very early in model development (a long time prior to the project), before implementing the flux condition it was found that 2-D Antarctic grounding line positions and migration were modeled very poorly with grid resolutions more than a few km (as found in MISMIP model intercomparisons).So we have no "with or without" pairs of simulations to compare, and cannot do them because "without" produces very poor results.The metrics 4,5 and 6 emerge from the model with the flux condition as an intrinsic part of the physics.We can say that the model compares well in MISMIP and MISMIP3D idealized intercomparison projects, in equilibrated grounding-line movements.
## -page 1566, lines 21-23: It would be interesting to learn more about the quantitative score for the other parameters as well.does this imply for the relative importance of a specific parameter?MM, again this runs into the challenge of avoiding excessive overlap between a sequence of 3 papers.But yes, we can see this being of value and will add it to our revised submission.The question of relative importance of parameters is dependent on which metric one is using.That's why we find it important to document the sensitivity of a number of ice chronology statistics to the ensemble parameters.But it should be emphasized that a parameter that has high impact on say some ice sheet characteristic (eg area of warm basal ice at LGM) might have minimal impact on the final score if there are no observational constraints that feed the scoring scheme closely related to that characteristic.
## -page 1566, lines 8-11: Could the authors speculate as to why the WAIS volume is less sensitive to the shelf flow parameter than the EAIS volume?(It seems that the parameters concerning shelf dynamics and shelf-ocean interaction generally havemore impact on WAIS than on EAIS.) Looking at back at earlier sensitivity studies (around some other baselines before the C1154 large ensemble was complete) this is not the case.Furthermore, if one considers the grounded RMS fits to present-day (ALBMAP) ice thickness, then the shelf flow parameter has a much bigger impact on WAIS (317 to 352 m grounded, and 374 to 414 m floating) then EAIS (204 to 213 m grounded) (even for the included sensitivity set).
We will add a short discussion of the need to examine many metrics/characteristics to the revised submission.
## -Fig.2: A zoom of the Ross region would be helpful since the melt patterns cannot be seen at this scale.
There is only a small degree of melt located in a small spot at the grounding line and the shelf edge.There is practically zero freeze-up.This can be seen in the map.A zoom would not provide any further insight.
## -In addition to figures 7, 8 and 9, could the authors include maps illustrating for instance the lateral boundaries for the baseline run and for ensemble members particularly close to/far from present-day observations of the Antarctic Ice Sheet?
We do not understand what is meant by lateral boundaries.We repeat that ensemble analysis is the focus of a separate submission and beyond the purpose of this submission.
## Technical corrections ## -Since there are so many parameter and variable names, they should be as short and self-explanatory as possible.As an example, vol0gw could be renamed something like V_WAIS etc The example of V_WAIS fails to capture the fact that the vol0gw is for grounded (g), present day (0), volume (vol) of wais (w).However a bit more explanation might be useful.So we've expanded the description of this variable as follows: for present-day WAIS (vol0gw; where volume=vol, present day=0, grounded=g, and WAIS=w) and for EAIS (vol0ge; EAIS=e), total grounded ice volume for the LGM C1155 (vol20g; LGM=20), the zonal position of the Ross shelf grounding line (RISg; RIS=Ross Ice Shelf, gl=grounding line) along the 81 â Ů ę S line of latitude, and the shelf areas for ROS and for RON-FIL.
## -It would be helpful if the authors could repeat the variable names instead of just using their symbol/abbreviation in the main text.
We respectfully disagree here as this negates the whole point of using abbreviations.
As stated above, the standard solution is a table of variable definitions, and we already provide this.As per response to reviewer #1, we have gone through the text with appropriate edits to ensure that variable usage far from initial introduction/context has its definition appropriately recalled.
## -Abbreviations which are only introduced in the main text should be repeated in the figure captions (e.g., GLZ, ACZ,...).
We have now done this for the more obscure SSM, GLZ, ACZ, and SFZ (for the first figure where usage appears), but do not see the need to do this for the much more commonly used or contextually obvious abbreviations such as WAIS and RON-FIL.
Remaining grammatical (missing word...) corrections have all been implemented.
Interactive comment on The Cryosphere Discuss., 7, 1533, 2013. C1156 For the specific question above, power relationships are more appropriate for transition between values that differ by more than one order of magnitude (otherwise linear variation in the controlling parameter (in this case Se) will be too heavily weighted towards large values.
## p. 1544.line7This model participated in the MISMIP3D.Could you please state how it compared with the other models.Results with a similar model are very reasonable for MISMIP3D, generally well within the range of the other models.In a related sensitivity study, Drouet et al. (2013) found that the use of the Schoof (2007) grounding-line flux parameterization in this C1150 type of model changes the time-dependent response somewhat, but not hugely, and on relatively small spatial and short time scales compared to those here.Drouet, A.S., D. Docquier, G. Durand, R. Hindmarsh, F. Pattyn, O. Gagliardini and T.## -The model description should be as complete as possible.Right now, several papers are required in order to gather the necessary information for the parameterizations.For the scope of this article, it would really help if all parameters were (re-)introduced in the manuscript, even if that means repeating some equations from prior publications (e.g., for the enhancement values.)Thispaper summarizes key features of the PSU model that is fully detailed in Pollard C1151