Thank you very much for taking the time to review our manuscript and providing us with the comments to improve the current state of the manuscript. Please find below a point-by-point answer to your comments.

Specific comments

1. There is a language mistake in the abstract in the sentence on l. 4-5.
   • The language mistake has been corrected.
2. I do not fully understand how the need to consider rapid changes in state distribution follows from the description of models 1-3 in the introduction as stated in the sentence on l. 68, starting with "Consequently, ...".
   • Thank you for pointing this out. We removed the “Consequently, …” from the sentence since the need to consider rapid changes arises from the application itself and not from the different surrogate modeling techniques.
3. In section 2.3 it would improve readability if it could be clearly stated what is being modeled as this would link the section better to the previous sections. While section 2.1 refers to "common assumption" and data that is commonly available and not available,  it is not fully clear how this links to the assumptions used in the current work. From reading the previous sections, one could guess that it is the four variables describing the stress state (the reduced stress tensor). Later, however, it appears that it is only S_hmin, S_Hmax and S_v which are modeled. The type of heterogeneity that is assumed or allowed in the model must also be stated.
   • Thank you for your comment. We added a clarification of what is model and they assumptions regarding the heterogeneity at the start of section 2.3. “Note that we focus on the construction of surrogate models for the S_hmin, S_Hmax, and S_v For better illustration of the general concepts, we assume that the material properties are homogeneous and isotropic within each layer. However, the presented concepts are not restricted to these assumptions. Further details regarding the model setup are listed in section 3.”
   • We also added an additional explanation how the data is used within the surrogate model construction, to better clarify the methodology. “The input for the GPR machine learning algorithm is the product of the basis functions and the training snapshots. The basis functions ψ provide the characteristic behavior of the model. The training snapshots are a controlled environment with known μ. This allows the GPR to derive the reduced coefficients θ_rb that complete the mapping between the input and output space since both are known for the training data.”
4. In section 2.2, I think the constitutive law for linear elasticity and the notation for the related rock parameters should be stated for completeness. Now, Poisson's ratio and Young's modulus are defined, but this is in the caption of Figure 3. I think it would improve readability to define these central parameters when the governing equations are introduced. 
   • The notations of the related rock parameters have been added as well as the constitutive relationships.
5. How the boundary conditions are defined is not fully clear to me. Based on what is written it seems to be Dirichlet conditions (with the indicated variability) at the top, but the conditions at the other boundaries are not clear to me. It would be good if the boundary conditions at each face of the domain could be precisely stated. If I understand correctly, the boundary conditions are chosen so that the principal axis of the stress align with the coordinate system. 
   • The definition of the boundary conditions has been added. “Throughout the study, we allow variations for the boundary conditions, the material properties, and the interface depths. For the displacement in the y-direction (eastern boundary), we apply a Dirichlet boundary condition with values varying between four to six meters, whereas the variations for the x-direction range from 0.2 m to 0.6 m (northern boundary). Different boundary conditions are needed to ensure reasonable variations of the stress magnitudes as observed in data records. The top boundary is assigned with a zero Neumann boundary condition, and all remaining boundaries are subjected to zero Dirichlet boundary conditions normal to the boundary (roller boundary conditions).”
6. I was missing the definition of the notation introduced in Fig. 2, but I realized it comes later (related to eq. (5)), to improve readability, the notation should be introduced when it is first used. 
   • The definition of the notation has been added to the caption of Fig. 2 to improve the readability.
7. The definition of the lambdas is missing related to eq (4). 
   • The lambdas denote the singular values, a corresponding clarification has been added.
8. When I read the section on the synthetic model, I got curious about how the methodology would perform considering layering which is not uniform in depth. It would be nice to comment on this. This also relates to my previous comment on heterogeneity.
   • We added a brief discussion about the extension to more complex geometrical models in the Discussion section. “The geometry considered in this study is simple, to focus on the surrogate model construction itself and its implications for geomechanical modeling. Both the non-intrusive and intrusive versions of the reduced basis method have been applied to geothermal real-case studies with a complex geometrical setup, demonstrating a similar performance as presented in this study (Degen and Cacace, 2021; Degen et al., 2021a, b, 2022a, b, c). This implies, that especially for fixed geometries, no degradation in the surrogate model quality or performance is expected. Also, considering variations for complex geometries is possible. However, one should keep in mind that the base assumption for the method is the existence of a low-dimensional parameter space. Hence, if too many geometrical parameters are varied at once, this assumption breaks down, and the approach will become inefficient.” 
9. In the pdf of the manuscript that I downloaded, figure 4 a) and b) and related figures are a bit blurry. It seems like the FE simulations and the NI-RB simulations have a very good agreement, but it is a bit difficult to see and due to the resolution it did not help to zoom in on the figure.
   • All images regarding the comparison of the FE and NI-RB simulation have been saved as vector graphics. The problem originated from the conversion of the program. Therefore, all images concerning this issue have been re-exported to improve the quality.
10. It is not clear to me how "Model Error" is defined, it would be good to have the exact definition. 
    • To calculate the model error, we use the mean squared error. An explanation has been added.
11. Based on the information provided, it is not fully clear to me what is shown in Figure 4 c) and d), it would be good if a definition of "sensitivity" was provided. It would also be good with a comment explaining the "equal" and "non-equal" "x- and y-strains". 
    • Sensitivity refers to the value of the sensitivity index. An explanation has been added in the captions of all figures containing the expression. Further, explanations are found in Section 2.4, a reference to this section has been inserted. Section 2.4 contains several papers explaining the concept of sensitivity indices in greater detail, which is beyond the scope of this paper. A clarification regarding the equal and non-equal strain scenarios has been added.
12. There is a mistake in the sentence starting on l. 386.
    • The language mistake has been corrected.
13. Caption is missing for Figure 9 d).
    • The caption for Figure 9d) has been added.
14. I do not understand the bias in the sensitivity analysis that according to the authors is caused by layer thickness and can be mitigated by choosing the same number of elements for each layer. While I can see that the results change if the same number of elements are chosen for each layer, why this is the case is not clear.  It is also not clear what "cells" refer to. If the same is meant by "cells" and "elements" it would be good to stick with one terminology.

• We added an explanation regarding the cause of the bias. “To elaborate, consider a model consisting of two layers, where the top layer has a thickness of 500 m and the base layer has a thickness of 200 m. If we discretize the model with a vertical resolution of 100 m, we require five elements for the first layer and two elements for the second layer. Calculating the sensitivity indices takes the variation of the stress distribution in each element into consideration. Consequently, the top layer will have more elements contributing and, because of that, likely a higher impact. By discretizing the model with an equal number of elements per layer (e.g., five elements for both the top and base layer), we can avoid this effect. For further information regarding potential biases, we refer to Degen and Wellmann (2024).”
  • Cells and elements refer indeed to the same thing. We unified the notation to avoid any misunderstandings.

1. The boundary conditions for the Case study presented in Section 4.6 should be stated.
   • We added a description for the boundary conditions.
2. In the discussion, it is commented (l. 474) that large shifts in the state can occur also for hydrological studies when the permeability contrast is large (e.g. pressure discontinuity if there are almost impermeable zones). This is correct, but it is not in agreement with what was stated in the introduction on l. 72. I think some information needs to be added in order to avoid confusion.
   • Thank you for pointing this out. This was indeed misleadingly formulated. We reformulated the sentence in the introduction and added further explanations. “Especially, the later part distinguishes the work significantly from previous studies in geothermal applications (e.g., Degen et al., 2022a), where only smooth state variable distributions of the pore pressure and/or temperature have been considered so far. Depending on the permeability contrast, the pore pressure exhibits rapid changes as well. However, these scenarios have not yet been investigated with respect to the non-intrusive reduced basis method.”
3. In the discussion, l. 460. it is referred to visible deviations between the reduced and full-order models, but I am not able to see this in the figure.
   • The deviations are quite small. We added a description of where to observe examples of these deviations.
4. It would be good to include a discussion on how this methodology can/is be combined with real data on stress magnitudes, e.g. from downhole experiments. A short comment is given in the sentence starting on l. 111 but it is not clear to me if or how this is done in the present work. 
   • We extended the discussion to include the potential incorporation of real data on the stress magnitudes in form of probabilistic uncertainty quantification methods.

Thank you very much for taking the time to review our manuscript and providing us with the comments to improve the current state of the manuscript. Please find below a point-by-point answer to your comments.

1. The manuscript does not justify the chosen threshold for identifying influential parameters, which may limit interpretability for other applications.
   • Unfortunately, there is no general way of determining the threshold way. We did provide a justification at the end of section 2.4. To further clarify the point, we extended the given explanation. This contains also a reference to a recent preprint, where the aspect of the threshold is discussed in detail. “There is no general way of determining this threshold. For the purpose of this study, we set it to 10^-2 since values below this threshold would be difficult to validate against typical in-situ stress measurement accuracies (Desroches et al., 2021b; Martin, 2007; Morawietz et al., 2020). his value is also in accordance with threshold values typically employed in literature (Cosenza et al., 2013; Degen et al., 2021a, b; Degen and Wellmann (2024); Sin et al., 2011; Tang et al., 2006; Vanrolleghem et al., 2015). The threshold is a dimensionless number, which is determined by the division of two variances, for further regarding the threshold itself and its determination, we refer to Degen and Wellmann (2024).”
2. The assumption of linear elasticity is appropriate for the presented examples but limits generalizability. A brief discussion of potential extensions to nonlinear constitutive behavior, such as elasto-plasticity, would be beneficial.

• We added a brief discussion about the extension to both nonlinear and coupled applications in the Discussion section. “As mentioned before, the non-intrusive reduced basis method has been developed in order to provide an efficient extension of the RB method for nonlinear applications (Hesthaven and Ubbiali, 2018). In contrast to this study, Degen et al. (2022a) consider fully coupled thermo-hydro-mechanical simulations, demonstrating that the current approach is also extendable to nonlinear and coupled relationships. The implications for nonlinear applications are further detailed in Degen et al. (2023) for different subsurface applications, including the highly nonlinear Richard's equations. The results of these studies demonstrate great promise also for nonlinear geomechanical applications, which is especially important when considering potential extensions to elasto-plasticity.”

1. Figure 2 could benefit from additional annotations explaining the transitions between steps in the NI-RB method. Table 1 lacks justification for the chosen parameter ranges.
   • We added additional annotations describing the transition between the various steps in Fig. 2.
   • The parameter ranges in Table 1 are mostly inspired by the latter case study. An explanation, along with the corresponding references, has been added to Table 1.
2. The manuscript lacks direct comparisons between the NI-RB method and alternative surrogate modeling techniques, such as Proper Orthogonal Decomposition with Galerkin projection or PINNs. Including this context would provide a better understanding of the NI-RB method's relative advantages and limitations.
   • We provided a detailed comparison of the NI-RB method with other surrogate modeling techniques in two recent publications and no conceptual differences for this comparison is expected in this study. We added an explanation, as well as the reference to previous studies.
3. The manuscript does not explore potential limitations, such as handling highly nonlinear systems or strong parameter couplings.
   • See answer for comment 2.
4. The global sensitivity results do not explicitly linked to decision-making processes, such as optimizing nuclear waste repository design.
   • Thank you for raising this point. We indeed do not link the sensitivity analyses to decision-making processes and rather focus on the conceptual aspects to introduce this novel surrogate modeling technique for the community that use 3D geomechanical-numerical models. The linkage to decision-making processes is a very interesting aspect, which is subject to future work and would involve a 3D case study. An explanation has been added.