Supplementary material for "Numerical issues in modeling ice sheet instabilities such as binge-purge type cyclic ice stream surging"

Running the models:

GSM: A full description of the GSM with the associated code archive and complete documentation will soon be submitted and then linked here. We are happy to provide the code upon individual request. The input files used specifically for "Numerical issues in modeling ice sheet instabilities such as binge-purge type cyclic ice stream surging" are described here. The GSM uses bed masks to set the sediment cover of each grid cell. The bed masks for all four resolutions (25 km, 12.5 km, 6.25 km, and 3.125 km) with an abrupt transition between soft sediment and hard bedrock are:

inputGSM01.bedMskLISsqHighRes25kmSymmetric.dat
inputGSM02.bedMskLISsqHighRes12p5kmSymmetric.dat
inputGSM03.bedMskLISsqHighRes6p25kmSymmetric.dat
inputGSM04.bedMskLISsqHighRes3p125kmSymmetric.dat

The GSM bed masks with a small (3.125 km) and wide (25 km) transition zone between soft sediment and hard bedrock are:

inputGSM05.bedMskLISsqHighRes3p125kmSymmetricSmallSmoothTransition.dat
inputGSM06.bedMskLISsqHighRes3p125kmSymmetricSmoothTransition.dat

Similar to the sediment cover, the GSM uses a topography mask to set the bed elevation. The topography files with a small (3.125 km) and wide (25 km) transition zone between the pseudo-Hudson Strait and Hudson Bay at 200 m below sea level, the surrounding at sea level, and the ocean grid cells at 500 m below sea level are:

inputGSM07.topo3p125kmSymmetricSmallSmooth.dat
inputGSM08.topo3p125kmSymmetricSmooth.dat

The GSM setup used here requires 9 input parameters. The parameters are read as parameter vectors from a parameter file. The GSM parameter vectors used in this paper are:

1.00 0.30 0.81 1.20 1.30 2.60 1.44 1.20 0.30 (parameter vector 00)
1.00 2.00 0.80 1.20 0.75 1.07 1.00 0.50 0.30 (parameter vector 01)
0.59 1.20 1.00 1.20 0.52 2.87 1.27 0.82 0.20 (parameter vector 02)
1.00 2.00 1.00 0.54 0.27 1.52 1.17 0.95 0.20 (parameter vector 03)
1.00 1.50 1.00 1.20 0.60 1.52 0.97 0.60 0.10 (parameter vector 04)
C_rmu, C_fslid, lapsr, PDDmelt, hpre, PrecRef, rTnorth, n_b

PISM: Instructions on how to install and run PISM can be found in the PISM online manual. The experiments in "Numerical issues in modeling ice sheet instabilities such as binge-purge type cyclic ice stream surging" were conducted with PISM v2.0.2. On a Unix/Linux-based system, the PISM base setup can be run by

mpiexec -n 8 pismr -Mx 120 -My 120 -Mz 60 -Lz 1e4 -y 2e5 \
        -Mbz 20 -Lbz 1000 \
        -bootstrap -i inputFile.nc -o outputfFile.nc \
        -stress_balance ssa+sia -stress_balance.sia.max_diffusivity 1e3 \
        -adapt_ratio 0.01 \
        -grid.registration center \
        -hydrology null \
        -extra_file extraOutputFile.nc -extra_vars thk,temppabase,velsurf_mag,velbar_mag,flux_mag,diffusivity,bmelt,taud_mag,tauc,tillphi,velbase_mag -extra_times \0:1000:2e5 \
        -ts_file extraTimeSeriesFile.nc -ts_times 0:100:2e5"

inputPISM01.pismInputFilesDist.py uses a parameter file (containing the parameter vectors) to create the inputFile.nc required in the command above. The PISM parameter vectors used in this paper are:

     0.54     19.81    205.25  7.68e-12    236.14  9.01e-09 (parameter vector 00)
     0.56     19.44    408.81  4.55e-12    232.60  9.45e-09 (parameter vector 01)
     0.60     21.13    237.37  2.65e-12    221.58  3.05e-09 (parameter vector 02)
     0.54     24.78    411.83  8.89e-12    239.09  5.88e-09 (parameter vector 03)
     0.61     26.85    221.41  3.42e-12    223.38  4.70e-09 (parameter vector 04)
     0.52     19.50    188.54  3.22e-12    236.63  7.92e-09 (parameter vector 05)
     0.81     25.69    173.93  3.82e-12    242.02  1.97e-09 (parameter vector 06)
     0.66     29.94    148.59  7.41e-12    224.49  9.31e-09 (parameter vector 07)
     0.55     15.25      92.53  5.63e-12    244.18  4.42e-09 (parameter vector 08)
     0.74     19.80    401.31  3.02e-12    228.42  3.58e-09 (parameter vector 09)
      soft,      hard,     Bmax,           Sb,      Tmin,            St

Analysis scripts to determine the surge characteristics (GSM and PISM)

The analysis scripts (four Python and one Bash script) are used to conduct the analysis described in "Numerical issues in modeling ice sheet instabilities such as binge-purge type cyclic ice stream surging".

1. analysis01.runPeaks.sh: Top-level bash script that sequentially runs all four python scripts.
2. analysis02.rPeakidSQts.py: Identifies the surges (peaks) and calculates the surge characteristics for every parameter vector and model setup.
3. analysis03.readStats.py: Takes the output of analysis02.rPeakidSQts.py and calculates the average (across all parameter vectors) change in surge characteristics between the base setup and every comparison model setup.
4. analysis04.rRMSEmeanBIASidSQts.py: Calculates the ice volume RMSE (Root Mean Square Error) and mean bias for every parameter vector and model setup (compared to the base setup)
5. analysis05.rmseMeanBias.py: Takes the output of analysis04.rRMSEmeanBIASidSQts.py and calculates the average (across all parameter vectors) ice volume RMSE and mean bias between the base setup and every comparison model setup.

Video 01: Glacial Systems Model (GSM) - surge onset

Shown are the ice sheet surface elevation (white/black contour lines), the Shallow Ice Approximation (SIA) part of the deformation work, the heating due to the Shallow Shelf Approximation (SSA) dynamics, the heat generated by basal shearing, the basal temperature with respect to the pressure melting point, the ice velocity at the base of the ice sheet in the x-direction (positive values represent eastward velocities), and the ice velocity at the base of the ice sheet in the y-direction (positive values represent northward flow) for parameter vector 1. The horizontal grid resolution is 3.125 km and the maximum model time step is 1 yr. The animation time step is 10 yr with a total run time of 200 yr. The animation is focused on the pseudo-Hudson Strait mouth with the ocean situated east of 450 km. The green and grey lines mark the warm-based area and the ice sheet margin, respectively. Note that the velocities in the GSM are defined on a staggered grid and the velocity fields in x- and y-direction are therefore shifted by half a grid cell length in x- and y-direction, respectively.

Video 02: Glacial Systems Model (GSM) - surge propagation and termination

Shown are the ice sheet surface elevation (white/black contour lines), the Shallow Ice Approximation (SIA) part of the deformation work, the heating due to the Shallow Shelf Approximation (SSA) dynamics, the heat generated by basal shearing, the basal temperature with respect to the pressure melting point, the ice velocity at the base of the ice sheet in the x-direction (positive values represent eastward velocities), and the ice velocity at the base of the ice sheet in the y-direction (positive values represent northward flow) for parameter vector 1. The horizontal grid resolution is 3.125 km and the maximum model time step is 1 yr. The animation time step is 10 yr with a total run time of 200 yr. The magenta line outlines the soft sediment pseudo-Hudson Bay and Hudson Strait area. The green and grey lines mark the warm-based area and the ice sheet margin, respectively. Note that the velocities in the GSM are defined on a staggered grid and the velocity fields in x- and y-direction are therefore shifted by half a grid cell length in x- and y-direction, respectively.

Video 03: Glacial Systems Model (GSM) - abrupt vs. 25 km wide transition between a soft and hard bed

Shown are the ice sheet surface elevation (white/black contour lines), the basal temperature with respect to the pressure melting point (left column), the ice velocity at the base of the ice sheet in the x-direction (center column, positive values represent eastward velocities), and the ice velocity at the base of the ice sheet in the y-direction (right column, positive values represent northward flow) for parameter vector 1. The upper row uses an abrupt transition between soft sediment and hard bedrock, whereas the bottom row has a 25 km wide transition zone. The horizontal grid resolution is 3.125 km and the maximum model time step is 1 yr. The animation time step is 10 yr with a total run time of 4 kyr. The magenta line outlines the soft sediment pseudo-Hudson Bay and Hudson Strait area. The hatched magenta area in the bottom row indicates the transition region surrounding the pseudo-Hudson Strait. The green and grey lines mark the warm-based area and the ice sheet margin, respectively. Note that the velocities in the GSM are defined on a staggered grid and the velocity fields in x- and y-direction are therefore shifted by half a grid cell length in x- and y-direction, respectively.

Video 04: Glacial Systems Model (GSM) - pseudo-Hudson Bay and Hudson Strait topography - two different slopes

Shown are the ice sheet surface elevation (white/black contour lines), the basal temperature with respect to the pressure melting point (left column), the ice velocity at the base of the ice sheet in the x-direction (center column, positive values represent eastward velocities), and the ice velocity at the base of the ice sheet in the y-direction (right column, positive values represent northward flow) for parameter vector 1. The upper row uses a 200 m deep pseudo-Hudson Bay and Hudson Strait topography with a 3.125 km transition zone (steep slope), whereas the bottom row has a 25 km wide transition zone (more gentle slope). The horizontal grid resolution is 3.125 km and the maximum model time step is 1 yr. The animation time step is 10 yr with a total run time of 4.5 kyr. The magenta line outlines the soft sediment pseudo-Hudson Bay and Hudson Strait area. The hatched magenta areas indicate the transition region (sediment and topography) surrounding the pseudo-Hudson Strait. The green and grey lines mark the warm-based area and the ice sheet margin, respectively. Note that the velocities in the GSM are defined on a staggered grid and the velocity fields in x- and y-direction are therefore shifted by half a grid cell length in x- and y-direction, respectively.

Video 05: Glacial Systems Model (GSM) - constant vs. resolution-dependent temperature ramp

Shown are the ice sheet surface elevation (white/black contour lines), the basal temperature with respect to the pressure melting point (left column), the ice velocity at the base of the ice sheet in the x-direction (center column, positive values represent eastward velocities), and the ice velocity at the base of the ice sheet in the y-direction (right column, positive values represent northward flow) for parameter vector 1. The upper row uses a constant temperature ramp (T_ramp = 0.0625°C), whereas the bottom row has a resolution-dependent temperature ramp (T_ramp = 0.5°C). T_exp is 28 in both cases. The horizontal grid resolution is 25 km and the maximum model time step is 1 yr. The animation time step is 1 kyr with a total run time of 200 kyr. The magenta line outlines the soft sediment pseudo-Hudson Bay and Hudson Strait area. The green and grey lines mark the warm-based area and the ice sheet margin, respectively. Note that the velocities in the GSM are defined on a staggered grid and the velocity fields in x- and y-direction are therefore shifted by half a grid cell length in x- and y-direction, respectively.

Video 06: Parallel Ice Sheet Model (PISM) - surge-affected area

Shown is the ice sheet surface elevation (white/black contour lines), the magnitude of the horizontal ice velocity at the base of the ice sheet, the basal temperature with respect to the pressure melting point, and the ice basal melt rate from energy conservation and subshelf melt for parameter vector 1 and the PISM base setup. The animation focuses on the inner 500 to 2500 km of the model domain. The horizontal grid resolution is 25 km and the maximum model time step is 1 yr. The animation time step is 1 kyr with a total run time of 200 kyr.

Video 07: Glacial Systems Model (GSM) - flat vs pseudo-Hudson Bay and Hudson Strait topography

Shown are the ice sheet surface elevation (white/black contour lines), the basal temperature with respect to the pressure melting point (left column), the ice velocity at the base of the ice sheet in the x-direction (center column, positive values represent eastward velocities), and the ice velocity at the base of the ice sheet in the y-direction (right column, positive values represent northward flow) for parameter vector 1. The upper row uses a flat topography, whereas the bottom row has a 200 m deep pseudo-Hudson Bay and Hudson Strait topography. The horizontal grid resolution is 3.125 km and the model time step is 1 yr. The animation time step is 10 yr with a total run time of 4 kyr. The magenta line outlines the soft sediment pseudo-Hudson Bay and Hudson Strait area. The hatched magenta area in the bottom row indicates the transition region surrounding the pseudo-Hudson Strait. The green and grey lines mark the warm-based area and the ice sheet margin, respectively. Note that the velocities in the GSM are defined on a staggered grid and the velocity fields in x- and y-direction are therefore shifted by half a grid cell length in x- and y-direction, respectively.

Video 08: Glacial Systems Model (GSM) - flat vs pseudo-Hudson Bay and Hudson Strait topography - basal heating

Shown are the ice sheet surface elevation (white/black contour lines), the Shallow Ice Approximation (SIA) part of the deformation work (left column), the heating due to the Shallow Shelf Approximation (SSA) dynamics (center column), and the heat generated by basal shearing (right column) for parameter vector 1. The upper row uses a flat topography, whereas the bottom row has a 200 m deep pseudo-Hudson Bay and Hudson Strait topography. The horizontal grid resolution is 3.125 km and the model time step is 1 yr. The animation time step is 10 yr with a total run time of 5 kyr. The magenta line outlines the soft sediment pseudo-Hudson Bay and Hudson Strait area. The hatched magenta area in the bottom row indicates the transition region surrounding the pseudo-Hudson Strait. The green and grey lines mark the warm-based area and the ice sheet margin, respectively.

Video 09: Parallel Ice Sheet Model (PISM) - horizontal grid resolution scaling

Shown are the ice sheet surface elevation (white/black contour lines), the basal temperature with respect to the pressure melting point (upper panels), and the magnitude of the horizontal ice velocity at the base of the ice sheet (lower panels) for parameter vector 1 using PISM. The animation focuses on the inner 500 to 2500 km of the model domain. The horizontal grid resolution is 50 km (left), 25 km (center), and 12.5 km (right). The model time step is 1 yr and the animation time step is 1 kyr. The model and animation run time is 200 kyr.