07 Sep 2021
07 Sep 2021
A comparison of the performance of depthintegrated icedynamics solvers
 ^{1}Complutense University of Madrid, Madrid, Spain
 ^{2}Geosciences Institute CSICUCM, Madrid, Spain
 ^{3}Potsdam Institute for Climate Impact Research, Potsdam, Germany
 ^{4}Climate and Global Dynamics Laboratory, National Center for Atmospheric Research, Boulder, CO 80305, USA
 ^{5}School of GeoSciences, University of Edinburgh, Edinburgh, UK
 ^{1}Complutense University of Madrid, Madrid, Spain
 ^{2}Geosciences Institute CSICUCM, Madrid, Spain
 ^{3}Potsdam Institute for Climate Impact Research, Potsdam, Germany
 ^{4}Climate and Global Dynamics Laboratory, National Center for Atmospheric Research, Boulder, CO 80305, USA
 ^{5}School of GeoSciences, University of Edinburgh, Edinburgh, UK
Abstract. In the last decade, the number of icesheet models has increased substantially, in line with the growth of the glaciological community. These models use solvers based on different approximations of ice dynamics. In particular, several depthintegrated dynamics approximations have emerged as fast solvers capable of resolving the relevant physics of ice sheets at the continen tal scale. However, the numerical stability of these schemes has not been studied systematically to evaluate their effectiveness in practice. Here we focus on three such solvers, the socalled Hybrid, L1L2SIA and DIVA solvers, as well as the wellknown SIA and SSA solvers as boundary cases. We investigate the numerical stability of these solvers as a function of grid resolution and the state of the ice sheet. Under simplified conditions with constant viscosity, the maximum stable timestep of the Hybrid solver, like the SIA solver, has a quadratic dependence on grid resolution. In contrast, the DIVA solver has a maximum timestep that is independent of resolution, like the SSA solver. Analysis indicates that the L1L2SIA solver should behave similarly, but in practice, the complexity of its implementation can make it difficult to maintain stability. In realistic simulations of the Greenland ice sheet with a nonlinear rheology, the DIVA and SSA solvers maintain superior numerical stability, while the SIA, Hybrid and L1L2SIA solvers show markedly poorer performance. At a grid resolution of ∆x = 4 km, the DIVA solver runs approximately 15 times faster than the Hybrid and L1L2SIA solvers. Our analysis shows that as resolution increases, the icedynamics solver can act as a bottleneck to model performance. The DIVA solver emerges as a clear outlier in terms of both model performance and its representation of the iceflow physics itself.
 Preprint
(1750 KB) 
Supplement
(54 KB)  BibTeX
 EndNote
Alexander Robinson et al.
Status: final response (author comments only)

RC1: 'Comment on tc2021239', Mauro Perego, 11 Oct 2021
In this paper the authors investigate the stability and performance of well known depthintegrated ice sheet models (SSA, SSA+SIA, L1L2, DIVA). The stability is studied theoretically on 1d problems with the assumption of uniform viscosity, and for explicit time integration schemes. The theoretical results are confirmed, for the most part, by numerical experiments. The authors also compare the stability and performance of these different models using the ice sheet codes CISM and Yelmo for modeling the Greenland ice sheet. Among these models, DIVA has better stability properties, and consequently performs better.
Overall the paper is well written and I think its contribution is very important. The theoretical and numerical stability results provided therein will certainly guide modelers and developers in choosing what depthintegrated models to use/implement.
I'm concerned about the discrepancy of their stability analysis and the results obtained using the L1L2 model implemented in CISM and Yelmo. I think it might have to do with the peculiar discretization chosen see detailed review below and I wonder whether they could get a better agreement if they changed the discretization used in the analysis.
I also suggest that the authors develop the theoretical analysis in an unbounded domain, which would make it simpler and cleaner.
Detailed review: Title: I think this paper is about stability more than performance, so I would suggest changing it into "A comparison of the stability and performance of ..."
 Abstract: please mention that the stability analysis is performed for an explicit time discretization scheme.
 page 1, line 20. Saying that the solution of Stokes problem is still infeasible is too strong of a statement. It would require a lot of resources, large clusters and parallel scalable implementations, but I would not say it's infeasible nowadays.
 Eq. (6) maybe it's worth pointing out that these equations implicitly define the velocity, as the viscosity depends on the velocity. Explicit formulas can be obtained. Also, it is possible to account for a sliding term in the SIA model by adjusting the SIA velocity with the term ρ g H ∇ s / β. I guess this is the formulation you use in the Greenland runs.
 Section 2.2: SSA solver relies on the hypothesis that the velocity is uniform in the vertical direction. So the velocity, at any zcoordinate is the same as the depthaveraged velocity. I think it would be clearer in equations (8), (9) and (10) to substitute the velocity and the basal velocity with the averaged velocity.
 Eq. (14). I would point out that the generalized integrals depend on the velocity (except from the n=1 case).
 Eq. (19). One can go directly to (17) to (19) byt taking the limit for \beta that goes to infinity, which gives in fact the noslip condition.
After eq. (19). The twosteps recipe to compute the velocity only holds for the case n=1, otherwise one needs to compute the vertical velocity in order to compute the viscosity needed for computing F_2 in eq. (19). Please explain how the problem is solved in the general case n >= 1.  page 10, line 1: For the stability analysis you can keep infinite domains. Just define x = n Δ x, with n being any integer (not just natural). I think this simplifies the analysis since you do not have to worry about boundaries, and in particular you can properly invert the circulant matrix arising from the discretization. It's understood that when solving the problem in practice you'll work on finite domains and have boundary effects.
 Eq. (43): Only the solution with the negative sign is acceptable.
 Section 3.2 I see a potential issue in the way the stability analysis is conducted. I would think that in most codes, when solving the thickness equation (59), the SSA and SIA velocities are discretized in the same way, whereas in the proposed analysis the SIA velocity term is expressed as a function of the thickness and treated differently. I think it would be best to consider the full velocity (u_{sia} + u_{ssa}) in (59) as a single object, take the derivative of H (u_{sia} + u_{ssa}), to get the advective and divergence terms as done in the DIVA case, and discretize them as done in the Diva case. In this way, the effective hybrid velocity should be the same as the depthaveraged hybrid velocity in (61).
 Section 3.3, Similarly, I would not treat separately the terms of the L1L2 depthaveraged velocity in eq. (71), and I would avoid reformulations in (72) but simply discretize the thickness equation as for the Diva method. This might be the source of discrepancy seen when solving the problem with the CISM and Yelmo discretizations.
 Eq. (72) It took me a long time to figure out how the equations have been derived. I would add some more steps to make it clearer.
 Section 3.4 I think the instability in the L1L2SIA model has been observed by BISICLES developers as well. To alleviate this, I believe BISICLES adopts a modification of the L1L2 method, consisting of avoiding the step in eq. (28), and taking the velocity to be uniform in the vertical direction and equal to the basal velocity. It might be worth mentioning this. (In your analysis, this modified L1L2 method would be indistinguishable from the SSA method, because you consider a uniform viscosity and n=1.)
 Section 4. Can you comment on the fact that the fitted exponent p goes from 1.5 to 2.8 instead of from 1 to 2 as in the theory. Does it depend on the Glen's law exponent n?
 Page 28, Line 22. There are two possible values for r, say r_{+} and r_{}. We have that r_{+} > 1, which implies, see (91), that the effect of a point x on the solution at y, grows exponentially with the distance between x and y. This is clearly not physical, so r_{} is the only physical solution. The two solutions exist because we are considering an infinite domain. Anyway, as I mentioned before I think that there is no need to do the analysis for a finite domain. That complicates the analysis in the Appendix without really adding any value.

RC2: 'Review comment of: A comparison of the performance of depthintegrated icedynamics solvers', Anonymous Referee #2, 15 Oct 2021
Review of “A comparison of performance of depthintegrated icedynamics solvers”
Summary:
The manuscript by A. Robinson and colleagues proposes an analysis of the performance of several depthintegrated stress balance approximations by combining a theoretical analysis of the performance with a comparison of results from several existing ice flow models that include them. They start with a 1d case and then compare the results for a more realistic simulation of the Greenland ice sheet. The paper is well written and usually clear to follow, though a few steps (listed below) could use some more details to facilitate following the derivation of the stability analysis. In the case of the L1L2SIA solver, the results are quite different between the stability analysis and the two Yelmo and CISM applications, but there is no clear explanation of this difference. I would be curious to learn what the authors think it the source of this discrepancy. Detailed comments below list relatively minor points, mostly to clarify the manuscript.
Detailed Comments:
p.1 l.6: Add that the simulations are done with finite differences and that you use con Neumann stability analysis.
p.1 l.18: thousands > maybe put a more exact number
p.1 l.19: kilometers > square kilometers
p.1 l.20: add that this is difficult to do at the high resolution required to accurately represent important processes
p.2 l.7: I don’t think the SIA is the most widely used approximation anymore (at least not used alone), it’s more historical and used to be the most common one
p.4 Eq.6: I was a bit confused by this integration, it would help to put a few more details on the steps to get to this form
p.5 Eq.9: Is that form assumed for all the approximations? If so, mention it here.
p.5 Eq.10: It is not clear how this integration is performed with the varying /mu, add indication of the main step/information needed, I could not figure out how to get to this result.
p.6 l.13: What are the implications of having a frozen bed with non zero basal stress? In which cases is that possible?
p.6 Eq.19: It might be simpler to say that in that case \beta >> 1 so \beta_eff converges toward 1/F_2
p.7 l.8: “in in” > “in its”
p.7 Eq.22: Add an explication about the derivation of this equation coming from the SIA and the parallel and normal components of \tau.
p.7 l.10: maybe add that this kind of analysis is designed for finite differences
p.9 Eq.32: Add where it comes from and the assumptions made.
p.9 l.33: Perturbation of what?
p.13 l.25: What is the conclusion?
p.13 l.21: Rewrite the entire equation given that assumptionp.17 l.9: What values are used in that range?
p.18 l .14: I would be a bit more nuanced that it is not really the case for sliding.
p.18 l.20: It does not really agree for \Delta_x < 10^3
p.18 l.33: Why could be the reasons why it is not stable as expected with more comprehensive models?
p.19 Fig.1: Add legend for the different lines on the figure. Caption: “in some other panels” > “In several panels”
p.20 l.8: How do you measure it?
p.20 l.11: How much does it adaptive timestep varies temporally?
p.26 l.1: Make sure that the results and data are actually archived
Alexander Robinson et al.
Alexander Robinson et al.
Viewed
HTML  XML  Total  Supplement  BibTeX  EndNote  

399  103  7  509  49  2  2 
 HTML: 399
 PDF: 103
 XML: 7
 Total: 509
 Supplement: 49
 BibTeX: 2
 EndNote: 2
Viewed (geographical distribution)
Country  #  Views  % 

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