tross.rst - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tross.rst (12024B)
---
1 .. include:: ../../global.txt
2
3 .. _sec-ross:
4
5 An SSA flow model for the Ross Ice Shelf in Antarctica
6 ------------------------------------------------------
7
8 As part of the EISMINT series of intercomparisons, MacAyeal and others :cite:`MacAyealetal`
9 successfully validated early-1990s ice shelf numerical models using velocity data for the
10 Ross ice shelf. The data were from the RIGGS survey :cite:`RIGGS2`, acquired in the period
11 1973--1978 and measured at a few hundred locations in a grid across the shelf. Substantial
12 modelling developments followed EISMINT-Ross, including inverse modeling to recover
13 depth-averaged viscosity :cite:`RommelaereMacAyeal` and parameter-sensitivity studies
14 :cite:`HumbertGreveHutter`. Previous PISM versions set up the EISMINT-Ross flow model and
15 performed the diagnostic computation, with RIGGS data for validation.
16
17 However, availability of rich new data sets for ice sheet modeling, including the ALBMAP
18 v1 :cite:`LeBrocqetal2010` ice sheet geometry, bedrock, and climate data set, and the
19 radar-derived (InSAR) MEaSUREs Antarctica Velocity Map :cite:`Rignotetal2011`, allows us to
20 use more complete, recent, and higher-resolution data for the same basic job. Furthermore
21 one can extend the diagnostic Ross ice shelf calculation both to other ice shelves around
22 Antarctica and to time-evolving ("prognostic") cases using the eigencalving
23 :cite:`Levermannetal2012` mechanisms.
24
25 The scripts in this subsection are found in directory ``examples/ross/``. In summary, the
26 script ``preprocess.py`` downloads easily-available data and builds a NetCDF input file
27 for PISM. For the diagnostic computation we document first, the script ``run_diag.sh`` (in
28 subdirectory ``examples/ross/diagnostic/``) runs PISM. The script ``plot.py`` shows a
29 comparison of observations and model results, as in :numref:`fig-rosspython`.
30
31 Preprocessing input data
32 ^^^^^^^^^^^^^^^^^^^^^^^^
33
34 NSIDC_ requires registration to download the `MEaSUREs InSAR-Based Antarctica Ice Velocity
35 Map <measures-ross_>`_; please register, log in, and get the *first*, 900 m version
36 of it (``antarctica_ice_velocity_900m.nc``) here:
37
38 https://n5eil01u.ecs.nsidc.org/MEASURES/NSIDC-0484.001/1996.01.01/
39
40 Put this file in ``examples/ross``.
41
42 The script ``preprocess.py`` then downloads ALBMAP using ``wget``; these two files total
43 around 350 Mb. Then it uses NCO_ to cut out the relevant portion of the grid and CDO_ to
44 conservatively-interpolate the high-resolution (900 m) velocity data onto the coarser (5
45 km) geometry grid used in ALBMAP. The script ``nc2cdo.py`` from directory ``util/``,
46 prepares the NetCDF file for the application of CDO, which requires complete projection
47 information. Run
48
49 .. code-block:: none
50
51 cd examples/ross/
52 ./preprocess.py
53
54 to download ALBMAP and finish the pre-processing.
55
56 The NetCDF file ``Ross_combined.nc`` produced by ``preprocess.py`` contains ice thickness,
57 bed elevations, surface temperature, net accumulation, as well as latitude and longitude
58 values. All of these are typical of ice sheet modelling data, both in diagnostic runs and
59 as needed to initialize and provide boundary conditions for prognostic (evolutionary)
60 runs; see below for the prognostic case with these data. The ``_combined`` file also has
61 variables ``u_ssa_bc`` and ``v_ssa_bc`` for the boundary values used in the (diagnostic
62 and prognostic) computation of velocity. They are used at all grounded locations and at
63 ice shelf cells that are immediate neighbors of grounded ice. The variable ``bc_mask``
64 specifies these locations. Finally the variables ``u_ssa_bc,v_ssa_bc``, which contain
65 observed values, are used after the run to compare to the computed interior velocities.
66
67 Diagnostic computation of ice shelf velocity
68 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
69
70 The diagnostic velocity computation bootstraps from ``Ross_combined.nc`` and performs a
71 short run; in the `211\times 211` grid case we demonstrate below, the key parts of the
72 PISM command are
73
74 .. code-block:: none
75
76 pismr -i ../Ross_combined.nc -bootstrap -Mx 211 -My 211 -Mz 3 -Lz 3000 -z_spacing equal \
77 -surface given -stress_balance ssa -energy none -no_mass -yield_stress constant -tauc 1e6 \
78 -pik -ssa_dirichlet_bc -y 1.0 -ssa_e 0.6 -ssafd_ksp_monitor
79
80 The computational grid here is the "native" `5` km data grid used in ALBMAP. Regarding the
81 options,
82
83 - The maximum thickness of the ice is `2766` m so we choose a height for the computational
84 box large enough to contain the ice (i.e. ``-Lz 3000``). Vertical grid resolution
85 is, however, unimportant in this case because we use the SSA stress balance only, and
86 the temperature set at bootstrapping suffices to determine the ice softness; thus the
87 options ``-Mz 3 -z_spacing equal``.
88
89 - Option ``-stress_balance ssa`` selects the SSA stress balance and turns off the SIA
90 stress balance computation, since our goal is to model the ice shelf. It also side-steps
91 a technical issue: PISM uses periodic boundary conditions at domain boundaries and most
92 fields in this setup are not periodic. Turning off SIA avoids operations such as
93 differencing surface elevation across the domain edges. For a more complete solution to
94 this technical issue see section :ref:`sec-jako` about a regional model using PISM's
95 "regional mode" ``pismr -regional`` and the option :opt:`-no_model_strip`.
96
97 - Option ``-y 1.0 -no_mass -energy none`` chooses a "diagnostic" run: in absence of
98 geometry evolution and stability restrictions of the energy balance model a
99 one-year-long run will be covered by exactly one time step.
100
101 - Option ``-pik`` is equivalent to ``-cfbc -part_grid -kill_icebergs -subgl`` in this
102 non-evolving example. Note that ``-kill_icebergs`` removes effectively-detached bits of
103 ice, especially in McMurdo sound area, so that the SSA problem is well-posed for the
104 grounded-ice-sheet-connected ice shelf.
105
106 - Option :opt:`-ssa_dirichlet_bc` forces the use of fields ``u_ssa_bc, v_ssa_bc, bc_mask``
107 described above. The field ``bc_mask`` is `1` at boundary condition locations, and `0`
108 elsewhere. For the prognostic runs below, the ice thickness is also fixed at boundary
109 condition locations, so as to prescribe ice flux as an ice shelf input.
110
111 - Options ``-yield_stress constant -tauc 1e6`` essentially just turn off the
112 grounded-ice evolving yield stress mechanism, which is inactive anyway, and force a high
113 resistance under grounded ice so it does not slide.
114
115 - Option ``-ssa_e 0.6`` is the single tuned parameter; this value gives good
116 correlation between observed and modeled velocity magnitudes.
117
118 - Option ``-ssafd_ksp_monitor`` provides feedback on the linear solver iterations
119 "underneath" the nonlinear (shear-thinning) SSA solver iteration.
120
121 There is no need to type in the above command; just run
122
123 .. code-block:: none
124
125 cd diagnostic/
126 ./run_diag.sh 2 211 0.6
127
128 Note ``run_diag.sh`` accepts three arguments: ``run_diag.sh N Mx E`` does a run
129 with ``N`` MPI processes, an ``Mx`` by ``Mx`` grid, and option
130 ``-ssa_e E``. The choices above give a run which only takes a few seconds, and it
131 produces output file ``diag_Mx211.nc``.
132
133 There are many reasonable choices for the effective softness of an ice shelf, as ice
134 density, temperature, and the presence of fractures all influence the effective softness.
135 Using an enhancement factor :opt:`-ssa_e 0.6` acknowledges that the physical justification
136 for tuning the ice softness is uncertain. One could instead use the temperature itself or
137 the ice density\ [#]_ as tuning parameters, and these are worthwhile experiments for the
138 interested PISM user.
139
140 The script ``plot.py`` takes PISM output such as ``diag_Mx211.nc`` to produce
141 :numref:`fig-rosspython`. The run shown in the figure used an enhancement factor of
142 `0.6` as above. The thin black line outlines the floating shelf, which is the actual
143 modeling domain here. To generate this figure yourself, run
144
145 .. code-block:: none
146
147 ../plot.py diag_Mx211.nc
148
149 .. figure:: figures/ross-results.png
150 :name: fig-rosspython
151
152 *Left*: Color is speed in m/a. Arrows are observed (white) and modeled (black)
153 velocities. *Right*: Comparison between modeled and observed speeds at points plotted
154 on the left.
155
156 Extending this example to other ice shelves
157 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
158
159 The SSA diagnostic solution described in this section can be easily applied to other ice
160 shelves in Antarctica, such as the Filchner-Ronne Ice Shelf modeled using PISM in
161 :cite:`AlbrechtLevermann2012`, for example.
162
163 Simply choose a different rectangular domain, within the area covered by the
164 whole-Antarctic data-sets used here, at the preprocessing stage. In particular you should
165 modify the lines "``ncks -O -d x1,439,649 -d y1,250,460 ...``" (for ALBMAP data) and
166 "``ncks -d x,2200,3700 -d y,3500,4700 ...``" (for MEaSUREs velocity data) in the
167 script ``examples/ross/preprocess.py``.
168
169 Prognostic modelling using eigencalving
170 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
171
172 Next we summarize how you can create an evolving-geometry model of the Ross ice shelf with
173 constant-in-time inflow across the fixed grounding line. See ``README.md`` and
174 ``run_prog.sh`` in ``examples/ross/prognostic/``. This example also demonstrates the
175 :opt:`-calving eigen_calving` model for a moving calving front :cite:`Levermannetal2012`.
176
177 Start by running ``preprocess.py`` in ``examples/ross/`` as described above. If
178 you have already done the diagnostic example above, then this stage is complete.
179
180 Then change to the ``prognostic/`` directory and run the default example:
181
182 .. code-block:: none
183
184 cd examples/ross/prognostic/
185 ./run_prog.sh 4 211 0.6 100
186
187 This 100 model year run on 4 processes and a 5 km grid took about forty minutes on a 2016
188 laptop. It starts with a bootstrapping stage which does a ``-y 1.0`` run, which generates
189 ``startfile_Mx211.nc``. It then re-initializes to start the prognostic run itself. See the
190 ``README.md`` for a bit more on the arguments taken by ``run_prog.sh`` and on viewing the
191 output files.
192
193 The PISM command done here is (essentially, and without showing diagnostic output choices)
194
195 .. code-block:: none
196
197 pismr -i startfile_Mx211.nc -surface given -stress_balance ssa \
198 -yield_stress constant -tauc 1e6 -pik -ssa_dirichlet_bc -ssa_e 0.6 \
199 -y 100 -o prog_Mx211_yr100.nc -o_size big \
200 -calving eigen_calving,thickness_calving -eigen_calving_K 1e17 \
201 -front_retreat_cfl -thickness_calving_threshold 50.0
202
203 Several of these options are different from those used in the diagnostic case. First,
204 while the command ``-pik`` is the same as before, now each part of its expansion, namely
205 ``-cfbc -part_grid -kill_icebergs -subgl``, is important. As the calving front evolves
206 (i.e. regardless of the calving law choices), option ``-part_grid`` moves the calving
207 front by one grid cell only when the cell is full of the ice flowing into it; see
208 :cite:`Albrechtetal2011`. The option ``-kill_icebergs`` is essential to maintain well-posedness
209 of the SSA velocity problem at each time step :cite:`Winkelmannetal2011`. See section
210 :ref:`sec-pism-pik`.
211
212 Option combination
213
214 .. code-block:: none
215
216 -calving eigen_calving,thickness_calving -eigen_calving_K 1e17 \
217 -front_retreat_cfl -thickness_calving_threshold 50.0
218
219 specifies that ice at the calving front will be removed if either a criterion on the
220 product of principal stresses is satisfied :cite:`Levermannetal2012`, namely ``eigen_calving``
221 with the given constant `K`, or if the ice thickness goes below the given threshold of 50
222 meters. See section :ref:`sec-calving`.
223
224 .. %FIXME Use evolving fracture density. See ``README.md``, ``preprocess_frac.py``, and
225 ``run_frac.sh`` in directory ``examples/ross/fracture_density/``. This example
226 demonstrates the fracture density transport model in :cite:`AlbrechtLevermann2012`.
227
228 .. rubric:: Footnotes
229
230 .. [#] High accumulation rates, cold firn with minimal compression, and basal freeze-on of
231 marine ice may all generate significant variation in shelf density.