Preprints
https://doi.org/10.5194/egusphere-2024-1093
© Author(s) 2024. This work is distributed under
the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
https://doi.org/10.5194/egusphere-2024-1093
© Author(s) 2024. This work is distributed under
the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
29 May 2024
| 29 May 2024
Status: this preprint is open for discussion and under review for Geoscientific Model Development (GMD).
Isogeometric analysis of the lithosphere under topographic loading: Igalith v1.0.0
Abstract. This paper presents methods from isogeometric finite element analysis for numerically solving problems in geoscience involving partial differential equations. In particular, we consider the numerical simulation of shells and plates in the context of isostasy. Earth's lithosphere is modeled as a thin elastic shell or plate floating on the asthenosphere and subject to topographic loads. We demonstrate the computational methods on the isostatic boundary value problem posed on selected geographic locations. For Europe, the computed lithospheric depression is compared with available Mohorovičić depth data. We also perform parameter identification for the effective elastic thickness of the lithosphere, the rock density, and the topographic load that are most plausible to explain the measured depths. An example of simulating the entire lithosphere of the Earth as a spherical shell using multi-patch isogeometric analysis is presented, which provides an alternative to spherical harmonics for solving partial differential equations on a spherical domain.
How to cite. Rosandi, R., Rosandi, Y., and Simeon, B.: Isogeometric analysis of the lithosphere under topographic loading: Igalith v1.0.0, EGUsphere [preprint], https://doi.org/10.5194/egusphere-2024-1093, 2024.
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 preprint. The responsibility to include appropriate place names lies with the authors.
Download & links
Download & links
- Preprint
(27597 KB) - Metadata XML
- BibTeX
- EndNote
Status: open (extended)
Comment types: AC – author | RC – referee | CC – community | EC – editor | CEC – chief editor
| : Report abuse
-
RC1: 'Comment on egusphere-2024-1093', Anonymous Referee #1, 30 Dec 2024
reply
This paper presents the applicability of Isogeometric Analysis for problems in geoscience. More particularly, it is herein considered the modelling of the lithosphere as a thin shell subject to gravitational and buoyancy forces exerted by the topographic features above the lithosphere and the upper layers of the athenosphere, respectively. The paper is overall very well written, easy to read, and the findings are accompanied/supported by numerical evidence. Nevertheless, I would like to have the remarks below addressed in an adequate manner before I can suggest publication of the paper:
- On page 2 the authors mention "Further features and capabilities of isogeometric analysis presented in this paper are the exact representation of curved domains, the coupling of multiple patches". The coupling of multiple patches is not a feature of IGA but rather a challenge of the method, from which standard finite element methods do not suffer. I would not enlist the coupling of multiple patches as a feature, but I would highlight that patch coupling is essential in IGA for problems of practical relevance,
- On page 5 the auhors provide the parametric description of the continuum using contravariant coordinates. They refer to these as curvilinear coordinates in general. Since the superscript is used, I would recommend calling these as contravariant and not simply curvilinear coordinates,
- On Section 2.1.1 Shell Configuration the authors expland in very detail in the variational formulation of the Kirchhoff-Love shell. These equations can be found in many resources, see for instance in Kiendl, Josef, et al. "Isogeometric shell analysis with Kirchhoff–Love elements." Computer methods in applied mechanics and engineering 198.49-52 (2009): 3902-3914. I would recommend reducing this part to the absolute minimum and just citing an appropriate resource,
- On page 7 the authors mention the following: "In the one-dimensional case with w = w(x), the bending term reduces further, which leads to a fourth-order differential equation for an Euler–Bernoulli beam when considering the strong formulation of the problem without boundary conditions". I think this is an oversimplification of what is in fact going on. The mathematical and mechanical effect on reducing a structural element from a 3d continuum to a beam introduces many additional mechanical effects (torsion, bending, twisting, etc.) that can not be understood as a reduction from a plate (or shell). I do not see a reason why this section is added and I would recommend removing it (see Bauer, A. M., et al. "Nonlinear isogeometric spatial Bernoulli beam." Computer Methods in Applied Mechanics and Engineering 303 (2016): 101-127. for more information),
- On page 7 the authors mention "We model the lithosphere as a thin elastic plate". Is this a thin elastic plate or a thin shell? The equations the authors laid down in Section 2.1.1 refer to a shell,
- On page 7 the authors refer to a "Kirchhoff-Love plate". To my knowledge "Kirchhoff-Love" is the designation for a shell, the corresponding plate is called just "Kirchhoff",
- On page 8 the authors write the following "Instead of working with the actual topographic elevation, we use a mass representation r obtained by taking the mass above the mid-surface of the lithosphere and normalizing it by some depth-independent reference density ϱr.We choose the reference density as the mean rock density from the current depth of the mid-surface to its initial depth and assume that it is homogeneous". This approach is not clear to me. I would like to ask the authors to add some reference if it is standard in the literature? I am not sure what this modelling implies,
- On page 9 the authors mention the following: "For the full Neumann problem, it can be shown using Korn’s inequality that the variational problem is well-posed, provided that A is a bounded domain with piecewise smooth boundary and the coefficient (ϱm −ϱr)g in the buoyancy term is bounded from below by a positive number.". Is there any reference for this statement? When the authors mention that something can be shown it is advisable to either add a reference or a proof, if there is no reference,
- On pages 9-10 the authors mention the following: "From a modeling point of view, the results may not reflect the physical reality since the Earth consists of different regimes and tectonic plates that interact with each other in a complex manner. Furthermore, due to the large scale of the simulation, the effects of flexural rigidity will not be visible. Nevertheless, we assume that the entire lithosphere can be modeled as a single spherical shell to showcase the capabilities of isogeometric analysis in numerical simulations on curved domains, especially on a spherical domain". Is this to be understood as an academic study without applicability to geoscience? I am not sure how to interpret this statement as it mentions that the results may not reflect the physical reality and that the study aims to showcase the capabilities of isogeometric analysis in numerical simulations on curved domains,
- On page 10 the authors mention the following: "to both parametrize the domain and construct finite element approximations
of solutions to the equations.". I would recommend rephrasing this to "to both parametrize the domain and construct finite element approximations of solutions to the corresponding partial differential equations." - On page 11 the authors mention the following: "In the following, we consider splines on the unit interval [0, 1] with an open knot sequence, i.e., the first and last knot values have multiplicity p+1. Then the knot sequence has the form". I would like to ask the authors to mention the interpolation effect that is enforced when the first and last know values have the prescribed multiplicity. Why are the authors restricting themselves to this choice and what is the implication of doing so?
- In the caption of Figure 5 the authors write "bivariate quadratic ...", I would recommend replacing "quadratic" with "biquadratic",
- Regarding Section 3.1.2 NURBS basis functions, I feel that it is spent quite a lot of real estate on detailing the B-Spline and NURBS basis functions, information that can be found easily in many other resources. Moreover, the authors do not cite these other resources in an adequate manner, but only Pielg and Tiller. I would like to ask the authors to add more references as needed,
- On page 16 the authors mention the following: "We can stack the control points on top of each other". I would recommend rephrasing this to "the control point coordinates are organized in vectors" That sound more appropriate for a scientific publication,
- On the same page the authors mention "There are various methods that can be employed to achieve this." when referring to NURBS multipatches. There are so many publications and studies on this matter, I would like to ask the authors to add adequate references. The reference list is herein quite poor,
- On the same page the authors mention the following: "by replacing shape functions at the boundary of each patch that coincides
with the boundary of another patch with interface functions that span over multiple patches". Could the authors clarify whether this is kind of a bending strip? See Kiendl et al. 2010 - On the same page the authors mention the following: "Another approach is to stitch shape functions at the patch interfaces together by imposing the C1 condition". Could the authors clarifys whether it is herein meant that the C^0 continuity is already assumed (meaning that the patches are conforming along their interface) and that the authors are just constraining the second row of control points in such a way as to achieve C^1 (or G^1) continuity? It is still unclear to me how the C^1 continuity is enforced (is this a Penalty, Lagrange Multipliers, Mortar, or any other approach?),
- On the same page the authors mention the following: "Aside from that, the computation of the null space
corresponding to the C1 constraint is generally a difficult task numerically.". Do the authors herein mean that the coupled system becomes non-convex? - Yet on the same page, the authors mention the following: "In this work, we will mainly consider planar domains that result from joining convex quadrilaterals along the sides. It has been shown that the class of bilinear G1 parametrizations is analysis-suitable, so that optimal convergence can be achieved in this setting". I understood that the authors would model the whole earth's lithosphere as a spherical shell. How can the domains be considered planar in such a case?
- On page 20 the authors mention the following: "It differs from the single-patch result due to the missing data outside of the simulation domain that are replaced by Neumann boundary conditions.". I think it is natural that the results naturally differ. Could the authors clarify whether the point of this comparison is to show that one can obtain satisfactory results also when using less DOFs? If not, could they please clarify what the exact point of this discussion is?
- On page 21 the authors mention the following: "When topographic load is the sought parameter, its initial value is set to 1km". Could the authors clarify why is the value of the topographic load measured in Km? Km should measure distances, right?
- I tried reproducing the results, but I stumbled upon some difficulties with the provided open-source repositories. I am able to run script script_central_java.m provided that all necessary third-party software is installed, but I am unable to execute script script_spherical_earth.m because I receive an error message at line 156 of file op_KL_bending_stress.m from GeoPDEs, namely: "Error using .^ Invalid data type. Argument must be numeric, char, or logical." I tried to debug the issue, but it is probably due to incompatibility of the version of the GeoPDEs library with the one used to produce the results in this study. Could the authors also upload your code also on GitHub and provide clear instructions on the versions of the third-party libraries needed to reproduce the results shown in this study?
- Otherwise, the results section is really great and I am impressed with the results shown in Figure 12. Great work, congratulations!
Citation: https://doi.org/10.5194/egusphere-2024-1093-RC1
Viewed
Total article views: 365 (including HTML, PDF, and XML)
Cumulative views and downloads
(calculated since 29 May 2024)
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 293 | 52 | 20 | 365 | 27 | 18 |
- HTML: 293
- PDF: 52
- XML: 20
- Total: 365
- BibTeX: 27
- EndNote: 18
Viewed (geographical distribution)
Total article views: 394 (including HTML, PDF, and XML)
Thereof 394 with geography defined
and 0 with unknown origin.
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
1
Latest update: 21 Jan 2025
RPTU Kaiserslautern-Landau
Yudi Rosandi
Universitas Padjadjaran
Bernd Simeon
RPTU Kaiserslautern-Landau
Short summary
We model Earth's lithosphere as a thin elastic shell and present numerical methods of isogeometric finite element analysis to simulate its deformation in isostatic equilibrium using technologies from computer-aided design. The simulations also serve as a basis for identifying parameters of the model that are most plausible to explain observed data. This research has been done to showcase the capabilities of isogeometric analysis in solving higher-order problems in geoscientific applications.
We model Earth's lithosphere as a thin elastic shell and present numerical methods of...