Introduction

Reconstructing the extent and behaviour of palaeo-ice sheets requires a multidisciplinary strategy that integrates a range of approaches, knowledge and information. Palaeo- ice-sheet fluctuations and dynamics have played a central role in the mechanisms of global climate change and large- scale oceanic and atmospheric reorganization. Assessments of the spatial extent and timing of Pleistocene glaciations on a global scale (Reference Denton and HughesDenton and Hughes, 1981; Reference PeltierPeltier, 1994) and within North America and Europe (e.g. Reference PorterPorter, 1984; Reference Ehlers and GibbardEhlers and Gibbard, 2004a, Reference Ehlers and Gibbardb), are therefore of critical importance to a variety of Earth sciences. However, while many reconstructions of ice-sheet extent and timing have utilized either numerical modelling or geological information, few have been constrained by linking numerically derived data with landform and sedimentary evidence in an explicit quantitative framework (e.g. Reference Kleman, Fastook and StroevenKleman and others, 2002).

Attempts to reconstruct ice sheets and their dynamics typically fall into one of two general strategies, both of which have advantages, limitations and involve numerous assumptions: (1) geological approaches, which use a variety of geomorphological and geochronological data through a set of erosional, transport and depositional rules to reconstruct the morphology and basal characteristics of an ice sheet at particular key snapshots, often coincident with maximum extent or a significant readvance (e.g. Reference Boulton, Smith, Jones and NewsomeBoulton and others, 1985; Reference PeltierPeltier, 1994; Reference Kleman and StroevenKleman and others, 1997; Ó Reference Ó Cofaigh, Pudswey, Dowdeswell and MorrisCofaigh and others, 2002), and (2) numerical modelling, which, with its internally explicit logic, consistently combines mathematical formulations of physical theory (be it simplified and reduced) and boundary conditions to reproduce the thickness and dynamics of ice sheets (e.g. Boulton and Payne, 1992; Fastook and Holmlund, 1994; Reference Siegert, Dowdeswell, Hald and SvendsenSiegert and others, 2001; Reference Takeda, Cox and PayneTakeda and others, 2002). Empirical data from a variety of sources, such as dated raised beaches, moraines, marine sediments and ice cores, may then be used to define the areal footprint or climatic forcing for numerical models (e.g. Reference Bintanja, van de Wal and OerlemansBintanja and others, 2005), and, as a result, the changing ice-sheet geometry and flow may be partially constrained using observational data.

While geological models provide insight into patterns of ice-sheet extent, advance and retreat, they are limited by the continuity, quality and availability of geomorphological evidence, which provides significant challenges when reconstructing a coherent and continuous glacial cycle (Reference Kleman, Hättestrand, Stroeven, Jansson, De Angelis, Borgstrom and KnightKleman and others, 2006). Ice-core oxygen and deuterium isotope records can provide temporal proxies for ice volume (e.g. Reference Anandakrishnan, Alley and WaddingtonAnandakrishnan and others, 1993) and palaeoclimate patterns, but, with the exception of promising recent work by Reference Clarke, Lhomme and MarshallClarke and others (2005) and Reference Lhomme, Clarke and MarshallLhomme and others (2005), do not provide insight into the spatial extent of ice. However, by driving the climatic component of numerical models by palaeoclimatic proxies derived from ice cores, ice-sheet growth and decay through the entire glacial cycle can be simulated. Yet, numerical ice-sheet models are plagued by significant methodological and constraining- data shortfalls concerning adequate parameterization, critical process inclusion, verification and the necessary climatic forcing. Thus, there is a critical need to integrate detailed field data with numerical ice-sheet models to better understand the linkages between glacial landform patterns and sediments and glaciologically plausible scenarios for their development (e.g. Reference ClarkClark, 1997; Reference Ehlers and GibbardEhlers and Gibbard, 2003; Reference HagdornHagdorn, 2003).

Until recently, the use of geomorphological and geological evidence as a calibration and verification tool for model predictions has been largely qualitative, using phrases such as optimum, close resemblance, compatible and acceptable to describe the correspondence between model predictions and field data (e.g. Reference Marshall and ClarkeMarshall and others, 2002; Reference Zweck and HuybrechtsZweck and Huybrechts, 2005). There have been successful attempts to quantitatively assess model output against geomorphological evidence, but these are limited in scope. For example, Reference Tarasov and PeltierTarasov and Peltier (2004) assess their suite of North American numerical ice-sheet simulations using geophysical inversions of relative sea-level data from specific sites. The offsets between geophysical evidence and numerical simulations were calculated using the root mean square (rms), and were used to differentiate between climatic forcing and to evaluate the deglacial evolution of the North American ice-sheet complex (Reference Tarasov and PeltierTarasov and Peltier, 2004).

Here we present and illustrate a set of analyses that provides a basis for a more comprehensive and systematic quantitative assessment framework that addresses significant limitations which have hampered qualitative studies and which are only likely to intensify as numerical modelling, remote sensing, field and dating techniques improve and become more sophisticated and data-rich. A Geographic Information System (GIS) toolbox is developed which is capable of quantitatively and systematically assessing output from a numerical ice-sheet model against geomorphological information, with the aim of ranking particular model simulations according to their level of correspondence with field data. In this approach, two distinct geomorphological types are used to constrain model output: (i) end moraines associated with two distinct advanced ice-sheet positions during the last glacial cycle, providing a control on marginal extent, and (ii) the direction of glacial lineations (and, by inference, the direction of palaeo-ice flow), providing a control on subglacial dynamics.

Methodology

The main aim of this study is to present a quantitative methodology that automatically identifies which of a suite of numerical ice-sheet experiments best matches field data in terms of a level of correspondence of modelled ice margins with end moraines and basal flow with glacial lineations (Fig. 1). Automated proximity and conformity analysis (APCA) (Reference Napieralski, Li and HarborNapieralski and others, 2006) is used to quantify the level of agreement between predicted and observed ice margin positions, and this is coupled with automated flow direction analysis (AFDA) (Reference Li, Napieralski, Harbor and HubbardLi and others, 2007) which compares modelled basal ice-flow directions with glacial lineation directions as a second measure of model performance. The combined approach is illustrated using output from a numerical model of the Fennoscandian ice sheet (FIS), which was used because the technique development work was part of a larger study of the dynamics and chronology of the FIS (Reference NapieralskiNapieralski, 2005).

Fig. 1. Overall approach taken for analysis. For this study, the analysis involves the comparison of marginal data, followed by the comparison of flow direction data (option I).

APCA (Reference Napieralski, Li and HarborNapieralski and others, 2006) is a GIS-based technique that utilizes a system of buffers and overlays to systematically assess the level of correspondence between linear or curvilinear features (Fig. 2). A series of buffers (zones) around the modelled ice-sheet margin position are superimposed on a separate buffer that is drawn around the position of a known moraine. The percentage overlap between the buffers indicates the proximity of the modelled ice-sheet margin to the moraine. Where the moraine and the modelled ice-sheet margin are close, a large percentage of the overlap occurs in the buffers close to the ice-sheet margin, and the simulation consequently receives a higher APCA score (Fig. 2a). Low APCA scores indicate that the predicted ice margin is distanced from the end moraines, and thus the model has overestimated or underestimated the ice-sheet margin position. In addition, the form of the APCA curve (a plot of the cumulative buffer area with distance from the ice margin (buffer number;Fig. 2b)) is a measure of the conformity of the moraine and the ice margin, such that it indicates whether the moraine trace and ice-sheet margin are parallel to each other or at an angle (Fig. 2b). For this study, 150 buffers, each 10 km in width (1500 km in total), were generated around a modelled ice margin, and one buffer with a width of 10 km was generated around end moraines.

Fig. 2. APCA uses a system of GIS-based buffering to determine the general proximity (a) and parallel conformity (b) between linear features. The area under the curve is used to determine which modelled output fits the empirical data best, based on the distance and angle between features (modified from Reference Napieralski, Li and HarborNapieralski and others, 2006).

AFDA (Reference Li, Napieralski, Harbor and HubbardLi and others, 2007) is used to compare modelled basal ice-flow directions with glacial lineation directions in order to rank the performance of different model experiments. A numerical ice-sheet model predicts basal flow conditions for each gridcell including basal temperatures (e.g. cold or warm bed), and ice-flow velocity and direction. In field mapping, glacial lineations are recorded at a wide variety of spatial scales (from striae to streamlined hills) based on where these features occur, and on maps often represent the presence of multiple features. In AFDA, the directional offset between an observed direction of lineation and the model-predicted ice-flow direction in each gridcell is denoted as a residual, which ranges from 0° to 180°, where 0° indicates perfect coincident, parallel flow, and 180° indicates that the model-predicted ice-flow direction is opposite to the glacial lineation direction. The resultant mean and variance of residuals for the whole lineation set is calculated to indicate the correspondence between model- predicted ice-flow direction and field-observed lineation. Applying this calculation to various time slices of model output, temporal variations of the level of correspondence between model prediction and field observations can be illustrated to evaluate the performance of different model experiments (Fig. 3). In addition, AFDA also reports the occurrence (area percentages) of model-predicted ice-free (no ice)/frozen-bed conditions where observed lineations exist, which is a critical component of ice-sheet model validation and indicates poor agreement between the model prediction and field observations (Reference Li, Napieralski, Harbor and HubbardLi and others, 2007). Where there are two or more sets of overriding lineations (two or more ice-flow directions recorded), AFDA treats these as separate lineation sets formed at different times, and compares them separately against modelled flow directions. Finally, instead of applying AFDA to all possible locations, in the approach presented here the analysis focuses on a subset of sites within the maximum extent of the ice sheet, as these are areas of high-quality field data and complex ice-flow history of particular glaciological interest.

Fig. 3. Steps in applying AFDA. (a) Field-based glacial lineations and model outputs. (b) Overlay model outputs and field evidence to produce a series of residual datasets for different time slices. (c) Plot resultant mean of residual values against their corresponding time slices to identify temporal patterns of correspondence between predicted directions and field observations. (d) Frequency analysis (rose diagram) of selected time slices (e.g. d and f) provides detailed information on the distribution of residuals across the area and can be used to evaluate the level of correspondence (from Reference Li, Napieralski, Harbor and HubbardLi and others, 2007).

Combining APCA and AFDA, variations in scores on both measures can be tracked through time to investigate how the correspondence between numerical model output and geomorphological data varies during the temporal evolution of an ice sheet. In our initial analysis we prioritized APCA scores of margin matching to bracket a subset of optimal simulations and subsequently used AFDA to assess the level of correspondence between basal flow and glacial lineations (Fig. 1, option I). Alternatively, one could use AFDA as the primary method to select simulations that match critical palaeo-ice-flow indicators, and then use APCA to determine which of these simulations best matches key marginal positions as recorded by moraines (Fig. 1, option II). Finally, it would also be possible to normalize, weight and combine APCA and AFDA scores for selected moraines and flow-sets, thereby providing the user with the ability to control and customize a particular hierarchy of varied geomorphological evidence to best constrain the model trajectory.

Ice-sheet model

Since the analysis presented here is not intended to be a detailed evaluation of a particular model, but rather aims to demonstrate a novel approach to evaluating model performance against geomorphologic data, we selected a readily available model and used standard assumptions and formulations. We used a three-dimensional (3-D) finite- difference flow model that is a modified version of that used by Hubbard (Reference Hubbard1999, Reference Hubbard2006) and Reference HubbardHubbard and others (2006) to reconstruct the Scottish and the Icelandic ice sheets at 1 and 2 km resolution, respectively. Its formulation, assumptions and implementation are described in detail in Hubbard (Reference Hubbard1999, Reference Hubbard2006) and Reference HubbardHubbard and others (2006) and the parameter settings used in this study are presented in Table 1. The thermomechanical flow component requires boundary conditions of subglacial topography, geothermal heat flux and surface air temperature. Apart from the computation of longitudinal stresses, the constitutive relation and numerics of the model are standard (Reference NyeNye, 1953; Reference MahaffyMahaffy, 1976; Reference HuybrechtsHuybrechts, 1986; Reference Marshall and ClarkeMarshall and Clarke, 1997; Reference Payne and BaldwinPayne and Baldwin, 1999; Reference Van der VeenVan der Veen, 1999). The model is thermally coupled and dynamics at the glacier bed are handled with a prescribed Weertman-type sliding relation (Reference WeertmanWeertman, 1964) triggered when basal temperatures approach the pressure-melting point. Isostatic adjustment is handled using an elastic-lithosphere/relaxed-asthenosphere approach which has gained favour over other alternatives (e.g. Reference Le Meur and HuybrechtsLe Meur and Huybrechts, 1996; Reference HagdornHagdorn, 2003). Longitudinal stress gradients are computed using an ice-stretching algorithm that has been validated against measurements at Haut Glacier d’Arolla, Switzerland (Reference HubbardHubbard, 2000), and become significant when considering ice dynamics over rough and steep topography, where sliding leads to fast flow, and at calving margins. Ice lost to calving is included as an additional mass-balance term and is related to the geometry and thickness of the calving front in the manner of Reference Payne and SugdenPayne and Sugden (1990). The model has been successfully validated against EISMINT (European Ice-Sheet Modelling Initiative) type 1 and 2 experiments in two and three dimensions (Reference Huybrechts and PayneHuybrechts and Payne, 1996).

Table 1. Parameters, values and units used in the ice-flow model

Climate is coupled to the upper surface boundary of the flow model using a simple elevation/mass-balance function in the manner of Reference HagdornHagdorn (2003) where the net mass balance at any grid node varies parabolically about the equilibrium-line altitude (ELA) and plateaus at a maximum mass-balance value nominally related to the mean annual precipitation at that node. Though crude, this formulation has a number of advantages in that it is straightforward to apply and manipulate with minimal parameterization, accounts for the maritime-to-continental climatic gradients that exist across Scandinavia, and allows for climate change to be easily implemented through manipulation of the regional ELA field. This regional ELA field is derived from data based on contemporary glaciers held in the World Glacier Inventory (http://nsidc.org/data/glacier_inventory/query.html) which were interpolated across the model domain using multivariate spatial regression. Past climate cooling is initiated through depression of the regional ELA field, which in this study is scaled to the Greenland Icecore Project (GRIP) oxygen isotope time series. Though this climate-driver is simplistic, by manipulating the regional ELA field and introducing spatial precipitation gradients to emulate aridity with the onset of cold conditions, it captures the broad characteristics of the FIS mass-balance distribution (Reference HagdornHagdorn, 2003).

The coupled mass-balance–flow model is applied at a 10 km resolution to a basal topography generated from a combination of ETOPO5 (bathymetry), GTOPO30, and Swedish and Norwegian national elevation data. A realistic geothermal heat-flux distribution is prescribed for the basal thermal boundary condition based on the work from borehole measurements collated by Reference Näslund, Jansson, Fastook, Johnson and AnderssonNäslund and others (2005). Using previous studies as a guideline (e.g. Reference HubbardHubbard, 1999, Reference Hubbard2006; Reference Marshall and ClarkeMarshall and others, 2002; Reference HagdornHagdorn, 2003; Reference Zweck and HuybrechtsZweck and Huybrechts, 2005; Reference HubbardHubbard and others, 2006), a subset of the input parameters was selected and perturbed in a suite of 40 ensemble experiments. The primary input variables selected include those affecting the level of topographic representation, isostatic response time, mass balance (maximum accumulation rate and elevation) and ice dynamics (flow enhancement and the sliding parameter) (Table 2).

Table 2. A summary of the parameters altered during the course of the study. Only a subset of model experiments is used to illustrate the steps taken to differentiate between model outputs based on levels of correspondence with field data

Fennoscandian Ice Sheet Field Data and Modelling

Study area

During the Pleistocene, most of northern Europe experienced repeated glaciations (Reference Mangerud, Jansen and LandvikMangerud and others, 1996; Reference Sejrup, Larsen, Landvik, King, Haflidason and NesjeSejrup and others, 2000). The presence of glacial landforms, such as well-defined end moraines and glacial lineations, provides evidence for the extent and flow patterns of the ice sheet during the last glaciation (Fig. 4). Most evidence indicates that the last glaciation initiated in the Scandinavian Mountains (Reference FredinFredin, 2002) and reached its maximum along its western flanks between 28 and 20 kyr BP (limited by the continental slope), its southern flanks between 20 and 18 kyr BP, and its eastern flanks between 18 and 16 kyr BP (Reference Boulton, Dongelmans, Punkari and BroadgateBoulton and others, 2001; Reference RinterknechtRinterknecht and others, 2004, Reference Rinterknecht2006). In addition, there is an abundance of field data indicating that the core of the FIS was frozen (occurrences of preserved relict landscapes that show little or no erosion;see Reference Kleman and StroevenKleman and Stroeven, 1997; Reference Fabel, Stroeven, Harbor, Kleman, Elmore and FinkFabel and others, 2002; Reference Hättestrand and StroevenHättestrand and Stroeven, 2002; Reference Stroeven, Fabel, Hättestrand and HarborStroeven and others, 2002), which should have resulted in relatively thick ice sheets (Reference Kleman and StroevenKleman and others, 1997; Reference HättestrandHättestrand, 1998; Reference Kleman and HättestrandKleman and Hättestrand, 1999). Several numerical model studies of the extent and timing of the FIS (e.g. Fastook and Holmlund, 1994; Reference Boulton, Dongelmans, Punkari and BroadgateBoulton and others, 2001; Reference Siegert, Dowdeswell, Hald and SvendsenSiegert and others, 2001; Reference Charbit, Ritz and RamsteinCharbit and others, 2002; Reference HagdornHagdorn, 2003) and geomorphological and stratigraphical-based reconstructions (e.g. Reference Boulton, Smith, Jones and NewsomeBoulton and others, 1985; Reference Kleman and StroevenKleman and others, 1997; Reference SvendsenSvendsen and others, 1999), while broadly consistent for ice-sheet extent (due to the presence of end moraines to indicate maximum ice-sheet extent), yield disparate results for other ice-sheet characteristics, such as ice thickness. For example, whereas till consolidation, isostatic data and relict landscapes all are consistent with the former presence of a thick ice sheet (2700–3000 m) (Reference Boulton, Smith, Jones and NewsomeBoulton and others, 1985; Reference Kleman and StroevenKleman and others, 1997; Reference Dowdeswell and SiegertDowdeswell and Siegert, 1999; Reference FjeldskaarFjeldskaar, 2000), some botanical and other geological data are more consistent with the former presence of a thin ice sheet (2100–2500 m) (Reference Nesje, Anda, Rye, Lien, Hole and BlikraNesje and others, 1987; Reference Nesje and DahlNesje and Dahl, 1992; Reference Brook, Nesje, Lehman, Raisbeck and YiouBrook and others, 1996; Reference Vorren and PlassenVorren and Plassen, 2002).

Fig. 4. Distribution of end moraines and lineations used in this project. Moraines 1–4 are of LGM age, and 5–7 are of YD age (from Reference Boulton, Dongelmans, Punkari and BroadgateBoulton and others, 2001). Glacial lineations from Kiruna, Sweden, (Reference Kleman and StroevenKleman and others, 1997) were used to compare simulated ice- flow direction with lineations. Within this study area, there is an abundance of cross-cutting lineations indicating a complex ice- flow history (see upper left box). Glacial lineations that share common direction have been separated and classified according to their physical characteristics (from Reference Kleman and StroevenKlema and others, 1997). Thus the lineations are divided into distinct flow fans (as shown in upper right box).

Results

Margin comparison (APCA)

An APCA score of 15 000 (150 buffers at 100%) indicates that 100% of the predicted ice-sheet margin position falls within 10 km (one buffer ring) of the position of the mapped end moraine, which at this resolution represents the perfect match. For each experiment, we report the highest APCA values for each moraine (Table 3) as a normalized value (0– 1). Virtually all of the ensemble experiments corresponded well to LGM1 for extended periods, which is no surprise given that this moraine marks the extreme western offshore limit before falling off to the Norwegian continental shelf slope (Table 3; Fig. 5a). Thus, any ice-sheet advance under lowered sea level is ultimately going to be pinned at this significant break in slope. In contrast, APCA scores for LGM2–4 varied widely between experiments, and rarely did all three achieve high APCA scores during the same experiment. Due to its north–south orientation in Denmark, LGM2 was not reproduced well by any of the ensemble experiments (Figs 5a and 6a). Rather, the ice margins were almost perpendicular to LGM2, signifying a challenge in duplicating its inconsistent shape and orientation relative to the local topography and main ice accumulation centres.

Fig. 5. Normalized APCA scores (level of agreement) for experiment 13, plotted over time for the LGM moraines (top) and the YD moraines (bottom). Bold sections indicate the age range of dates (e.g. Reference Tschudi, Ivy-Ochs, Schluchter, Kubik and RainioTschudi and others, 2000; Reference RinterknechtRinterknecht and others, 2004) for each moraine. Results show which moraines were reproduced the best and the length of time this correspondence occurred. For example, best correspondence occurs for LGM1 and YD5, but while trying to attain a LGM4 configuration, the modelled ice-sheet extent overshoots the target (peaking twice in APCA score while advancing and retreating past the moraine).

Table 3. Summary of APCA scores from selected model simulations reporting the highest for each of the seven moraines for the model experiments. The scores are normalized, so that a score of 0.00 indicates a relative mismatch between modelled ice-sheet extent and a given moraine location and a score of 1.00 indicates a perfect correspondence between model prediction and field data. The APCA scores for each experiment were extracted only from time slices during which end-moraine development has been estimated to have occurred (from Reference Boulton, Dongelmans, Punkari and BroadgateBoulton and others, 2001). See Figure 4 for location of the moraines

The ensemble experiments frequently underestimated or overestimated the southern and southeastern extent of the ice sheet (Figs 5a and 6a). Between 26 and 22 kyr BP the predicted margin advanced over the mapped moraine position, and between 19 and 14 kyr BP the modelled ice extent retreated over the mapped moraine position for many of the runs presented here. This situation is caused by mass balance being overestimated across the majority of the ice sheet during LGM conditions and reflects the fact that the feedback between atmospheric frigidity and moisture content and thus the onset of extreme aridity (restricting mass balance) is not included in the model. Hence, the ice sheet advances and continues to expand rapidly, with ELA lowering up to the LGM, and does not experience the period of marginal stability which would occur as the ice sheet effectively starved itself of precipitation. However, the ice sheet is also highly sensitive, to the point of being intrinsically unstable and to some extent ‘self-governing’ in these far-field marginal zones characterized by shallow gradients, basal sliding and fast flow which are located up to 1000km away from the controlling topographic accumulation centres. Hence these experiments also reveal a strong element of decoupling at the fringes of the topographic model space.

Fig. 6. Modelled ice-sheet geometries with ice surface elevation contours corresponding to (a) LGM moraines 1–4 and (b) YD moraines 5–7.

YD moraines proved relatively difficult to reproduce in the experiments (Table 3), except for YD5. The model was able to generate ice margins that were moderately close to YD6, but they had relatively poor conformity, and, as a result, APCA scores were average throughout most of the model simulations. YD7 was the most difficult moraine for the model to reproduce;as the ice sheet experienced deglaciation, the modelled southeastern margin fails to retreat much further than northern mainland Europe, thus never reaching the location of the major YD moraines in southern Finland (see Figs 5b and 6b; Table 3). This may be attributed to the generalized mass-balance scheme used in this model, which overestimates precipitation rates over the central/eastern portion of the ice sheet, causing the ice sheet to deglaciate at a much slower rate. On the few occasions the model did agree well with this moraine (e.g. runs 21, 30), there was also a significant change in the LGM configuration so that the level of agreement for LGM3 and LGM4 moraines was poor. As stated, these issues are primarily a reflection of the application of an overly simplistic mass-balance parameterization based on modern-day distributions of precipitation and ELA and forced with a linearly scaled GRIP isotope record to invoke climate change. The result is a model that is handicapped when it comes to reproducing the complex temporal and spatial interplay between temperature, precipitation, mass balance and the concomitant marginal ice-sheet fluctuations across the model domain. What is most important here is that the APCA results quantitatively highlight the limitations of the model.

To evaluate overall model performance, APCA scores were examined for time periods consistent with estimated dates of moraine development (based primarily on Reference Boulton, Dongelmans, Punkari and BroadgateBoulton and others, 2001). Scores from the seven moraines (Table 3) were totalled to provide an overall level of correspondence between modelled ice extent and major moraines (Fig. 7). This allowed the selection of four ‘best fit’ runs (29, 30, 32 and 33; Fig. 7) that were then further evaluated with AFDA. Compared to run 29, run 30 has the basal sliding parameter increased, run 32 has enhanced thermodynamic coupling and run 33 has a reduced elevation-mass-balance feedback.

Fig. 7. Total APCA scores for the subset of model experiments used in this study, based on the normalized APCA scores for all seven moraines. The graph illustrates which experiments best agreed with all seven moraines and illustrates the relative ability of the model to reproduce LGM or YD moraines. The optimum scoring experiments were selected for AFDA (runs 29–33).

Flow direction analysis (AFDA)

Mean residuals between predicted and observed flow directions for the four flow sets were calculated from 105 to 10 kyr (1.0 kyr interval) for the four ‘best fit’ runs with the highest total APCA scores. The temporal variation of mean residuals reflects changing correspondence between model and field data throughout the 95 kyr modelled glacial cycle, due to alterations in ice-flow directions as a result of changing ice- sheet geometry (such as the location of ice divides) and basal thermal regime during stages of growth and decay. There were distinct periods when model output displayed relatively high levels of correspondence with each flow set (Fig. 8). For example, output corresponded well with flow set 38 at 75, 46 and 12 kyr, and with flow sets 24, 31 and 34 (which have similar flow directions) at 86, 66, 38 and 21 kyr. This potentially corroborates the Reference Kleman and StroevenKleman and others (1997) estimate that deglaciation fan 38 formed before the LGM, perhaps during marine oxygen isotope stage 5, and synchronous fans 31 and 34 potentially formed during the LGM.

Fig. 8. Temporal variations of calculated resultant mean and variance of residual values (a (38), c (31)) between predicted directions and flow datasets and occurrence of ice-free and frozen-bed conditions (b (38), d (31)) for flow datasets 31 and 38 from 105 yr to 10 kyr BP (flow fans 24 and 34 fluctuate in relative harmony with fan 31) (figure from Reference Li, Napieralski, Harbor and HubbardLi and others, 2007).

However, because the flow sets were treated as synchronous events when running AFDA, set 24 typically did not score well. Reference Kleman and StroevenKleman and others (1997) classify it as a time- transgressive deglacial system developed incrementally during a westward shift in ice-divide position. The age difference between proximal and distal parts of this swarm may be of the order 2–3 kyr. The non-synchronous character of this swarm introduces an error difficult to quantify when comparing with the true time-slice flow pattern generated by the ice-sheet model. If a particular flow set develops in a time-transgressive manner, then it is unlikely that the mean residual approach would capture this effectively. However, mapping patterns of residuals in a sequence of time-steps has the potential to reveal such a pattern of landform development.

AFDA results can be particularly useful in that they highlight complex changes in basal thermal conditions and flow directions during each model experiment. For example, the frequency distributions of mean residuals for flow sets 31 and 38 during run 30 at 75 kyr and 21 kyr show considerable spread in residual values (Fig. 9). At 75 kyr, the mean residual was relatively low for flow set 38 (27.4°), and high for flow set 31 (135.9°). At 21 kyr, in contrast, modelled ice- flow directions agree better with field data for flow set 31 (10.2° vs 101.0°). In addition, the agreement between model prediction and observation can be further assessed by examining the variance of residuals at different time slices (Fig. 9). For example, the variance of flow set 31 at 21 kyr is 0.72% and indicates a reasonable fit. However, the variance alone cannot be used to judge if there is a reasonable fit. For example, the flow set 31 at 75 kyr also has a low variance (0.98%). It does not mean there is a good fit in this time slice because the mean residual is 135.9°. Overall, variance is influenced by the size of the flow set but also reflects the presence of primary and secondary flow directions.

Fig. 9. Frequency distributions (rose diagrams) of residual values from run 30 at time slices 21 kyr (a, c) and 75 kyr (b, d) for flow sets 31 (a, b) and 38 (c, d) (from Reference Li, Napieralski, Harbor and HubbardLi and others, 2007).

Discussion

Efforts to integrate geomorphological data with numerical ice-sheet models have been limited by an inability to systematically measure levels of correspondence between model output and field data. The methodology presented here provides a first step in systematic, quantitative comparison of modelled timing, extent and general flow patterns against suites of field data. APCA provides the capability to quantify levels of correspondence between curvilinear features, such as predicted margin extent and end moraines, and to reveal the timing and duration of good correspondence. AFDA compares modelled ice-flow directions with glacial lineations to identify the temporal pattern of the offset between simulated and observed ice-flow directions. With a methodology that combines these two techniques, an ensemble suite of numerical experiments can be used to quantitatively assess models across parameter space in terms of their performance as measured by correspondence with geomorphological datasets.

In palaeo-ice-sheet modelling, where the degrees of freedom may outweigh the available constraints, instances of equifinality are likely, as widely differing combinations of input parameters and forcing may generate configurations that generally correspond with the available field evidence for ice extent and flow directions. In such circumstances, the problem is underdetermined and it will be hard to differentiate between models in terms of a general qualitative comparison with field data. However, the approach used here provides an objective, quantitative measure of correspondence. This allows for a systematic and thorough exploration of model performance across the range of viable parameter values. This has the potential to yield new insights into ice-sheet behaviour, and may also allow for the identification of key areas for differentiating between alternative ice-sheet models that can then guide new field data collection. For example, in this study experiments 30 and 32 had similar high levels of overall agreement with the seven moraines (Fig. 7;Table 3), despite differences in their basal boundary conditions (Table 1). Under these circumstances, the contrasting but complex interplay at key locations of modelled cold and warm basal ice which effectively protects and erodes the substrate respectively, could be tested in the field by recourse to predicted complex cosmogenic nuclide exposure histories at these sites.

Although discrimination between modelling runs can be accomplished in several ways, here margins were compared first, followed by ice-flow directions. Each moraine can be compared individually against model output, but caution should be taken when using individual moraines to discriminate between model experiments. For example, virtually all experiments yield a good correspondence with LGM1 and YD5 moraines, since the model is fairly insensitive to the substantial topographic relief of the western margin. Using only the APCA scores from these areas to assess the performance of the model would produce misleading results because, despite obtaining high APCA scores at these moraines during experiment 32, the model was unable to accurately reproduce LGM3 and YD7 moraines. In contrast, summing the APCA scores across all seven moraines generates a total score that indicates a more representative level of agreement for each experiment.

Alternatively, APCA scores could be weighted to match a segment of the ice sheet more prominently than other locations (e.g. if certain locations have better chronological constraints than others). For example, in examining the dynamics of the eastern margin of the FIS, a weighted value could be applied to the eastern moraines to influence the overall correspondence with field data. More importantly, the success of APCA is dependent on the distribution of end moraines, which needs to reflect the non-synchronous growth and decay behaviour of an ice sheet by using well-dated end moraines evenly distributed throughout the study area.

It is critical, however, to evaluate the distribution of residuals in more detail because the mean value may be a product of a small or large standard deviation of offset values, a bimodal distribution of values or an abundance of predicted frozen bed zones or ice-free areas. An experiment that predicts frozen bed conditions or no ice when a given flow set indicates some sliding should have occurred receives a poor AFDA score because it has failed to reproduce formative conditions for that particular set of flow traces and period of time. It is important to examine the residual distributions underlying the mean values to fully understand how accurately the model is reproducing the direction of glacial lineations.

Future applications of this method for this ice sheet could include the use of other flow sets, including various flow fan types, or other specific types of field evidence, such as eskers, drumlins or till fabric analyses rather than lineations. The inclusion of cosmogenic nuclide or radiocarbon dates (i.e. of moraines) could provide temporal constraints to further refine the range of acceptable model parameters, and thereby improve the accuracy of model simulations. Furthermore, in the particular example used here, the physics and spatial resolution of the ice-sheet model could also be improved. For example, the model is based on a 10km gridcell resolution and thus fails to capture critical subgrid topographic interactions, especially in the initial stages of mountain inception and the subsequent reduced flux to lower elevations through deep-cut valleys and fjords, which are less than 10km wide. Additionally, because the effect of subglacial hydrostatic pressures on ice flow is not included, the model may misrepresent the occurrence of ice streams, and because they exert an important control on the 3-D geometry of the ice sheet, by implication, the thermal regime. Thus, increasing the modelling resolution coupled with developments in the model physics would enhance these relationships and improve model efficacy.

It is important to note that APCA and AFDA were designed as general spatial analysis tools. Although in this paper they are illustrated in an application to the FIS, they can be applied to other ice-sheets to calibrate and validate model output with field evidence. Models that correspond well with field evidence can then be used to evaluate the dynamics of ice sheets during ice-sheet inception, growth and decay, including the stages of ice-sheet advance, streaming events and evolution of the subglacial thermal regime and topographic limits on ice-sheet elevation. Evaluating alternate models for ice sheets is particularly important given that CLIMAP estimates of the thickness and volume of global ice sheets during the LGM (Denton and Hughes, 1981) have been questioned and new reconstructions have been recommended (Reference Mix, Bard and SchneiderMix and others, 2001). To achieve the best possible results, any new model reconstruction should be constrained to geomorphological and chronological evidence, in order to more accurately reproduce ice-sheet thickness and volume. Reconstructions that are consistent with field evidence will provide more reliable data as input to global climate models which, in return, will be better equipped to reconstruct the environmental conditions during glacial inception, advance, maximum and decay.

Conclusions

We have developed two GIS-based techniques, APCA and AFDA, to quantify the level of correspondence between ice- sheet modelling runs and field data. Successful application of APCA requires the availability of well-defined, well-dated and evenly distributed evidence for ice-marginal positions. In the work presented here, this evidence consisted of LGM and YD moraines with radiocarbon and cosmogenic dates for the Fennoscandian ice sheet. AFDA can be used to assess ice-sheet model output by revealing discrete periods when simulated ice-flow directions agreed and conflicted with field evidence for former ice-flow directions (lineations) (i.e. Reference Kleman and StroevenKleman and others, 1997), but does assume a synchronous deposition of glacial lineations. Together, APCA and AFDA allow a user to assess modelling runs based on geomorphological data in terms of the full parameter space of model uncertainties and sensitivities. The methodology used in this study is not limited to the FIS, but rather can be used to evaluate numerical models that simulate any palaeo-ice sheet. APCA could also be used to compare simulations of valley glacier extents with sets of recessional moraines and trimlines, and simulations of present-day glaciers and ice sheets against current and historic ice-marginal positions. Those interested in applying or adapting this approach can consult Reference Napieralski, Li and HarborNapieralski and others (2006) and Reference Li, Napieralski, Harbor and HubbardLi and others (2007) for details on each method.