tindex.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
---
tindex.rst (17958B)
---
1 .. include:: ../../global.txt
2
3 .. _sec-jako:
4
5 Example: A regional model of the Jakobshavn outlet glacier in Greenland
6 =======================================================================
7
8 Jakobshavn Isbrae is a fast-flowing outlet glacier in western Greenland that drains
9 approximately 7% of the area of the Greenland ice sheet. It experienced a large
10 acceleration following the loss of its floating tongue in the 1990s
11 :cite:`JoughinAbdalatiFahnestock`, an event which seems to have been driven by warmer
12 ocean temperatures :cite:`Hollandetal2008`. Because it is thick, has a steep surface
13 slope, has a deep trough in its bedrock topography (Figure :numref:`fig-jako-basin-topg`),
14 and has a thick layer of low-viscosity temperate ice at its base :cite:`Luethietal2009`,
15 this ice flow is different from the ice streams in West Antarctica or Northeast Greenland
16 :cite:`TrufferEchelmeyer`.
17
18 This section describes how to build a PISM regional model of this outlet glacier
19 :cite:`DellaGiustina2011` using scripts from ``examples/jako/``. The same strategy should
20 work for other outlet glaciers. We also demonstrate the PISM regional mode ``pismr
21 -regional``, and Python `drainage-basin-delineation tools <regional-tools_>`_ which can be
22 downloaded from the PISM source code website. Such regional models allow modest-size
23 computers to run high resolution models [#]_ and large ensembles. Regional analysis is
24 justified if detailed data is available for the region.
25
26 The geometric data used here is the SeaRISE :cite:`Bindschadler2013SeaRISE` 1 km dataset
27 for the whole Greenland ice sheet. It contains bedrock topography from recent CReSIS radar
28 in the Jakobshavn area. We also use the SeaRISE 5 km data set which has climatic mass
29 balance from the Greenland-region climate model RACMO :cite:`Ettemaetal2009`.
30
31 A regional ice flow model generally needs ice flow and stress boundary conditions. For
32 this we use a 5 km grid, whole ice sheet, spun-up model state from PISM, described in
33 :ref:`sec-start` of this *Manual*. You can download the large NetCDF result from
34 the PISM website, or you can generate it by running a script from :ref:`sec-start`.
35
36 .. figure:: figures/jako-mask-topg.png
37 :name: fig-jako-basin-topg
38
39 A ``regional-tools`` script computes a drainage basin mask from the surface DEM (left;
40 Modis background) and from a user-identified terminus rectangle (blue). The regional
41 model can exploit high-resolution bedrock elevations inland from Jakobshavn fjord
42 (right; meters asl).
43
44 Get the drainage basin delineation tool
45 ---------------------------------------
46
47 The drainage basin tool ``regional-tools`` is at |pism-regional-tools-url|.
48 Get it using ``git`` and set it up as directed in its ``README.md``. Then come back to the
49 ``examples/jako/`` directory and link the script. Here is the quick summary:
50
51 .. code-block:: bash
52
53 cd ~/usr/local/ # the location you want
54 git clone https://github.com/pism/regional-tools.git
55 cd regional-tools/
56 python setup.py install # may add "sudo" or "--user"
57 cd PISM/examples/jako/
58 ln -s ~/usr/local/regional-tools/pism_regional.py . # symbolic link to tool
59
60
61 Preprocess the data and get the whole ice sheet model file
62 ----------------------------------------------------------
63
64 Script ``preprocess.sh`` downloads and cleans the 1 km SeaRISE data, an 80 Mb file called
65 ``Greenland1km.nc``. [#]_ The script also downloads the SeaRISE 5 km data set
66 ``Greenland_5km_v1.1.nc``, which contains the RACMO surface mass balance field (not
67 present in the 1 km data set). If you have already run the example in :ref:`sec-start`
68 then you already have this file and you can link to it to avoid downloading:
69
70 .. code-block:: none
71
72 ln -s ../std-greenland/Greenland_5km_v1.1.nc
73
74 The same script also preprocesses a pre-computed 5 km grid PISM model result
75 ``g5km_gridseq.nc`` for the whole ice sheet. This provides the boundary conditions, and
76 the thermodynamical initial condition, for the regional flow model we are building. If you
77 have already generated it by running the script in section :ref:`sec-gridseq` then link
78 to it,
79
80 .. code-block:: none
81
82 ln -s ../std-greenland/g5km_gridseq.nc
83
84 Otherwise running ``preprocess.sh`` will download it. Because it is about 0.6 Gb this may
85 take some time.
86
87 So now let's actual run the preprocessing script:
88
89 .. code-block:: none
90
91 ./preprocess.sh
92
93 Files ``gr1km.nc``, ``g5km_climate.nc``, and ``g5km_bc.nc`` will appear. These can be
94 examined in the usual ways, for example:
95
96 .. code-block:: none
97
98 ncdump -h gr1km.nc | less # read metadata
99 ncview gr1km.nc # view fields
100
101 The boundary condition file ``g5km_bc.nc`` contains thermodynamical spun-up variables
102 (``enthalpy,bmelt,bwat``) and boundary values for the sliding velocity
103 (``u_ssa_bc,v_ssa_bc``); these have been extracted from ``g5km_gridseq.nc``.
104
105 None of the above actions is specific to Jakobshavn, though all are specific to Greenland.
106 If your goal is to build a regional model of another outlet glacier in Greenland, then you
107 may be able to use ``preprocess.sh`` as is. The SeaRISE 1 km data set has recent CReSIS
108 bed topography data only for the vicinity of the Jakobshavn outlet, however, and it is
109 otherwise just BEDMAP. Because outlet glacier flows are bed-topography-dominated,
110 additional bed elevation data should be sought.
111
112 Identify the drainage basin for the modeled outlet glacier
113 ----------------------------------------------------------
114
115 Here we are going to extract a "drainage basin mask" from the surface elevation data (DEM)
116 in ``gr1km.nc``. The goal is to determine, in part, the locations outside of the drainage
117 basin where boundary conditions taken from the precomputed whole ice sheet run can be
118 applied to modeling the outlet glacier flow itself.
119
120 The basin mask is determined by the gradient flow of the surface elevation. Thus
121 generating the mask uses a highly-simplified ice dynamics model (namely: ice flows down
122 the surface gradient). Once we have the mask, we will apply the full PISM model in the
123 basin interior marked by the mask. Outside the basin mask we will apply simplified models
124 or use the whole ice sheet results as boundary conditions.
125
126 The script ``pism_regional.py`` computes the drainage basin mask based on a user choice of
127 a "terminus rectangle"; see Figure :numref:`fig-jako-basin-topg`. There are two ways to
128 use this script:
129
130 - To use the graphical user interface (GUI) mode.
131
132 Run
133
134 .. code-block:: none
135
136 python pism_regional.py
137
138 Select ``gr1km.nc`` to open. Once the topographic map appears in the Figure window, you
139 may zoom enough to see the general outlet glacier area. Then select the button "Select
140 terminus rectangle". Use the mouse to select a small rectangle around the Jakobshavn
141 terminus (calving front), or around the terminus of another glacier if you want to model
142 that. Once you have a highlighted rectangle, select a "border width" of at least 50
143 cells. [#]_ Then click "Compute the drainage basin mask." Because this is a large data
144 set there will be some delay. (Multi-core users will see that an automatic parallel
145 computation is done.) Finally click "Save the drainage basin mask" and save with your
146 preferred name; we will assume it is called ``jakomask.nc``. Then quit.
147
148 - To use the command-line interface.
149
150 The command-line interface of ``pism_regional.py`` allows one to re-create the mask
151 without changing the terminus rectangle choice. (It also avoids the slowness of the GUI
152 mode for large data sets.) In fact, for repeatability, we will assume you have used this
153 command to calculate the drainage basin:
154
155 .. code-block:: none
156
157 python pism_regional.py -i gr1km.nc -o jakomask.nc -x 360,382 -y 1135,1176 -b 50
158
159 This call generates the red region in :numref:`fig-jako-basin-topg`. Options ``-x A,B -y
160 C,D`` identify the grid index ranges of the terminus rectangle, and option ``-b`` sets
161 the border width. To see more script options, run with ``--help``.
162
163 Cut out the computational domain for the regional model
164 -------------------------------------------------------
165
166 We still need to "cut out" from the whole ice sheet geometry data ``gr1km.nc`` the
167 computational domain for the regional model. The climate data file ``g5km_climate.nc`` and
168 the boundary condition file ``g5km_bc.nc`` do not need this action because PISM's coupling
169 and SSA boundary condition codes already handle interpolation and/or subsampling for such
170 data.
171
172 You may have noticed that the text output from running ``pism_regional.py`` included a
173 cutout command which uses ``ncks`` from the NCO tools. This command also appears as a
174 global attribute of ``jakomask.nc``:
175
176 .. code-block:: none
177
178 ncdump -h jakomask.nc | grep cutout
179
180 Copy and run the command that appears, something like
181
182 .. code-block:: none
183
184 ncks -d x,299,918 -d y,970,1394 gr1km.nc jako.nc
185
186 This command is also applied to the mask file; note the option ``-A`` for "append":
187
188 .. code-block:: none
189
190 ncks -A -d x,299,918 -d y,970,1394 jakomask.nc jako.nc
191
192 Now look at ``jako.nc``, for example with "``ncview -minmax all jako.nc``". This file is
193 the full geometry data ready for a regional model. The field ``ftt_mask`` identifies the
194 drainage basin, outside of which we will use simplified time-independent boundary
195 conditions. Specifically, outside of the ``ftt_mask`` area, but within the computational
196 domain defined by the extent of ``jako.nc``, we will essentially keep the initial
197 thickness. Inside the ``ftt_mask`` area all fields will evolve normally.
198
199 Quick start
200 -----------
201
202 The previous steps starting with the command "``./preprocess.sh``" above, then using the
203 command-line version of ``pism_regional.py``, and then doing the ``ncks`` cut-out steps,
204 are all accomplished in one script,
205
206 .. code-block:: none
207
208 ./quickjakosetup.sh
209
210 Running this takes about a minute on a fast laptop, assuming data files are already
211 downloaded.
212
213 Spinning-up the regional model on a 5 km grid
214 ----------------------------------------------
215
216 To run the PISM regional model we will need to know the number of grid points in the 1 km
217 grid in ``jako.nc``. Do this:
218
219 .. code-block:: none
220
221 ncdump -h jako.nc |head
222 netcdf jako {
223 dimensions:
224 y = 425 ;
225 x = 620 ;
226 ...
227
228 The grid has spacing of 1 km, so our computational domain is a 620 km by 425 km rectangle.
229 A 2 km resolution, century-scale model run is easily achievable on a desktop or laptop
230 computer, and that is our goal below. A lower 5 km resolution spin-up run, matching the
231 resolution of the 5 km whole ice sheet state computed earlier, is also achievable on a
232 small computer; we do that first.
233
234 The boundary condition fields in ``g5km_bc.nc``, from the whole ice sheet model result
235 ``g5km_gridseq.nc``, may or may not, depending on modeller intent, be spun-up adequately
236 for the purposes of the regional model. For instance, the intention may be to study
237 equilibrium states with model settings special to the region. Here, however we assume that
238 some regional spin-up is needed, if for no other reason that the geometry used here (from
239 the SeaRISE 1km data set) differs from that in the whole ice sheet model state.
240
241 We will get first an equilibrium 5 km regional model, and then do a century run of a 2 km
242 model based on that. While determining "equilibrium" requires a decision, of course, a
243 standard satisfied here is that the ice volume in the region changes by less than 0.1
244 percent in the final 100 model years. See ``ice_volume_glacierized`` in ``ts_spunjako_0.nc``
245 below.
246
247 The 5 km grid [#]_ uses ``-Mx 125 -My 86``. So now we do a basic run using 4 MPI
248 processes:
249
250 .. code-block:: none
251
252 ./spinup.sh 4 125 86 &> out.spin5km &
253
254 .. This takes 4.5055 proc-hours on bueler-gazelle
255
256 You can read the ``stdout`` log file while it runs: "``less out.spin5km``". The run takes
257 about 4.4 processor-hours on a 2016 laptop. It produces three files which can be viewed
258 (e.g. with ``ncview``): ``spunjako_0.nc``, ``ts_spunjako_0.nc``, and ``ex_spunjako_0.nc``.
259 Some more comments on this run are appropriate:
260
261 - Generally the regridding techniques used at the start of this spin-up run are
262 recommended for regional modeling. Read the actual run command by
263
264 .. code-block:: none
265
266 PISM_DO=echo ./spinup.sh 4 125 86 | less
267
268 - We use ``-i jako.nc -bootstrap``, so we get to choose our grid, and (as usual in PISM
269 with ``-bootstrap``) the fields are interpolated to our grid.
270
271 - A modestly-fine vertical grid with 20 m spacing is chosen, but even finer is
272 recommended, especially to resolve the temperate ice layer in these outlet glaciers.
273
274 - There is an option :opt:`-no_model_strip` ``10`` asking ``pismr -regional`` to put a 10
275 km strip around edge of the computational domain. This strip is entirely outside of the
276 drainage basin defined by ``ftt_mask``. In this strip the thermodynamical spun-up
277 variables ``bmelt,tillwat,enthalpy,litho_temp`` from ``g5km_bc.nc`` are held fixed and
278 used as boundary conditions for the conservation of energy model. A key part of putting
279 these boundary conditions into the model strip are the options
280
281 .. code-block:: none
282
283 -regrid_file g5km_bc.nc -regrid_vars bmelt,tillwat,enthalpy,litho_temp,vel_ssa_bc
284
285 - Dirichlet boundary conditions ``u_ssa_bc,v_ssa_bc`` are also regridded from
286 ``g5km_bc.nc`` for the sliding SSA stress balance, and the option ``-ssa_dirichlet_bc``
287 then uses them during the run. The SSA equations are solved as usual except in the
288 ``no_model_strip`` where these Dirichlet boundary conditions are used. Note that the
289 velocity tangent to the north and south edges of the computational domain is
290 significantly nonzero, which motivates this usage.
291
292 - The calving front of the glacier is handled by the following command-line options:
293
294 .. code-block:: none
295
296 -front_retreat_file jako.nc -pik
297
298 This choice uses the present-day ice extent, defined by SeaRISE data in
299 ``Greenland1km.nc``, to determine the location of the calving front (see
300 :ref:`sec-prescribed-retreat`). Recalling that ``-pik`` includes ``-cfbc``, we are
301 applying a PIK mechanism for the stress boundary condition at the calving front. The
302 other PIK mechanisms are largely inactive because prescribing the maximum ice extent,
303 but they should do no harm (see section :ref:`sec-pism-pik`).
304
305 .. figure:: figures/jako-csurf.png
306 :name: fig-jako-csurf
307
308 Left: modeled surface speed at the end of a 2 km grid, 100 model year, steady
309 present-day climate run. Right: observed surface speed, an average of four winter
310 velocity maps (2000,2006--2008) derived from RADARSAT data, as included in the SeaRISE
311 5 km data set :cite:`Joughinetal2010`, for the same region. Scales are in meters per
312 year.
313
314 Century run on a 2 km grid
315 --------------------------
316
317 Now that we have a spun-up state, here is a 100 model year run on a 2 km grid with a 10 m
318 grid in the vertical:
319
320 .. code-block:: none
321
322 ./century.sh 4 311 213 spunjako_0.nc &> out.2km_100a &
323
324 This run requires at least 6 GB of memory, and it takes about 16 processor-hours.
325
326 It produces a file ``jakofine_short.nc`` almost immediately and then restarts from it
327 because we need to regrid fields from the end of the previous 5 km regional run (in
328 ``spunjako_0.nc``) and then to "go back" and regrid the SSA boundary conditions from the 5
329 km whole ice sheet results ``g5km_bc.nc``. At the end of the run the final file
330 ``jakofine.nc`` is produced. Also there is a time-series file ``ts_jakofine.nc`` with
331 monthly scalar time-series and a spatial time-dependent file ``ex_jakofine.nc``. The
332 surface speed at the end of this run is shown in :numref:`fig-jako-csurf`, with a
333 comparison to observations.
334
335 Over this 100 year period the flow appears to be relatively steady state. Though this is
336 not surprising because the climate forcing and boundary conditions are time-independent, a
337 longer run reveals ongoing speed variability associated to subglacially-driven sliding
338 cyclicity; compare :cite:`vanPeltOerlemans2012`.
339
340 The ice dynamics parameters chosen in ``spinup.sh`` and ``century.sh``, especially the
341 combination
342
343 .. code-block:: none
344
345 -topg_to_phi 15.0,40.0,-300.0,700.0 -till_effective_fraction_overburden 0.02 \
346 -pseudo_plastic -pseudo_plastic_q 0.25 -tauc_slippery_grounding_lines
347
348 are a topic for a parameter study (compare :cite:`AschwandenAdalgeirsdottirKhroulev`) or a
349 study of their relation to inverse modeling results (e.g. :cite:`Habermannetal2013`).
350
351 Plotting results
352 ----------------
353
354 :numref:`fig-jako-csurf` was generated using pypismtools_, NCO_ and CDO_. Do
355
356 .. code-block:: none
357
358 ncpdq -a time,z,y,x spunjako_0.nc jako5km.nc
359 nc2cdo.py jako5km.nc
360 cdo remapbil,jako5km.nc Greenland_5km_v1.1.nc Greenland_5km_v1.1_jako.nc # FIXME: if fails, proceed?
361 ncap2 -O -s "velsurf_mag=surfvelmag*1.;" Greenland_5km_v1.1_jako.nc \
362 Greenland_5km_v1.1_jako.nc
363 basemap-plot.py -v velsurf_mag --singlerow -o jako-velsurf_mag.png jakofine.nc \
364 Greenland_5km_v1.1_jako.nc
365
366 To choose a colormap ``foo.cpt`` add option ``--colormap foo.cpt`` in the last command.
367 For this example ``PyPISMTools/colormaps/Full_saturation_spectrum_CCW.cpt`` was used.
368
369 .. rubric:: Footnotes
370
371 .. [#] PISM can also do 1 km runs for the whole Greenland ice sheet; see this `news item
372 <https://pism-docs.org/wiki/doku.php?id=news:first1km>`_.
373
374 .. [#] If this file is already present then no actual download occurs, and preprocessing
375 proceeds. Thus: Do not worry about download time if you need to preprocess again.
376 The same comment applies to other downloaded files.
377
378 .. [#] This recommendation is somewhat Jakobshavn-specific. We want our model to have an
379 ice-free down flow (western) boundary on the resulting computational domain for the
380 modeled region.
381
382 .. [#] Calculate ``620/5 + 1`` and ``425/5 + 1``, for example.