the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 01 Oct 2025
A Fully Implicit Second Order Method for Viscous Free Surface Stokes Flow – Application to Glacier Simulations
Abstract. We present a fully implicit finite element method for viscous free surface Stokes problems which are common in glaciology and geodynamics. Standard fully implicit iterative solvers diverge for such problems when using large time steps. We address this by incorporating a stabilization term that vanishes upon convergence. This approach enables the use of large time steps, thereby significantly improving the efficiency of simulations. Furthermore, we leverage the fully implicit solver to develop second-order time-stepping schemes that mitigate the high errors associated with first-order methods for large time steps. The methods are tested and evaluated on model domains relevant to geodynamics and glaciology.
-
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
(1607 KB)
-
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
- Preprint
(1607 KB) - Metadata XML
- BibTeX
- EndNote
- Final revised paper
Journal article(s) based on this preprint
Interactive discussion
Status: closed
-
CEC1: 'No compliance with the policy of the journal', Juan Antonio Añel, 07 Oct 2025
-
AC1: 'Reply on CEC1', Josefin Ahlkrona, 29 Oct 2025
Dear Executive Editor,
Thank you for noticing this. We have now moved away from github and created a zenodo-repository for the edited version of the Biceps code used in this paper: https://zenodo.org/records/17371013
All files needed to run both the Biceps and Elmer-simulations are found at https://zenodo.org/records/17475634 where we have now made sure all relevant data is available, as well as bash-scripts that reproduce all convergence plots. The Elmer-code can as before be found at https://zenodo.org/records/7892181
We hope that we now comply with the rules of GMD. We will update the code availability section accordingly after potential review.
Sincerely,
Josefin Ahlkrona, André Löfgren and Clara HenryCitation: https://doi.org/10.5194/egusphere-2025-4359-AC1
-
AC1: 'Reply on CEC1', Josefin Ahlkrona, 29 Oct 2025
-
RC1: 'Comment on egusphere-2025-4359', Anonymous Referee #1, 15 Nov 2025
The article presents an ad-hoc-stabilized scheme for implicit time-evolution of the viscous free-surface Stokes equations. The idea relies on a relaxed explicit Euler discretization of the transport theorem that ties the discrete solution at two consecutive time steps. The main technical concern pertains to the introduction of a relaxation parameter controlling the size of the time-step appearing in the regularization. Its role is not discussed in full detail.
Since this parameter effectively reduces the time step size in the stabilization (but only in the stabilization and not in the general time-marching scheme) it introduces a mismatch in the timescales (and time grids) involved in the discretization. Thus, if the value of this parameter is not one, then it may introduce an inconsistency and lead to loss of convergence. The inconsistency, however, may only become dominant for much smaller values of the time step than the ones feasible in the article's simulations. The information provided in the text is not nearly enough to determine if this indeed would be the case, or if some additional measures are taken to avoid the issue. As a side comment, this "reduced time step" could potentially be interpreted as an intermediate stage (as in a Runge-Kutta method). This technique is relatively common for stiff problems, but the numerical solutions at these intermediate steps should not be given any physical meaning (just like the stages of an R-K method are not). In that sense the value of the parameter theta could be taken to be of the form 1/n where n is the number of intermediate stages to be discarded, interpreting only every "n-th" step as a physical solution. Unfortunately, the manuscript does not discuss this parameter with enough details to determine if this interpretation or idea is correct.
Nevertheless, the authors provide numerical evidence and sufficient comparisons that seem to indicate stability and reliability of the technique in some realistic situations of simplified (2D) glaciar flows.
The article is in general well written and, modulo the aforementioned discussion of the relaxation parameter, and the correction of several details scattered through the manuscript (see files attached), the referee can recommend the article for publication.
Annotations on the manuscript (see the attached file)
1.- limit
2.- Illustrate explicitly in the figure the domain \omega
3.- either "the boundary terms vanish" or "the boundary term vanishes"
4.- Illustrate the region \omega explicitly in figure 1
5.- "At the [free or top] surface... or even better, "on \Gamma_s"
6.- Indicate explicitly "on \Gamma_s" in equation (5)
7.- "on"
8.- Indicate explicitly "on \Gamma_b" in equation (8)
9.- Are all the curly brackets really necessary?
10.- This expression and several related ones later on seem to be missing the terms corresponding to the y coordinate. In the numerical experiments at the end the problem is considered in only 2D, but the rest of the discussion pertains to the full 3D problem, and so the "y-terms" must be included or an explicit comment should be made at this point regarding why only a 2D slice is considered. Compare for example with eq (14), where the "y-terms" do appear and then eq. (15) when they disappear back.
11.- See previous remark
12.- See remark 10
13.- W^{1,2}
14.- The line numbers in the algorithm are missing
15.- This statement is purely speculative. Either provide a proof/stringer argument or supress the sentence.
16.- do
17.- "... on weak form, recalling that $\partial_t(\rho\mathbf g \cdot v) = 0$, which leads to"
18.- The first equality seems unnecessary
19.- The stabilization seems to consist of connecting the terms at times k and k+1 through an explicit Euler discretization. In particular, the substitution of eq (17) into the right hand side of (1). The Euler discretization comes with the addition of a factor $\theta$ controlling the size of the time step. This is referred to as "user defined parameter". However, unless theta is identically equal to one, this would result on an inconsistent discretization and therefore would lead to failure of convergence. The role of the parameter needs to be discussed in more detail, taking into account the loss of consistency that it introduces.
20.- Just as in comment 15, this statement is purely speculative and it must be supported strongly or dismissed. In what specific sense do these solutions "mimic" the ones arising from an implicit discretization? Rather than this speculative statement it could be said that "numerical experiments seem to indicate that this results in more stable discretizations" simply pointing to the experimental evidence without speculating.
21.- This line belongs in the previous paragraph
22.- Line numbers missing
23.- This expression is missing the "y-terms" (see 10)
24.- indeed
25.- missing "y-terms" (see 10)
26.- This is not a distance, but it is a difference
27.- better factoring out $\Delta t$
28.- "...side of (15) and the expression above..."
29.- " , leading to: "
30.- this sentence should be right after the equation
31.- "is substituted into the Stokes equations"
32.- "A key observation is that, if the iteration converges,"
33.- this is not a distance, in this case "difference" can be used
34.- mark
35.- this is not a distance, but in this case "displacement" can be used
36.- two
37.- from
38.- "y-terms" missing, see (10)
39.- Is this the domain in Figure 3? If so, it should be mentioned
40.- This is not an error, unless you have an error estimator, it is simply the convergence tolerance.
41.- This figure is not referenced explicitly in the text. It should.
42.- This is not measuring an error, so the term "accuracy" is incorrect. Do you mean that this is used to determine convergence of the iteration?
43.- are
44.- this is not an error, it is the iteration tolerance
45.- Once again, this is not an error, unless there is an error estimator that is not mentioned anywhere
46.- "stabilization to be stable" is there such a thing as an unstable stabilization? Then it is no stabilization at all
47.- Do you mean "convergence of the second order method"?
48.- types
49.- remove comma
50.- Except for the first three points, all other errors are large to the point of making the numerical solution completely meaningless, so making a distinction between the intermediate and large step sizes seems a moot point
51.- Forward Euler
52.- This is not an error, it is simply the relative change between iterations
53.- Iteration tolerance
54.- convergence or iteration tolerance
55.- stabilization
56.- This point is quite important and is somewhat glossed over very lightly. It should be discussed in more detail, or at least references to a deeper discussion should be provided
57.- The notation is not clear. Is this a gaussian random variable and the first argument of octave( , ) is the mean and the second the standard deviation? Is it a different type of random variable? If so, what distribution does it follow and what do the arguments control?
58.- There is a double comma
59.- ms^{-2}
60.- convergence tolerance
61.- This is not an error, is simply relative change. An iteration can converge to the wrong solution, thus even id the relative change is 0 the error might be very large. Since there is no error estimator invovled, the term error should be avoided.
62.- scheme
63.- keep
64. - If I recall correctly there was an issue with the use of github. Should this be updated?-
AC2: 'Reply on RC1', Josefin Ahlkrona, 13 Feb 2026
The comment was uploaded in the form of a supplement: https://egusphere.copernicus.org/preprints/2025/egusphere-2025-4359/egusphere-2025-4359-AC2-supplement.pdf
-
AC2: 'Reply on RC1', Josefin Ahlkrona, 13 Feb 2026
-
RC2: 'Comment on egusphere-2025-4359', Boris Kaus, 16 Jan 2026
This is a very interesting and well-written manuscript that improves my earlier free surface stabilisation (FSSA) algorithm to allow for significantly larger timesteps in glacier and geodynamics codes that have a free surface.
The authors demonstrate that a purely implicit timestepping algorithm doesn’t give much advantage compared to an explicit scheme with FSSA if a first-order timestepping scheme is employed. This surprised me somewhat (I had thought otherwise), but is rigorously shown here. In addition, they show that using a second-order timestepping scheme can be made second-order by using a modified FSSA algorithm.
This is initially described for a simple linear viscous setup, which is expanded to a more realistic ice flow simulation in section 7.
Overall, this is a very nice contribution that will be helpful for both the glaciology and geodynamics communities. I think it can be accepted with minor revisions, and I just have a few minor remarks:
- You seem to use Picard fixed-point iterations. If you have a good initial guess, a faster way to converge nonlinear iterations is Newton iterations, which should work particularly well for moderately nonlinear problems such as powerlaw viscous materials with an exponent of n=3. This would potentially also help the convergence if the implicit timesteps, even though it requires deriving the appropriate Jacobian. Is this not commonly done in the glaciology community, or are there other reasons you stick with Picard?
- eq. 29: Can you give a function or algorithm on how “octave_?” is computed? The way this is written is a bit puzzling to me, but this might be because I am not from glaciology...
Minor typos:
l. 299: kg m^{-3}
Fig. 8 caption: “stabilisation”
Citation: https://doi.org/10.5194/egusphere-2025-4359-RC2 - AC3: 'Reply on RC2', Josefin Ahlkrona, 13 Feb 2026
Interactive discussion
Status: closed
-
CEC1: 'No compliance with the policy of the journal', Juan Antonio Añel, 07 Oct 2025
Dear authors,
Unfortunately, after checking your manuscript, it has come to our attention that it does not comply with our "Code and Data Policy".
https://www.geoscientific-model-development.net/policies/code_and_data_policy.html
You have archived the Biceps code on GitHub. However, GitHub is not a suitable repository for scientific publication. GitHub itself instructs authors to use other long-term archival and publishing alternatives, such as Zenodo. Therefore, the current situation with your manuscript is irregular. Please, publish your code in one of the appropriate repositories and reply to this comment with the relevant information (link and a permanent identifier for it (e.g. DOI)) as soon as possible, as we can not accept manuscripts in Discussions that do not comply with our policy. Also, please include the relevant primary input/output data.Also, you must include a modified 'Code and Data Availability' section in a potentially reviewed manuscript, containing the information of the new repositories.
I must note that if you do not fix this problem, we cannot continue with the peer-review process or accept your manuscript for publication in our journal.
Juan A. Añel
Geosci. Model Dev. Executive EditorCitation: https://doi.org/10.5194/egusphere-2025-4359-CEC1 -
AC1: 'Reply on CEC1', Josefin Ahlkrona, 29 Oct 2025
Dear Executive Editor,
Thank you for noticing this. We have now moved away from github and created a zenodo-repository for the edited version of the Biceps code used in this paper: https://zenodo.org/records/17371013
All files needed to run both the Biceps and Elmer-simulations are found at https://zenodo.org/records/17475634 where we have now made sure all relevant data is available, as well as bash-scripts that reproduce all convergence plots. The Elmer-code can as before be found at https://zenodo.org/records/7892181
We hope that we now comply with the rules of GMD. We will update the code availability section accordingly after potential review.
Sincerely,
Josefin Ahlkrona, André Löfgren and Clara HenryCitation: https://doi.org/10.5194/egusphere-2025-4359-AC1
-
AC1: 'Reply on CEC1', Josefin Ahlkrona, 29 Oct 2025
-
RC1: 'Comment on egusphere-2025-4359', Anonymous Referee #1, 15 Nov 2025
The article presents an ad-hoc-stabilized scheme for implicit time-evolution of the viscous free-surface Stokes equations. The idea relies on a relaxed explicit Euler discretization of the transport theorem that ties the discrete solution at two consecutive time steps. The main technical concern pertains to the introduction of a relaxation parameter controlling the size of the time-step appearing in the regularization. Its role is not discussed in full detail.
Since this parameter effectively reduces the time step size in the stabilization (but only in the stabilization and not in the general time-marching scheme) it introduces a mismatch in the timescales (and time grids) involved in the discretization. Thus, if the value of this parameter is not one, then it may introduce an inconsistency and lead to loss of convergence. The inconsistency, however, may only become dominant for much smaller values of the time step than the ones feasible in the article's simulations. The information provided in the text is not nearly enough to determine if this indeed would be the case, or if some additional measures are taken to avoid the issue. As a side comment, this "reduced time step" could potentially be interpreted as an intermediate stage (as in a Runge-Kutta method). This technique is relatively common for stiff problems, but the numerical solutions at these intermediate steps should not be given any physical meaning (just like the stages of an R-K method are not). In that sense the value of the parameter theta could be taken to be of the form 1/n where n is the number of intermediate stages to be discarded, interpreting only every "n-th" step as a physical solution. Unfortunately, the manuscript does not discuss this parameter with enough details to determine if this interpretation or idea is correct.
Nevertheless, the authors provide numerical evidence and sufficient comparisons that seem to indicate stability and reliability of the technique in some realistic situations of simplified (2D) glaciar flows.
The article is in general well written and, modulo the aforementioned discussion of the relaxation parameter, and the correction of several details scattered through the manuscript (see files attached), the referee can recommend the article for publication.
Annotations on the manuscript (see the attached file)
1.- limit
2.- Illustrate explicitly in the figure the domain \omega
3.- either "the boundary terms vanish" or "the boundary term vanishes"
4.- Illustrate the region \omega explicitly in figure 1
5.- "At the [free or top] surface... or even better, "on \Gamma_s"
6.- Indicate explicitly "on \Gamma_s" in equation (5)
7.- "on"
8.- Indicate explicitly "on \Gamma_b" in equation (8)
9.- Are all the curly brackets really necessary?
10.- This expression and several related ones later on seem to be missing the terms corresponding to the y coordinate. In the numerical experiments at the end the problem is considered in only 2D, but the rest of the discussion pertains to the full 3D problem, and so the "y-terms" must be included or an explicit comment should be made at this point regarding why only a 2D slice is considered. Compare for example with eq (14), where the "y-terms" do appear and then eq. (15) when they disappear back.
11.- See previous remark
12.- See remark 10
13.- W^{1,2}
14.- The line numbers in the algorithm are missing
15.- This statement is purely speculative. Either provide a proof/stringer argument or supress the sentence.
16.- do
17.- "... on weak form, recalling that $\partial_t(\rho\mathbf g \cdot v) = 0$, which leads to"
18.- The first equality seems unnecessary
19.- The stabilization seems to consist of connecting the terms at times k and k+1 through an explicit Euler discretization. In particular, the substitution of eq (17) into the right hand side of (1). The Euler discretization comes with the addition of a factor $\theta$ controlling the size of the time step. This is referred to as "user defined parameter". However, unless theta is identically equal to one, this would result on an inconsistent discretization and therefore would lead to failure of convergence. The role of the parameter needs to be discussed in more detail, taking into account the loss of consistency that it introduces.
20.- Just as in comment 15, this statement is purely speculative and it must be supported strongly or dismissed. In what specific sense do these solutions "mimic" the ones arising from an implicit discretization? Rather than this speculative statement it could be said that "numerical experiments seem to indicate that this results in more stable discretizations" simply pointing to the experimental evidence without speculating.
21.- This line belongs in the previous paragraph
22.- Line numbers missing
23.- This expression is missing the "y-terms" (see 10)
24.- indeed
25.- missing "y-terms" (see 10)
26.- This is not a distance, but it is a difference
27.- better factoring out $\Delta t$
28.- "...side of (15) and the expression above..."
29.- " , leading to: "
30.- this sentence should be right after the equation
31.- "is substituted into the Stokes equations"
32.- "A key observation is that, if the iteration converges,"
33.- this is not a distance, in this case "difference" can be used
34.- mark
35.- this is not a distance, but in this case "displacement" can be used
36.- two
37.- from
38.- "y-terms" missing, see (10)
39.- Is this the domain in Figure 3? If so, it should be mentioned
40.- This is not an error, unless you have an error estimator, it is simply the convergence tolerance.
41.- This figure is not referenced explicitly in the text. It should.
42.- This is not measuring an error, so the term "accuracy" is incorrect. Do you mean that this is used to determine convergence of the iteration?
43.- are
44.- this is not an error, it is the iteration tolerance
45.- Once again, this is not an error, unless there is an error estimator that is not mentioned anywhere
46.- "stabilization to be stable" is there such a thing as an unstable stabilization? Then it is no stabilization at all
47.- Do you mean "convergence of the second order method"?
48.- types
49.- remove comma
50.- Except for the first three points, all other errors are large to the point of making the numerical solution completely meaningless, so making a distinction between the intermediate and large step sizes seems a moot point
51.- Forward Euler
52.- This is not an error, it is simply the relative change between iterations
53.- Iteration tolerance
54.- convergence or iteration tolerance
55.- stabilization
56.- This point is quite important and is somewhat glossed over very lightly. It should be discussed in more detail, or at least references to a deeper discussion should be provided
57.- The notation is not clear. Is this a gaussian random variable and the first argument of octave( , ) is the mean and the second the standard deviation? Is it a different type of random variable? If so, what distribution does it follow and what do the arguments control?
58.- There is a double comma
59.- ms^{-2}
60.- convergence tolerance
61.- This is not an error, is simply relative change. An iteration can converge to the wrong solution, thus even id the relative change is 0 the error might be very large. Since there is no error estimator invovled, the term error should be avoided.
62.- scheme
63.- keep
64. - If I recall correctly there was an issue with the use of github. Should this be updated?-
AC2: 'Reply on RC1', Josefin Ahlkrona, 13 Feb 2026
The comment was uploaded in the form of a supplement: https://egusphere.copernicus.org/preprints/2025/egusphere-2025-4359/egusphere-2025-4359-AC2-supplement.pdf
-
AC2: 'Reply on RC1', Josefin Ahlkrona, 13 Feb 2026
-
RC2: 'Comment on egusphere-2025-4359', Boris Kaus, 16 Jan 2026
This is a very interesting and well-written manuscript that improves my earlier free surface stabilisation (FSSA) algorithm to allow for significantly larger timesteps in glacier and geodynamics codes that have a free surface.
The authors demonstrate that a purely implicit timestepping algorithm doesn’t give much advantage compared to an explicit scheme with FSSA if a first-order timestepping scheme is employed. This surprised me somewhat (I had thought otherwise), but is rigorously shown here. In addition, they show that using a second-order timestepping scheme can be made second-order by using a modified FSSA algorithm.
This is initially described for a simple linear viscous setup, which is expanded to a more realistic ice flow simulation in section 7.
Overall, this is a very nice contribution that will be helpful for both the glaciology and geodynamics communities. I think it can be accepted with minor revisions, and I just have a few minor remarks:
- You seem to use Picard fixed-point iterations. If you have a good initial guess, a faster way to converge nonlinear iterations is Newton iterations, which should work particularly well for moderately nonlinear problems such as powerlaw viscous materials with an exponent of n=3. This would potentially also help the convergence if the implicit timesteps, even though it requires deriving the appropriate Jacobian. Is this not commonly done in the glaciology community, or are there other reasons you stick with Picard?
- eq. 29: Can you give a function or algorithm on how “octave_?” is computed? The way this is written is a bit puzzling to me, but this might be because I am not from glaciology...
Minor typos:
l. 299: kg m^{-3}
Fig. 8 caption: “stabilisation”
Citation: https://doi.org/10.5194/egusphere-2025-4359-RC2 - AC3: 'Reply on RC2', Josefin Ahlkrona, 13 Feb 2026
Peer review completion
Journal article(s) based on this preprint
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 458 | 204 | 37 | 699 | 24 | 25 |
- HTML: 458
- PDF: 204
- XML: 37
- Total: 699
- BibTeX: 24
- EndNote: 25
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
Josefin Ahlkrona
A. Clara J. Henry
André Löfgren
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
- Preprint
(1607 KB) - Metadata XML
Dear authors,
Unfortunately, after checking your manuscript, it has come to our attention that it does not comply with our "Code and Data Policy".
https://www.geoscientific-model-development.net/policies/code_and_data_policy.html
You have archived the Biceps code on GitHub. However, GitHub is not a suitable repository for scientific publication. GitHub itself instructs authors to use other long-term archival and publishing alternatives, such as Zenodo. Therefore, the current situation with your manuscript is irregular. Please, publish your code in one of the appropriate repositories and reply to this comment with the relevant information (link and a permanent identifier for it (e.g. DOI)) as soon as possible, as we can not accept manuscripts in Discussions that do not comply with our policy. Also, please include the relevant primary input/output data.
Also, you must include a modified 'Code and Data Availability' section in a potentially reviewed manuscript, containing the information of the new repositories.
I must note that if you do not fix this problem, we cannot continue with the peer-review process or accept your manuscript for publication in our journal.
Juan A. Añel
Geosci. Model Dev. Executive Editor