ICESHEET 1 . 0 : A program to produce paleo-ice sheet reconstructions with minimal assumptions

We describe a program that produces paleo-ice sheet reconstructions using an assumption of steady state, perfectly plastic ice flow behaviour. It incorporates three input parameters: ice margin, basal shear stress and basal topography. Though it is unlikely that paleo-ice sheets were ever in complete steady-state conditions, this method can produce an ice sheet without relying on complicated and unconstrained parameters such as climate and ice dynamics. This makes it advantageous 5 to use in glacial-isostatic adjustment ice sheet modelling, which are often used as input parameters in global climate modelling simulations. We test this program by applying it to the modern Greenland Ice Sheet and Last Glacial Maximum Barents Sea ice sheet and demonstrate the optimal parameters that balance computational time and accuracy.


Introduction
Reconstructing past ice sheets is a complex task, due to the large number of parameters that can affect their growth and retreat.For example, Tarasov et al. (2012) presented a glacial systems model that contained 39 parameters that could be tuned, which included climatology, Earth rheology, ice physics and margin chronology.Many of these parameters are poorly constrained by available observations.In particular, past climate is often parameterized based on ice core data from Greenland and Antarctica, or reconstructions from speleothems that are located far from where the ice sheets existed.
Since past climatic parameters are generally only well characterized in areas outside of where paleo-ice sheets existed, ice sheet reconstructions that are independently determined using evidence of glacial-isostatic adjustment (GIA) are often used in paleo-climate simulations (e.g.Braconnot et al., 2007Braconnot et al., , 2012)).One of the most commonly used GIA based reconstructions of glaciation is the ICE-xG series (e.g.Peltier, 2004;Peltier et al., 2015).They produce configurations of ice sheets that minimize the misfits of geodetic and relative sea level data, with limited regard to the physical realism of the ice sheet itself.Another commonly used reconstruction is the ANU model (e.g.Lambeck et al., 2010), which was developed using an assumed peak ice elevation at the centre of ice sheets, and using a parabolic ice profile to the margins.In their formulation, each flowline ray is allowed to have different basal shear stress values, but is less flexible in regards to the direction of the flowline, and spatial variability in basal shear stress along it.
The program presented in this paper produces a physically realistic ice sheet reconstructions taking into account changes in basal shear stress and topography, while being simple enough that it does not depend on numerous parameters with large uncertainties.The goal of this program is to provide an compromise between the GIA-only ice sheet re-constructions that have limited or no physics applied to their construction, and the full glacial systems models that demand considerable computational resources.The reconstructions are based on the assumption of perfectly plastic, steady-state ice conditions.It allows for the rapid determination of paleoice sheet configurations, which is desirable when matching observations of GIA.We present an example application of this program to the Barents Sea Ice Sheet, a relatively shortlived portion of the Eurasian Ice Sheet complex, by trying to match an existing GIA based model.We also apply the model to the contemporary Greenland ice sheet to provide an indication of how well the model is capable of reconstructing a known ice sheet geometry.Ultimately, the goal would be to reconstruct, in a timestepped fashion, the entire history of an ice sheet complex.In this case, the basal topography is relatively well determined (since there is no existing ice), and the basal shear stress can be established to a certain extent by the surficial geology and geomorphology.The ice topography and basal shear stress are determined through time using external evidence, such as the nature of GIA.An example of this is presented for the western Laurentide Ice Sheet by Gowan et al. (2016).

Theory
The reconstructions produced by the ICESHEET program are based on the assumption that ice rheology adheres to perfectly plastic, steady-state conditions (i.e.ignoring lateral shear stresses, and assuming that the ice surface is not dynamically changing).The two-dimensional form of this theory was derived by Nye (1952), and neglects variability in topography and longitudinal changes in stress.In this equation, the ice surface gradient is directly related to the strength of the ice-bed interface, or basal shear stress.The basal shear stress is related to a number of factors, including basal geology, sediment thickness and strength, hydrology, temperature and bed roughness.
The ice surface elevation is E, s is the distance along ice flowline profile, τ o is the shear stress at the base of the ice sheet, which balances the driving stress, ρ i is the density of ice, g is the gravity at the Earth's surface, and H is the ice thickness.If the distance from the ice sheet margin to the centre of the ice sheet is known, then the thickness along the profile between the two points can be calculated using the following formula (Cuffey and Paterson, 2010).
In this equation, L is the distance between the margin and centre of the ice sheet, and x is the distance from the centre.Though this equation is simple, it can be used to make a rough estimate of the thickness of ice sheets, neglecting basal topography (Cuffey and Paterson, 2010).Equation ( 2) was used to create the ANU ice sheet reconstructions (i.e.Lambeck et al., 1998Lambeck et al., , 2006Lambeck et al., , 2010)).The weakness of using this equation is that the centre of the ice sheet has to be assumed a priori.It also does not take into account changing basal shear stress conditions or changes in topography.
In order to overcome problems with spatial changes in basal topography and shear stress, in addition to the uncertainties in the location of the ice sheet centre, Reeh (1982) and Fisher et al. (1985) presented expanded version of Eq. ( 1) that allows for changes in the direction of the flowline.The equation becomes the following partial differential equation.
The coordinate system is set up so that x points towards the centre of the ice sheet, and y is parallel to the margin.Presented in the notation used by Reeh (1982), Eq. ( 1) is substituted into the left side of Eq. ( 3) with the ice thickness represented in terms of ice surface elevation and basal topography elevation B, and substituting in a characteristic thickness, The above equation describes the change in ice thickness over an arbitrary surface.This partial differential equation can be solved by the method of characteristics (Kamke, 1965).The x and y partial derivatives in Eq. ( 4) are substituted by p = ∂E/∂x and q = ∂E/∂y, then rearranged in terms of p.
The solution to the partial differential equation then becomes three ordinary differential equations that are solved simultaneously, using the method of characteristics (Reeh, 1982).
Equation ( 6) gives the direction of local maximum steepness, while the other two equations describe how the elevation changes spatially in the x direction.Fisher et al. (1985) expanded Eq. ( 8) to allow for changes in basal shear stress (in terms of the characteristic thickness, H f ).
These equations are solved by numerical integration to determine the course and gradient of an ice flowline.In the next subsection, we note some of the improvements to the original methodology, including adjustments to the base topography with realistic GIA, dealing with margins that are in marine environments, automatic determination of ice sheet saddles, and adjusting for the presence of nunataks.
It is important to note that assuming perfectly plastic, steady-state conditions for the ice sheet is not accurate in areas where the ice sheet was highly dynamic, or where lateral shear stress was an important factor.Due to this, the output basal shear stress is unlikely to reflect the true basal shear stress in those areas.

Algorithm to reconstruct ice sheets
In order to solve Eqs. ( 6)-( 8), initial values for E, y and q are required.Starting the calculation at the margin is convenient from the perspective of reducing a priori assumptions on ice distribution, though it leads to a singularity because the ice thickness is zero (E = B).Consequently, the value of E at the margin must be set to be a nominal value (in the sample problems presented in this study, 1 m).Although the actual thickness of ice near the margin may be as high as tens of metres, the choice of starting value will not have a large effect on the final model.For instance, the distance from the margin required in Eq. ( 2) to reach 10 m from a starting value of 1 m, and a low basal shear stress value (5 kPa) is 90 m, substantially smaller than the uncertainty in the margin location for paleo-ice sheets (Clark et al., 2012;Gowan, 2013;Hughes et al., 2016).For simplicity, the value of q is defined to be zero at the margin.This can be justified because near the margin the value of term H f /(E − B) will dominate Eq. ( 5) in the defined coordinate system.
The ice sheet reconstruction is calculated in a piece-wise manner (see Fig. 1 for an illustration of the steps involved).The ice flowline calculation is initiated at intervals along the margin, which are user defined.The flowline calculation proceeds until it reaches a particular elevation (a user defined contour interval), at which point the program checks to see if any flowlines cross over, or if a saddle point in the ice sheet has been reached.A sequential list of the modelling steps is given below.
1.All parameters (ice sheet margin, shear stress map, topography map) are converted from geographical coordinates to a Cartesian coordinate system prior to the execution of the program.
2. Estimates of the basal shear stress for the area of interest are read into the program.The shear stress values must be adjusted for each time epoch to produce an appropriate ice sheet configuration.
3. The basal topography data for the area of interest are read in.For the first iteration of ice sheet model development, it uses modern topography or topography adjusted for changes in global mean sea level (in practice, it has limited impact on the final reconstruction, i.e. < 100 m near the edge of the ice sheet and much less than that in the interior, even with predominantly marine based ice sheets).In subsequent iterations, the topography is adjusted for glacial-isostatic adjustment, to take into account the fact that the ice sheet will deform the Earth, and that the ice sheets will cause changes to sea level.The modified topography is calculated before running the ice sheet program.In the Barents Sea Ice Sheet sample problem, we use the CALSEA program to calculate GIA (Nakada and Lambeck, 1987;Lambeck et al., 2003).CALSEA computes glacial-isostatic adjustment using a spherically symmetric Earth, with a Maxwell rheology mantle and elastic lithosphere, using the PREM model (Dziewonski and Anderson, 1981) for other Earth model parameters.In includes time evolving shorelines and rotational feedback.
4. The program reads in the margin, and defines locations along the perimeter where the flowline calculation initiates.The minimum distance along the margin between where flowline calculation is initiated is user-defined.
The program defines the initial direction of flow to be perpendicular to the margin, away from the centre of the ice sheet.
5. The margin is set to have an initial ice thickness of 1 m.If the margin is located where the topography is below sea level, it is assumed that the margin corresponds to the grounding line of the ice sheet.A conservative estimate of the thickness of ice at this point is set to , where ρ seawater is the density of sea water and ρ ice is the density of ice, which is the thickness of ice corresponding to the equivalent mass of the water column at that point.There is a check to make sure that the ice surface slope between adjacent points on the boundary is not too steep for the given Step 7 Step up from 10 to 20 m x Step 11

Eliminate crossovers
Step basal shear stress values.If it is, the ice thickness at the point with the lower elevation is increased.This check is only done where B < 0.
6.The calculation of ice elevation contours is a recursive process.If the contour crosses over itself (signifying a saddle on the surface of the ice sheet), the contour polygon is split, and the calculation is continued as separate polygons (see step 12).
7. The program searches for points on the contour that are below the next contour elevation.The elevation may be above the next contour elevation along the margin, or if a point coincides with a nunatuk (see step 14).
It then calculates the flowline by numerical integration of Eqs. ( 6)-( 8), using the Runge-Kutta method (Press, 1992).When it reaches the next contour elevation, the calculation stops.
8. If the flowline calculation cannot reach the next contour elevation, which happens when the topography is too high (H → 0, or E < B), the point is flagged and not included in the next contour (Fig. 1).
9. If the flowline direction changes sufficiently so that q ≥ H f /(E − B) (i.e.p approaches zero), the local coordinate system is rotated so that p is in the direction of maximum flow.
10.If the calculated flowline goes outside the last calculated contour polygon, it is flagged and the point is not included in the next contour.This happens when the ice surface is near its peak height.This can also happen in areas where there is a sudden change in topography or basal shear stress, which causes a deflection in the flowline direction (Fig. 1).
11.After the flowlines are calculated for each applicable point along the polygon, the program checks to see if any of the calculated flowlines cross over.Offending crossovers are eliminated using a motorcycle algorithm (e.g.Vigneron and Yan, 2014).The eliminated flowlines are flagged and not included in the next contour (Fig. 1).
12. At this point, an initial polygon of the next elevation contour can be constructed.This is checked to ensure that it is a simple polygon (i.e. a polygon that does not cross over itself).If it is not, then the program breaks it into several polygons, and determines whether they represent domes (ice gradient is increasing towards the centre of the polygon) or saddles (the ice gradient is decreasing towards the centre of the polygon).Where a saddle is identified, it is determined to have reached its peak elevation and is eliminated from subsequent calculations (Fig. 1).
13.The ice elevation and thickness for all points on a valid polygon (including flagged points) are written to file.
14.The polygon is resampled using the user-defined distance interval.There is also a check using Eq. ( 2) to estimate the distance to the next contour.If the difference in estimated distance between adjacent points is greater than the user defined distance threshold, additional points are included.This process excludes flagged points, and may incorporate basal topographic highs, where flowline calculation will not be initiated (Fig. 1).
This process is repeated for each time interval of interest.After calculation of the ice reconstruction, the calculated elevation values are averaged into a grid to be used as input for a GIA calculation program.The grid is created using a continuous smoothing algorithm, which is part of Generic Mapping Tools (Smith and Wessel, 1990).
3 Sample reconstruction -Barents Sea Ice Sheet

Setup
The Barents Sea Ice Sheet was predominantly marine based, and likely formed by the merging of isolated ice caps over Svalbard, Franz Josef Land, Novaya Zemlya and the Scandinavia Ice Sheet (Ingólfsson and Landvik, 2013).The hypothesis explaining the glaciation of the entire Barents Sea is that GIA warped the floor of the Barents Sea upwards, favouring the formation of grounded ice.At the Last Glacial Maximum (LGM) (about 20 ka), the ice sheet covered the entire continental shelf region west of Novaya Zemlya (Ingólfsson and Landvik, 2013).The extent was probably limited in the Kara Sea east of Novaya Zemlya, compared to the mid-Weichselian (45-55 ka) glaciation.At the LGM, the ice thickness was likely greatest to the east of Svalbard, on the basis of the pattern of paleo-sea level reconstructions (Lambeck, 1995).
In this sample problem, the ice sheet extent is taken as the "most likely" configuration at 20 ka from the DATED project (Hughes et al., 2016).Since the Barents Sea Ice Sheet merged with the Scandinavian Ice Sheet at the LGM, the margin is cut off far enough south so that the northern part of the Scandinavian Ice Sheet is sufficiently represented.The basal topography used in this problem is from IBCAO (Jakobsson et al., 2012).The basal topography of Svalbard takes into account the thickness of modern ice cover.There is no published information on the thickness of ice on Novaya Zemlya, so we use contemporary ice surface topography.The basal shear stress was initially parameterized on the basis of topography and bedrock geology.The values were adjusted in order to produce an ice thickness distribution that is similar to the GIA based ANU model (Lambeck, 1995;Lambeck et al., 2006Lambeck et al., , 2010)).Exact matching of ice thickness in the sample problem to the ANU model was not attempted, since it is of low resolution, and has a different margin configuration to that of Hughes et al. (2016).Specifically, it is less extensive along the Bear Island Trough.In order to approximate the ice thickness from the ANU model, the basal shear stress was set to be high along the northern part of the ice sheet, and relatively low in the southern Barents Sea.Both the topography and basal shear stress values are sampled at 5 km (Fig. 2).This purpose of this test is to demonstrate that GIA has an impact on the ice sheet reconstruction.This test only in- cludes the Barents Sea Ice Sheet for the calculation of GIA.
In a full glacial reconstruction (e.g.Gowan et al., 2016), it is necessary to include the effects of far field ice sheets, and realistic ice sheet growth and decay.

Resolution test
In order to test the optimal parameters for producing ice sheet reconstructions, a series of tests with different distance and contour intervals were performed, the results can be found in Table 1.This test involved using modern topography minus the approximate 133 m reduction in global mean sea level at 20 ka (Fig. 2, Lambeck et al., 2014).The shear stress and basal topography values are shown in Fig. 2. Figure 3 shows how changes in basal shear stress and basal topography affect the modelled ice sheet.The spacing between contours is greater in areas of low basal topography and shear stress, which replicates ice flow from areas of high to low basal topography, and around barriers that resist ice flow.The program execution time largely depends on the chosen sampling interval along the contour polygons (Table .1).The reference ice sheet configuration used a distance interval of 1 km, and a contour interval of 10 m (Fig. 4).Unsurprisingly, considering the 5 km resolution grid, all tests using distance intervals 5 km or less produced nearly identical reconstructions, as they captured the details of the grids.Using a contour interval of 20 m gives almost the same result as as 10 m, with diminishing accuracy when increased above this, without significant reductions in execution time.The optimal parameters for matching the reference configuration and fast execution time are a 5 km spacing and 20 m contour interval (Table 1).Increasing the distance parameter decreases the execution time, but is unable to match the reference reconstruction, particularly in the mountainous regions of Svalbard and Scandinavia.There is a tendency towards overestimating the ice thickness when the initiation distance is larger than 5 km (Fig. 4).During the initial phases of GIA based ice model development, it may be prudent to decrease the resolution of the grids to quickly determine an estimate of basal shear stress, then increase the resolution when refinement is necessary.

GIA test
When an ice sheet grows, the basal topography is modified by GIA, which will significantly impact the Barents Sea Ice Sheet example.Therefore, in order to obtain an accu-rate characterization of the ice sheet surface topography and thickness, it is necessary to re-run the program with the modified basal topography.The Earth model used in this sample problem is spherically symmetric and includes a 90 km thick elastic lithosphere, 4 × 10 20 Pa s upper mantle viscosity and 10 22 Pa s lower mantle, which is in the range of best fitting models for this region (Lambeck et al., 2010).The distance interval used is 5 km and contour interval is 20 m.Since there is a viscous component of the response, the ice sheet is allowed to grow linearly from 30 ka (when glaciation in the Barents Sea is presumed to be similar to present, Mangerud et al., 1998) to 20 ka, then linearly decrease back to present levels at 10 ka.After the first iteration of GIA, the ice sheet contribution to global mean sea level is subtracted to determine the Earth deformation.When combined with the actual global mean sea level at this time (−133 m), it should give a reasonable estimate of local basal topography.
The results show that one iteration of GIA has a significant effect on ice sheet reconstruction, and in this case increases the total volume by about 5.8 % (Fig. 5).In addition, since the basal topography becomes more depressed towards the centre of the ice sheet relative to the initial reconstruction, the reconstructed ice surface topography is lower and has a more gentle gradient.A second iteration of GIA had only a minor effect on the reconstructed ice sheet (0.4 % increase in volume from the first iteration).
Additional tests by Gowan (2014) for the full deglacial Laurentide Ice Sheet showed that there is only a weak dependence on reconstructed ice volume and the Earth model used to compute GIA.For three layer (lithosphere, upper mantle, lower mantle) Earth models, the ice volume varied most with changes in lower mantle viscosity at LGM extent, but the difference was less than 0.5 % (though smaller ice sheets will have less dependence on the lower mantle).Towards the end of deglaciation, there was more dependence on upper mantle viscosity, but again, the volume difference was less than 0.5 %.Though the volume was close to the same, there were slight differences in the distribution of ice, though not by more than 100 m in extreme cases.Therefore, the recommendation when creating an ice sheet model is to include at least one iteration of GIA, but the chosen Earth model is not as important.continuity based inversion of the contemporary ice thickness (Morlighem et al., 2014).Reeh (1982) reconstructed the Greenland Ice Sheet reasonably well using the methodology explained earlier using a constant basal shear stress of 90 kPa.Since ICESHEET can have spatially variable basal shear stress and account for variable topography, it is possible to refine this.Advances in remote sensing over the last 30 years also allow a more accurate comparison to contemporary topography.The goal of this example is to determine the misfit between the ICESHEET reconstructed ice surface topography and the contemporary ice sheet using a methodology analogous to the reconstruction of a paleo-ice sheet.The input grounded ice margin and basal topography data come from the IceBridge BedMachine Greenland, Version 2 data set (Morlighem et al., 2014(Morlighem et al., , 2015)).The basal shear stress value domains were designed the same way as a paleo-ice sheet would be constructed.The domains were constructed purely on the basis of basal topography (Fig. 6), since information on basal geology is limited.They were predominantly divided into areas of rugged topography (i.e.mountainous regions), flat-lying areas, and fjords.There intentionally was no attempt to divide it on the basis of modern ice flow patterns, given that it may not be possible to deduce them for a paleo-ice sheet.The shear stress values in the domains were adjusted iteratively in order to try to match the observed ice surface topography.In a paleo-ice sheet, it will not be possible to know what the ice surface topography was a priori.In that case, other sources of data (i.e.GIA) must be used as the basis for the reconstruction.

Results
The resulting reconstruction is shown in Fig. 6.For comparison purposes, the ice sheet is averaged into a 25 km grid.The reconstructed ice sheet surface topography has an average difference of −37 ± 2 m (within 200 m of the true topography for most of the ice sheet).The largest errors (> 400 m) occur in places where there are narrow ice streams near the edge of the ice sheet, which could not be parameterized using the coarse resolution shear stress domains.In general, the shear stress values are highest in the mountainous regions in southeastern Greenland.The basal shear stress is lowest in the centre of the ice sheet, likely reflecting the flat-lying basal topography.Direct inversions for basal shear stress have only been performed for some of the ice streams in (e.g.Sergienko et al., 2014;Shapero et al., 2016).In the study by Sergienko et al. (2014), the basal shear stress exhibited a banded pattern, alternating between low (< 50 kPa) to high (> 150 kPa) values over spatial ranges of 5-20 km.Shapero et al. (2016) found that the basal shear stress directly under fast flowing ice streams was almost negligible, but at the sides it could exceed 375 kPa.If averaged over a larger area, these values are consistent with the 100-200 kPa values in our reconstruction (Fig. 6).The resolution test was also performed with the Greenland simulation (Table 1).In this sample, the 5 km distance interval, 20 m contour interval does not perform quite as well as in the Barents Sea example.This is a result of having a larger area of mountainous terrain.Still, less than 3 % of the elements are greater than 100 m different from the reference reconstruction.If the area of focus is predominantly mountainous, it may be prudent to decrease the distance interval.The volume of the Greenland Ice Sheet, taken directly from the data set by Morlighem et al. (2014)

Figure 2 .
Figure 2. Basal topography used in the resolution test, which is modern topography minus the 133 m drop in global mean sea level at 20 ka.Also shown in brown is the 20 ka ice margin (Hughes et al., 2016) and the location of places described in the text.Sv -Svalbard.FJL -Franz Josef Land.NZ -Novaya Zemlya.BRT -Bear Island Trough.(b) Basal shear stress values used in the example in this paper.

Figure 3 .
Figure 3. Example from central-western Svalbard of how spatial changes in basal topography and basal shear stress affect the reconstructed ice surface topography of the ice sheet (see text).The contour interval is 100 m in the figure, though this sample was calculated with a 5 km spacing and 20 m contour interval.The dark black lines are the modern-day coastlines.

Figure 5 .
Figure 4. (a) Reference ice sheet reconstruction with a 1 km spacing and 10 m contour interval using the topography and basal shear stress in Fig. 2. (b) The difference between a model calculated with a 20 km spacing and 20 m contour interval and the reference reconstruction shown in (a).This demonstrates that the lower resolution tends to overestimate the ice surface elevation in mountainous regions.

Figure 6 .
Figure 6.Sample reconstruction of the contemporary Greenland Ice Sheet.(a) Modern topography.The brown line is the current grounded ice margin, the black lines are the modern day coastlines.(b) Reconstructed topography.(c) Final iterated basal shear stress domains and values used for the reconstruction.(d) Difference between the observed topography and reconstructed topography.

Geosci. Model Dev., 9, 1673-1682, 2016 www.geosci-model-dev.net/9/1673/2016/
12 Schematic illustrating the steps in calculating the ice sheet, illustrating steps 7, 8, 10, 11, 12 and 14 in Sect.2.2.The black lines indicate the initial elevation contour, blue lines indicate calculated flowlines, red lines indicate the next elevation contour, black circles indicate flowline initiation points, unfilled circles indicate added initiation points for the next elevation contour, crosses indicate flagged points that are not included in the next elevation contour.

Table 1 .
(Sambridge et al., 2009)n test Execution time on Terrawulf III(Sambridge et al., 2009), Dual Intel Xeon X5650 at 2.66 GHz running OpenSuse 13.2.Compiled with ifort 15 with -O2 flag.b Percent of 0.5 • longitude by 0.25 • elements that are > 100 m different from the reference model (out of 23 205 total elements for the Barents Sea, and 21 901 for Greenland).
(Gowan et al., 2016) km 3 .From Table1, the reconstructed volume is within 5 % of this value, except in the lowest resolution tests.ICESHEET 1.0 is a program that can quickly create reconstructions of paleo-ice sheets, with a given margin configuration and estimated basal shear stress.We have provided two proof of concept examples showing reconstructions of the modern Greenland Ice Sheet and the Barents Sea Ice Sheet at the LGM.It is recommended that at least one iteration of GIA is included to best characterize the thickness and ice surface topography.It is also recommended (if a 5 km basal topography grid is used) to use a flowline spacing interval of 5 km and contour interval of 20 m for optimal calculation speed.This program has been used to create a full late glacial GIA based ice sheet reconstruction of the western Laurentide ice sheet(Gowan et al., 2016).It is ideal for producing ice sheet reconstructions that have minimal input assumptions, but are glaciologically plausible.A suite of ice sheet reconstructions through a glacial cycle could be used as independent inputs for climate and ice sheet dynamics modelling.