Submitted as: model description paper 12 Jan 2021
Submitted as: model description paper  12 Jan 2021
icepack: a new glacier flow modeling package in Python, version 1.0
 ^{1}Polar Science Center, Applied Physics Laboratory, University of Washington, Seattle, WA, USA
 ^{2}Department of Earth and Space Sciences, University of Washington, Seattle, WA, USA
 ^{1}Polar Science Center, Applied Physics Laboratory, University of Washington, Seattle, WA, USA
 ^{2}Department of Earth and Space Sciences, University of Washington, Seattle, WA, USA
Abstract. We introduce a new software package called icepack for modeling the flow of glaciers and ice sheets. Icepack is built on the finite element modeling library Firedrake, which uses the domainspecific language UFL for describing weak forms of partial differential equations. The diagnostic models in icepack are formulated through action principles that are specified in UFL. The components of each action functional can be substituted for different forms of the user’s choosing, which makes it easy to experiment with the model physics. The action functional itself can be used to define a solver convergence criterion that is independent of the mesh and requires little tuning on the part of the user. Icepack includes the 2D shallow ice and shallow stream models. We have also defined a 3D hybrid model based on spectral semidiscretization of the BlatterPattyn equations. Finally, icepack includes a GaussNewton solver for inverse problems that runs substantially faster than the BFGS method often used in the glaciological literature. The overall design philosophy of icepack is to be as usable as possible for a wide a swath of the glaciological community, including both experts and novices in computational science.
Daniel Shapero et al.
Status: final response (author comments only)

RC1: 'Comment on gmd2020419', Douglas Brinkerhoff, 09 Feb 2021
 AC1: 'Reply on RC1', Daniel Shapero, 07 Apr 2021

RC2: 'Comment on gmd2020419', Anonymous Referee #2, 19 Feb 2021
Please see the attachment containing my detailed review.
 AC2: 'Reply on RC2', Daniel Shapero, 07 Apr 2021
 AC4: 'Reply on RC2', Daniel Shapero, 07 Apr 2021

CC1: 'Community comment on gmd2020419', Ed Bueler, 19 Feb 2021
Summary.
Icepack is an important new model, and a description paper is appropriate. Four features of Icepack stand out: its use of Firedrake/Python, the flexible actionprinciple design of its stressbalance solver module, its fromthestart attention to data assimilation and inverse modeling, and its design as a modeling environment and language instead of a readytorun model. All of these choices are addressed appropriately in this manuscript. Readers are very likely to try the model, which fullfills a major purpose of a description paper at this earlyish development stage.
However, this description paper can be improved in three significant ways:
A. Greater attention to the meaning of Icepack as a *timedependent and geometryevolving model*, and to related defaults.
B. Avoidance of *bad linguistic habits* inherited (mostly) from the ice sheet modeling literature.
C. As currently laid out, the paper treats the reader mostly as a potential codeveloper of the model, while Icepack's *effectiveness as a simulation tool for science* is muddled.
Fully addressing concern A would require major code extensions, which are not my intention in commenting on this. Rather, the manuscript should make the reader aware of which evolution aspects are wellhandled by current Icepack and which are in future development.
The above leading concerns will be addressed in more detail below, in a list which addresses specific line numbers. The associated potential improvement "A","B","C" is listed when appropriate.
Linebyline comments.13, AC: Presumably all software packages can be effective in the hands of experts, and here it is suggested that Icepack standsout as better for nonexperts. My question is whether it will be effective for nonexperts interested in announced goals 1 and 4 (among those listed on lines 1518). This is hard to believe given how Icepack seems to be designed around, and the manuscript focused on, goals 2 and 3. Consider the reader who wants to simulate glacier extent for some years into the future for a mountain glacier or Greenland, who can supply a simulated climate (atmosphere/ocean), a bed topography, and a current geometry in some data files. Does this paper convince me that they need *less* expert knowledge to use Icepack for that purpose than other existing ice sheet models? Not yet. Of course a full usage answer occurs in online tutorials and examples, not just the manuscript. Nonetheless the absence of attention to transport equation boundary conditions, and to how mass/energy surface inputs are handled, gives me the impression that such a reader is on a long codevelopment path with the Icepack authors, requiring the development of much expert knowledge before first useful results. Said a different way, expert knowledge is required for any software that does not have an aggressive scheme for putting reasonable defaults into the hands of novices. To paraphrase a recent Firedrake paper [Farrell et al 2020], it is a mechanismvspolicy concern. Icepack has a library of mechanisms, but an absence of apparent policy means expert knowledge is needed to recover usability for realscience applications. I rolled my eyes at the implications of the sentence on line 13, and these concerns remained after reading the whole paper.
2139: I also think Python+Firedrake+PETSc is the most promising environment to build a new ice flow simulation library/model. I'm on board!
4446, B: Let me vote to *not* maintain the "diagnostic"/"prognostic" linguistic tribalism. The world calls these equations "conservation of momentum" and "conservation of mass". (Indeed one should remind the reader that the latter describes thickness evolution because of how glacier models normally parameterize fluid geometry.) An "also known as the 'prognostic equation' [cite]" is appropriate, but another paper using this tribal language will cause yet more students to need to unlearn silly language in order to read the mainstream fluids and numerical PDE literature.
4647: Here! Here! It is a infinitedimensional DAE! Good. Many readers will be unfamiliar with the concept; I cite [Ascher & Petzold 1998] for that but there may be better references.
56, equation (1) and nearby, AB: Two concerns. First, it is later acknowledged that this is not really an advection equation (lines 376377), so one does not need to call it that here either. The SIA is not some weird alternative universe of weakwilled modelers, it is what *all models should produce* in the large. That model, the only clearlyunderstood coupled model, makes equation (1) a diffusion, as noted. Surely calling this a "transport equation" or even "thickness transport equation" is adequate. One then pointsout that q=hu is one way to parameterize flux, and that u comes from a coupled equation *driven by \nabla s* (even in the Stokes case). It might be acknowledged that (1) does not have a PDE "type" in the classical sense. (It is a DAE, after all.) Second, equation (1) holds with what boundary conditions? This manuscript maintains the tradition of pretending not to notice. To quote [Schoof & Hewitt 2013], "A sometimes weakly perceived point in glaciology is that the model above is in fact a freeboundary problem ...", and *this applies to (1) regardless of the stress balance model*. Is this paper just going to pretend boundary conditions for the main, and first, equation don't exist? (Or pretend that all glaciers end in cliffs of sufficient height so that the singular change of coordinates (16) is no problem?) This manuscript could even cite existing wellposedness literature for (1), especially [Calvo et al 2003] and [Jouvet & Bueler 2012], which confirm the meaningfulness of the freeboundary problem formulation and the resultant simultaneous Dirichlet (h=0) and Neumann (q=0) boundary conditions at grounded margins. And then say what Icepack will when solving (1).
68, B: It is not clear why the conservation of momentum model for a slow fluid should ever be called a "momentum transport" model. There is no d/dt for momentum! Surely "conservation of momentum equation" or "stress balance equation" suffice.
74, B: No, the BlatterPattyn (BP) equation is not a "hybrid". This word is perfectly descriptive when categorizing ice sheet momentum models, and I am not sure who decided to mangle the language in this way. (It precedes this manuscript.) There are principled models, of which SIA, SSA, and BP are three, arising from small parameter arguments. There were horizontal hybrids (Ritz) and there still are vertical hybrids (Pollard, Bueler, Winkelmann). Hybrids, by definition, combine more than one principled model in some (less principled!) manner. If you are trying to say that BP is a shallow model, but lessshallow than SIA and SSA ... that's true! And if you are saying that BP can balance stresses in a plug or a shear flow, that's true too. Say those things.
101102, B: Regarding the sentence "The action principle can be viewed as a consequence of the Onsager reciprocity relations for systems near to equilibrium ...": I am trying to imagine a reader for whom reading that would be a useful learning experience regarding the action principle, as opposed to an irritating one.
108: I do not understand the intended meaning of "vastly more convenient numerically" here. The cost of evaluating the objective, versus a merit function like the square of the residual norm, is the same. The cost of evaluating a residual is the same whether or not it happens to be a gradient of an objective. Perhaps what is meant is more like: "Minimizing a convex functional provides better convergence guarantees than solving a general system of nonlinear equations ...", which is certainly true. At least a forward reference to section 4.2 would help explain the phrase.
115, equation (5), B: Now the bad linguistic habits are being inherited from the FE world. This is not "mass" in any sense other than that any matrix with entries a_ij = \int psi^i psi^j dx is always called a "mass matrix" in the FE world. And the reader may not have that FE bad habit already, in which case confusion ensues because (5) has nothing to do with mass conservation! My suggestion for this term, implicitly acknowledging how it arises in the SIA, is to call it "localization".
115124: The action combining (5) and (6) is already convex and coercive *in L^2*. The reason to add the penalty term is (presumably) because the FE space choices want to work in H^1. (Adding (7) with ell>0 makes J coercive in H^1.) This penalized form is a perfectly reasonable idea, and it makes the SIA behave more like other stress balances, *and* it is one of those unprincipled things one does to get it all to work properly ... if you add another kludge you could even call yourself a "hybrid".
126: This author's name is Bueler not Beuler. (The latter is a PETSc ts_type, so it is easy to get confused.)
137140: There is a double negative in equations (12) and (14) which is not used in (6) and (8), respectively. Does this reflect anything important?
148149: It is of course true that for general boundary conditions the "shallow stream equations do not have a simple analytical solution." What is probably meant here, however, is that there is not a wellknown, exactlysolvable, basalfrictionincluded boundary value problem suitable for testing. But that's not true, and indeed the *very first theoretical paper in glaciology* provides one, namely [Böðvarsson 1955]! See the full story, and a derived marine flowline exact solution suitable for testing, in [Bueler 2014]. The method of manufactured solutions applies, of course, but at loss of clarity on the meaning and magnitude of the resulting numerical errors, and with great danger of testing the wrong parts of a (nonlinear) system phase space, and with loss of the history of mathematical glaciology.
158: I guess the phrase "individual fastflowing glacier" arises because of the sense that the SSA does not handle slow flow very well, which is true, and that fast ice is separated by slow flow. Nonetheless, I would replace "features of an individual fastflowing glacier" by "features of fastflowing glaciers".
170171: I am not clear on "a tensor product basis of Lagrange finite elements in the vertical and higherorder polynomials in a single vertical layer". Should the first "vertical" actually be "horizontal"? The novice reader not already used to FE stuff could use a figure or more words here, I suspect.
175, equation (16), A: From here forward, in the manuscript, the lateral boundary condition for equation (1) is even less clear because this change of variables takes the lateral free margin and blows it up. (In the grounded ice case, and where "blow up" is used in the mathematical sense.) With what consequences for simulating moving margins? Section 4.1 does not address this (i.e. how the boundary condition for the transport equation behaves as h>0; a boundary condition is not mentioned). I believe, on the other hand, that the SIA groundedmargin example in section 5.2 does not use (16), which as described applies only applying to "the hybrid flow model", but section 5.2 does not mention the margin anyway. I am concerned that positive consequences of using (16) are emphasized while negative consequences are not mentioned, and that Icepack has no path around these negative consequences.
215: For clarity to readers not already thinking about Firedrake's extruded meshes, I suggest adding a parenthetical: "linear or quadratic elements on (horizontal) triangles".
216: Instead of "Lagrange interpolating polynomial basis", which readers accustomed to numerical analysis terminology (outside FE) might associate to the Lagrange method of finding interpolating polynomials, it is probably clearer to just say "Lagrange elements". Perhaps emphasize the nodal placement? (Legendre nodes, concentrated near endpoints, versus the equallyspaced nodes of Lagrange elements.) Orthogonality is not the only reason why Legendre is better; node location implies better approximation properties.
227: "Smoke test" is an unfamiliar idiom to me. As in testing machinery or testing electronics or smokingout bugs or whatamIsmoking? Pity the nonnative speaker.
227: "most minimalistic" > "minimal"
234235: I have received the advice "write about what you have done, not what you haven't", and I think it applies here. Using two layers in this way assumes one polythermal structure and will not allow transition to the other. No need to consider or mention it.
253289: This chunk of text is important, and it describes one of the best aspects of Icepack's use of Python. I think it should be put at the start of subsection 2.3, which will give clarity that when equations (22) and (23) appear, they are primarily an example of submodel plugability.
283289: I'm on board with keyword arguments!
298299: Suggested replacement for the mushy phrase: "describing realistic glacier flows, but it is not sufficient by itself" > "parameterizing fluidity".
300301, C: I seriously doubt that more than one user in 100 will want to substitute their own meltfraction dependence, because no one has the data for it. So: "likely" > "possible". (This is a clear case where good defaults are more important than modularity.)
297305, C: This paragraph reads, to me, like an argument that expert knowledge *will* be required for the Icepack user. ("users must calculate the fluidity field themselves") That is, unless Icepack includes both defaults and higherlevel simulation drivers which allow the novice user to never know about these possibilities. More broadly, what does Icepack design offer, when it comes to avoiding explicit choices of all submodels, on the way to a first scientificallyuseful model for the beginner?
318320: Mushy sentences. Needed?
345, C: The sentence "The key classes that users interact with are flow models and solvers" is right at the heart of my concern C. Indeed, after developing PISM with much (developer group) concern about flow models and solvers, we spent the next 10 years talking to users about their data formats, parameterstudy schema, and surface process models (supra, sub, and calving front). Barely had time to work on flow models and solvers ... until I quit daytoday PISM development. The byfar largest population of scientists care about how flowing ice interacts with their climates. They use ice flow models on the assumption that ice flow modelers know how to model ice flow! This "key classes" sentence is not describing users, it is describing codevelopers.
368: "Practicing ... may be unfamiliar with" is unnecessary. Just say what you have to say.
362381, A: Three concerns about section 4.1. First, as already noted, transport equation (1) has boundary values and they are again silent here. Second, there are multiple transport problems in glacier modeling: equation (1), the advectiondominated enthalpy equation, the ice shelf damage transport equation, the age equation, ... etc. Is this paragraph covering them all, or just equation (1)?
368373, AC: The third concern, regarding this paragraph, is that Icepack apparently has not yet setup effective adaptive timestepping. (I have not gotten a runtime error from a modern explicit ODE solver for a long time! Have you? Some runs are slow cause they are stiff, indeed.) This should be addressed/acknowledged. In particular, assuming a perfectlyimplemented *implicit* solver, how does \Delta t scale with \Delta x? If the problem were really advective, so that the goal is to model influences as they travel at the characteristic speed, then the answer would be \Delta t = O(\Delta x) for accuracy, even with your implicit scheme. (Better than O(\Delta x^2), yes indeed ... but then explicit steppers would be fine ... but (1) is not an advection ... we must think more.) Though the problem is actually some diffusive/advective mix, some defined scaling of \Delta t with \Delta x is still needed for accuracy. This can come from an adaptive timestepper, or be designed from scratch. Then there is the matter of the user's data time scale for surface mass balance, etc., which implies "events" in the adaptive timestepper (when data is read). In fact, the actual questions a glacier model must answer, regarding time steps, are these: How does timestepping change under spatial grid refinement? Does it converge in spacetime refinement? Is it robust over realistic geometries and inputs? What is the user interface? This paragraph convinces me that Icepack is not yet there, which should not be completely buried.
363: "equation equation"
365, A: The "In the interest of simplicity" phrase here tells me the authors simply have not thoughtthrough the timeevolution of glacier geometry and the needed boundary conditions. (You have a scheme for unconditionalstably generating glacier surfaces which conserves mass in the presence of ablation? Then don't hide it!) Please take the problem seriously: Address how you maintain reasonable margin shape and positive thickness in an implicit scheme. It is o.k. to admit that time steps cannot be arbitrarily long in your scheme, to get convergence; the "conditional/unconditional stable" language is an artifact of linear PDE theory, and your problem, taken seriously, is superdupernonlinear.
376377, B: This sentence is a very good argument for *not* calling (1) an advection equation.
383444, A: Sections 4.2 and 4.3 convince me that when the Icepack team takes on aspects of model design they care about then it comes out very well! These sections suggest the momentum balance ("diagnostic" ... grrr) solver has great defaults. Likewise with sections 4.4 and 4.5. (Now for serious attention to timestepping, mass conservation, and the user's surface mass and energyflux data, and clarity on Icepack's TODO list.)
459462: GaussNewton is a good choice for this purpose, as pointed out for exactly this purpose, inversion of glacier models, by [Habermann et al 2012]. Which should be cited.
463466: Symbolic derivatives. Again I am on board with the benefits of a Firedrakebased tool chain.
469, A: Here "terminus boundary condition" refers to the momentum balance equation. There is, as far as I can tell, no regard for the "terminus boundary condition" needed when solving equation (1), i.e. the mass conservation terminus boundary condition, even in the case of a cliff, much less a grounded margin.
474487: Extruded meshes. On board with Firedrake. But if the model claims to have a solution of (1) then there would be some mechanism for addressing regions which become icefree within a time step? No mention thereof.
493515, A: Section 4.5 describes a nice solution to a real problem. In fact, can one consider this kind of fit to the glacier profile at a grounded margin? Not obvious how to proceed, I agree, but note that the thickness field h of a glacier has the similar continuousbutnondifferentiable character as the pressure on the ocean calving front.
517522: Rich problemaware preconditioners and solvers. On board with Firedrake.
523, A: This honest, clear sentence is not at all expected by the reader who remembers lines 1519 in the Introduction. Please let the reader know earlier.
533586: The MISMIP+ experiment is part of the 2008onward tradition of running ice sheet models in rectangular boxes (MISMIP, EISMINTHOM, ...). This is fine in a field where one also has some laboratory fluid models, i.e. actual stuff, which fit in boxes. But this tradition has a distressing impact on new model development, which is to modularize around the choices one makes between intercomparisonspecified boxes and their associated boundary conditions. (Yes, that is what Icepack looks like right now.)
588597, A: This example would be more impressive with hybrid physics, right?, but Icepack is not there yet. (Or else it would have been applied here.) Readers of GMD should be taken to be serious people. Tell them the score, and I don't mean relative to what other models are capable of. What needs to improve?
598: Regarding "computationally cheap enough ... in a matter of minutes on a desktop": I suppose the excuse for this sentence is that other people get away with writing such stuff? The question is how run time scales with mesh resolution. The reader can handle a plot, if you want to generate one.
604647: Nice examples! Ice shelf physics is a scope where Icepack is a convincing choice for a research project.
662663: Regarding "The person learning ... is largely absent from the discussion", have you looked at the PISM User's Manual? (Start from page 1.) This opinion ranges from disputable to insulting, but mostly reflects not paying too much attention to unpublished prior literature (i.e. online manuals). On exactly that topic, the online Icepack documentation is excellent.
666713: Saying these HCI principles out loud is a very worthwhile aspect of the manuscript.
677678, C: There is a *big* gap between "user interface in an interpreted language" and a "program that can only run in batch mode". Cbased PETSc programs like the PISM ice sheet model, and many other commandline programs, live in the interior of this gap. For example, please consider these three HCI principles as they apply to commandline git. It is neither "a user interface in an interpreted program" nor a "batch mode" only thing. (Most commandline tools are not!) Git has a steep learning curve, because it is a DAG modeling language, but such commandline design can hit your principles too. Indeed, a good antidote to your false dichotomy is [Brown et al 2014], and addressing UI points made there would increase the credibility of section 6. (You'll see that Icepack is definitely a Brownapproved design librarywise, but doing science is not equal to designing a library API either.) Ultimately any science application of Icepack will be a map from inputs (data) to outputs (simulations or inversions), and the usability of that map (e.g. in ensemble simulations) is different from the developability of it in an interactive environment. Python is a great environment for experimentation, and for ice sheet modeling, and progressive evaluation is a desirable principle, but an interactive Python session is not the only alternative to 1990s climatemodel design. Presumably Icepack usage is intended to progress from allinteractive mode to parallel production/ensemble runs in batch systems anyway? Address that?
681683, C: Regarding "The API ... one step at a time": There is a *huge* amount of expert knowledge required to do science with an "ice sheet model" which does not have a policy for doing things "one step at a time". Imagine a paper that proposed a newandbetter WRF model but said something like this about the demands on the user!
682683, C: Also, this "user" is a codeveloper.
693713, C: This basic point about abstraction gradient is good. But solver/discretization components/choices is not the only such gradient. From talking to a lot of novice glacier modelers, I assert the key abstraction gradient, which an ice sheet model must finesse, regards ice flow physics. (What aspects of the full, coupled, nonlinear dynamics are the next ones that the user's constructed Icepack model should/can handle for the science goal? How to buildin that physics without unnecessary parameters? How to generate intermediate results which reveal which processes are missing, re the science goal?) Thinking on this stuff is where Icepack developers could make their next real progress. That means deemphasizing solver components/discretizations in the user's view, primarily by setting aggressive defaults, even as the developers must get solvers right (which is the strength of Icepack).
721: "phyics"
724: Having reviewed quite a few ice sheet modeling papers, every single one claimed something about its usability. Less mush, more answers please.
725: I'm on board with (2), and wish for more careful analysis of Icepack's contribution here, but "quantifying" in (1) is not what you are doing in this manuscript. (Nor, probably, do you have the ability to do it.)
732740, AC: If I saw "4. robust timeevolution and climate interaction tools", or similar, on this list then I would be more of a believer in Icepack's future. Will you be able to break out of stressbalanceandinversemodelsolverplayground mode, and start answering some of the science questions enumerated in the Introduction?
ReferencesAscher, U. and Petzold, L. (1998). Computer Methods for Ordinary Differential Equations and DifferentialAlgebraic Equations, SIAM, Philadelphia.
Böðvarsson G (1955) On the flow of icesheets and glaciers. Jökull,
5, 1–8 (find this toolittleread paper at: https://github.com/bueler/bodmarine/blob/master/Bodvardsson1955.pdf)Brown, J., Knepley, M. G., & Smith, B. F. (2014). Runtime extensibility and librarization of simulation software. Computing in Science & Engineering, 17(1), 3845.
Bueler, E. (2014). An exact solution for a steady, flowline marine ice sheet. Journal of Glaciology, 60(224), 11171125.
Calvo, N., Díaz, J. I., Durany, J., Schiavi, E., & Vázquez, C. (2003). On a doubly nonlinear parabolic obstacle problem modelling ice sheet dynamics. SIAM Journal on applied mathematics, 63(2), 683707.
Farrell, P. E., Kirby, R. C., & MarchenaMenendez, J. (2020). Irksome: Automating RungeKutta timestepping for finite element methods. arXiv preprint arXiv:2006.16282.
Habermann, M., Maxwell, D., & Truffer, M. (2012). Reconstruction of basal properties in ice sheets using iterative inverse methods. Journal of Glaciology, 58(210), 795808.
Jouvet, G., & Bueler, E. (2012). Steady, shallow ice sheets as obstacle problems: wellposedness and finite element approximation. SIAM Journal on Applied Mathematics, 72(4), 12921314.
Schoof, C., & Hewitt, I. (2013). Icesheet dynamics. Annual Review of Fluid Mechanics, 45, 217239.
 AC3: 'Reply on CC1', Daniel Shapero, 07 Apr 2021
Daniel Shapero et al.
Model code and software
icepack: glacier flow modeling with the finite element method in Python Daniel Shapero, Jessica Badgeley, and Andrew Hoffman https://doi.org/10.5281/zenodo.4318150
Daniel Shapero et al.
Viewed
HTML  XML  Total  BibTeX  EndNote  

342  148  17  507  5  4 
 HTML: 342
 PDF: 148
 XML: 17
 Total: 507
 BibTeX: 5
 EndNote: 4
Viewed (geographical distribution)
Country  #  Views  % 

Total:  0 
HTML:  0 
PDF:  0 
XML:  0 
 1