Love number computation within the Ice-sheet and Sea-level System Model (ISSM v4.24)

Caron, Lambert; Ivins, Erik; Larour, Eric; Adhikari, Surendra; Metivier, Laurent

doi:https://doi.org/10.5194/egusphere-2024-3414

Preprints

https://doi.org/10.5194/egusphere-2024-3414

Preprints

03 Feb 2025

| 03 Feb 2025

Love number computation within the Ice-sheet and Sea-level System Model (ISSM v4.24)

Lambert Caron, Erik Ivins, Eric Larour, Surendra Adhikari, and Laurent Metivier

Abstract. The Love number solver presented here is a new capability within the Ice-sheet and Sea-level System Model (ISSM) for computing the solid-Earth response to tidal forcing and surface mass loading. This new capability allows solving zero-frequency free oscillation equations of motion decomposed into the well-known y_i system and enables high wave number computations with spherical harmonic truncation degree of ~10,000 and above. It facilitates capturing the high-resolution response of the solid Earth to a step-function forcing in terms of gravity potential changes, vertical and horizontal bedrock displacement, and polar motion. The model incorporates compressible isotropic elasticity and three forms of linear viscoelasticity for mantle rheology: the Maxwell, Burgers, and Extended Burgers Materials (EBM). We detail our approach to the paralellization and numerical optimization of the solver, and report the accuracy of our results with respect to community benchmark solutions. Our main motivation is to facilitate simulations of a coupled system of surface mass transport (e.g., changes in polar ice sheets and sea level) and solid Earth models at kilometer-scale lateral resolutions with numerical efficiency that supports the exploration of large model ensembles and uncertainty quantification.

Received: 02 Nov 2024 – Discussion started: 03 Feb 2025

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.

Download & links

Preprint (PDF, 1618 KB)

Supplement (92 KB)

Download & links

Lambert Caron, Erik Ivins, Eric Larour, Surendra Adhikari, and Laurent Metivier

Status: final response (author comments only)

RC1:
'Comment on egusphere-2024-3414', Anonymous Referee #1, 17 Mar 2025

General comments
This manuscript describes a new capability for evaluating viscoelastic Love Numbers that has been built within the ISSM coupled ice sheet modeling framework. The new functionality allows to model the solid Earth response to surface loads up to very short wavelengths, and supports different rheological laws including the EBM transient model. The manuscript is very well written, it is technically sound and discusses in detail a comprehensive set of benchmarks with community reference solutions for incompressible LNs in order to validate the numerical approach. The authors also devoted considerable effort in improving the code efficiency, in view of its possible applications to ensemble modeling or Bayesian inversion frameworks, hence I believe that this new tool will be of great relevance for the GIA community. I have just one major point that in my opinion should be addressed before recommending publication, which is illustrated below.
Specific comments
My main (and only) concern about the manuscript is about the stability of the approach discussed in the paper in the case of compressible models. A compressible self-gravitating model is known to be affected by gravitational instabilities both in the elastic and viscoelastic regimes (see e.g. Hanyk et al., 1999 or Vermeersen and Mitrovica, 2000) which, from a mathematical point of view, manifest themselves as singularities in the Laplace-transformed solution on the real positive axis. Within the classic framework of viscoelastic normal modes those correspond to exponentially growing modes, which are usually dropped since their timescale is generally much larger than those of interest for GIA, and only stable modes associated with singularities on the real negative axis are kept to evaluate the time-domain solution. However, if the time-domain LNs are computed according to eq (35), there is effectively no way to filter out those modes (since the individual modes are not identified at all, which is one the main advantages of the proposed method), and therefore it is possible that the Laplace solution is sampled near a singular point, which would break the inversion scheme. I do not have any specific suggestion about how this problem can be circumvented, but I think that it is an important aspect that would merit some discussion.
Minor points:
Page 9, equation 17: For a tidal forcing, the gravity potential is computed without considering the direct term (i.e. on the r.h.s. ‘k’ is used instead of ‘1+k’). Is it to be compliant with the GRACE convention for the definition of potential, in which the direct term is not included? If this is the case, I suggest making this point explicit.
Page 17, lines 335-340: the “pseudo-fluid Love numbers” and the “degree n_vmax where the viscous part of the LNs becomes negligible” are first mentioned here, but the meaning of those concepts is made explicit only in section 5 below. Maybe it could be beneficial to the reader if a few more words are given here about those two points, anticipating that they will be discussed in more detail in section 5.
Technical corrections
Page 3, line 75: Missing period after the reference to Kierulf et al.
Page 4, line 92: maybe βtⁿ should read βt^1/n, or n ≅ 3 should be n ≅ 1/3 ?
Page 7, line 122: “an mantle” -> “a mantle” (or “a solid mantle”?)
Page 12, line 123: my understanding is that rⁿ should be r^-n
Page 15, line 290: “This methods” -> “This method”
Page 21, line 357: “theses” -> “these”
Page 25, lines 439-440: “point to with” -> “point to”
Please note that some entries in the bibliography are missing the DOI.

Citation: https://doi.org/10.5194/egusphere-2024-3414-RC1
- AC1: 'Reply on RC1', Lambert Caron, 04 Aug 2025
  
  Reviewer 1:
  This manuscript describes a new capability for evaluating viscoelastic Love Numbers that has been built within the ISSM coupled ice sheet modeling framework. The new functionality allows to model the solid Earth response to surface loads up to very short wavelengths, and supports different rheological laws including the EBM transient model. The manuscript is very well written, it is technically sound and discusses in detail a comprehensive set of benchmarks with community reference solutions for incompressible LNs in order to validate the numerical approach. The authors also devoted considerable effort in improving the code efficiency, in view of its possible applications to ensemble modeling or Bayesian inversion frameworks, hence I believe that this new tool will be of great relevance for the GIA community. I have just one major point that in my opinion should be addressed before recommending publication, which is illustrated below.
  Specific comments
  My main (and only) concern about the manuscript is about the stability of the approach discussed in the paper in the case of compressible models. A compressible self-gravitating model is known to be affected by gravitational instabilities both in the elastic and viscoelastic regimes (see e.g. Hanyk et al., 1999 or Vermeersen and Mitrovica, 2000) which, from a mathematical point of view, manifest themselves as singularities in the Laplace-transformed solution on the real positive axis. Within the classic framework of viscoelastic normal modes those correspond to exponentially growing modes, which are usually dropped since their timescale is generally much larger than those of interest for GIA, and only stable modes associated with singularities on the real negative axis are kept to evaluate the time-domain solution. However, if the time-domain LNs are computed according to eq (35), there is effectively no way to filter out those modes (since the individual modes are not identified at all, which is one the main advantages of the proposed method), and therefore it is possible that the Laplace solution is sampled near a singular point, which would break the inversion scheme. I do not have any specific suggestion about how this problem can be circumvented, but I think that it is an important aspect that would merit some discussion.
  [Authors’ response]: We thank the reviewer for bringing up this topic with relevant references. We agree that this is an important point and have added a new paragraph (see section 3.5) to discuss this potential downside of the Post-Widder transform. Our discussion builds on Hanyk et al and Vermeersen and Mitrovica to suggest mitigation strategies for the user.
  Minor points:
  Page 9, equation 17: For a tidal forcing, the gravity potential is computed without considering the direct term (i.e. on the r.h.s. ‘k’ is used instead of ‘1+k’). Is it to be compliant with the GRACE convention for the definition of potential, in which the direct term is not included? If this is the case, I suggest making this point explicit.
  [Authors' response]: This is not related to GRACE, but to the definition of the system. We assume that the geoid is based on the gravity potential of the Earth, reflecting both changes of mass in the interior and at the surface (hence why 1 is included in the surface loading case). Since tidal forcing originates from the Moon, Sun or other celestial bodies (except for the pole tide case), including the direct term would reflect the gravity field of the combined system instead of the geoid. For the rotational feedback (aka the pole tide) and for GRACE applications, we have a dedicated process that uses the PMTF and degree-2 tidal love numbers to compute the resulting deformation. In the end, the specifics depends on the intended application by the user, the scope of this paper is to provide them with the values of k. We have revised the paragraph following eq 17 to explain this choice.
  Page 17, lines 335-340: the “pseudo-fluid Love numbers” and the “degree n_vmax where the viscous part of the LNs becomes negligible” are first mentioned here, but the meaning of those concepts is made explicit only in section 5 below. Maybe it could be beneficial to the reader if a few more words are given here about those two points, anticipating that they will be discussed in more detail in section 5.
  [Authors' response]: We have added a few sentences in section 4 as suggested, and referenced the appropriate section for more information.
  Technical corrections
  Page 3, line 75: Missing period after the reference to Kierulf et al.
  Page 4, line 92: maybe βtⁿ should read βt^1/n, or n ≅ 3 should be n ≅ 1/3 ?
  Page 7, line 122: “an mantle” -> “a mantle” (or “a solid mantle”?)
  Page 12, line 123: my understanding is that rⁿ should be r^-n
  Page 15, line 290: “This methods” -> “This method”
  Page 21, line 357: “theses” -> “these”
  Page 25, lines 439-440: “point to with” -> “point to”
  Please note that some entries in the bibliography are missing the DOI.
  [Authors' response]: All of these points have been corrected. We thank the reviewer for a thorough examination of the manuscript.
  
  Citation: https://doi.org/10.5194/egusphere-2024-3414-AC1
RC2:
'Comment on egusphere-2024-3414', Anonymous Referee #2, 01 Jul 2025

The comment was uploaded in the form of a supplement: https://egusphere.copernicus.org/preprints/2025/egusphere-2024-3414/egusphere-2024-3414-RC2-supplement.pdf

Citation: https://doi.org/10.5194/egusphere-2024-3414-RC2
- AC2: 'Reply on RC2', Lambert Caron, 04 Aug 2025
  
  Reviewer 2:
  Review: Love number computation within the Ice-sheet and Sea-level System Model (ISSM V4.24) – Caron et al.
  Summary:
  The manuscript of Caron et al. presents the theory for computing Love numbers for Maxwell, Burgers, and Extended Burgers rheologies. In this discussion they present a general overview of the parallelization strategy and demonstrate their optimization strategy. Furthermore, they compare their Love numbers to those of Spada et al. (2011), a benchmarking paper in the GIA community, and find good agreement for a Maxwell rheology. This body of work and the associated code is an important contribution to GIA modeling and coupled GIA/ice sheet modeling.
  The code may also have utility for other deformation problems given suitable modifications. Overall, this manuscript is well written and in fairly good shape. Below I outline a few major items and several minor items that would help improve the quality of the manuscript. None of these should take very long to implement and thus, I suggest this manuscript be accepted following minor revisions.
  [Authors' response]: We thank the reviewer for their thorough examination of both the content and presentation of the manuscript. Their comments have significantly improved the clarity and quality of the paper.
  Major requests:
  • Abstract: Please add the something about the Love number approach being for radial 1D (i.e., spherical symmetric earth models) so the limitations of this approach is more clear to the general audience. You could also add a sentence on the benefits of this approach (e.g., that it is computationally quick compared to codes with 3D viscosity structures) to balance out the negative.
  [Authors' response]: We have clarified the abstract as requested.
  • The dot used to show the divergence is at the bottom “.” as opposed to the middle “∙”. Please fix this in equation 2, 18. The dot product in equation 12 is at the bottom “.” and should be in the middle “∙”
  [Authors' response]: Corrected.
  • Please make sure all variables and super-/subscripts are defined within the text. Table 1 is a helpful guide, but it doesn’t contain everything and it makes for an unpleasant read to hunt for variable definitions that could easily be stated in the text. For example: eq 11-14: SH degree n and SH order m are not explained nor are the bounds for their summations; eq 19 – what is E; eq 29 – explain 𝜌# in text please; etc
  [Authors' response]: We have modified the text throughout the manuscript to label variables as they are introduced.
  • Line 166-171: There is a disagreement in the ordering of the text and yi in comparison to the equations 11-14. For a savvy reader this will not be a problem, but for others this could be confusing. I suggest first presenting equations 11, 12, 13, and 14 in order. So, introduce the full displacement vector field, the full traction vector field, etc. Then introduce “yi” and state the order of the vector elements following the convention of Alterman et al.
  [Authors' response]: We moved this sentence after the y_i equations, as suggested, and clarified explicitly which one related to each physical variable therein.
  • Equation 17 – Please state what N(\theta,\phi,t) is (i.e., how is it different from N^{m}_{n}, the
  normalization factor or the number of layers N_{L}) and also consider a different variable to avoid
  confusion
  [Authors' response]: Noted, we changed the number of layers to N_L to \mathcal{L}, the spherical harmonics Normalization N^{m}_{n} to \mathcal{N}, and kept N for the notation of the geoid, as standard in the field.
  • Is there a standard you followed for using ~, ‘, and f_{i} for the different factors in 38. It would be helpful if there was and if you could please explain it.
  [Authors' response]: We have added a sentence (section 5.1.1) to make our approach to notations explicit: $f_{1..11}$ refer to the precomputed terms in order in which the appear in eq. (39), terms with a prime symbol $'$ are used for non-dimensionalized quantities (e.g. $\mu'$ is the non-dimensionalized shear modulus) while terms with a tilde $~$ are used for other scaling constants.
  • Figure 7: Please add a bounding box to left figure to show the location of the right figure.
  [Authors' response]: Done, this is also referenced in the caption now.
  • Figure 8: Please add something to the caption describing where the h love numbers are plotting.
  [Authors' response]: We have amended the caption to indicate that the curve for h tends to be identical to that of k (red line).
  • Multi-panel figures: Please choose a consistent method for referencing multi panel figures. For example, Figure 8 uses (a) and (b), but none of the other figures use letters.
  [Authors' response]: We have removed the letters from Figure 8 to make all figures consistent.
  • Multi-panel Figures: Please ensure there is a description of each panel in the caption. For example, Figure 2 does not explain both plots. Second please modify the captions so that it is easier to find where the descriptions of each panel is. Bold and italic font might be helpful, particularly when using “Right”, “Top Row”, etc.
  [Authors' response]: We revised most of the captions in this way, highlighting text pointing to rows and columns in bold. For figure 1 and 3, there was little to add to the information provided by the figure frame titles so we opted not to.
  
  • Section 5.2: As described, it is unclear how experiment 1 & 2 are changing the mesh discretization. Both are described as “increasing the radial grid resolution ….”. The term layer is also used, but in this context, it is unclear if the authors are referring to the layers shown in Figure 1 or if they mean a spherical shell. A figure showing the change to the mesh for each experiment would be helpful.
  [Authors' response]: We mean neither the layers themselves nor a spherical shell. For the sake of performing the propagation of the yi system within each layer (from bottom to top boundary), we perform the iterative procedure explained in section 3.3, updating the system after each radial step Δr. This subdivision creates a 1D grid from r_{inner} to the surface. The question asked in comparing experiments 1 and 2 is twofold: 1) for a given number of grid points (equivalent resources), is it better that Δr be uniform across layers, or that each layer have the same number of propagation steps? 2) how fine should this grid be, in order to solve for surface Love numbers accurately, before yielding diminishing returns?
  We have added a short paragraph at the start of section 5.2 to better explain what get across these points.
  
  Minor requests:
  • Line 28: “… model disparities related to solid Earth structure and mantle Rheologies.” Please add “1D” to the solid Earth structure to be clearer.
  [Authors' response]: Added “radial”.
  • Line 75: missing a period at the end of sentence “ … Kierulf et al., 2022).”
  [Authors' response]: Corrected.
  • Line 79: Table 4 is introduced before Table 1, 2, and 3. Consider revising the text of table locations/labels so things appear in the order in which they are discussed.
  [Authors' response]: Done.
  • Line 83: Since you are discussing changes in Earth’s shape would \dot{J}_2 be more appropriate than just J_2?
  [Authors' response]: Corrected.
  • Line 103: “… Heaviside function forcing …” -> “… Heaviside forcing function …”, maybe easier to read?
  [Authors' response]: Agreed, changed.
  • Line 161: “… and material …” -> “… and then material …”
  [Authors' response]: Changed to “and the material isotropic.”.
  • Line 205: What is meant by “resolved”?
  [Authors' response]: Changed to “solved”.
  • Equations 24, 25, 26: The 0 on the right-hand side should be written as a vector and not a scalar.
  [Authors' response]: Fixed.
  • Table 3: Please add units to your times. Currently one must assume they are in seconds like the precision you report in the table caption.
  [Authors' response]: Changed “Computation time (s)” to “Computation time (seconds)” at the beginning of the caption to make this unit harder to miss. We would prefer to not include it in each of the table cells in order to avoid clutter.
  • Line 440: delete “with”
  [Authors' response]: Done.
  
  Citation: https://doi.org/10.5194/egusphere-2024-3414-AC2

Lambert Caron, Erik Ivins, Eric Larour, Surendra Adhikari, and Laurent Metivier

Supplement

https://doi.org/10.5194/egusphere-2024-3414-supplement

Lambert Caron, Erik Ivins, Eric Larour, Surendra Adhikari, and Laurent Metivier

Viewed

Total article views: 694 (including HTML, PDF, and XML)

HTML	PDF	XML	Total	Supplement	BibTeX	EndNote
378	298	18	694	42	34	29

HTML: 378
PDF: 298
XML: 18
Total: 694
Supplement: 42
BibTeX: 34
EndNote: 29

Views and downloads (calculated since 03 Feb 2025)

Month	HTML	PDF	XML	Total
Feb 2025	52	28	3	83
Mar 2025	52	16	2	70
Apr 2025	27	164	5	196
May 2025	22	12	1	35
Jun 2025	20	10	2	32
Jul 2025	37	28	3	68
Aug 2025	86	35	2	123
Sep 2025	82	5	0	87

Cumulative views and downloads (calculated since 03 Feb 2025)

Month	HTML	PDF	XML	Total
Feb 2025	52	28	3	83
Mar 2025	52	16	2	70
Apr 2025	27	164	5	196
May 2025	22	12	1	35
Jun 2025	20	10	2	32
Jul 2025	37	28	3	68
Aug 2025	86	35	2	123
Sep 2025	82	5	0	87

Viewed (geographical distribution)

Total article views: 682 (including HTML, PDF, and XML) Thereof 682 with geography defined and 0 with unknown origin.

Country	#	Views	%

Latest update: 11 Sep 2025

Short summary

Presented here is a new model of the solid-Earth response to tides and mass changes in ice sheets, oceans, and groundwater, in of terms of gravity change and bedrock motion. The model is capable simulating mantle deformation including elasticity, transient and steady-state viscous flow. We detail our approach to numerical optimization, and report the accuracy of results with respect to community benchmarks. The resulting coupled system features kilometer-scale resolution and fast computation.


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