the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 30 May 2025
Exploiting Physics-Based Machine Learning to Quantify Geodynamic Effects – Insights from the Alpine Region
Abstract. Geodynamical processes are important to understand and assess the evolution of the Earth system as well as its natural resources. Given the wide range of characteristic spatial and temporal scales of geodynamic processes, their analysis routinely relies on computer-assisted numerical simulations. To provide reliable predictions such simulations need to consider a wide range of potential input parameters, material properties as they vary in space and time, in order to address associated uncertainties. To obtain any quantifiable measure of these uncertainties is challenging both because of the high computational cost of the forward simulation and because data is typically limited to direct observations at the near surface and for the present day state. To account for both of these challenges, we present how to construct efficient and reliable surrogate models that are applicable to a wide range of geodynamic problems using a physics-based machine learning method. In this study, we apply our approach to the case study of the Alpine region, as a natural example for a complex geodynamic setting where several subduction slabs as imaged by tomographic methods interact below a heterogeneous lithosphere. We specifically develop surrogates for two sets of observables, topography and surface velocity, to provide models that can be used in probabilistic frameworks to validate the underlying model structure and parametrization. We additionally construct models for the deeper crustal and mantle domains of the model, to improve the system understanding. For this last family of models, we highlight different construction methods to develop models to either allow evaluations in the entirety of the 3D model or only at specific depth intervals.
Competing interests: At least one of the (co-)authors is a member of the editorial board of Geoscientific Model Development. The peer-review process was guided by an independent editor, and the authors also have 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
(14281 KB)
-
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
- Preprint
(14281 KB) - Metadata XML
- BibTeX
- EndNote
- Final revised paper
Journal article(s) based on this preprint
Interactive discussion
Status: closed
-
RC1: 'Comment on egusphere-2025-1925', Sergio Zlotnik, 16 Jun 2025
-
AC1: 'Reply on RC1', Denise Degen, 17 Jul 2025
Thank you very much for taking the time to review our paper and helping us to improve its current state. Below, we outline how we will address the comments in the manuscript itself.
- Section 2.1: is the forward model linear? or the viscosity \eta is a function of velocities \bf{u}?This is very relevant to the discussion of section 4.1.1 (errors for eta are smaller that for rho)
- We solve for a linear forward Stokes Flow problem, that is, we do not consider any functional behavior of the effective viscosity parameter either in terms of solid velocities nor strain rates. We will add this information to section 2.1 in the revised manuscript.
- Section 2.2: Are the snapshots centered before the svd?Number of snapshots? Size of the test set? This is explained later, but maybe it can be repeated in 2.2.
- Reviewer 1 is right, the snapshots are centered before executing the svd. We will add all required information to section 2.2 in the revised manuscript.
- Section 4.1.1: How local and global errors are measured?It is stated that errors are of order of 1m (tol 1E-4)although in Fig 2 errors are 1 to 5%. How should I interpret these results?
- Global errors are calculated as the RMSE of the back-projected solution. We will add the exact formula to the revised paper accordingly. The problem with global errors is that, locally, the error can exceed these values. Therefore, we additionally investigated the local error behavior (by taking the difference between the full- and reduced order solution) for designated solutions. This we discussed in Fig. 2, where we demonstrate how the maximum local difference is an order of magnitude higher than the global error. We will add a more detailed description to the revised paper to better highlight the differences between these two error metrics and their importance for the evaluation of the approximation quality.
- I could not follow the discussion on the viscosity parameterization. Are there 6 or 2 viscosity parameters?
- The forward model has six viscosity parameters. However, during the construction of the surrogate model, we decided to combine some layers in terms of their viscosity parameterization, and consider only the viscosity of the lithosphere (including slabs) and asthenosphere as free parameters. We will add some clarifying text to the revised manuscript.
- "However, as we will detail in subsection 4.3, this surrogate model is currently not constructible because ofchallenges related to the data set." why?
- The reason behind our statement is that, during the construction of the surrogate model for the case where we combined both density and viscosity, we did find that some of the full order forward solutions showcased velocities and topographies that vary by several orders of magnitude with respect to the other family of solutions. Worth to mention is that these were solutions to fully-solved forward simulations (that is we did not report any crash or premature ending of the forward model) and that they appear to follow a random pattern within the parameter space. The presence of random variations, that is, solutions that do not follow the expected physical behavior, hindered to construct a surrogate model. While addressing the comments received during the first round of review, community post and reviewer#2, we did carry out an in-depth revision of the workflow to extract the results from the full order solutions to construct the surrogate model and we realized that the appearance of these random „unexpected“ solutions was due to an error in the shell script we use to identify the timestep from where to extract the full order solutions (for those cases it was picking an earlier timestep before the system was fully isostatically balance). We would like to stress here that we double check all three scenarios discussed in the manuscript and we found that this error only occurred for scenario 3, and only for few realizations within this scenario, being the reason, we did not spot the problem before. Given that we identified the source of the problem, we can now carry out the construction of the surrogate model for this last scenario as well, which we will integrate in the revised version of the manuscript.
- The use of POD is fairly clear. But for reproducibility it would be necessary that the configuration of the NN is included in the manuscript or in the public data.
- The main hyperparameter of the NN are listed in the tables of Figure 2, additionally we provided a data folder containing the scripts for the surrogate model construction of all surrogate models used in the paper. In the “SurrogateModelScripts” folder is always a file which contains the configuration of the NN. To better highlight this, we will add a description (including references) to the revised manuscript.
Citation: https://doi.org/10.5194/egusphere-2025-1925-AC1 - Section 2.1: is the forward model linear? or the viscosity \eta is a function of velocities \bf{u}?This is very relevant to the discussion of section 4.1.1 (errors for eta are smaller that for rho)
-
AC1: 'Reply on RC1', Denise Degen, 17 Jul 2025
-
CC1: 'ASPECT vs. LaMEM results as used in egusphere-2025-1925', Boris Kaus, 08 Jul 2025
Degen et al. (2025) use geodynamic forward simulations of a reasonably simply 3D model setup of the Alps to train surrogate models. In some cases, they obtain very large surface velocities in the forward models which they attribute to numerical instabilities in the employed software package LaMEM. A comparison with the software package ASPECT gives smaller surface velocities, which they interpret as further confirmation of this conclusion.
In the attached review, I reproduce the results of Degen et al. (2025) and demonstrate that they made a number of major mistakes in their forward models:
- The initial ASPECT model setup is flipped upside down, and they thus simulate the deformation of an overturned lithosphere at 400 km depth, rather than the setup employed in the LaMEM models. A visual inspection of some of the ASPECT results shows "sticky air" bubbles rising through the mantle.
- Moreover, in several cases they appear to have used LaMEM simulation results that were not yet in isostatic balance, which gives very large surface velocities and small topographies, Running the same simulations further in time until they are in isostatic balance results in much smaller velocities.
Whereas the general idea of using geodynamic models to train surrogate models remains an interesting one, the forward models used in the current manuscript need to be redone.
In the current form, the manuscript should not be accepted.
-
AC2: 'Reply on CC1', Denise Degen, 17 Jul 2025
Thank you for your in-depth analysis of our manuscript and, above all, for taking time to look into the forward models presented (especially those „problematic“ ones). All comments raised did help us to identify the source of the problem we encountered in the original study, which we are now ready to correct and address in the revised manuscript.
Reading through the lists of comments raised, we identified two main points, for which we provide in what follows our replies. As a general remark, we noticed that all criticism relates to the forward model solutions. While we acknowledge the points raised, which as stated above, helped us in improving the final application presented in the study, we should point out that the main subject of the manuscript is to describe a novel approach to construct physics-preserving surrogate models, which, for the case at hands, can be applied to long-term geodynamics problems to aid probabilistic (inversion) analysis, including uncertainty quantification. This said, the methodology we presented is only indirectly, if at all, linked to the full-order forward model used to construct interpolation functions for the final mapping.
- Point 1 - The initial ASPECT model setup is flipped upside down, and they thus simulate the deformation of an overturned lithosphere at 400 km depth, rather than the setup employed in the LaMEM models. A visual inspection of some of the ASPECT results shows "sticky air" bubbles rising through the mantle.
- We acknowledge the comment. Indeed, the model was „flipped upside down“. We have corrected it and will integrate the new results in the revised manuscript.
- Point 2 - Moreover, in several cases they appear to have used LaMEM simulation results that were not yet in isostatic balance, which gives very large surface velocities and small topographies, Running the same simulations further in time until they are in isostatic balance results in much smaller velocities.
- Thank you for looking at the simulations and for pointing out the potential source of the inconsistent behavior for those specific forward models. Indeed, that was the case. In order to automatize the extraction of the timestep results from each forward simulation, we relied on a shell script. The targeted timestep coincides with a computation time of 0.27 Myr, which, given also the analysis done in the community post, corresponds to a time where the system is in isostatic balance. This said, we did find an error in the script, leading to the extraction of the results at 0.027 and 0.0027 Myr for a few simulations. Worth to note here is that we have double-checked all cases (also those where we varied single parameters) and found that the error only happened randomly, for very few cases, exactly those for which we spotted unexpected behavior in the final solution (out of isostatic balance). Given that all simulations did run to the end, without any error, and given the random appearance in the parameter space of these unexpected results, we did not detect this error prior to the initial submission. Meanwhile, we have already corrected this error in the shell script, which enabled us to extract the right timestep for all simulations and to allow the construction of a surrogate model for also the last scenario. We will update the manuscript with these new results in the revised version accordingly.
Citation: https://doi.org/10.5194/egusphere-2025-1925-AC2 - Point 1 - The initial ASPECT model setup is flipped upside down, and they thus simulate the deformation of an overturned lithosphere at 400 km depth, rather than the setup employed in the LaMEM models. A visual inspection of some of the ASPECT results shows "sticky air" bubbles rising through the mantle.
-
RC2: 'Comment on egusphere-2025-1925', Anonymous Referee #2, 10 Jul 2025
-
AC3: 'Reply on RC2', Denise Degen, 17 Jul 2025
We would like to thank reviewer 2 for his/her comments.
We would like to repeat, as in the response to the community comment, that the main objective of our study is rather on describing a novel methodology to construct surrogate models for geodynamic simulations and to allow for probabilistic inverse studies and not on the respective geodynamic simulation methods that the surrogate model will be applied to. However, the comments from reviewer 2 are about aspects related to some of the forward simulations performed with the open-source software LaMEM and only applied to the last scenario presented and discussed, scenario 3. Thus, the comments of reviewer 2 exclusively focus on aspects that make up only a small portion of the entire work. We thank reviewer 2 for these comments and want to elaborate with some detail to those: As a matter of fact, the methodology presented in the manuscript is not dependent on the forward dynamic solver adopted. The problems encountered while constructing the surrogate for this last scenario do not impact the scientific merit of the results in principle, and we therefore are convinced that those problems do not „downgrade“ the surrogate model quality per se. As detailed in our answer, we confirm that there was an error in the extraction of the timesteps and we acknowledge all comments raised from reviewer 2, whose expertise with the specific software has helped improving parts of the applications presented. Having identified this error in the script we used to extract the corrected timesteps, and having realized that the error only affected a few simulations for the last scenario, we are now able to construct a surrogate model for this last case in the revised manuscript. Let us disclaim here that the problems encountered, wrong timestep extraction, are solely related to the combined scenario 3. In addition, all simulations, including those for which we did extract the results prematurely, did run without any problem and all data has been already stored on an internal server. Therefore, there is no need to redo any forward computation. We only need to correct the extraction, which we did, in order to finalize the analysis for the last scenario, which we will integrate in the revised version of the manuscript.
Please find below a description of how we will address the comments in the paper.
- The development of surrogate models represents a noteworthy and valuable contribution. However, it is important to express concerns regarding the reliability of the underlying forward models. Several simulations terminated prematurely due to a crash in a third-party library (details are provided in Appendix A). These failed simulations were not correctly identified and were instead interpreted as successful, albeit exhibiting “instabilities” such as low topography and high vertical surface velocities. In reality, these anomalies are attributable to the premature termination of the simulation, which prevented the system from reaching isostatic equilibrium. Substituting the crashing solver component with an alternative resolves the issue entirely.
- First of all, we would like to clarify that all our simulations did not „terminate prematurely“ and we did not face any problem with any specific „third-party library“. The source of the problem stems from an error in the automatic extraction of the results through a shell script. All results were supposed to be extracted at 0.27 Myr, a time when the system is in isostatic balance. However, the error in the shell script leads to an occasional extraction after 0.027 and 0.0027 Myr. Because of the random nature of this occurrence, we did not spot the problem before. Therefore, we cannot agree with the conclusive remark from reviewer#2 stating that the source of errors stems from solvers failure. As already stated in our general comment above, all results from the forward simulations have been already stored, so there is no need to carry out additional „full-order simulations“. A general note of caution here: the problem rarely lies on the library (MUMPS), rather on the way it is used. MUMPS is a well-known and widely used distributed-memory parallel solver (the same apply to superLU_DIST). The reason for us to prefer to rely on MUMPS instead of a more frontal solver as superLU_DIST stems mainly from the former being designed for a relatively broader range of matrices and relying on a dynamic pivoting over hybrid parallelization. This said, it is relatively hard for us to pin-point a likely potential reason for the crashing simulations experienced by reviewer#2: we had too few details on the configuration, architecture, library versions and similar aspect on his/her local installation (a configuration log would help us in this regard) to be able to reproduce the crashing simulation.
- Although the affected simulations are explicitly reported only for the combined training set, it cannot be definitively ruled out that similar problems may exist in the other sets. Therefore, it is strongly recommended that all forward simulations be recomputed to ensure consistency and reliability. As this task exceeds the scope of a standard major revision, the manuscript should be reconsidered for resubmission after the complete regeneration of the training datasets.
- We double checked the extraction regarding all data sets and can confirm that the issue is solely related to the combined scenario 3, where we changed both the density and viscosity. Given that the error is not related to the forward solution, but rather to a wrong extraction of the results, no further forward simulations are required to generate the training datasets for the last scenario. The training dataset for scenario 1 and 2 does not require any modification either.
- To reduce the computational cost associated with regenerating these datasets, an optimized set of input parameters and solver configurations may be employed. Appendix B summarizes a series of such optimizations, which together yield a reduction in run time by more than two orders of magnitude.
- Thank you for the suggestions. However, as stated above, all simulation results for all scenarios, including scenario 3, have been previously stored. Therefore, there is no need to rerun any simulations, but rather to redo the extraction step. A revised version of the data set will be provided in the revised manuscript.
- Rather than including the datasets themselves, the supplementary material should provide scripts used to extract relevant information from the simulation outputs. It is also recommended to include automation scripts for input file generation and forward model execution. Appendix C presents an example Python implementation of this functionality.
- As pointed out in the general comment, the main objective of the paper is to leverage surrogate model construction for geodynamical simulations to enable probabilistic inverse methods in this field. Therefore, the data folder concerning the surrogate models was constructed with this main aim in mind. To ease the surrogate model construction, the simulations results are provided in form of data sets. Our decision was motivated mainly by the wish to enhance ease, from the user side, to reproduce the surrogate model construction and its applications (the main subject of the manuscript) without involving any installation of the forward solver used for generating the training set. Nevertheless, to ease the reproducibility of the forward model, we will include the mentioned shell script in the data folder containing the information about the forward model.
- Detailed comments
- These details comments relate to specific aspects that will be considered and changed, if not completely removed in the revised manuscript.
Citation: https://doi.org/10.5194/egusphere-2025-1925-AC3 -
RC3: 'Reply on AC3', Anonymous Referee #2, 21 Jul 2025
“First of all, we would like to clarify that all our simulations did not „terminate prematurely“ and we did not face any problem with any specific „third-party library“.
To allow an independent assessment of the correctness of this statement, please upload the log files of all forward simulations.
In addition to the log files, please provide the complete software infrastructure used to accomplish the following tasks:
A) Preparation of the input files for each realization
B) Submission of each realization using a workload manager (e.g., slurm)
C) Error checking for forward model runs (this part is absolutely critical)
D) Extraction of input data for surrogate models from forward simulation outputs
According to GMD policy, all of these critical software components must be made publicly available.
“The source of the problem stems from an error in the automatic extraction of the results through a shell script. All results were supposed to be extracted at 0.27 Myr, a time when the system is in isostatic balance. However, the error in the shell script leads to an occasional extraction after 0.027 and 0.0027 Myr.”
This error appears rather puzzling, especially considering that your simulations do not output data at any of the mentioned time steps (0.27, 0.027, or 0.0027 Myr). The closest timestamps I found in the output logs are: 2.796436e-01, 2.964359e-02, and 2.843589e-03. Furthermore, there is no apparent pattern in the problematic realization numbers (82, 89, and 95), which suggests the issue occurs randomly. Please provide a detailed explanation of this bug, along with the actual extraction script.
Why was 0.27 Myr chosen? This duration is not sufficient to reach isostatic balance. Based on your typical model results, I estimate that at least 0.5 Myr is required. To remain on the safe side, 1 Myr should be used. Alternatively, I suggest implementing a custom termination criterion in LaMEM based on a surface velocity threshold. Therefore, I strongly recommend that all forward models be recalculated for a longer simulation period (at least 0.5 Myr, preferably 1 Myr).
Additionally, I recommend increasing the spatial resolution, as some geological units appear to be underrepresented. Ideally, a resolution sensitivity test should be included.
There is no need to upload gigabytes of extracted data, as this can easily be regenerated. Of course, uploading it remains optional.
In any case, the complete software infrastructure related to the forward modeling process must be shared in accordance with GMD policy.Citation: https://doi.org/10.5194/egusphere-2025-1925-RC3
- The development of surrogate models represents a noteworthy and valuable contribution. However, it is important to express concerns regarding the reliability of the underlying forward models. Several simulations terminated prematurely due to a crash in a third-party library (details are provided in Appendix A). These failed simulations were not correctly identified and were instead interpreted as successful, albeit exhibiting “instabilities” such as low topography and high vertical surface velocities. In reality, these anomalies are attributable to the premature termination of the simulation, which prevented the system from reaching isostatic equilibrium. Substituting the crashing solver component with an alternative resolves the issue entirely.
-
AC3: 'Reply on RC2', Denise Degen, 17 Jul 2025
-
AC4: 'Comment on egusphere-2025-1925 Regarding RC3', Denise Degen, 01 Aug 2025
We would like to thank reviewer 2 for his/her comments.
On a general note, we would like to stress again the focus of the study, that is, to showcase a workflow to construct surrogate models that conform to the actual physics at play in order to enable setting up and solving inverse problems as applied to typical geodynamic simulations.
This said, we notice that, unfortunately, the comments raised by reviewer 2 again do focus on a particular aspect related to a subset of forward full-order geodynamic simulations, but do not address the other scientific aspects of the study. While we already stated that, thanks to the community post as well as the first round of comments raised by reviewer 2, we have been able to identify the source of errors for some forward simulations (due to a wrong timestep selection during the postprocessing stage), we would like to clarify that the error made relates to a minor part of the study (a single sub-set of the simulations) and should not be considered as downgrading the scientific merit of the study, which stems from the surrogate model construction and its usability for geodynamic applications. As stated in our first round of replies, we are ready to provide a revised manuscript where the surrogate model will be constructed also for the final subset of simulations, and we hope that with the new addition, the message and merit of the study will be clarified.
In the attached document, we provide a point-by-point answer to all comments raised during the second round of posting.
Interactive discussion
Status: closed
-
RC1: 'Comment on egusphere-2025-1925', Sergio Zlotnik, 16 Jun 2025
These are comments to the manuscript "Exploiting Physics-Based Machine Learning to Quantify Geodynamic Effects – Insights from the Alpine Region" by Denise Degen and co authors submitted to GMD.
The manuscript propose the use of numerical tools based on model order reduction (MOR) to assess the effect of a set of (geodynamic) parameters on the elevation and surface velocity of the Alpine region.
The use or MOR tools allows for a very fast solution of a problem that would be much more expensive with traditional numerical techniques. Therefore, opening the door for inverse problems and studies of the effect of combined parameters that would be impractical otherwise.The combination of techniques proposed (POD+NN) has some advantages with respect to other reduced order methods as it does not require building discrete operators in the evaluation. On the other hand it is more difficult to assess the precision of the surrogate.
The manuscript is very well written and properly structured. It is easy to follow. I think it deserved publication subject to some minor comments that I believe will improve its clarity.
The comments follow.
Section 2.1
is the forward model linear? or the viscosity \eta is a function of velocities \bf{u}?
This is very relevant to the discussion of section 4.1.1
(errors for eta are smaller that for rho)
Section 2.2
Are the snapshots centered before the svd?
Number of snapshots? Size of the test set? This is explained later, but maybe it can be repeated in 2.2.
Section 4.1.1
How local and global errors are measured?
It is stated that errors are of order of 1m (tol 1E-4)
although in Fig 2 errors are 1 to 5%. How should I interpret these results?I could not follow the discussion on the viscosity parameterization. Are there 6 or 2 viscosity parameters?
"However, as we will detail in subsection 4.3, this surrogate model is currently not constructible because of
challenges related to the data set." why?The use of POD is fairly clear. But for reproducibility it would be necessary that the configuration of the NN is included in the manuscript or in the public data.
Citation: https://doi.org/10.5194/egusphere-2025-1925-RC1 -
AC1: 'Reply on RC1', Denise Degen, 17 Jul 2025
Thank you very much for taking the time to review our paper and helping us to improve its current state. Below, we outline how we will address the comments in the manuscript itself.
- Section 2.1: is the forward model linear? or the viscosity \eta is a function of velocities \bf{u}?This is very relevant to the discussion of section 4.1.1 (errors for eta are smaller that for rho)
- We solve for a linear forward Stokes Flow problem, that is, we do not consider any functional behavior of the effective viscosity parameter either in terms of solid velocities nor strain rates. We will add this information to section 2.1 in the revised manuscript.
- Section 2.2: Are the snapshots centered before the svd?Number of snapshots? Size of the test set? This is explained later, but maybe it can be repeated in 2.2.
- Reviewer 1 is right, the snapshots are centered before executing the svd. We will add all required information to section 2.2 in the revised manuscript.
- Section 4.1.1: How local and global errors are measured?It is stated that errors are of order of 1m (tol 1E-4)although in Fig 2 errors are 1 to 5%. How should I interpret these results?
- Global errors are calculated as the RMSE of the back-projected solution. We will add the exact formula to the revised paper accordingly. The problem with global errors is that, locally, the error can exceed these values. Therefore, we additionally investigated the local error behavior (by taking the difference between the full- and reduced order solution) for designated solutions. This we discussed in Fig. 2, where we demonstrate how the maximum local difference is an order of magnitude higher than the global error. We will add a more detailed description to the revised paper to better highlight the differences between these two error metrics and their importance for the evaluation of the approximation quality.
- I could not follow the discussion on the viscosity parameterization. Are there 6 or 2 viscosity parameters?
- The forward model has six viscosity parameters. However, during the construction of the surrogate model, we decided to combine some layers in terms of their viscosity parameterization, and consider only the viscosity of the lithosphere (including slabs) and asthenosphere as free parameters. We will add some clarifying text to the revised manuscript.
- "However, as we will detail in subsection 4.3, this surrogate model is currently not constructible because ofchallenges related to the data set." why?
- The reason behind our statement is that, during the construction of the surrogate model for the case where we combined both density and viscosity, we did find that some of the full order forward solutions showcased velocities and topographies that vary by several orders of magnitude with respect to the other family of solutions. Worth to mention is that these were solutions to fully-solved forward simulations (that is we did not report any crash or premature ending of the forward model) and that they appear to follow a random pattern within the parameter space. The presence of random variations, that is, solutions that do not follow the expected physical behavior, hindered to construct a surrogate model. While addressing the comments received during the first round of review, community post and reviewer#2, we did carry out an in-depth revision of the workflow to extract the results from the full order solutions to construct the surrogate model and we realized that the appearance of these random „unexpected“ solutions was due to an error in the shell script we use to identify the timestep from where to extract the full order solutions (for those cases it was picking an earlier timestep before the system was fully isostatically balance). We would like to stress here that we double check all three scenarios discussed in the manuscript and we found that this error only occurred for scenario 3, and only for few realizations within this scenario, being the reason, we did not spot the problem before. Given that we identified the source of the problem, we can now carry out the construction of the surrogate model for this last scenario as well, which we will integrate in the revised version of the manuscript.
- The use of POD is fairly clear. But for reproducibility it would be necessary that the configuration of the NN is included in the manuscript or in the public data.
- The main hyperparameter of the NN are listed in the tables of Figure 2, additionally we provided a data folder containing the scripts for the surrogate model construction of all surrogate models used in the paper. In the “SurrogateModelScripts” folder is always a file which contains the configuration of the NN. To better highlight this, we will add a description (including references) to the revised manuscript.
Citation: https://doi.org/10.5194/egusphere-2025-1925-AC1 - Section 2.1: is the forward model linear? or the viscosity \eta is a function of velocities \bf{u}?This is very relevant to the discussion of section 4.1.1 (errors for eta are smaller that for rho)
-
AC1: 'Reply on RC1', Denise Degen, 17 Jul 2025
-
CC1: 'ASPECT vs. LaMEM results as used in egusphere-2025-1925', Boris Kaus, 08 Jul 2025
Degen et al. (2025) use geodynamic forward simulations of a reasonably simply 3D model setup of the Alps to train surrogate models. In some cases, they obtain very large surface velocities in the forward models which they attribute to numerical instabilities in the employed software package LaMEM. A comparison with the software package ASPECT gives smaller surface velocities, which they interpret as further confirmation of this conclusion.
In the attached review, I reproduce the results of Degen et al. (2025) and demonstrate that they made a number of major mistakes in their forward models:
- The initial ASPECT model setup is flipped upside down, and they thus simulate the deformation of an overturned lithosphere at 400 km depth, rather than the setup employed in the LaMEM models. A visual inspection of some of the ASPECT results shows "sticky air" bubbles rising through the mantle.
- Moreover, in several cases they appear to have used LaMEM simulation results that were not yet in isostatic balance, which gives very large surface velocities and small topographies, Running the same simulations further in time until they are in isostatic balance results in much smaller velocities.
Whereas the general idea of using geodynamic models to train surrogate models remains an interesting one, the forward models used in the current manuscript need to be redone.
In the current form, the manuscript should not be accepted.
-
AC2: 'Reply on CC1', Denise Degen, 17 Jul 2025
Thank you for your in-depth analysis of our manuscript and, above all, for taking time to look into the forward models presented (especially those „problematic“ ones). All comments raised did help us to identify the source of the problem we encountered in the original study, which we are now ready to correct and address in the revised manuscript.
Reading through the lists of comments raised, we identified two main points, for which we provide in what follows our replies. As a general remark, we noticed that all criticism relates to the forward model solutions. While we acknowledge the points raised, which as stated above, helped us in improving the final application presented in the study, we should point out that the main subject of the manuscript is to describe a novel approach to construct physics-preserving surrogate models, which, for the case at hands, can be applied to long-term geodynamics problems to aid probabilistic (inversion) analysis, including uncertainty quantification. This said, the methodology we presented is only indirectly, if at all, linked to the full-order forward model used to construct interpolation functions for the final mapping.
- Point 1 - The initial ASPECT model setup is flipped upside down, and they thus simulate the deformation of an overturned lithosphere at 400 km depth, rather than the setup employed in the LaMEM models. A visual inspection of some of the ASPECT results shows "sticky air" bubbles rising through the mantle.
- We acknowledge the comment. Indeed, the model was „flipped upside down“. We have corrected it and will integrate the new results in the revised manuscript.
- Point 2 - Moreover, in several cases they appear to have used LaMEM simulation results that were not yet in isostatic balance, which gives very large surface velocities and small topographies, Running the same simulations further in time until they are in isostatic balance results in much smaller velocities.
- Thank you for looking at the simulations and for pointing out the potential source of the inconsistent behavior for those specific forward models. Indeed, that was the case. In order to automatize the extraction of the timestep results from each forward simulation, we relied on a shell script. The targeted timestep coincides with a computation time of 0.27 Myr, which, given also the analysis done in the community post, corresponds to a time where the system is in isostatic balance. This said, we did find an error in the script, leading to the extraction of the results at 0.027 and 0.0027 Myr for a few simulations. Worth to note here is that we have double-checked all cases (also those where we varied single parameters) and found that the error only happened randomly, for very few cases, exactly those for which we spotted unexpected behavior in the final solution (out of isostatic balance). Given that all simulations did run to the end, without any error, and given the random appearance in the parameter space of these unexpected results, we did not detect this error prior to the initial submission. Meanwhile, we have already corrected this error in the shell script, which enabled us to extract the right timestep for all simulations and to allow the construction of a surrogate model for also the last scenario. We will update the manuscript with these new results in the revised version accordingly.
Citation: https://doi.org/10.5194/egusphere-2025-1925-AC2 - Point 1 - The initial ASPECT model setup is flipped upside down, and they thus simulate the deformation of an overturned lithosphere at 400 km depth, rather than the setup employed in the LaMEM models. A visual inspection of some of the ASPECT results shows "sticky air" bubbles rising through the mantle.
-
RC2: 'Comment on egusphere-2025-1925', Anonymous Referee #2, 10 Jul 2025
-
AC3: 'Reply on RC2', Denise Degen, 17 Jul 2025
We would like to thank reviewer 2 for his/her comments.
We would like to repeat, as in the response to the community comment, that the main objective of our study is rather on describing a novel methodology to construct surrogate models for geodynamic simulations and to allow for probabilistic inverse studies and not on the respective geodynamic simulation methods that the surrogate model will be applied to. However, the comments from reviewer 2 are about aspects related to some of the forward simulations performed with the open-source software LaMEM and only applied to the last scenario presented and discussed, scenario 3. Thus, the comments of reviewer 2 exclusively focus on aspects that make up only a small portion of the entire work. We thank reviewer 2 for these comments and want to elaborate with some detail to those: As a matter of fact, the methodology presented in the manuscript is not dependent on the forward dynamic solver adopted. The problems encountered while constructing the surrogate for this last scenario do not impact the scientific merit of the results in principle, and we therefore are convinced that those problems do not „downgrade“ the surrogate model quality per se. As detailed in our answer, we confirm that there was an error in the extraction of the timesteps and we acknowledge all comments raised from reviewer 2, whose expertise with the specific software has helped improving parts of the applications presented. Having identified this error in the script we used to extract the corrected timesteps, and having realized that the error only affected a few simulations for the last scenario, we are now able to construct a surrogate model for this last case in the revised manuscript. Let us disclaim here that the problems encountered, wrong timestep extraction, are solely related to the combined scenario 3. In addition, all simulations, including those for which we did extract the results prematurely, did run without any problem and all data has been already stored on an internal server. Therefore, there is no need to redo any forward computation. We only need to correct the extraction, which we did, in order to finalize the analysis for the last scenario, which we will integrate in the revised version of the manuscript.
Please find below a description of how we will address the comments in the paper.
- The development of surrogate models represents a noteworthy and valuable contribution. However, it is important to express concerns regarding the reliability of the underlying forward models. Several simulations terminated prematurely due to a crash in a third-party library (details are provided in Appendix A). These failed simulations were not correctly identified and were instead interpreted as successful, albeit exhibiting “instabilities” such as low topography and high vertical surface velocities. In reality, these anomalies are attributable to the premature termination of the simulation, which prevented the system from reaching isostatic equilibrium. Substituting the crashing solver component with an alternative resolves the issue entirely.
- First of all, we would like to clarify that all our simulations did not „terminate prematurely“ and we did not face any problem with any specific „third-party library“. The source of the problem stems from an error in the automatic extraction of the results through a shell script. All results were supposed to be extracted at 0.27 Myr, a time when the system is in isostatic balance. However, the error in the shell script leads to an occasional extraction after 0.027 and 0.0027 Myr. Because of the random nature of this occurrence, we did not spot the problem before. Therefore, we cannot agree with the conclusive remark from reviewer#2 stating that the source of errors stems from solvers failure. As already stated in our general comment above, all results from the forward simulations have been already stored, so there is no need to carry out additional „full-order simulations“. A general note of caution here: the problem rarely lies on the library (MUMPS), rather on the way it is used. MUMPS is a well-known and widely used distributed-memory parallel solver (the same apply to superLU_DIST). The reason for us to prefer to rely on MUMPS instead of a more frontal solver as superLU_DIST stems mainly from the former being designed for a relatively broader range of matrices and relying on a dynamic pivoting over hybrid parallelization. This said, it is relatively hard for us to pin-point a likely potential reason for the crashing simulations experienced by reviewer#2: we had too few details on the configuration, architecture, library versions and similar aspect on his/her local installation (a configuration log would help us in this regard) to be able to reproduce the crashing simulation.
- Although the affected simulations are explicitly reported only for the combined training set, it cannot be definitively ruled out that similar problems may exist in the other sets. Therefore, it is strongly recommended that all forward simulations be recomputed to ensure consistency and reliability. As this task exceeds the scope of a standard major revision, the manuscript should be reconsidered for resubmission after the complete regeneration of the training datasets.
- We double checked the extraction regarding all data sets and can confirm that the issue is solely related to the combined scenario 3, where we changed both the density and viscosity. Given that the error is not related to the forward solution, but rather to a wrong extraction of the results, no further forward simulations are required to generate the training datasets for the last scenario. The training dataset for scenario 1 and 2 does not require any modification either.
- To reduce the computational cost associated with regenerating these datasets, an optimized set of input parameters and solver configurations may be employed. Appendix B summarizes a series of such optimizations, which together yield a reduction in run time by more than two orders of magnitude.
- Thank you for the suggestions. However, as stated above, all simulation results for all scenarios, including scenario 3, have been previously stored. Therefore, there is no need to rerun any simulations, but rather to redo the extraction step. A revised version of the data set will be provided in the revised manuscript.
- Rather than including the datasets themselves, the supplementary material should provide scripts used to extract relevant information from the simulation outputs. It is also recommended to include automation scripts for input file generation and forward model execution. Appendix C presents an example Python implementation of this functionality.
- As pointed out in the general comment, the main objective of the paper is to leverage surrogate model construction for geodynamical simulations to enable probabilistic inverse methods in this field. Therefore, the data folder concerning the surrogate models was constructed with this main aim in mind. To ease the surrogate model construction, the simulations results are provided in form of data sets. Our decision was motivated mainly by the wish to enhance ease, from the user side, to reproduce the surrogate model construction and its applications (the main subject of the manuscript) without involving any installation of the forward solver used for generating the training set. Nevertheless, to ease the reproducibility of the forward model, we will include the mentioned shell script in the data folder containing the information about the forward model.
- Detailed comments
- These details comments relate to specific aspects that will be considered and changed, if not completely removed in the revised manuscript.
Citation: https://doi.org/10.5194/egusphere-2025-1925-AC3 -
RC3: 'Reply on AC3', Anonymous Referee #2, 21 Jul 2025
“First of all, we would like to clarify that all our simulations did not „terminate prematurely“ and we did not face any problem with any specific „third-party library“.
To allow an independent assessment of the correctness of this statement, please upload the log files of all forward simulations.
In addition to the log files, please provide the complete software infrastructure used to accomplish the following tasks:
A) Preparation of the input files for each realization
B) Submission of each realization using a workload manager (e.g., slurm)
C) Error checking for forward model runs (this part is absolutely critical)
D) Extraction of input data for surrogate models from forward simulation outputs
According to GMD policy, all of these critical software components must be made publicly available.
“The source of the problem stems from an error in the automatic extraction of the results through a shell script. All results were supposed to be extracted at 0.27 Myr, a time when the system is in isostatic balance. However, the error in the shell script leads to an occasional extraction after 0.027 and 0.0027 Myr.”
This error appears rather puzzling, especially considering that your simulations do not output data at any of the mentioned time steps (0.27, 0.027, or 0.0027 Myr). The closest timestamps I found in the output logs are: 2.796436e-01, 2.964359e-02, and 2.843589e-03. Furthermore, there is no apparent pattern in the problematic realization numbers (82, 89, and 95), which suggests the issue occurs randomly. Please provide a detailed explanation of this bug, along with the actual extraction script.
Why was 0.27 Myr chosen? This duration is not sufficient to reach isostatic balance. Based on your typical model results, I estimate that at least 0.5 Myr is required. To remain on the safe side, 1 Myr should be used. Alternatively, I suggest implementing a custom termination criterion in LaMEM based on a surface velocity threshold. Therefore, I strongly recommend that all forward models be recalculated for a longer simulation period (at least 0.5 Myr, preferably 1 Myr).
Additionally, I recommend increasing the spatial resolution, as some geological units appear to be underrepresented. Ideally, a resolution sensitivity test should be included.
There is no need to upload gigabytes of extracted data, as this can easily be regenerated. Of course, uploading it remains optional.
In any case, the complete software infrastructure related to the forward modeling process must be shared in accordance with GMD policy.Citation: https://doi.org/10.5194/egusphere-2025-1925-RC3
- The development of surrogate models represents a noteworthy and valuable contribution. However, it is important to express concerns regarding the reliability of the underlying forward models. Several simulations terminated prematurely due to a crash in a third-party library (details are provided in Appendix A). These failed simulations were not correctly identified and were instead interpreted as successful, albeit exhibiting “instabilities” such as low topography and high vertical surface velocities. In reality, these anomalies are attributable to the premature termination of the simulation, which prevented the system from reaching isostatic equilibrium. Substituting the crashing solver component with an alternative resolves the issue entirely.
-
AC3: 'Reply on RC2', Denise Degen, 17 Jul 2025
-
AC4: 'Comment on egusphere-2025-1925 Regarding RC3', Denise Degen, 01 Aug 2025
We would like to thank reviewer 2 for his/her comments.
On a general note, we would like to stress again the focus of the study, that is, to showcase a workflow to construct surrogate models that conform to the actual physics at play in order to enable setting up and solving inverse problems as applied to typical geodynamic simulations.
This said, we notice that, unfortunately, the comments raised by reviewer 2 again do focus on a particular aspect related to a subset of forward full-order geodynamic simulations, but do not address the other scientific aspects of the study. While we already stated that, thanks to the community post as well as the first round of comments raised by reviewer 2, we have been able to identify the source of errors for some forward simulations (due to a wrong timestep selection during the postprocessing stage), we would like to clarify that the error made relates to a minor part of the study (a single sub-set of the simulations) and should not be considered as downgrading the scientific merit of the study, which stems from the surrogate model construction and its usability for geodynamic applications. As stated in our first round of replies, we are ready to provide a revised manuscript where the surrogate model will be constructed also for the final subset of simulations, and we hope that with the new addition, the message and merit of the study will be clarified.
In the attached document, we provide a point-by-point answer to all comments raised during the second round of posting.
Peer review completion
Journal article(s) based on this preprint
Data sets
Non-Intrusive Reduced Basis Method - Case Study of the Alpine Region Denise Degen, Ajay Kumar, Magdalena Scheck-Wenderoth, Mauro Cacace https://doi.org/10.5281/zenodo.14755256
LaMEM and ASPECT input and data files corresponding to Exploiting Physics-Based Machine Learning to Quantify Geodynamic Effects – Insights from the Alpine Region Ajay Kumar https://doi.org/10.5281/zenodo.15478977
Model code and software
Non-Intrusive Reduced Basis Method - Case Study of the Alpine Region Denise Degen, Ajay Kumar, Magdalena Scheck-Wenderoth, Mauro Cacace https://doi.org/10.5281/zenodo.14755256
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 2,154 | 110 | 31 | 2,295 | 33 | 45 |
- HTML: 2,154
- PDF: 110
- XML: 31
- Total: 2,295
- BibTeX: 33
- EndNote: 45
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
Ajay Kumar
Magdalena Scheck-Wenderoth
Mauro Cacace
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
- Preprint
(14281 KB) - Metadata XML
These are comments to the manuscript "Exploiting Physics-Based Machine Learning to Quantify Geodynamic Effects – Insights from the Alpine Region" by Denise Degen and co authors submitted to GMD.
The manuscript propose the use of numerical tools based on model order reduction (MOR) to assess the effect of a set of (geodynamic) parameters on the elevation and surface velocity of the Alpine region.
The use or MOR tools allows for a very fast solution of a problem that would be much more expensive with traditional numerical techniques. Therefore, opening the door for inverse problems and studies of the effect of combined parameters that would be impractical otherwise.
The combination of techniques proposed (POD+NN) has some advantages with respect to other reduced order methods as it does not require building discrete operators in the evaluation. On the other hand it is more difficult to assess the precision of the surrogate.
The manuscript is very well written and properly structured. It is easy to follow. I think it deserved publication subject to some minor comments that I believe will improve its clarity.
The comments follow.
Section 2.1
is the forward model linear? or the viscosity \eta is a function of velocities \bf{u}?
This is very relevant to the discussion of section 4.1.1
(errors for eta are smaller that for rho)
Section 2.2
Are the snapshots centered before the svd?
Number of snapshots? Size of the test set? This is explained later, but maybe it can be repeated in 2.2.
Section 4.1.1
How local and global errors are measured?
It is stated that errors are of order of 1m (tol 1E-4)
although in Fig 2 errors are 1 to 5%. How should I interpret these results?
I could not follow the discussion on the viscosity parameterization. Are there 6 or 2 viscosity parameters?
"However, as we will detail in subsection 4.3, this surrogate model is currently not constructible because of
challenges related to the data set." why?
The use of POD is fairly clear. But for reproducibility it would be necessary that the configuration of the NN is included in the manuscript or in the public data.