tparameter-study.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
---
tparameter-study.rst (8125B)
---
1 .. include:: ../../global.txt
2
3 .. _sec-paramstudy:
4
5 An ice dynamics parameter study
6 -------------------------------
7
8 The readers of this manual should not assume the PISM authors know all the correct
9 parameters for describing ice flow. While PISM must have *default* values of all
10 parameters, to help users get started,\ [#]_ it has more than two hundred user-configurable
11 parameters. The goal in this manual is to help the reader adjust them to their desired
12 values. While "correct" values may never be known, or may not exist, examining the
13 behavior of the model as it depends on parameters is both a nontrivial and an essential
14 task.
15
16 For some parameters used by PISM, changing their values within their ranges of
17 experimental uncertainty is unlikely to affect model results in any important manner (e.g.
18 ``constants.sea_water.density``). For others, however, for instance for the exponent in
19 the basal sliding law, changing the value is highly-significant to model results, as we'll
20 see in this subsection. This is also a parameter which is very uncertain given current
21 glaciological understanding :cite:`CuffeyPaterson`.
22
23 To illustrate a parameter study in this Manual we restrict consideration to just two
24 important parameters for ice dynamics,
25
26 - `q=` :config:`basal_resistance.pseudo_plastic.q`: exponent used in
27 the sliding law which relates basal sliding velocity to basal shear stress in the SSA
28 stress balance; see subsection :ref:`sec-basestrength` for more on this parameter, and
29 - `e=` :config:`stress_balance.sia.enhancement_factor`: values larger
30 than one give flow "enhancement" by making the ice deform more easily in shear than is
31 determined by the standard flow law :cite:`LliboutryDuval1985`, :cite:`PatersonBudd`;
32 applied only in the SIA stress balance; see section :ref:`sec-rheology` for more on this
33 parameter.
34
35 By varying these parameters over full intervals of values, say `0.1\le q \le 1.0`
36 and `1 \le e \le 6`, we could explore a two-dimensional parameter space. But of
37 course each `(q,e)` pair needs a full computation, so we can only sample this
38 two-dimensional space. Furthermore we must specify a concrete run for each parameter pair.
39 In this case we choose to run for 1000 model years, in every case initializing from the
40 stored state ``g10km_gridseq.nc`` generated in the previous section :ref:`sec-gridseq`.
41
42 The next script, which is ``param.sh`` in ``examples/std-greenland/``, gets values
43 `q \in \{0.1,0.5,1.0\}` and `e\in\{1,3,6\}` in a double ``for``-loop. It
44 generates a run-script for each `(q,e)` pair. For each parameter setting it calls
45 ``spinup.sh``, with the environment variable ``PISM_DO=echo`` so that ``spinup.sh`` simply
46 outputs the run command. This run command is then redirected into an appropriately-named
47 ``.sh`` script file:
48
49 .. literalinclude:: scripts/run-5-study.sh
50 :language: bash
51
52 Notice that, because the stored state ``g10km_gridseq.nc`` used `q=0.5` and
53 `e=3`, one of these runs simply continues with no change in the physics.
54
55 To set up and run the parameter study, without making a mess from all the generated files,
56 do:
57
58 .. literalinclude:: scripts/run-5-setup.sh
59 :language: bash
60 :lines: 3-
61
62 The result of the last command is to generate nine run scripts,
63
64 .. code-block:: none
65
66 p10km_0.1_1.sh p10km_0.1_3.sh p10km_0.1_6.sh
67 p10km_0.5_1.sh p10km_0.5_3.sh p10km_0.5_6.sh
68 p10km_1.0_1.sh p10km_1.0_3.sh p10km_1.0_6.sh
69
70 The reader should inspect a few of these scripts. They are all very similar, of course,
71 but, for instance, the ``p10km_0.1_1.sh`` script uses options ``-pseudo_plastic_q 0.1``
72 and ``-sia_e 1``.
73
74 .. figure:: figures/ivol-param.png
75 :name: fig-ivolparamstudy
76
77 Time series of ice volume ``ice_volume_glacierized`` from nine runs in our parameter study
78 example, with parameter choices `(q,e)` given.
79
80 We have not yet run PISM, but only asked one script to create nine others. We now have the
81 option of running them sequentially or in parallel. Each script itself does a parallel
82 run, over the ``NN=4`` processes specified by ``param.sh`` when generating the run
83 scripts. If you have `4 \times 9 = 36` cores available then you can do the runs
84 fully in parallel (this is ``runparallel.sh`` in ``examples/std-greenland/``):
85
86 .. literalinclude:: scripts/run-5-parallel.sh
87 :language: bash
88
89 Otherwise you should do them in sequence (this is ``runsequential.sh`` in
90 ``examples/std-greenland/``):
91
92 .. literalinclude:: scripts/run-5-serial.sh
93 :language: bash
94
95 On the same old 2012-era 4 core laptop, ``runsequential.sh`` took a total of just under 7
96 hours to complete the whole parameter study. The runs with `q=0.1` (the more
97 "plastic" end of the basal sliding spectrum) took up to four times longer than the
98 `q=0.5` and `q=1.0` runs. Roughly speaking, values of `q` which are
99 close to zero imply a subglacial till model with a true yield stress, and the result is
100 that even small changes in overall ice sheet state (geometry, energy, ...) will cause
101 *some* location to exceed its yield stress and suddenly change flow regime. This will
102 shorten the time steps. By contrast, the `e` value is much less significant in
103 determining run times.
104
105 .. figure:: figures/p10km-comparison.png
106 :name: fig-paramstudy
107
108 Surface speed ``velsurf_mag`` from a 10 km grid parameter study. Right-most subfigure
109 is observed data from ``Greenland_5km_v1.1.nc``. Top row: `q=0.1` and
110 `e=1,3,6` (left-to-right). Middle row: `q=0.5`. Bottom row: `q=1.0`.
111 All subfigures have common color scale (velocity m/year), as shown in the right-most
112 figure, with 100 m/year contour shown in all cases (solid black).
113
114 On a supercomputer, the ``runparallel.sh`` script generally should be modified to submit
115 jobs to the scheduler. See example scripts ``advanced/paramspawn.sh`` and
116 ``advanced/paramsubmit.sh`` for a parameter study that does this. (But see your system
117 administrator if you don't know what a "job scheduler" is!) Of course, if you have a
118 supercomputer then you can redo this parameter study on a 5 km grid.
119
120 Results from these runs are seen in :numref:`fig-ivolparamstudy` and
121 :numref:`fig-paramstudy`. In the former we see that the `(0.5,3)` run simply continues the
122 previous initialization run. In some other graphs we see abrupt initial changes, caused by
123 abrupt parameter change, e.g. when the basal sliding becomes much more plastic (`q=0.1`).
124 In all cases with `e=1` the flow slows and the sheet grows in volume as discharge
125 decreases, while in all cases with `e=6` the flow accelerates and the sheet shrinks in
126 volume as discharge increases.
127
128 In :numref:`fig-paramstudy` we can compare the surface speed model results to
129 observations. Roughly speaking, the ice softness parameter `e` has effects seen
130 most-clearly by comparing the interior of the ice sheet; scan left-to-right for the
131 `e=1,3,6` subfigures. The basal sliding exponent `q` has effects seen
132 most-clearly by comparing flow along the very steep margin, especially in the southern
133 half of the ice sheet; scan top-to-bottom for `q=0.1,0.5,1.0`, going from
134 nearly-plastic at top to linear at bottom.
135
136 From such figures we can make an informal assessment and comparison of the results, but
137 objective assessment is important. Example objective functionals include:
138
139 #. compute the integral of the square (or other power) of the difference between the model
140 and observed surface velocity :cite:`AschwandenAdalgeirsdottirKhroulev`, or
141 #. compute the model-observed differences between the histogram of the number of cells
142 with a given surface speed :cite:`BKAJS`.
143
144 Note that these functionals are measuring the effects of changing a small number of
145 parameters, namely two parameters in the current study. So-called "inversion" might use
146 the same objective functionals but with a much larger parameter space. Inversion is
147 therefore capable of achieving much smaller objective measures :cite:`Habermannetal2013`,
148 :cite:`Larouretal2012`, :cite:`Priceetal2011`, though at the cost of less understanding,
149 perhaps, of the meaning of the optimal parameter values.
150
151 .. rubric:: Footnotes
152
153 .. [#] See :ref:`sec-parameter-list` for the full list.