the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 25 Jun 2025
MinSIA v1: a lightweight and efficient implementation of the shallow ice approximation
Abstract. Simulations of ice flow have recently been boosted to an unprecedented numerical performance by machine learning techniques. This paper aims at keeping classical numerics competitive in this field. It introduces a new numerical scheme for the shallow ice approximation. Key features are a semi-implicit time-stepping scheme in combination with a dynamic smoothing of the nonlinearity in the slope-dependence of the flow velocity. As a first step, the software MinSIA presented here provides a lightweight implementation of the new scheme in MATLAB. An implementation in Python is under development. MinSIA is designed for simulations with several million nodes on standard desktop PCs and allows for spatial resolutions of 25 m or even finer. The numerical scheme performs particularly well for heavily glaciated topographies with moderately inclined ice surfaces. In turn, the advantage of the scheme decreases slightly for alpine topographies with steep walls during phases of moderate glaciation.
Competing interests: The author is a member of the editorial board of Geoscientific Model Development. The peer-review process was guided by an independent editor, and the author also has no other competing interests to declare.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors. Views expressed in the text are those of the authors and do not necessarily reflect the views of the publisher.-
Notice on discussion status
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
-
Preprint
(13098 KB)
-
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
- Preprint
(13098 KB) - Metadata XML
- BibTeX
- EndNote
- Final revised paper
Journal article(s) based on this preprint
Interactive discussion
Status: closed
-
RC1: 'Comment on egusphere-2025-2242', Daniel Moreno-Parada, 31 Jul 2025
-
AC1: 'Reply on RC1', Stefan Hergarten, 26 Aug 2025
Dear Daniel Moreno-Parada,
thank you very much for you review! I appreciate your comments very much in times when most of the reviewers comment only on the way how the introduction is written. On the other hand, it looks as if you want a paper different from what I want, which makes it a bit complicated. Since you made several comments, I uploaded my preliminary repsonses as an attachment.
Best,
Stefan
-
AC1: 'Reply on RC1', Stefan Hergarten, 26 Aug 2025
-
RC2: 'Comment on egusphere-2025-2242', Thomas Zwinger, 29 Aug 2025
Please find my comments on the presented manuscripts in the attached PDF
Kind Regards,
Thomas Zwinger
-
AC2: 'Reply on RC2', Stefan Hergarten, 03 Sep 2025
Dear Thomas Zwinger,
thank you very much for your thorough review!
For me, the central point is the "applicability" of the SIA to fine grids since it defines whether or not the entire exercise makes sense. Therefore, I want to discuss this point here, while I uploaded the detailed responses as a separate document.
Your statements about the SIA at fine resolutions confuse me. One could even interpret them in such a way that the SIA became increasingly wrong by itself if we use finer grids. I guess that this is not what you want to say.
The mathematical analyses, typically using a vertical/horizontal scaling factor epsilon refer to the situation that we expand or compress the domain (and the topography) horizontally. If we compress it and make the valleys narrower, horizontal velocity gradients increase and the SIA becomes worse. So far clear, but we cannot identify the spatial resolution of the grid with the horizontal length scale in the denominator of epsilon. If we cover the same topography with a finer grid (without introducing additional roughness), nothing happens.
For a given topography, the SIA is good if the ice thickness is small (what the SIA is made for). It becomes worse with increasing thickness and finally becomes better again when the ice overflows the small valleys. If the valleys are not wide enough in relation to their depth, horizontal stresses become relevant in nature. Then the SIA overestimates the velocity in the middle of the valley and underestimates it at the sides. This effect is basically independent of the spatial resolution for fine grids, but the accuracy does not automatically become better if we decrease the resolution so that we cannot resolve the valleys any more. Therefore, I disagree to your statement that an SIA appropriate mesh resolution would be > 1000 m. From my point of view, the typical test examples include regions (valleys) where the SIA is sufficient and regions where it is inaccurate.
As a central point, however, I think that we must distinguish properly between the physical deficiencies of the SIA (narrow valleys, strong sliding, ...) and numerical issues. The SIA leads to some kind of nonlinear diffusion equation (although the linearized diffusivity is 3D instead of D, see detailed responses). Diffusion equations dampen spatial undulations exponentially. If L is the wavelength of the undulation, the respective decay (damping) constant is D/(2*pi*L)^2. So short-wavelength undulations are damped very rapidly by the SIA itself. Tentatively, I would say that the assumption of hydrostatic pressure even makes damping too fast (since it affects the entire vertical column), but this is practically not relevant.
Numerically, however, it is opposite. Explicit schemes cannot capture the fast damping, but overshoot and generate growing oscillations from small undulations. The hydrostatic approximation makes these oscillations even worse. But in any case, this is a deficiency of the numerical schemes and not a deficiency of of the SIA itself, as I would read from your comments.
Summarizing the first point, I strongly disagree to your points that the SIA becomes worse with increasing spatial resolution (still thinking that this is not what you wanted to say) and that the numerical issues with the SIA are a deficiency of the SIA itself and not of the numerical scheme. Therefore, I still believe that the new numerical scheme brings the solution closer to the exact solution of the SIA (with all its problems in narrow valleys and if sliding is strong).
There will probably be an extension by horizontal stresses in the future, but I do not know so far what this extension will look like.
As you also raise an issue about the correction for steep slopes in your major points, let me briefly comment on this here. The equations for the hydrostatic pressure and for the horizontal shear stress are the same for the original shallow water equations and for the SIA. To my knowledge, this is not coincidence, but arises because it makes no difference up to this point whether we start from Navier--Stokes or Stokes and what the rheology looks like. As far as I know, Hutter (1983, doi 0.1007/978-94-015-1167-4) already proposed a version for an inclined bed.
The central idea (beyond zero pressure at the free surface) behind this extension is that the velocity is parallel to the bed (not necessarily horizontal), which reduces the hydrostatic pressure compared to the original shallow water equations. It is good to know that this version is consistent with the slab solution (I only have the 2010 lecture notes of Ralf Greve and not the textbook, but it should be similar). I wrote a short proof of this, as requested, in the detailled responses.
In essence, I still think that the version my correction is better for steep slopes (ice surface almost parallel to the bed), while the original form may be better for glacier fronts with an almost horizontal bed (although it does not comply with the condition defined by Hutter, 1983). So I found it the best option to leave it to the users which version to take.
Best regards,
Stefan-
RC3: 'Reply on AC2', Thomas Zwinger, 05 Sep 2025
In order to react to the statement of the author in the rebuttal letter (quote): "Your statements about the SIA at fine resolutions confuse me. One could even interpret them in such a way that the SIA became increasingly wrong by itself if we use finer grids. I guess that this is not what you want to say." let me break down the essential point of my critics to avoid such guessing: My point rather was that the missing physics in SIA needed to resolve topographic features presented on length-scales well below multiples of the ice thickness cannot be overcome by running it on ever finer grids. Even if one is numerically able to produce a stable solution to SIA on such fine resolved topographies/grids, we cannot expect it to be physically correct. I think that statement is well backed by the literature I cite. Hence, from a physical point of view, I do not see the added value of 25 m resolution grids on rugged terrains like the Alps in combination to the chosen approximation explained in the manuscript. In my opinion the presented article also does not allow the reader to quantify the errors introduced by such a setup, as it only shows variations of runs that all are performed with the model itself.
Citation: https://doi.org/10.5194/egusphere-2025-2242-RC3 -
AC3: 'Reply on RC3', Stefan Hergarten, 10 Sep 2025
Dear Thomas Zwinger, thanks for your rapid response! Since your main concern did not make sense to me the way I interpreted it, the clarification is definitely helpful. I never questioned the deficiency that the SIA overestimates velocities in the center of valleys and underestimates velocities at the sides as long as the valleys are not sufficiently wide. Best regards, Stefan
Citation: https://doi.org/10.5194/egusphere-2025-2242-AC3
-
AC3: 'Reply on RC3', Stefan Hergarten, 10 Sep 2025
-
RC3: 'Reply on AC2', Thomas Zwinger, 05 Sep 2025
-
AC2: 'Reply on RC2', Stefan Hergarten, 03 Sep 2025
-
EC1: 'Comment on egusphere-2025-2242', Ludovic Räss, 08 Sep 2025
Dear authors,
The manuscript draft has its merits but would require major revisions to carefully address the concerns raised by the reviewers.
In particular, the revised version should make sure to clearly stress the benefits and limitations of the presented approach, and make sure to perform all relevant benchmarks in the regime or configuration that is most valid for the SIA framework.
Also, note that the high-resolution could also be formulated in grid points rather than in meters, which could further highlight that being able to solve SIA with many DoFs is not only useful when aiming at small grid spacing but would also allow to investigate larger domains.
Regarding convergence, SIA is a PDE and as such it should be solved with mesh convergence. Specifically, it would be relevant to provide mesh convergence experiment and report about the findings during revision. Also, running SIA on resolved vs filtered coarse topography may not show much discrepancy in terms of resolving physics, but may be undeniably better in terms of achieving converged results.
Finally, if the authors are not willing to provide a (less performant) Python version to support open-source code usage, I would remove the related bits from the manuscript's abstract.
According to the above, I encourage the authors to submit a major revision of their work for further assessment. Thanks
Ludovic Räss
GMD Topical EditorCitation: https://doi.org/10.5194/egusphere-2025-2242-EC1
Interactive discussion
Status: closed
-
RC1: 'Comment on egusphere-2025-2242', Daniel Moreno-Parada, 31 Jul 2025
Dear author,
Please, find attached my report document.
Best,
Daniel
-
AC1: 'Reply on RC1', Stefan Hergarten, 26 Aug 2025
Dear Daniel Moreno-Parada,
thank you very much for you review! I appreciate your comments very much in times when most of the reviewers comment only on the way how the introduction is written. On the other hand, it looks as if you want a paper different from what I want, which makes it a bit complicated. Since you made several comments, I uploaded my preliminary repsonses as an attachment.
Best,
Stefan
-
AC1: 'Reply on RC1', Stefan Hergarten, 26 Aug 2025
-
RC2: 'Comment on egusphere-2025-2242', Thomas Zwinger, 29 Aug 2025
Please find my comments on the presented manuscripts in the attached PDF
Kind Regards,
Thomas Zwinger
-
AC2: 'Reply on RC2', Stefan Hergarten, 03 Sep 2025
Dear Thomas Zwinger,
thank you very much for your thorough review!
For me, the central point is the "applicability" of the SIA to fine grids since it defines whether or not the entire exercise makes sense. Therefore, I want to discuss this point here, while I uploaded the detailed responses as a separate document.
Your statements about the SIA at fine resolutions confuse me. One could even interpret them in such a way that the SIA became increasingly wrong by itself if we use finer grids. I guess that this is not what you want to say.
The mathematical analyses, typically using a vertical/horizontal scaling factor epsilon refer to the situation that we expand or compress the domain (and the topography) horizontally. If we compress it and make the valleys narrower, horizontal velocity gradients increase and the SIA becomes worse. So far clear, but we cannot identify the spatial resolution of the grid with the horizontal length scale in the denominator of epsilon. If we cover the same topography with a finer grid (without introducing additional roughness), nothing happens.
For a given topography, the SIA is good if the ice thickness is small (what the SIA is made for). It becomes worse with increasing thickness and finally becomes better again when the ice overflows the small valleys. If the valleys are not wide enough in relation to their depth, horizontal stresses become relevant in nature. Then the SIA overestimates the velocity in the middle of the valley and underestimates it at the sides. This effect is basically independent of the spatial resolution for fine grids, but the accuracy does not automatically become better if we decrease the resolution so that we cannot resolve the valleys any more. Therefore, I disagree to your statement that an SIA appropriate mesh resolution would be > 1000 m. From my point of view, the typical test examples include regions (valleys) where the SIA is sufficient and regions where it is inaccurate.
As a central point, however, I think that we must distinguish properly between the physical deficiencies of the SIA (narrow valleys, strong sliding, ...) and numerical issues. The SIA leads to some kind of nonlinear diffusion equation (although the linearized diffusivity is 3D instead of D, see detailed responses). Diffusion equations dampen spatial undulations exponentially. If L is the wavelength of the undulation, the respective decay (damping) constant is D/(2*pi*L)^2. So short-wavelength undulations are damped very rapidly by the SIA itself. Tentatively, I would say that the assumption of hydrostatic pressure even makes damping too fast (since it affects the entire vertical column), but this is practically not relevant.
Numerically, however, it is opposite. Explicit schemes cannot capture the fast damping, but overshoot and generate growing oscillations from small undulations. The hydrostatic approximation makes these oscillations even worse. But in any case, this is a deficiency of the numerical schemes and not a deficiency of of the SIA itself, as I would read from your comments.
Summarizing the first point, I strongly disagree to your points that the SIA becomes worse with increasing spatial resolution (still thinking that this is not what you wanted to say) and that the numerical issues with the SIA are a deficiency of the SIA itself and not of the numerical scheme. Therefore, I still believe that the new numerical scheme brings the solution closer to the exact solution of the SIA (with all its problems in narrow valleys and if sliding is strong).
There will probably be an extension by horizontal stresses in the future, but I do not know so far what this extension will look like.
As you also raise an issue about the correction for steep slopes in your major points, let me briefly comment on this here. The equations for the hydrostatic pressure and for the horizontal shear stress are the same for the original shallow water equations and for the SIA. To my knowledge, this is not coincidence, but arises because it makes no difference up to this point whether we start from Navier--Stokes or Stokes and what the rheology looks like. As far as I know, Hutter (1983, doi 0.1007/978-94-015-1167-4) already proposed a version for an inclined bed.
The central idea (beyond zero pressure at the free surface) behind this extension is that the velocity is parallel to the bed (not necessarily horizontal), which reduces the hydrostatic pressure compared to the original shallow water equations. It is good to know that this version is consistent with the slab solution (I only have the 2010 lecture notes of Ralf Greve and not the textbook, but it should be similar). I wrote a short proof of this, as requested, in the detailled responses.
In essence, I still think that the version my correction is better for steep slopes (ice surface almost parallel to the bed), while the original form may be better for glacier fronts with an almost horizontal bed (although it does not comply with the condition defined by Hutter, 1983). So I found it the best option to leave it to the users which version to take.
Best regards,
Stefan-
RC3: 'Reply on AC2', Thomas Zwinger, 05 Sep 2025
In order to react to the statement of the author in the rebuttal letter (quote): "Your statements about the SIA at fine resolutions confuse me. One could even interpret them in such a way that the SIA became increasingly wrong by itself if we use finer grids. I guess that this is not what you want to say." let me break down the essential point of my critics to avoid such guessing: My point rather was that the missing physics in SIA needed to resolve topographic features presented on length-scales well below multiples of the ice thickness cannot be overcome by running it on ever finer grids. Even if one is numerically able to produce a stable solution to SIA on such fine resolved topographies/grids, we cannot expect it to be physically correct. I think that statement is well backed by the literature I cite. Hence, from a physical point of view, I do not see the added value of 25 m resolution grids on rugged terrains like the Alps in combination to the chosen approximation explained in the manuscript. In my opinion the presented article also does not allow the reader to quantify the errors introduced by such a setup, as it only shows variations of runs that all are performed with the model itself.
Citation: https://doi.org/10.5194/egusphere-2025-2242-RC3 -
AC3: 'Reply on RC3', Stefan Hergarten, 10 Sep 2025
Dear Thomas Zwinger, thanks for your rapid response! Since your main concern did not make sense to me the way I interpreted it, the clarification is definitely helpful. I never questioned the deficiency that the SIA overestimates velocities in the center of valleys and underestimates velocities at the sides as long as the valleys are not sufficiently wide. Best regards, Stefan
Citation: https://doi.org/10.5194/egusphere-2025-2242-AC3
-
AC3: 'Reply on RC3', Stefan Hergarten, 10 Sep 2025
-
RC3: 'Reply on AC2', Thomas Zwinger, 05 Sep 2025
-
AC2: 'Reply on RC2', Stefan Hergarten, 03 Sep 2025
-
EC1: 'Comment on egusphere-2025-2242', Ludovic Räss, 08 Sep 2025
Dear authors,
The manuscript draft has its merits but would require major revisions to carefully address the concerns raised by the reviewers.
In particular, the revised version should make sure to clearly stress the benefits and limitations of the presented approach, and make sure to perform all relevant benchmarks in the regime or configuration that is most valid for the SIA framework.
Also, note that the high-resolution could also be formulated in grid points rather than in meters, which could further highlight that being able to solve SIA with many DoFs is not only useful when aiming at small grid spacing but would also allow to investigate larger domains.
Regarding convergence, SIA is a PDE and as such it should be solved with mesh convergence. Specifically, it would be relevant to provide mesh convergence experiment and report about the findings during revision. Also, running SIA on resolved vs filtered coarse topography may not show much discrepancy in terms of resolving physics, but may be undeniably better in terms of achieving converged results.
Finally, if the authors are not willing to provide a (less performant) Python version to support open-source code usage, I would remove the related bits from the manuscript's abstract.
According to the above, I encourage the authors to submit a major revision of their work for further assessment. Thanks
Ludovic Räss
GMD Topical EditorCitation: https://doi.org/10.5194/egusphere-2025-2242-EC1
Peer review completion
Journal article(s) based on this preprint
Model code and software
MinSIA v1: a lightweight and efficient implementation of the shallow ice approximation Stefan Hergarten https://doi.org/10.5281/zenodo.15362846
Video supplement
MInSIA examples Stefan Hergarten http://hergarten.at/minsia/examples/index.php
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 1,925 | 234 | 32 | 2,191 | 36 | 63 |
- HTML: 1,925
- PDF: 234
- XML: 32
- Total: 2,191
- BibTeX: 36
- EndNote: 63
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
Stefan Hergarten
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
- Preprint
(13098 KB) - Metadata XML
Dear author,
Please, find attached my report document.
Best,
Daniel