Towards a GIS assessment of numerical ice-sheet model performance using geomorphological data

Abstract

A major difficulty in assimilating geomorphological information with ice-sheet models is the lack of a consistent methodology to systematically compare model output and field data. As an initial step in establishing a quantitative comparison methodology, automated proximity and conformity analysis (APCA) and automated flow direction analysis (AFDA) have been developed to assess the level of correspondence between modelled ice extent and ice-marginal features such as end moraines, as well as between modelled basal flow directions and palaeo-flow direction indicators, such as glacial lineations. To illustrate the potential of such an approach, an ensemble suite of 40 numerical simulations of the Fennoscandian ice sheet were compared to end moraines of the Last Glacial Maximum and the Younger Dryas and to glacial lineations in northern Sweden using APCA and AFDA. Model experiments evaluated in this manner were ranked according to level of correspondence. Such an approach holds considerable promise for optimizing the parameter space and coherence of ice-flow models by automated, quantitative assessment of multiple ensemble experiments against a database of geological or glaciological evidence.

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 (Denton and Hughes, 1981; Peltier, 1994) and within North America and Europe (e.g. Porter, 1984; Ehlers and Gibbard, 2004a, b), 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. Kleman 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. Boulton and others, 1985; Peltier, 1994; Kleman and others, 1997; Ó Cofaigh 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; Siegert and others, 2001; Takeda 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. Bintanja 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 (Kleman and others, 2006). Ice-core oxygen and deuterium isotope records can provide temporal proxies for ice volume (e.g. Anandakrishnan and others, 1993) and palaeoclimate patterns, but, with the exception of promising recent work by Clarke and others (2005) and Lhomme 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. Clark, 1997; Ehlers and Gibbard, 2003; Hagdorn, 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. Marshall and others, 2002; Zweck and Huybrechts, 2005). There have been successful attempts to quantitatively assess model output against geomorphological evidence, but these are limited in scope. For example, Tarasov 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 (Tarasov 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) (Napieralski 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) (Li 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 (Napieralski, 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 (Napieralski 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 Napieralski and others, 2006).

AFDA (Li 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 (Li 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 Li 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 (1999, 2006) and Hubbard 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 (1999, 2006) and Hubbard 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 (Nye, 1953; Mahaffy, 1976; Huybrechts, 1986; Marshall and Clarke, 1997; Payne and Baldwin, 1999; Van der Veen, 1999). The model is thermally coupled and dynamics at the glacier bed are handled with a prescribed Weertman-type sliding relation (Weertman, 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. Le Meur and Huybrechts, 1996; Hagdorn, 2003). Longitudinal stress gradients are computed using an ice-stretching algorithm that has been validated against measurements at Haut Glacier d’Arolla, Switzerland (Hubbard, 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 Payne 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 (Huybrechts 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 Hagdorn (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 (Hagdorn, 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 Näslund and others (2005). Using previous studies as a guideline (e.g. Hubbard, 1999, 2006; Marshall and others, 2002; Hagdorn, 2003; Zweck and Huybrechts, 2005; Hubbard 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

Margin matching and flow-direction analysis

The comparison of modelled output and field data was made using ESRI’s ArcGIS and Arc/Info software and was carried out in two steps. First, the APCA level of correspondence between observed and predicted ice-marginal extent was determined for each 1000 year time-step of model runs. Modelled ice-sheet extent was compared with field data for maximum extent during the Weichselian (for simplification, despite age differences, denoted here as Last Glacial Maximum (LGM)) and the Younger Dryas (YD) stade. The output was ranked according to the level of agreement with field data, and for the model simulations that corresponded best with end-moraine suites, the second step in the analysis was a determination of the level of correspondence between observed and predicted ice-flow directions using AFDA (Fig. 1, option I).

For AFDA, specific parts of the study region were selected for the determination of residuals between modelled and field-based ice-flow directions. The criteria for the selection of these locations included (i) areas with a complex ice-flow history as shown by cross-cutting features (this is generally true for areas which varied in proximity to consecutive ice sheet divide locations), (ii) areas of specific glaciological interest (e.g. locations near major ice streams) and (iii) areas with an abundance of high-quality field data. Glacial lineations indicating ice-flow directions were grouped according to flow sets or fans (defined by Boulton and Clark, 1990a, b), and in this study flow sets by Kleman and others (1997) were digitized, designated direction values (in degrees) and then rasterized to the same resolution as model output (10 km).

Fennoscandian Ice Sheet Field Data and Modelling

Study area

During the Pleistocene, most of northern Europe experienced repeated glaciations (Mangerud and others, 1996; Sejrup 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 (Fredin, 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 (Boulton and others, 2001; Rinterknecht and others, 2004, 2006). 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 Kleman and Stroeven, 1997; Fabel and others, 2002; Hättestrand and Stroeven, 2002; Stroeven and others, 2002), which should have resulted in relatively thick ice sheets (Kleman and others, 1997; Hättestrand, 1998; Kleman and Hättestrand, 1999). Several numerical model studies of the extent and timing of the FIS (e.g. Fastook and Holmlund, 1994; Boulton and others, 2001; Siegert and others, 2001; Charbit and others, 2002; Hagdorn, 2003) and geomorphological and stratigraphical-based reconstructions (e.g. Boulton and others, 1985; Kleman and others, 1997; Svendsen 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) (Boulton and others, 1985; Kleman and others, 1997; Dowdeswell and Siegert, 1999; Fjeldskaar, 2000), some botanical and other geological data are more consistent with the former presence of a thin ice sheet (2100–2500 m) (Nesje and others, 1987; Nesje and Dahl, 1992; Brook and others, 1996; Vorren 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 Boulton and others, 2001). Glacial lineations from Kiruna, Sweden, (Kleman 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 Klema and others, 1997). Thus the lineations are divided into distinct flow fans (as shown in upper right box).

Field data

Seven end moraines from LGM and YD age were digitized from previous reconstructions of the FIS (primarily from Kleman and others, 1997; Boulton and others, 2001) (Fig. 4). Four marginal positions were selected to represent the LGM, and three to represent the YD. These moraines were selected because they are relatively well defined, represent different sectors of the ice sheet, are reasonably well dated (Tschudi and others, 2000; Rinterknecht and others, 2004, 2006) and thus provide two distinct markers for ice-marginal position during the last 26 000 years of ice-sheet history. Glacial lineations, defined here as landforms that are produced parallel to the local basal ice flow, such as drumlins, crag- and-tails, flutes, glaciotectonic folds and striae and sedi- mentological imprints such as till fabrics, were digitized and grouped into flow fans (from Kleman and others, 1997).

Although four separate areas were chosen for ice-flow direction analysis as part of our larger study, we use the Kiruna area here as a case study of the methodology. Flow fans within this study area were extracted and lineations within each fan that were outside the study area were ignored. The four flow sets are: 24 (a deglaciation fan), 31 and 34 (synchronous fans from the LGM) and 38 (a deglaciation fan from marine oxygen isotope stage 5) (flow-set numbering scheme from Kleman and others, 1997; Fig. 4). Drumlin characteristics and age (Hättestrand and others, 2004) and cosmogenic isotope dating of deglacial landscapes (Stroeven and others, 2003) broadly confirm these age relationships. The number of lineations within the flow sets ranged from 313 (set 38) to 18 (set 34), and had a predominant flow direction towards either the north (sets 31, 34), northeast (set 24) or southeast (set 38). Each lineation within a particular flow fan was designated a direction value (in degrees) and then rasterized to correspond with modelling output. Of importance to the technique is that the glacial lineations were given a direction based on assumptions made on general ice-flow patterns.

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. Tschudi and others, 2000; Rinterknecht 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 Boulton 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 Boulton 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 Kleman 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 Li and others, 2007).

However, because the flow sets were treated as synchronous events when running AFDA, set 24 typically did not score well. Kleman 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 Li 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 (Mix 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. Kleman 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 Napieralski and others (2006) and Li and others (2007) for details on each method.