The impact of model resolution on the simulated Holocene retreat of the southwestern Greenland ice sheet using the Ice Sheet System Model (ISSM)

Geologic archives constraining the variability of the Greenland ice sheet (GrIS) during the Holocene provide targets for ice sheet models to test sensitivities to variations in past climate and model formulation. Even as data–model comparisons are becoming more common, many models simulating the behavior of the GrIS during the past rely on meshes with coarse horizontal resolutions (≥10 km). In this study, we explore the impact of model resolution on the simulated nature of retreat across southwestern Greenland during the Holocene. Four simulations are performed using the Ice Sheet System Model (ISSM): three that use a uniform mesh and horizontal mesh resolutions of 20, 10, and 5 km, and one that uses a nonuniform mesh with a resolution ranging from 2 to 15 km. We find that the simulated retreat can vary significantly between models with different horizontal resolutions based on how well the bed topography is resolved. In areas of low topographic relief, the horizontal resolution plays a negligible role in simulated differences in retreat, with each model instead responding similarly to retreat driven by surface mass balance (SMB). Conversely, in areas where the bed topography is complex and high in relief, such as fjords, the lower-resolution models (10 and 20 km) simulate unrealistic retreat that occurs as ice surface lowering intersects bumps in the bed topography that would otherwise be resolved as troughs using the higher-resolution grids. Our results highlight the important role that high-resolution grids play in simulating retreat in areas of complex bed topography, but also suggest that models using nonuniform grids can save computational resources through coarsening the mesh in areas of noncomplex bed topography where the SMB predominantly drives retreat. Additionally, these results emphasize that care must be taken with ice sheet models when tuning model parameters to match reconstructed margins, particularly for lower-resolution models in regions where complex bed topography is poorly resolved.

1 Introduction

As the cryosphere community continues to make strides in understanding processes that govern variability of the present-day ice sheets, geologic proxies constraining past ice sheet change provide important clues as to how ice sheets may have responded to past climate change (Alley et al., 2010). Decades of research have led to the development of high-resolution geologic reconstructions that detail the spatial pattern and rate of retreat of the Greenland ice sheet (GrIS) over the last deglaciation as it evolved towards its present-day geometry (Weidick, 1968; Bennike and Bjorck, 2002; Young and Briner, 2015).

Southwestern Greenland is an area that experienced a large reduction in ice sheet extent. The ice margin retreated on the order of 150 km inland from the present-day coastline in response to warming during the early and middle Holocene (Briner et al., 2016). This landscape is punctuated by widely traceable moraine sequences (Weidick, 1968; Ten Brink and Weidick, 1974) that extend nearly 600 km throughout western Greenland and provide a constraint on the past retreat pattern of the GrIS in this region; the chronology of these moraines continues to be refined (Weidick et al., 2004, 2012; Young et al., 2013; Larsen et al., 2015; Lesnek and Briner, 2018). This history provides a benchmark for ice sheet model–data comparisons that will further enhance our understanding of the processes that influenced GrIS variability during the past while at the same time helping to highlight deficiencies in existing model frameworks.

Currently, ice sheet models simulating the evolution of the GrIS are focused either on long-term spin-ups over a glacial cycle (Huybrechts, 2002; Applegate, 2012) or its evolution during the last deglaciation (Tarasov and Peltier, 2002; Simpson et al., 2009; Lecavalier et al., 2014, 2017; Buizert et al., 2018). While many of these studies were primarily concerned with capturing the overall mass changes of the GrIS, one lineage of studies incorporated datasets of past GrIS change to develop a data-constrained model of its evolution over the last deglaciation. This was achieved by pairing an ice sheet model with a glacial isostatic adjustment and relative sea level model (Tarasov and Peltier, 2002; Simpson et al., 2009; Lecavalier et al., 2014, 2017). By incorporating data constraining the location of the GrIS beyond the present-day coastline, its vertical extent through time (i.e., ice thinning records), and records of relative sea level, the studies by Lecavalier et al. (2014, 2017) represent the most comprehensive model of GrIS change during the last deglaciation, with the results recently being compared against geologic archives of ice margin change (Larsen et al., 2015; Young and Briner, 2015; Sinclair et al, 2016).

While the models of Lecavalier et al. (2014, 2017) capture the timing and retreat pattern associated with the deglaciation in many locations well, large mismatches occur, particularly in southwestern Greenland and in areas where fast flow may have dominated (Young and Briner, 2015; Sinclair et al., 2016). The climatology used to force ice sheet models through time remains a primary source of uncertainty, and great strides have been made to improve our understanding of past climate history in Greenland through improved reconstructions of temperature (e.g., Kobashi et al., 2017; Lecavalier et al., 2017) and methods involving data assimilation of paleoclimate proxies with climate model output (Hakim et al., 2016; Buizert et al., 2018). Although recent experiments have investigated sensitivities to model formulation (Zekollari et al., 2017) and horizontal resolution over past climates (Zekollari et al., 2017; Seguinot et al., 2016; Golledge et al., 2012), testing the sensitivity of simulated ice retreat to the ice flow dynamics model (i.e., the level of complexity in its numerical approximations) and to model resolution, both in time and space, still remains an important area of research.

With regards to model setup, the use of coarse model resolutions (≥10 km grid spacing) might explain some of the model–data discrepancy (Larsen et al., 2015; Young and Briner, 2015; Sinclair et al., 2016). Driven by how well models resolve subglacial topography, simulations of the present-day GrIS have shown an important dependence on model resolution for accurately simulating ice flow across Greenland (Greve and Herzfeld, 2013; Aschwanden et al., 2016). Dependence on model resolution also extends to modeling future ice mass loss, where higher-resolution models simulate more mass loss than models with lower resolutions (Greve and Herzfeld, 2013). Although some work has focused on model resolution and its impact on simulated mass flux from the GrIS for in the present and the future, how model resolution affects simulated retreat in paleo-ice-sheet modeling studies is not well constrained. Prior work demonstrates that low-resolution grids limit a model's ability to capture features such as ice streams and marine-terminating outlet glaciers (Aschwanden et al., 2016), which might be on the order of a few kilometers in width. Additionally, in land-terminating portions of an ice sheet, low model resolution may lead to large jumps in the snowline, which ultimately can lead to large advances or retreats in the ice margin on the order of the model resolution (Young and Briner, 2015; Sinclair et al., 2016), and therefore limit the model's ability to capture smaller-scale ice marginal fluctuations (i.e., km scale).

In this study, we present results from regional ice sheet modeling experiments in southwestern Greenland during the Holocene using the three-dimensional thermomechanical Ice Sheet System Model (ISSM). We build on earlier efforts that focused on this ice sheet sector (e.g., Van Tatenhove et al., 1995). Ice model resolution is the primary target for assessment here, with four separate simulations being run, each with its own horizontal resolution ranging from 20 to 2 km. In this study we do not attempt to obtain a perfect match between the simulated model retreat and that derived from the geologic reconstructions, as that requires further sensitivity studies that are not the current motivation for this work. Instead, since model resolution is a constraint that is typically not explored when studying the past due to the computational cost, in this study we aim to determine whether increased model resolution is worth the computation time for simulating past ice sheet retreat.

2 Model description and setup

2.1 Ice sheet model

We use the Ice Sheet System Model v4.13 (ISSM; Larour et al., 2012), a finite-element, thermomechanical ice sheet model. We choose the higher-order approximation of Blatter (1995) and Pattyn (2003), hereafter referred to as BP, to solve the momentum balance equations. Although recent work has used the higher-order approximation in simulations over past time periods (Zekollari et al., 2017), this ice flow approximation is still rarely used when simulating over paleoclimate timescales. We use this approximation, however, as our choice is based upon representing the past dynamics of the ice sheet history as best as possible even though computational time is increased over conventional paleoclimate ice sheet models using the more common shallow ice approximation (SIA; Hutter, 1983).

The model domain for this study (Fig. 1) focuses on the southwestern region of Greenland, where geologic proxies detail Holocene ice retreat from the present-day coastline (Weidick, 1968; Ten Brink and Weidick, 1974). By constraining our domain to southwestern Greenland, the number of mesh elements within the model can be minimized when compared to modeling the entire GrIS, thus reducing the computational load. The model domain extends from the present-day coastline to the ice sheet divide. The southern and northern borders of the domain coincide with areas of minimal north-to-south across-boundary flow based upon present-day ice surface velocities from Rignot and Mouginot (2012). The associated boundary conditions used to drive the model are discussed in Sect. 2.6.

Figure 1 Present-day interferometric synthetic aperture radar (InSar) ice surface velocities from Rignot and Mouginot (2012) for the Greenland ice sheet. The regional model domain is highlighted in black. Marked locations correspond to (1) Jakobshavn Isbræ, (2) Kangiata Nunâta Sermia (KNS), and (3) Nuuk.

2.2 Domain discretization

Typically, prior paleoclimate ice sheet modeling efforts across Greenland have used uniform meshes with a horizontal resolution of 20 km (Simpson, 2009; Lecavalier, 2014); a more recent model used a 10 km horizontal mesh resolution (Buizert et al., 2018). For the following experiments, the first three models are generated using a uniform triangular grid with horizontal resolutions of 20, 10, and 5 km. The fourth model (herein referred to as nonuniform) relies on anisotropic mesh adaptation, whereby the element size varies as a function of the bed topography (see Fig. S1 in the Supplement for visualization of mesh resolutions). The maximum horizontal mesh resolution is 15 km where gradients in the bed topography are smooth (primarily the interior bed over the domain) and becomes progressively finer in areas of high relief, with a minimum horizontal resolution of 2 km (mainly in fjord regions). The bed topography for each model is taken from BedMachine Greenland v3 (Morlighem et al., 2017) and is initialized with present-day ice surface elevation from the Greenland Ice Mapping Project digital elevation model of Howat et al. (2014). In Fig. 2, the corresponding bed height is shown for each model detailing the associated differences based on horizontal grid resolution. Generally, the bed topography is captured better using the higher-resolution mesh, with the nonuniform mesh (Fig. 2a) being able to best resolve valleys along the present-day coastline. The 5 km mesh captures the same general topographic features as the nonuniform mesh, albeit with less detail. At 10 km, individual valleys become unresolved, particularly around Nuuk and Jakobshavn Isbræ (see Fig. 1 for locations). The 20 km model fails to capture any topographic features that would hold glacier outlets.

Figure 2 Associated bed topography maps for the nonuniform high-resolution mesh (a), uniform 5 km mesh (b), uniform 10 km mesh (c), and uniform 20 km mesh (d).

2.3 Surface mass balance (SMB)

We use the positive degree day method outlined in Tarasov and Peltier (1999) to construct the necessary accumulation and ablation history used to drive our ice sheet model during the past from monthly mean temperature and precipitation fields. The spatial monthly mean surface air temperature and precipitation climatology spanning the period 1980–2010 is taken from Box et al. (2013). The surface air temperatures are then scaled based upon isotopic variations in the Greenland Ice Core Project (GRIP) δ¹⁸O record (Dansgaard et al., 1993) as follows:

$$\begin{array}{}\text{(1)}& \mathit{\delta }\mathrm{\Delta }T\left(t\right)=\mathit{\delta }\left({}^{\mathrm{18}}\mathrm{O}\left(t\right)+\mathrm{34.83}\right),\end{array}$$

where d=2.4 ^∘C ‰⁻¹ (Huybrechts, 2002). Anomalies from Eq. (1) are applied to the present-day climatology to create a temperature forcing back through time. Precipitation rate changes 7.3 % for every 1 ^∘C of temperature change derived in Eq. (1) (Huybrechts, 2002). For the positive degree day scheme, snow melts first (0.006 m  ^∘C⁻¹ day⁻¹) followed by bare ice (0.0083 m  ^∘C⁻¹ day⁻¹) with allocation for the formation of superimposed ice included (see supplemental information in Le Morzadec et al., 2015). The temperature forcing is adjusted throughout the run using a lapse rate correction of 5 ^∘C km⁻¹ (Abe-Ouchi et al., 2007) to account for changes in ice surface height throughout the simulation, while elevation-dependent desertification is included (Budd and Smith, 1981) to ensure reduction in precipitation by a factor of 2 for every kilometer change in ice sheet surface elevation. Further details regarding the positive degree day and accumulation scheme implemented within the ISSM can be found in Le Morzadec et al. (2015).

2.4 Thermal model and basal drag

The thermal evolution of the ice is captured using an enthalpy formulation described in Aschwanden et al. (2012), which includes formulations for both temperate and cold ice. Transient surface air temperatures are imposed at the ice surface, while geothermal heat flux (from Shapiro and Ritzwoller, 2004) is applied at the base. The model contains five vertical layers, with spacing between layers decreasing modestly towards the base. To simulate the vertical distribution of temperature within the ice sheet, we rely on quadratic finite elements (i.e., P1 × P2) along the z axis as a means for our vertical interpolation. Details of the implementation and description of these higher-order vertical finite elements can be found in Cuzzone et al. (2018). Through using higher-order finite elements as a means for vertical interpolation, this method allows the ice sheet model to capture sharp thermal vertical gradients particularly at the bed, which is an improvement over conventional methods using a linear vertical interpolation, despite having fewer vertical layers. This ultimately limits the necessity for a large number of vertical layers in our ice model and therefore decreases computational load.

To capture spatial variations in sliding, the spatially varying basal drag coefficient (k) in Eq. (2) is derived using inverse methods (Morlighem et al., 2010; Larour et al., 2012), providing the best match between modeled and interferometric synthetic aperture radar (InSAR) surface velocities (Rignot and Mouginot, 2012). This is performed independently for each model resolution.

$$\begin{array}{}\text{(2)}& {\mathit{\tau }}_{\mathrm{b}}=-k^{\mathrm{2}}N{v}_{\mathrm{b}},\end{array}$$

where the τ_b represents the basal stress, N represents effective pressure, and v_b represents magnitude of the basal velocity.

Since the drag coefficient (k) derived using this methodology is constrained to modern day, we adopt an approach based upon Hindmarsh and LeMeur (2001) and Greve (2005) to construct a spatially varying temperature-dependent scaling parameter (λ) as a function of time.

$$\begin{array}{}\text{(3)}& {\mathit{\lambda }}_{\mathrm{t}}=e^{\left(T_{\mathrm{b}\left(\mathrm{modern}\right)}-T_{\mathrm{b}\left(\mathrm{t}\right)}\right)/\mathit{\alpha }},\end{array}$$

where T_b(modern) is the basal temperature relative to pressure melting derived from a thermal steady-state computation for modern day (Seroussi et al., 2013), T_b(t) is the basal temperature relative to pressure melting at time t, and α is a constant scaling factor (^∘C) often referred to as the sub-melt parameter (Hindmarsh and Le Meur, 2001). For these simulations α is set equal to 5. This number was chosen as it allows for a Last Glacial Maximum (LGM) GrIS simulated ice volume that is consistent with other ice sheet models that restrict ice extent to only present-day land (Applegate et al., 2012; Robinson et al., 2011). It is noted that values for this parameter lack a theoretical basis (Hindmarsh and Le Meur, 2001) and are often set to a value that prevents numerical instabilities from arising. Lastly, we scale the spatially varying basal drag coefficient (k) as a function of λ, with maximum values capped at 300 to limit numerical instabilities that may arise from unreasonably large numbers:

$${\mathit{\tau }}_{\mathrm{b}}=-{\mathit{\lambda }}_{\mathrm{t}}min\left(\mathrm{300},k^{\mathrm{2}}\right)N{v}_{\mathrm{b}}.$$

For this approach, the basal stress (τ_b) increases as the basal temperatures decrease relative to present day, with virtually no sliding occurring for high values of k. Conversely, the basal stress τ_b decreases as basal temperatures increase, with high sliding for low values of k. Lastly, the ice hardness, B, is temperature dependent following the rate factors given in Cuffey and Paterson (2010, p. 75). We initialize B by solving for a present-day thermal steady state (Cuzzone et al., 2018), while during forward runs B evolves transiently through time.

2.5 Experimental setup and boundary conditions for the regional domain

We impose Dirichlet boundary conditions for the southern, northern, and ice divide boundaries, while the flux at the ice front in unconstrained. To create the necessary transient boundary conditions (ice thickness, temperature, and velocity), we perform a continental-scale GrIS simulation from the LGM (21 500 years ago) to present day. This continental-scale simulation uses the BP ice flow approximation and is performed on a five-layer nonuniform mesh ranging in horizontal resolution from 3 km in areas of high present-day surface velocities to 20 km over the ice interior. It is performed using forcings and parameterizations similar to the regional model, as described in Sect. 2.3, 2.4, and 2.5.

For the regional model, we initialize the model with present-day geometry and run a relaxation centered at 12 000 years ago, applying the appropriate interior ice boundary conditions of ice thickness, ice temperature, and the x and y component of ice velocity from the continental-scale GrIS simulation. This time period is chosen as the ice margin over southwestern Greenland was near or at the present-day coastline with the margin remaining stable during this interval (Young and Briner, 2015). The relaxation simulates 20 000 years until the ice volume is in equilibrium. From here, the models are run transiently to present day. Since the ISSM currently does not have the capability of modeling solid earth viscoelastic deformation transiently, we include an offline time-dependent forcing that accounts for changes in relative sea level from glacial isostatic adjustment (Caron et al., 2018), which modifies the land area available for glaciation and impacts the presence of floating ice. While grounding line migration is simulated in these experiments, calving and submarine melting of floating ice are not included.

2.6 Present-day thermal steady-state ice surface velocities

The thermal steady-state ice surface velocities for the present day are shown for each individual model (Fig. 3). Generally, representation of faster ice flow along the coast improves with increasing resolution (i.e., increasing RMSE for lower-resolution models compared to Rignot and Mouginot, 2012). This is primarily attributed to an improved representation of subglacial topography and ice thickness in the more highly resolved models (Aschwanden et al., 2016). As many of the outlet glaciers along this margin have troughs that are on the order of a few kilometers in width, the lower-resolution models (10 and 20 km) do not fully resolve the fast-flowing ice streams of Jakobshavn Isbræ and outlets to its north. Outlet glaciers in the southern portion of the domain near Kangiata Nunâta Sermia (KNS) are also less well resolved in the lower-resolution models, although the general swath of higher velocities is captured well for most fast-flowing areas of the ice sheet when compared to the observations (Rignot and Mouginot, 2012). It is noted that the nonuniform mesh represents these faster flow features best when compared to observations in most regions due to its high resolution. Accordingly, the associated mass flux at the ice margin is representative of these differences in model resolution, with the 10 and 20 km models having an approximately 25 % increase in mass flux (GT yr⁻¹) compared to the observations, while the 5 km and nonuniform mesh have an approximately 5 % to 9 % increase in mass flux compared to observations.

Figure 3 Present-day steady-state ice surface velocities for each individual model and their differences from observations (Rignot and Mouginot, 2012, shown in Fig. 1).

3 Results

3.1 Relaxed state at 12 ka

The models are relaxed for 20 000 years using a constant climate corresponding to 12 ka (Fig. 4a). The four models simulate decreasing ice volume with decreasing model resolution; the 20 km model simulates approximately 6 % less total ice volume than the nonuniform model. Ice surface velocities for the relaxed states (Fig. 5) depict the role of model horizontal resolution in capturing fjords and narrow outlets close to the model domain edge (i.e., present-day coastline). Generally, the two higher-resolution models (nonuniform and 5 km) capture narrow, fast flow in these outlets, whereas the lower-resolution models simulate a more diffuse pattern in ice surface velocities. This is primarily the case in areas where the bed topography is better resolved in the higher-resolution models and, therefore, confines the flow to narrow outlets. Ice velocities are reduced in the higher-resolution models for areas where the low bed topography that channels ice flow is interrupted by bumps and depressions in the bed. These features become less resolved in the lower-resolution models, with the 20 km model simulating much higher ice surface velocities in the Nuuk and Jakobshavn areas. Consistently, ice mass flux (in GT yr⁻¹) along the ice front (at the present-day coastline) is 34 % and 14 % higher than the nonuniform model for the 20 and 10 km model, respectively, which is the primary driver for lower simulated ice volumes for the relaxed 12 ka state.

Figure 4 (a) Ice volume evolution for the 12 ka constant climate relaxation. (b) Transient ice volume evolution for the simulations from 12 ka to present day.

Figure 5 Relaxed ice surface velocities at 12 ka for each model and differences from the nonuniform model.

3.2 Simulated ice volume (12 ka to present day)

The ice volume evolution for each model is shown in Fig. 4b. Generally, the 20 and 10 km models simulate the lowest present-day ice volumes, but they also begin at 12 ka with lower ice volumes than the higher-resolution models. Each model follows a similar trend, with ice volume loss occurring between 12 and 1 ka, followed by an increase in ice volume to present day.

3.3 Large-scale simulated retreat

Figure 6 shows the simulated extent at 11.2, 10.5, and 9.5 ka for each model (middle and bottom row). All four models generally show ice retreat from the coastline occurring between 11.5 and 11.2 ka in the northern portion of the model domain, whereas the ice margin experiences little to no retreat farther south. Despite differences in horizontal resolution, all models show a similar magnitude and pattern of retreat, with higher-resolution models depicting details in the ice margin similar in scale and sinuosity to the mapped pattern of moraines (see Fig. 6a for mapped moraines). Similarity between the magnitude and pattern of retreat also occurs at 10.5 and 9.5 ka amongst all models. In contrast to the northern portion of the model domain, the southern portion features a simulated retreat that varies widely based upon model resolution. For example, at 10.5 ka the higher-resolution models (5 km and nonuniform) exhibit little retreat from the coastline, whereas the 10 and 20 km models show upwards of 50–60 km of retreat. Differences in retreat between the higher- and lower-resolution models are further seen at 9.5 ka. Over land-terminating portions in the southern area, the modeled ice margin retreats similarly within all models; however, in the fjord regions (e.g., inland of Nuuk), only the 10 and 20 km models show ice margin retreat (of ∼50–70 km), whereas the higher-resolution models exhibit no ice margin retreat.

Figure 6 (a) Mapped moraines and the existing chronology of ice retreat over the northern and southern portion of our domain. The mapped moraines and corresponding ages of retreat were taken from Lesnek and Briner (2018). (b, c) Simulated ice sheet margin for the different model resolutions shown over locations in the northern (b) and southern (c) domain. The present-day ice thickness is shown, derived from Morlighem et al. (2017) and Howat et al. (2014).

3.4 Simulated retreat (12 ka to present day) – along flow lines

To better illustrate simulated ice margin behavior through time, we analyze ice retreat along five specific flow lines (A through E; Fig. 7) across the domain. In Fig. 7, ice retreat along flow line A is shown for each model. All models show a similar trend with the highest retreat rate (upwards of 50–100 m yr⁻¹) occurring between approximately 11.5 and 10 ka. Between ∼10 ka and the present, all simulated ice margins generally reside within 10 km of the present-day ice margin. The retreat history simulated by the 10 and 20 km models exhibits a relatively stable ice margin for much of this period, whereas the higher-resolution models (i.e., nonuniform and 5 km) depict an ice margin that is characterized by higher-frequency variability on the order of 5–8 km. The retreat history along flow lines B and C is consistent in timing and pattern to flow line A (shown in Figs. S2 and S3).

Figure 7 (a) Simulated retreat along flow line A for each model. (b) Ice surface profiles shown at 500-year intervals (blue: older; green: younger; the red line indicates the simulated present-day ice surface profile), with the underlying bed topography (filled black area). The red tick mark on the x axis denotes the present-day ice margin (Howat et al., 2014). Readers should refer to this figure for locations of flow lines used in this study.

Differences in the bed topography between the four models reflect model resolution, with higher-resolution models capturing topography closer to reality (Fig. 7). Nevertheless, the bed topography along the flow line A is similar among the different models owing to the low-relief topography in this region. The low-angle ice surface responds to surface melt similarly among the four models along flow line A (see Fig. S2 for surface temperature and SMB along flow lines) and is likely why the retreat history is similar. Along flow line B (Fig. S3), bed topography in all models exhibits increasing elevation into the ice sheet interior. Whereas the nonuniform and 5 km models capture a trough between 30 and 120 km along the flow line, this feature is subtle in the 10 km model and nonexistent in the 20 km model. Similar to flow line A, however, the simulated retreat in flow line B is similar amongst the different models both in rate and magnitude of retreat. Generally, the ice margins exhibit retreat forced primarily through surface lowering in response to negative SMB (Fig. S2) because the ice margin retreats similarly through areas of varying bed topography. At approximately 150 km along the flow line, the ice retreats into an area with higher bed topography, and the ice surface profiles become steeper, thereby raising the equilibrium line altitude (ELA) and stabilizing the ice margin. Along flow line C (Fig. S4), the bed topography for the nonuniform and 5 km model is characterized by low elevations at the beginning and end of the flow line, with a peak in elevation in the middle, while the 10 and 20 km models exhibit a more consistent bed elevation along the flow line. The upward slope of the bed for the nonuniform and 5 km model tends to slow the retreat of the ice margin along flow line C during the early Holocene, although differences between the 10 and 20 km models are not dramatic. After 8 ka, ice retreat stabilizes in all models, similar to flow lines A and B, although the nonuniform and 5 km models exhibit marginal ice fluctuations on the order of 5–8 km.

In the southern portion of the domain, the fjords within the Nuuk region dominate the landscape. The ice margin in flow line D (Fig. 8) remains fixed for the entire simulation in the nonuniform and 5 km model simulations. In contrast, both the 10 and 20 km models depict retreat at ∼9.5 ka, after which the 10 km model quickly stabilizes and the 20 km model exhibits variability up to present day. All models fail to simulate a present-day ice margin that comes close to today's observed ice margin (Fig. 8). The lower-resolution models simulate retreat in this region on the order of 30 to 50 km, which is controlled primarily by the bed topography. In reality, a trough extends much of the distance along this flow line, which is captured well by both the nonuniform and 5 km mesh, where depths reach ∼500 m below sea level. Consequently, the nonuniform and 5 km models are better able to capture the stress balance and mass transport, as they simulate more realistic ice flow and delivery of ice mass to the margin in this region. In the 10 km model, surface lowering intersects a bed bump that is above sea level at approximately 40 km along flow line D. Inland of this rise in the bed, the 10 km model bed contains a shallow trough, which is capable of sustaining the ice margin throughout the remainder of the simulation. The 20 km model lacks any clear trough and instead captures a significant rise in the bed topography at 70 km along the flow line, where the other models resolve a trough. As the ice surface lowers along flow line D, it becomes increasingly influenced by this bed feature. Due to the upward slope and horizontal top of this bed feature, the margin varies in response to the high-frequency climate variability during the Holocene (Fig. S2).

Figure 8 Same as in Fig. 7 but for flow line D.

Flow line E (Fig. S5) follows a narrow and shallow trough south of flow line D. This shallow trough is only captured completely in the nonuniform model, although the 5 km model captures low topography along the ice margin. In the 10 and 20 km models there is no indication of a trough and instead the bed topography is high and generally flat. Similar to flow line D, downwasting via negative SMB (Fig. S2) drives ice retreat in the 10 and 20 km models. Because the nonuniform mesh captures a trough along flow line E, delivery of ice mass to the margin continues through the Holocene, stabilizing the ice margin position despite surface lowering through negative SMB.

4 Discussion

We find that model resolution plays a negligible role in the simulated ice margin history in the northern portion of our domain along flow lines A, B, and C. In the southern domain along flow lines D and E, however, model resolution plays a large role in the simulated ice margin history.

4.1 Retreat within the northern domain

In the northern portion of our model domain, geologic archives indicate retreat from the coast occurred between 12 and 11 ka (Kelly et al., 2015; Young and Briner, 2015; Fig. 6a), which is generally consistent with our simulations regardless of resolution. The subsequent retreat in all models towards the present-day ice margin is also generally consistent with the geologic reconstructions of ice margin retreat in this region (Fig. 6b and c). Van Tatenhove et al. (1996) provide one of the earliest ice sheet model–data comparisons for this region (around flow line C) during the last deglaciation. Van Tatenhove et al. (1996) compared three different ice sheet models ranging in resolution from 20 to 40 km to ice margin reconstructions constrained by radiocarbon ages from the region and indicated that model resolution played little role in the inter-model retreat differences. Van Tatenhove et al. (1996) pointed to the strong governing role of SMB in this region with little influence from ice streams. Likewise, simulations from the 20 km resolution model of Lecavalier et al. (2014) show reasonable agreement in the retreat across this region when compared to geologic reconstructions. Our results indicate that the bed topography in this region is well represented among the different models, despite their differences in resolution, and thus simulated ice margin history faithfully responds to SMB forcing and is not complicated by ice flow adjustments to underlying topography. However, one feature that stands out in the higher-resolution models (5 km and nonuniform) is the presence of high-frequency ice marginal fluctuations on the order of 5–8 km (Figs. 7, S3, S4). The geologic record indicates that small-scale marginal fluctuations are likely responsible for the moraine record and seem to be related to high-frequency variability in temperature (e.g., Young et al., 2013). Thus, models capable of capturing small-scale fluctuations in the ice margin history are valuable for comparing with geologic constraints of past ice sheet change. The inability of lower-resolution models to capture these features has been highlighted in previous work (Van Tatenhove et al., 1996; Larsen et al., 2015) and hampers data–model comparisons.

4.2 Retreat within the southern domain

In the southern portion of the model domain, fjord systems provide a different bed setting than in the north, presenting a significant challenge for modeling ice margin change. Deep and narrow troughs up to 500 m below sea level and 3–5 km wide seemingly played an important role in governing ice margin retreat. Many geologic archives that constrain past ice margin variability in this region (Sinclair et al., 2016) reveal rapid deglaciation from the present-day coastline to near the present-day margin at ∼10 ka (e.g., Larsen et al., 2014; Fig. 6a). None of our experiments match the geologic observations.

Our simulated ice retreat is highly dependent on model resolution in this area because the different models represent the bed topography quite differently. For example, only the nonuniform and 5 km models capture the deep fjords, whereas the 10 and 20 km models have unrealistic bed features that end up driving retreat (Figs. 8 and S5). Our simulations do not include calving or submarine melting, and, therefore, each model's simulated ice surface responds similarly to negative SMB. However, the ability of the high-resolution models to resolve the narrow and deep fjords allows the ice margin to persist, as the stress balance and mass transport is well captured. Since the fjords are not well represented in the low-resolution models, there is lower delivery of ice mass to the margin, and the simulated retreat is driven as the ice surface lowers and intersects elevated bump artifacts in the bed topography. While none of the models capture the timing or amount of retreat accurately, the high-resolution models in this case perform the worst, capturing negligible retreat. The rapid ice margin recession recorded by the geologic reconstructions in this marine-dominated region probably highlight the influence of calving and enhanced submarine melting of floating ice, neither of which are included in our model simulations. The lack of submarine melting, in particular, may lead to the model–data mismatch; available evidence (e.g., Dyke et al., 2014) supports influence of the warm Irminger Current during the early Holocene, which likely penetrated fjords up to the ice margin.

4.3 Different drivers of retreat

There are stark differences in processes affecting retreat of the land-dominated ice margin (i.e., SMB in the northern section of the model domain) and the marine-dominated ice margin (i.e., calving and submarine melt in the southern portion of the domain). These different drivers of ice margin change also affect different sectors of the contemporary GrIS (Sole et al., 2008; Straneo and Heimbach, 2013). Our results highlight that in areas of simple, low-relief bed topography, SMB drives the simulated retreat with few differences existing between models of varying spatial resolution. Therefore, efforts that attempt to match geologic reconstructions will be better served by focusing on representing SMB as accurately as possible. Conversely, in areas with complex, high-relief bed topography, such as in fjord settings, models that are unable to capture the deep and narrow troughs may unreasonably simulate retreat (see Åkesson et al., 2018). Anisotropic mesh capabilities play an important role in allowing a model to adjust its resolution spatially while using computer time efficiently.

For low-resolution models, care must be taken when attempting to capture the reconstructed retreat in areas of complex bed topography. For example, in order to satisfy relative sea level records used to constrain an ice sheet model of the GrIS, Lecavalier at al. (2014) artificially increased middle Holocene temperatures used to drive the ice sheet model. Although this resulted in a simulated ice margin history consistent with available geologic records, it is noted that such external forcings may drive unphysical retreat in areas of complex bed topography that may otherwise have been driven by ice dynamics. Another consideration is the regional setting presented here. Since these experiments are focused on the southwestern GrIS, an area that may be relatively topographically uniform, we expect the results of the marine-influenced region in the southern part of our domain to be most relevant for other portions of the GrIS, in particular in eastern Greenland where fjords dominate the landscape (Morlighem et al., 2017).

4.4 Model limitations

When simulating the retreat of the southwestern GrIS, our choice of climate forcing, using the GRIP δ¹⁸O record (Dansgaard et al., 1993), follows what has been a cornerstone in forcing Greenland ice sheet modeling over the paleoclimate record (Huybrechts, 2002; Greve et al., 2011; Applegate et al., 2012). This approach has been adjusted in Tarasov and Peltier (2002), Simpson et al. (2009) and Lecavalier et al. (2014) by synthetically increasing Holocene temperatures, with more recent simulations of the deglaciation of Greenland making use of more recent temperature proxy reconstructions that are better constrained throughout Greenland (Lecavalier et al., 2017; Buizert et al., 2018). Nevertheless, using a single, scaled paleoclimate record from Summit ignores the more likely history of a spatiotemporally variable climate history spanning the Holocene around Greenland (cf. Vinther et al., 2009). In any case, since the traditional approach (i.e., the GRIP δ¹⁸O scaling) assumes that the spatial variability in temperature and seasonality remains fixed to modern day, our results cannot fully reconcile how changes in the magnitude of warming and spatial variation of that warming affects our results. Additionally, the scaling of the basal friction coefficient introduces some uncertainties particularly when considering the temperature forcing throughout the simulation. Our method for scaling the basal friction coefficient through time follows a common approach used in many modeling studies over paleoclimate timescales (Hindmarsh and LeMeur, 2001; Greve ,2005). For these simulations, the evolution of the basal temperatures through time depends on the surface temperature forcing which is derived from the GRIP δ¹⁸O scaling. Therefore, changes to the surface temperature forcing can impact the evolution of the basal temperatures over time, which ultimately affects the ice sliding following this approach. This model limitation falls within the bounds of current ice sheet modeling efforts, whereby a lack of physically based basal sliding parameterizations exist. Despite this limitation, the conclusions presented here remain unaffected.

Although these simulations have no reasonable representation of calving, the results do indicate that models with a resolution of 10 km or greater would be likely unable to address calving processes in fjords, as typical fjord width is ≤10 km. In ice sheet models, calving is often related to water depth, considering past changes in eustatic and relative sea level (Huybrechts, 2002; Simpson et al., 2009; Lecavalier et al., 2014). Although the high-resolution models presented here do capture the narrow fjords, implementation of a calving scheme would currently be computationally intensive. One possibility for future work would be to force the model with high submarine basal melt rates as a proxy for calving, as done in Åkesson et al. (2018). Submarine melt has been shown to be an important mechanism driving both contemporary ice mass loss (Rignot, 2010) and past GrIS variability on glacial and interglacial timescales (Bradley et al., 2018; Tabone et al., 2018). Although few constraints do exist detailing past variations in ocean temperature for the Labrador Sea (Winsor et al., 2012; Gibb et al., 2015) and Disko Bay (Jennings et al., 2006), applying submarine melt rates to marine termini throughout our model domain would not be possible without significant uncertainty. Additionally, in applying basal melting to floating ice, it is uncertain whether 2 km resolution would be sufficient to accurately capture grounding line migration, as recent research (Seroussi and Morlighem, 2018) suggests that resolutions of 1 km or higher are often necessary to match present-day fluctuations.

5 Conclusions

We investigated how ice sheet model resolution influences the simulated Holocene retreat of the southwestern GrIS using the ISSM. Our focus on the southwestern GrIS is driven by two factors: first, the regional approach allowed for modeling at a high resolution (for the nonuniform and 5 km mesh) while minimizing computational costs that would increase significantly while running a GrIS-wide simulation. Second, the southwestern GrIS is an area where geologic archives indicate the ice sheet underwent large-scale and relatively well-known retreat during the Holocene.

The results presented here indicate that model resolution has a selective influence in the simulated retreat over southwestern Greenland during the Holocene. In areas where the bed topography is relatively simple, low-relief, and free from marine influence, model resolution plays an insignificant role in influencing the pattern and rate of retreat. Here, models with different resolutions respond similarly to SMB-driven retreat. On the other hand, in areas with complex and high-relief bed topography, such as deep troughs and fjords, the low-resolution models lead to unrealistic retreat. As all models in these simulations only respond to surface melt and, therefore, ice surface lowering (and no mass loss via calving or submarine melt), the low-resolution models (10 and 20 km) simulate ice retreat driven purely as a consequence of incorrectly capturing the bed geometry. As one example, ice surface lowering in these models intersects bed bumps that would otherwise be resolved as a trough in higher-resolution models.

Our results imply that computational resources can be saved when modeling certain portions of the GrIS. Conversely, the results also highlight the importance of model resolution in areas of complex topography. Ice sheet models using a nonuniform mesh can adapt grids to fit these constraints while using computation time efficiently. However, for models using uniform fine-scale meshing, resolving such features becomes computationally difficult, especially over long paleoclimate timescales. As ice sheet models sometimes rely on the geologic record for validation, care must be taken in evaluating model–data misfits. In areas of complex topography, over-tuning of model parameters or climatology may occur in low-resolution models that seek to match the reconstructed margin. We suggest that increased model resolution is critical in regions dominated by fjords (e.g., southeastern Greenland).

Future work with the ISSM will focus on using a model that has a lower resolution in areas driven mainly by SMB and a higher resolution in areas influenced by dynamical ice processes using nonuniform mesh capabilities. Future work will also seek to evaluate the sensitivity of using improved climate forcings (Hakim et al., 2016; Buizert et al., 2018), better representations of ice dynamics (calving and submarine melt), and more quantitative comparisons to improved ice margin reconstructions.