Introduction
In order to make predictions about the future behaviour of marine ice sheets, ice-sheet models need to adequately represent grounding-line motion. Fixed-grid models have been demonstrated to give inconsistent results when spatial resolution is varied (Reference Vieli and PayneVieli and Payne, 2005). Recent studies have shown that modelled grounding-line behaviour is convergent, but very high resolution is needed to achieve convergence (Reference Durand, Gagliardini, de Fleurian, Zwinger and Le MeurDurand and others, 2009; Reference Gladstone, Payne and CornfordGladstone and others, 2010b). Attempts to parameterize the grounding line have met with varying degrees of success (Reference Pattyn, Huyghe, De Brabander and De SmedtPattyn and others, 2006; Reference Pollard and DeContoPollard and DeConto, 2009; Gladstone and others, 2010b). However, grounding-line modelling typically has been tested with idealized simulations from, or similar to, the Marine Ice-Sheet Model Intercomparison Project (MISMIP; full project information can be found at http://homepages.ulb.ac.be/fpattyn/mismip/), which impose high basal drag and no ice-shelf buttressing and may represent scenarios that are more demanding than many real marine ice-sheet outlet glaciers of interest.
In the current study the self-consistency of modelled steady-state grounding-line positions is investigated over a range of idealized scenarios with parameter values relevant to a variety of real-world settings. A flowline ice-sheet model, described below, is used.
In the first set of experiments we carry out perturbed parameter ensembles of model simulations for each of the 24 different grounding-line parameterizations (GLPs) of Reference Gladstone, Payne and CornfordGladstone and others (2010b) at a fixed resolution (by GLP we refer to a special treatment of the gridcell containing the grounding line, such as modification of the basal drag coefficient in that gridcell). The rate factor, net accumulation, basal drag coefficient and bedrock slope are all varied. The aims of the ensembles are to investigate how grounding-line errors vary in response to varying parameters in response to choice of GLP.
Motivated by the results of the perturbed parameter ensembles, a second set of experiments is carried out in which the size of grounding-line errors is quantified as a function of resolution for varying basal drag coefficient and channel width (here a lateral drag parameterization for flowline models is used to represent buttressing). The aim of these experiments is to quantify the resolution required to simulate grounding-line positions to a specified accuracy (in this case 1 km) over a range of scenarios that includes parameter values appropriate to marine ice-sheet outlet glaciers featuring ice-shelf buttressing.
Finally, we discuss the relevance of these simulations to application of ice-sheet models to marine ice sheets such as Pine Island Glacier (PIG), West Antarctica.
Method: Model and Experiments
All the simulations presented here are carried out using the fixed-grid ice-stream ice-shelf (FGSTSF) model of Reference Gladstone, Lee, Vieli and PayneGladstone and others (2010a). This is identical to the FGSTSF model of Reference Vieli and PayneVieli and Payne (2005) except that the higher-order piecewise parabolic method (PPM) is used for thickness evolution (for a description of the PPM method, see Reference Gladstone, Lee, Vieli and PayneGladstone and others, 2010a). It is a vertically integrated (vertical shear is not represented) flowline model using the ‘shelfy-stream’ equation.
In addition, the second set of experiments (in which resolution is varied) is carried out using a parameterization for buttressing in which the ice sheet is assumed to flow within a channel of constant width W with a ‘no slip’ condition at the side walls along the length of the model (Reference Van der Veen and WhillansVan der Veen and Whillans, 1996; Reference Vieli and PayneVieli and Payne, 2003).
The force balance is given by
where H is the ice thickness, x is distance from the ice divide, ѵ is the effective viscosity, u is the depth- and width-averaged velocity, A is the temperature-dependent rate factor, β 2 is a basal drag coefficient, g is the acceleration due to gravity, ρ is the ice density and s is the surface height. The second term on the left-hand side is removed for the first set of experiments. ѵ is given by
The flow-law parameter, n, is set to 3 for all simulations in the current study.
The left-hand boundary of the domain represents the ice divide and has a zero-velocity boundary condition. The right-hand boundary represents the calving front of the floating ice shelf, and a force-balance boundary condition is used. See Reference Gladstone, Lee, Vieli and PayneGladstone and others (2010a) for a full model description.
Grounding-line parameterizations
24 GLPs are used, as presented by Reference Gladstone, Payne and CornfordGladstone and others (2010b). They are based on choice of an interpolation function for thickness across the gridcell containing the grounding line, which is used along with the flotation condition to determine grounding-line position at subgridscale precision. The six thickness interpolation schemes are summarized in Table 1.
Table 1. Summary of interpolation schemes in the GLPs used in this study
A choice is then made as to how the thickness function is used to apply a correction to the forcing terms: basal drag and gravitational driving stress. These forcing correction schemes are summarized in Table 2.
Table 2. Summary of forcing modifications in the GLPs used in this study
Linear basal drag scaling refers to use of the interpolated grounding-line position to scale the basal drag coefficient according to the proportion of the cell that is grounded.
‘Profile scaling’ refers to use of the full thickness function rather than simply grounding-line position in calculating a correction to the forcing terms. See Reference Gladstone, Payne and CornfordGladstone and others (2010b) for a full description of the GLPs.
GLPs are referred to by combining the abbreviations in Tables 1 and 2 (e.g. LI_B1 uses linear interpolation for thickness across the grounding line and linear scaling of the basal drag term in the gridcell containing the grounding line).
Experimental set-up
The experimental set-up is similar to MISMIP experiments 1 and 2, but with greater variation of parameters. In all experiments the bedrock is downsloping towards the ocean and linear. The net accumulation rate is uniform at 0.3 ma–1 except where stated otherwise. All simulations were spun up from a uniform ice sheet of 100m thickness.
Given that fixed-grid models erroneously exhibit a range of steady-state grounding-line positions (Reference Gladstone, Lee, Vieli and PayneGladstone and others, 2010a), we have carried out pairs of advance and retreat simulations. Advance simulations comprise spin-up to steady state with constant forcing. Given the uniform 100m initial thickness condition, this results in advance of the grounding line towards steady state in all simulations presented in the current study. The retreat simulation corresponding to any given advance simulation comprises spin-up close to steady state under enhanced forcing followed by a forcing reset and adjustment to a new steady state. This is intended to cause the grounding line to retreat towards its final steady state, which occurs in all simulations in the current study except where stated otherwise. The enhanced forcing differs for the two sets of experiments and is described at the start of each section. For each set of experiments, different methods for applying enhanced forcing were tested (modifying rate factor, accumulation and combinations of both), with the outcome that the method of forcing modification has no impact on the final grounding-line position so long as the forcing enhancement is sufficient that significant retreat occurs after the forcing reset. There are strong theoretical grounds (Reference SchoofSchoof, 2007) to believe that the advance and retreat simulations should reach the same steady-state grounding-line positions.
The above describes the set-up common to both sets of experiments in the present study. Details specific to each set of experiments are described at the start of the relevant section.
Perturbed Parameter Ensembles
Reference Gladstone, Payne and CornfordGladstone and others (2010b) demonstrated modest sensitivity of grounding-line errors to choice of GLP within a MISMIP-like set-up. Here we explore sensitivity of grounding-line errors to different parameter values and different GLP choices, including combined parameter and GLP choices. This allows a more general reassessment of the relative performance of GLPs and a first assessment of the implications of a non-MISMIP-like set-up for grounding-line errors.
Errors in steady-state grounding-line position, and the relative performance of the different grounding-line parameterizations (GLPs) presented by Reference Gladstone, Payne and CornfordGladstone and others (2010b), are now investigated over a range of values for four different inputs. The inputs comprise two internal model parameters (rate factor A and drag coefficient β 2), a forcing term (net accumulation a) and a domain set-up parameter (bedrock gradient b’). A and β 2 are varied widely between barely plausible upper and lower limits (Table 3). The bedrock gradient and net surface accumulation are varied not to extremes but rather to ensure a spread of steady-state grounding-line positions that lie comfortably within the model domain (i.e. within limits that allow simulations to run successfully to completion). These four inputs and their upper and lower limits (Table 3) define the input space for the ensemble experiments. Note that the range of values for A corresponds approximately to a temperature range of zero down to –258C.
Table 3. Parameter value ranges used in the ensemble experiments
The domain size is 2112 km from ice divide (left boundary of domain) to ice front (right boundary of domain). The gridpoint spacing Δx is 2.4 km and the time-step Δt is 0.2 years. The implementation of advance and retreat simulations is identical to that described in Reference Gladstone, Payne and CornfordGladstone and others (2010b, appendix A), with the forcing enhancement for the retreat experiments based on accumulation and rate factor, the only difference being that the total run length is greater in the current study. The analyses are based on steady-state grounding-line positions achieved through very long simulations. All advance experiments were run for 100 ka, and all retreat experiments were run for 200 ka (100 ka initial phase, 10 ka of gradual forcing reset and 90 ka to final steady state). Long run times were needed to ensure steady state for the whole range of inputs. In particular it takes longer to reach steady state when the basal drag is high and when bedrock slope is shallow.
The bedrock profile, linearly downsloping in all experiments, is given by
where b is the bedrock height relative to sea level, b’ is the magnitude of the prescribed bedrock gradient and all distances are in m. The form of Eqn (3) is chosen so as to prevent the steady-state grounding-line position leaving the domain in the case of the shallower bedrock profiles.
Ensembles of simulations are carried out to investigate sensitivity to the inputs and to compare the GLPs. Two different sampling techniques are used: One at a Time (OAT) and Latin Hypercube Sampling (LHS). These sampling techniques, described in the Appendix, are used to generate an OAT ensemble and an LHS ensemble, both of which are run for all 24 GLPs.
Errors are assessed here by the metric ‘retreat minus advance’ (RMA; Reference Gladstone, Payne and CornfordGladstone and others, 2010b). RMA is a quantification of the size of the region of locally stable grounding-line positions (Reference Gladstone, Lee, Vieli and PayneGladstone and others, 2010a) and is defined as
where x gr is the final grounding-line position from a retreat experiment and x ga is the final grounding-line position from the corresponding advance experiment. It follows from the theoretical work of Reference SchoofSchoof (2007) that RMA should be zero in these experiments, and it has been shown (Reference Gladstone, Lee, Vieli and PayneGladstone and others, 2010b) that RMA approaches zero as resolution increases for the ice-sheet model presented here. It is worth noting that x gr≥x ga for all simulations in the current study.
RMA for all 24 GLPs from the OAT experiments is shown in Figure 1. RMA shows some dependence on rate factor and bedrock slope, but a much greater dependence on basal drag coefficient. Higher basal drag, shallower bedrock and higher rate factors (i.e. warmer ice) all increase grounding-line errors. The GLPs that incorporate profile scaling in the correction to the basal drag term in the gridcell containing the grounding line give generally smaller errors, but the spread of errors due to choice of GLP is small compared with the spread of errors due to varying the basal drag coefficient. Accumulation has little discernible impact on grounding-line errors. There is a 10–20km variation in grounding-line errors associated with varying rate factor, comparable to the spread associated with choosing different GLPs. The variations in bedrock slope cause a larger variation in grounding-line errors, ranging from a few km for steeper slopes up to 20–30km for shallow slopes. But increasing basal drag coefficient causes errors from much less than 1 km for low basal drag (β 2≤3.2×109 Pa sm–1) up to 100 km for high basal drag.
Fig. 1. The RMA metric from the OAT simulations shown against (a) rate factor, (b) basal drag coefficient (β 2), (c) bedrock gradient and (d) accumulation. Line types separate the different forcing corrections that constitute the GLPs (see Table 2; Reference Gladstone, Lee, Vieli and PayneGladstone and others, 2010b): black lines indicate the default parameterization (linear in basal drag) whereas grey lines indicate profile scaling of basal drag; solid lines indicate no correction to the gravitational driving stress term whereas dashed lines indicate profile scaling. Note that the tendency of the steady-state grounding line to lie at or near a gridpoint (Reference Gladstone, Lee, Vieli and PayneGladstone and others, 2010b) means that several of the 24 GLPs approximately overlie each other in this figure.
The OAT ensembles show how grounding-line errors and GLP performance respond to varying one input at a time, but they cannot capture potentially nonlinear responses to varying more than one input. The LHS ensembles allow investigation of grounding-line errors and relative GLP performance over the full input space, varying all inputs together.
The RMA metric for the simplest GLP, LI_B1 (see Tables 1 and 2), is shown in Figure 2 for the full LHS ensemble. This shows similar patterns to the OAT simulations when plotted against each input separately, but with some scatter as the combined impacts of all inputs are now shown.
Fig. 2. The RMA metric for the linear interpolation GLP (LI_B1; Reference Gladstone, Payne and CornfordGladstone and others, 2010b) from the LHS simulations shown against (a) rate factor, (b) basal drag coefficient (β 2), (c) bedrock gradient and (d) accumulation. Crosses mark individual simulations from the LHS sample.
Basal drag is the only input for which the RMA values from the OAT simulations follow the same trend and are consistently of the same order of magnitude as those for the LHS simulations. This indicates that basal drag dominates not only the other inputs individually but also combinations of inputs (i.e. there is no strong nonlinear response in terms of grounding-line errors to varying multiple parameters).
Although it can be seen that errors increase with shallowing bedrock slope, the trend for increasing errors with increasing rate factor, seen in the OAT results (Fig. 1), cannot be seen in Figure 2. This does not mean that such a trend is not present, but if present it is hidden at least to visual inspection by the scatter induced by varying β 2.
The LHS ensembles for the other GLPs (not shown) all show similar behaviour to LI, though the size of the errors varies a small amount. The GLP H2_GB2 shows slightly less spread in steady-state grounding-line positions and also generally smaller values for the RMA metrics, consistent with the single input comparison shown by Reference Gladstone, Payne and CornfordGladstone and others (2010b).
The relative performance of the different GLPs measured by the RMA metric over the LHS ensemble is summarized in Figure 3 in which the GLPs are ranked according to RMA for each member of the ensemble. There is a clear difference between the orderings for low basal drag (in which the GLPs with profile-based forcing corrections perform relatively poorly) and the orderings for medium to high basal drag (in which the GLPs with profile-based forcing corrections do well). From a practical viewpoint, given that errors are much less than a gridcell at low basal drag values for all GLPs, the ordering of GLPs is not important for low basal drag values. Hence the GLPs involving profile-based forcing term adjustments can be said to offer better overall performance.
Fig. 3. Relative GLP performance over input space. The unlabelled x- and y-axes are basal drag coefficient and rate factor, respectively, both on log scales identical to those of the x-axes of Figure 2b and a, respectively (i.e. rate factor increases logarithmically from 9.8×10–26 Pa–3 s–1 to 5.4×10–24 Pa–3 s–1, and basal drag coefficient increases logarithmically from 3.156×108 Pa sm–1 to 3.156× 1011 Pa sm–1). The RMA metric is used to rank the simulations over GLPs for each member of the LHS ensemble. The rank is shown linearly in greyscale from black (best-performing GLP for the given input) to white (worst-performing GLP). The GLP labels are those of Reference Gladstone, Payne and CornfordGladstone and others (2010b) (see also Tables 1 and 2). The GLPs are arranged in approximate order of increasing complexity of thickness interpolation function from top to bottom, and in order of increasing complexity of forcing parameterization (where profile-based correction of forcing terms is considered more complex than linear or no correction) from left to right.
In terms of the different thickness profiles, it is noteworthy that linear interpolation (LI) performs relatively poorly, along with the higher-order thickness profiles, linear extrapolation (LE) and cubic interpolation (CI), which use upstream and downstream gradients as well as thickness values. The second-order harmonic-mean-based interpolation (H2) and the interpolation of Reference Pattyn, Huyghe, De Brabander and De SmedtPattyn and others (2006) offer the best performance.
Inputs other than β 2 were not found to be important in ranking the GLPs. The relationship between ranking and rate factor is far weaker than that between ranking and basal drag (Fig. 3). Similarly to rate factor, bed gradient and accumulation do not have a strong impact on the ranking of GLPs (not shown). However, the input space used here is not a superset of all real-world inputs. Rate factor and drag coefficient are not expected to vary beyond the range of inputs in the current study, but bedrock gradient and accumulation might. The sensitivity of errors to accumulation appears negligible, but this is not the case with bedrock gradient. Errors increase as bed gradient shallows for all GLPs in the current study. Also, the impact of bedrock features of shorter wavelength than the domain has not been investigated in this study.
In summary, the recommendations of Reference Gladstone, Payne and CornfordGladstone and others (2010b) are confirmed in this exploration of grounding-line errors over input space for multiple GLPs: while use of a GLP is important (Reference Gladstone, Lee, Vieli and PayneGladstone and others 2010a), choice of which GLP is of secondary importance (at least out of the GLPs presented here). In particular, the impact of basal drag coefficient on grounding-line errors dominates choice of GLP and also dominates other inputs and combinations of inputs.
Grounding-Line Errors and Resolution
Reference Gladstone, Payne and CornfordGladstone and others (2010b) demonstrated convergence of steady-state grounding-line position for the different GLPs given fixed inputs. In the previous section, grounding-line errors at 2.4 km resolution were shown to be highly sensitive to basal drag. We postulate that the high grounding-line errors for high basal drag are due to the step change in basal drag across the grounding line from nonzero to zero. In addition to this step change being problematic to represent in a conventional finite-difference-based fixed-grid model, there is also, at least in simulations without a GLP, a step change in the forcing regime when the grounding line advances or retreats by Δx due to the step change in area over which the basal drag force is applied. This forcing step change is mitigated by use of a moving grid model that tracks the grounding line (Reference Vieli and PayneVieli and Payne, 2005). Ideally a GLP would resolve this problem by allowing the grounding-line position to move smoothly through the gridcell, but in practice the GLPs used here show step-like behaviour (Reference Gladstone, Payne and CornfordGladstone and others, 2010b), hence grounding-line error is sensitive to basal drag.
Having established that grounding-line errors increase as basal drag increases, we now speculate that grounding errors decrease as ice-shelf buttressing or ice-stream lateral drag increases. The reason for this is that back-stress from ice-shelf buttressing, due to the lateral drag from non-slip side walls along the length of the floating ice shelf, imposes a force in the upstream direction that diminishes as basal drag increases across the grounding line, effectively reducing the step change in forcing across the grounding line.
In this section, grounding-line errors are quantified as a function of resolution for two sets of simulations. In the first set of simulations, 11 different values of β 2 are used, chosen on a log scale from 7.2×108 Pa sm–1 (i.e. the MISMIP value ×10–2) to 7.2×1011 Pa sm–1 (i.e. the MISMIP value ×10). Buttressing is not included in these simulations. In the second set of simulations a lateral drag parameterization is used (Reference Van der Veen and WhillansVan der Veen and Whillans, 1996; Reference Vieli and PayneVieli and Payne, 2003). Seven experiments are carried out in which the basal drag coefficient β 2 is fixed at 1010 Pa sm–1 and the parameterized channel width W takes the values 100, 200, 400, 800, 1600 and 3200 km and infinity (no lateral drag). The simulations are run at six different resolutions from 0.3 to 9.6 km, with a factor two change between each resolution. Note that 100 km is already wide for an ice stream. The current study aims to investigate the range of behaviours from no buttressing to significant buttressing. 100 km was chosen as the narrow end of the range of values for W because a much narrower channel width leads to grounding-line migration beyond the model domain for the current idealized set-up.
The domain is as in MISMIP experiments 1 and 2, but extended in the seawards direction to be 2400 km long. The bedrock height relative to sea level b(x) is given (in m) by
Rate factor is 2.15×10–25 Pa–3 s–1 (as in step 5 of MISMIP experiment 1). The net accumulation is 0.3 ma–1, enhanced to 0.5 ma–1 for the initial phase of the retreat simulations. The retreat simulations have an initial phase of 50 ka and a total run length of 100 ka. The advance simulations have a run length of 70 ka.
The linear interpolation GLP LI_B1 is used. This GLP is chosen because it is the simplest to implement, even in a two-dimensional ice-sheet model, and in terms of performance there is not a large difference between GLPs.
Steady state was established through visual inspection of grounding-line evolution plots (not shown), and all simulations reached steady state except for the four highest basal drag coefficient values (i.e. values >4.5×1010 Pa sm–1). Simulations using these four values are omitted from the rest of this section. This does not affect our analysis or conclusions since the high basal drag values are less relevant to our scenarios of interest, marine ice-sheet outlet glaciers, and because the ‘shelfy-stream’ model used here becomes less justifiable for high basal drag values when the neglected vertical shear becomes significant.
The steady-state grounding-line positions are shown against resolution in Figure 4 both for experiments with varying basal drag (Fig. 4a) and varying parameterized channel width (Fig. 4b). As expected, the advance and retreat simulations appear to converge towards the same steady-state grounding-line position for all variations of basal drag and parameterized channel width presented here. Both advance and retreat simulations appear to overestimate the grounding-line position. The cause of this is not fully understood, though it is hypothesized to result from the GLP. Reference Gladstone, Lee, Vieli and PayneGladstone and others (2010a) found that, without a GLP, advance simulations tended to massively underestimate the steady-state grounding-line position, whereas using the same GLP as the current study gave advance simulations that overestimate the grounding-line position.
Fig. 4. Steady-state grounding-line positions against resolution are shown in greyscale for (a)W=1(no lateral drag) and basal drag coefficient values from β 2 = 7.2×108 Pa sm–1 (black lines) through 1.4×109, 2.8×109, 5.7×109, 1.1×1010, 2.3×1010 and 4.5×1010 (light grey lines); and (b) β 2=1010 Pa sm–1 and parameterized channel width values from W= 100 km (black lines) through 200 km, 400 km, 800 km, 1600 km, 3200 km and infinity (light grey lines). Both advance (solid lines) and retreat (dashed lines) simulations are shown.
The previous section showed that grounding-line errors are higher for higher basal drag values. Figure 4a confirms that this finding is consistent over a range of resolutions and also demonstrates that convergence is poorer at higher basal drag values.
As buttressing is increased by reducing the parameterized channel width (Fig. 4b), the steady-state grounding-line position moves seawards, errors decrease and convergence with resolution improves.
In order to quantify the relationship between grounding-line errors and resolution at a finer scale than can be ascertained from Figure 4, we define and plot an error estimate ε. It is assumed (Reference Gladstone, Lee, Vieli and PayneGladstone and others, 2010a) that the most accurate solution for a given set of inputs is given by the highest resolution the model has been run at. The steady-state grounding-line error estimate ε gs is therefore defined as
where 
is the steady-state grounding-line position from a simulation at resolution Δx, and 
is the steady-state grounding-line position from a simulation run at maximum resolution (i.e. Δx = 300m in the current study). We set a target error of ε = 1km as being an acceptable error for grounding-line position in real-world applications to marine ice-sheet outlet glaciers, where grounding-line movement of tens of km can occur on a timescale of decades (Reference Joughin, Smith and HollandJoughin and others, 2010).
The error estimate ε is shown against resolution in Figure 5. Plotting ε instead of actual grounding-line position allows a log scale to be used. Convergence is apparent for all basal drag values (Fig. 5a) and for all parameterized channel widths (Fig. 5b) down to the target error of 1 km. With no buttressing, the target error is only achievable with basal drag coefficient <2.8×109 Pa sm–1 and at resolution finer than Δx = 1.2 km (Fig. 5a). With 100 km parameterized channel width and basal drag coefficient 1010 Pa sm–1, the target error is achievable at resolution Δx = 2.4km (Fig. 5b).
Fig. 5. The error estimate ε against resolution is shown in greyscale for (a) W=∞ (no lateral drag) and basal drag coefficient values from β 2 = 7.2×108 Pa sm–1 (black lines) through 1.4×109, 2.8×109, 5.7×109, 1.1×1010, 2.3×1010 and 4.5×1010 Pa sm–1 (light grey lines); and (b) β 2=1010 Pa sm–1 and parameterized channel width values from W= 100 km (black lines) through 200 km, 400 km, 800 km, 1600 km, 3200 km and infinity (light grey lines). Both advance (solid lines) and retreat (dashed lines) simulations are shown.
These experiments support the hypothesis that both increasing buttressing and decreasing basal drag reduce grounding-line errors.
Conclusions
We have assessed model performance for a variety of GLPs over a region of input space and specifically for varying basal drag and parameterized channel width over a range of resolutions, focusing on the known problem of poor convergence of grounding-line behaviour. The metrics used, RMA and ε, are measures of self-consistency of steady-state grounding-line behaviour and should be viewed as a model verification step rather than validation, which can be done only through comparison against observational data. The model and GLPs presented here have already been verified against a semi-analytic benchmark (Reference SchoofSchoof, 2007; Reference Gladstone, Payne and CornfordGladstone and others, 2010b) for a fixed set of inputs.
Reference Gladstone, Payne and CornfordGladstone and others (2010b) showed that the GLPs with profile-based forcing corrections performed well for a fixed set of inputs. Our findings support the ranking of GLPs found by Reference Gladstone, Payne and CornfordGladstone and others (2010b) over a wider input space, except for very low values of basal drag, in which case the errors for all GLPs are sufficiently small not to be of concern. We find that the spread in grounding-line errors due to choice of GLP is small compared with the spread induced by varying inputs, especially basal drag coefficient. The implementation of a GLP is essential for the model presented here, but the choice of GLP is secondary to the model setup. In other words, the scenario the model is to be used for has a much larger impact on resolution requirements than does GLP choice.
Our results show that numerical errors related to grounding-line behaviour are smaller in certain situations relevant to real-world outlet glaciers than in MISMIP experiments. In other words, MISMIP is too harsh a test for some real-world applications. More specifically, PIG in West Antarctica, predicted to cause several cm of sea-level rise over the next century (Reference Joughin, Smith and HollandJoughin and others, 2010), is relevant to the current study. PIG flows into an embayment in which buttressing occurs, and the basal drag coefficient under the ice stream and near the grounding line is less than the value of 1010 Pa sm–1 used in the buttressing experiments in the current study (Reference Vieli and PayneVieli and Payne, 2003). Our results suggest that PIG could be modelled adequately (from the perspective of self-consistent grounding-line behaviour) with a resolution of 1 km or coarser, a resolution that current ice-sheet models are already capable of achieving, at least for individual outlet glaciers (e.g. Reference Joughin, Smith and HollandJoughin and others, 2010).
However, the current study omits investigation of the impact of bedrock undulations of wavelength less than the model domain, which are thought to cause ice plains and temporary stability (Reference Joughin, Smith and HollandJoughin and others, 2010), and it omits ice-shelf melting, which is high especially for PIG (Reference Jenkins, Vaughan, Jacobs, Hellmer and KeysJenkins and others, 1997). A recent modelling study of PIG achieved a resolution of Δx = 160m in the vicinity of the grounding line, yet the solution still exhibited strong resolution dependence (Reference Joughin, Smith and HollandJoughin and others, 2010). We speculate that the ice-shelf basal melt parameterization of Reference Joughin, Smith and HollandJoughin and others (2010), which imposes shelf melt rates of over 100ma–1 next to the grounding line, adds to the convergence problems for grounding-line modelling.
The current study is based on experiments using a one-dimensional ‘shelfy-stream’ model (Reference Gladstone, Lee, Vieli and PayneGladstone and others, 2010a) with a linear basal drag law. The dynamical response of the system might vary for different model formulations, though qualitatively other fixed-grid (and, to a lesser extent, adaptive-grid) models are also challenged by convergence of grounding-line behaviour (Reference Durand, Gagliardini, de Fleurian, Zwinger and Le MeurDurand and others, 2009). We recommend that convergence of grounding-line behaviour needs to be demonstrated by any model being used to make statements or predictions about marine ice sheets and that the convergence experiments should be directly relevant in terms of model inputs/parameters to the physical regime being studied. In particular, shelf melt and bedrock undulations may affect convergence of modelled grounding-line behaviour and these effects have not yet been studied.
Acknowledgements
We thank the two anonymous reviewers for helpful comments which led to improvements in the clarity of the manuscript. This work was supported by the UK’s National Centre for Earth Observation cryosphere and climate themes, the UK Natural Environment Research Council (NERC) Joint Weather and Climate Research Programme, and the project ‘Understanding contemporary change in the West Antarctic ice sheet’ funded by NERC standard grants NE/E006256/1 and NE/E006108/1.
Appendix: Sampling Input Space
Sampling is the process of selecting a number of points within input space, i.e. selecting a number of parameter combinations to use in a simulation. The current study uses two sampling methods: One At a Time (OAT) sampling and Latin Hypercube Sampling (LHS).
The OAT ensembles are intended to provide a limited but easily interpretable exploration of input space. A default value is chosen for each of the four inputs, and one input is varied at a time while the other inputs are held at default value.
Given that some of the inputs vary over several orders of magnitude and that we are interested in learning about behaviour at all orders of magnitude, a sampling distribution that is uniform over the inputs would provide a sample that is unjustifiably sparse in the lower orders of magnitude. Hence we define a distribution for sampling that is uniform over the logs of the inputs. If P is an input then its default value P d is given by the log centre of the range of allowed values (see minimum and maximum values, Table 3), given by
where P max and P min are the upper and lower limits of the parameter, respectively. Each input is varied between its minimum and maximum values while the other inputs remain at their default values. Nine values are chosen for each input, spanning the range from minimum to maximum linearly in the log of P. Specifically, if i is an integer from 1 to 9,
where Pi is the ith value for P. This gives a total of 4 (inputs)×9 (values) = 36 members of the ensemble (i.e. 36 simulations to be carried out).
While OAT sampling gives easily interpretable outputs, it does not allow for quantification of possible nonlinear interactions between inputs. The LHS (e.g. Reference Santner, Williams and NotzSantner and others, 2003, p. 127) ensemble is intended to provide a more thorough exploration of input space. A full description of LHS is not given in the current study, but note that LHS provides a more space-filling design than random sampling while still preserving the underlying distribution. In practice this means that, for a sufficiently large sample, LHS allows all combinations of parameter variations to be studied, including possible nonlinear response to parameter combinations that would not be detectable using OAT sampling. A sample size of 100 has been used for the LHS.
