An SSA flow model for the Ross Ice Shelf in Antarctica¶
As part of the EISMINT series of intercomparisons, MacAyeal and others  successfully validated early-1990s ice shelf numerical models using velocity data for the Ross ice shelf. The data were from the RIGGS survey , acquired in the period 1973–1978 and measured at a few hundred locations in a grid across the shelf. Substantial modelling developments followed EISMINT-Ross, including inverse modeling to recover depth-averaged viscosity  and parameter-sensitivity studies . Previous PISM versions set up the EISMINT-Ross flow model and performed the diagnostic computation, with RIGGS data for validation.
However, availability of rich new data sets for ice sheet modeling, including the ALBMAP v1  ice sheet geometry, bedrock, and climate data set, and the radar-derived (InSAR) MEaSUREs Antarctica Velocity Map , allows us to use more complete, recent, and higher-resolution data for the same basic job. Furthermore one can extend the diagnostic Ross ice shelf calculation both to other ice shelves around Antarctica and to time-evolving (“prognostic”) cases using the eigencalving  mechanisms.
The scripts in this subsection are found in directory
examples/ross/. In summary, the
preprocess.py downloads easily-available data and builds a NetCDF input file
for PISM. For the diagnostic computation we document first, the script
examples/ross/diagnostic/) runs PISM. The script
plot.py shows a
comparison of observations and model results, as in Fig. 41.
Preprocessing input data¶
NSIDC requires registration to download the MEaSUREs InSAR-Based Antarctica Ice Velocity
Map; please register, log in, and get the first, 900 m version
of it (
Put this file in
preprocess.py then downloads ALBMAP using
wget; these two files total
around 350 Mb. Then it uses NCO to cut out the relevant portion of the grid and CDO to
conservatively-interpolate the high-resolution (900 m) velocity data onto the coarser (5
km) geometry grid used in ALBMAP. The script
nc2cdo.py from directory
prepares the NetCDF file for the application of CDO, which requires complete projection
cd examples/ross/ ./preprocess.py
to download ALBMAP and finish the pre-processing.
The NetCDF file
Ross_combined.nc produced by
preprocess.py contains ice thickness,
bed elevations, surface temperature, net accumulation, as well as latitude and longitude
values. All of these are typical of ice sheet modelling data, both in diagnostic runs and
as needed to initialize and provide boundary conditions for prognostic (evolutionary)
runs; see below for the prognostic case with these data. The
_combined file also has
v_bc for the boundary values used in the (diagnostic
and prognostic) computation of velocity. They are used at all grounded locations and at
ice shelf cells that are immediate neighbors of grounded ice. The variable
specifies these locations. Finally the variables
u_bc,v_bc, which contain
observed values, are used after the run to compare to the computed interior velocities.
Diagnostic computation of ice shelf velocity¶
The diagnostic velocity computation bootstraps from
Ross_combined.nc and performs a
short run; in the \(211\times 211\) grid case we demonstrate below, the key parts of the
PISM command are
pismr -i ../Ross_combined.nc -bootstrap -Mx 211 -My 211 -Mz 3 -Lz 3000 -z_spacing equal \ -surface given -stress_balance ssa -energy none -no_mass -yield_stress constant -tauc 1e6 \ -pik -ssa_dirichlet_bc -y 1.0 -ssa_e 0.6 -ssafd_ksp_monitor
The computational grid here is the “native” \(5\) km data grid used in ALBMAP. Regarding the options,
The maximum thickness of the ice is \(2766\) m so we choose a height for the computational box large enough to contain the ice (i.e.
-Lz 3000). Vertical grid resolution is, however, unimportant in this case because we use the SSA stress balance only, and the temperature set at bootstrapping suffices to determine the ice softness; thus the options
-Mz 3 -z_spacing equal.
-stress_balance ssaselects the SSA stress balance and turns off the SIA stress balance computation, since our goal is to model the ice shelf. It also side-steps a technical issue: PISM uses periodic boundary conditions at domain boundaries and most fields in this setup are not periodic. Turning off SIA avoids operations such as differencing surface elevation across the domain edges. For a more complete solution to this technical issue see section Example: A regional model of the Jakobshavn outlet glacier in Greenland about a regional model using PISM’s “regional mode”
pismr -regionaland the option
-y 1.0 -no_mass -energy nonechooses a “diagnostic” run: in absence of geometry evolution and stability restrictions of the energy balance model a one-year-long run will be covered by exactly one time step.
-pikis equivalent to
-cfbc -part_grid -kill_icebergs -subglin this non-evolving example. Note that
-kill_icebergsremoves effectively-detached bits of ice, especially in McMurdo sound area, so that the SSA problem is well-posed for the grounded-ice-sheet-connected ice shelf.
-ssa_dirichlet_bcforces the use of fields
u_bc, v_bc, vel_bc_maskdescribed above. The field
vel_bc_maskis \(1\) at boundary condition locations, and \(0\) elsewhere. For the prognostic runs below, the ice thickness is also fixed at boundary condition locations, so as to prescribe ice flux as an ice shelf input.
-yield_stress constant -tauc 1e6essentially just turn off the grounded-ice evolving yield stress mechanism, which is inactive anyway, and force a high resistance under grounded ice so it does not slide.
-ssa_e 0.6is the single tuned parameter; this value gives good correlation between observed and modeled velocity magnitudes.
-ssafd_ksp_monitorprovides feedback on the linear solver iterations “underneath” the nonlinear (shear-thinning) SSA solver iteration.
There is no need to type in the above command; just run
cd diagnostic/ ./run_diag.sh 2 211 0.6
run_diag.sh accepts three arguments:
run_diag.sh N Mx E does a run
N MPI processes, an
Mx grid, and option
-ssa_e E. The choices above give a run which only takes a few seconds, and it
produces output file
There are many reasonable choices for the effective softness of an ice shelf, as ice
density, temperature, and the presence of fractures all influence the effective softness.
Using an enhancement factor
-ssa_e 0.6 acknowledges that the physical justification
for tuning the ice softness is uncertain. One could instead use the temperature itself or
the ice density1 as tuning parameters, and these are worthwhile experiments for the
interested PISM user.
plot.py takes PISM output such as
diag_Mx211.nc to produce
Fig. 41. The run shown in the figure used an enhancement factor of
\(0.6\) as above. The thin black line outlines the floating shelf, which is the actual
modeling domain here. To generate this figure yourself, run
Extending this example to other ice shelves¶
The SSA diagnostic solution described in this section can be easily applied to other ice shelves in Antarctica, such as the Filchner-Ronne Ice Shelf modeled using PISM in , for example.
Simply choose a different rectangular domain, within the area covered by the
whole-Antarctic data-sets used here, at the preprocessing stage. In particular you should
modify the lines “
ncks -O -d x1,439,649 -d y1,250,460 ...” (for ALBMAP data) and
ncks -d x,2200,3700 -d y,3500,4700 ...” (for MEaSUREs velocity data) in the
Prognostic modelling using eigencalving¶
Next we summarize how you can create an evolving-geometry model of the Ross ice shelf with
constant-in-time inflow across the fixed grounding line. See
examples/ross/prognostic/. This example also demonstrates the
-calving eigen_calving model for a moving calving front .
Start by running
examples/ross/ as described above. If
you have already done the diagnostic example above, then this stage is complete.
Then change to the
prognostic/ directory and run the default example:
cd examples/ross/prognostic/ ./run_prog.sh 4 211 0.6 100
This 100 model year run on 4 processes and a 5 km grid took about forty minutes on a 2016
laptop. It starts with a bootstrapping stage which does a
-y 1.0 run, which generates
startfile_Mx211.nc. It then re-initializes to start the prognostic run itself. See the
README.md for a bit more on the arguments taken by
run_prog.sh and on viewing the
The PISM command done here is (essentially, and without showing diagnostic output choices)
pismr -i startfile_Mx211.nc -surface given -stress_balance ssa \ -yield_stress constant -tauc 1e6 -pik -ssa_dirichlet_bc -ssa_e 0.6 \ -y 100 -o prog_Mx211_yr100.nc -o_size big \ -calving eigen_calving,thickness_calving -eigen_calving_K 1e17 \ -front_retreat_cfl -thickness_calving_threshold 50.0
Several of these options are different from those used in the diagnostic case. First,
while the command
-pik is the same as before, now each part of its expansion, namely
-cfbc -part_grid -kill_icebergs -subgl, is important. As the calving front evolves
(i.e. regardless of the calving law choices), option
-part_grid moves the calving
front by one grid cell only when the cell is full of the ice flowing into it; see
. The option
-kill_icebergs is essential to maintain well-posedness
of the SSA velocity problem at each time step . See section
PIK options for marine ice sheets.
-calving eigen_calving,thickness_calving -eigen_calving_K 1e17 \ -front_retreat_cfl -thickness_calving_threshold 50.0
specifies that ice at the calving front will be removed if either a criterion on the
product of principal stresses is satisfied , namely
with the given constant \(K\), or if the ice thickness goes below the given threshold of 50
meters. See section Calving and front retreat.
High accumulation rates, cold firn with minimal compression, and basal freeze-on of marine ice may all generate significant variation in shelf density.