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 error​s are measured?It is stated that errors are o​f 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.

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.

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.