We sincerely thank the reviewer for the thorough evaluation and valuable feedback on our manuscript, and for the encouraging words. We have carefully addressed every point raised and made significant changes and additions to the manuscript. 

Reviewer 1 (Q1): “I would encourage the authors to double check for typos in equations before the final revision.” 

Response (Q1): We sincerely thank the reviewer for the positive feedback and the helpful suggestions; we have carefully re-checked the equations for any potential typos and inconsistencies.

Reviewer 1 (Q2): “[...]  I imagine different instances of ice have inherent variability in parameters such as A, E and nu, and I am curious whether there is enough information in the literature to constrain a reasonable range of these parameters in natural systems. Do they range over orders of magnitude? Are there any physically based ways of estimating them? This is largely beyond the scope of this paper, but a little bit of information in the discussion might be nice."

Response (Q2):  We thank the reviewer for raising this important point. To the best of our knowledge, currently there is no standard calibration protocol for triaxial ice rheology. Moreover, any calibration program, and the calibrated values, will depend on the assumed constitutive model.  For instance, assuming Maxwell vs. Kelvin configuration leads to different calibration protocols and potentially different calibrate parameters.  On the other hand, the existing experimental datasets on ice mechanics are often exploratory in nature, and not designed for calibrating specific models. With that being said, we tried to resort to a semi-systematic calibration procedure (as explained in the paper) which allowed us to directly calibrate some of the important parameters, such as those for Glen–Nye's law.

Action (Q2): The point regarding the reasonable ranges is very valid. We will add the following paragraph to the Discussion section:

“The parameters in Table 1 are calibrated for isotropic, freshwater polycrystalline ice. Under this assumption the literature indicates a comparatively narrow band for elastic constants, with Sinha (1989) reporting E ≈ 8.93-9.53GPa and ν ≈ 0.308–0.311 over 0 to -50°C. The Maxwell configuration excludes the reversible rate-dependent deformation (primary creep), which is known to lead to lower Young’s modulus values in comparison to commonly reported values. For the Glen–Arrhenius set {A,n,Q}, reported values show bounded variability with n ≈ 2-4, Q ≈ 60-120 kJ⋅mol^-1 while A spans orders of magnitude across studies due to microstructure and differences in calibration windows (Glen, 1955; Barnes et al., 1971; Zeitz et al., 2020). Within this context, our calibrated {A,n,Q}lie within established intervals and are consistent with ranges reported by Durham et al. (1983) for polycrystalline ice.”

Reviewer 1 (Q3): “[...] a 1km thick glacier might experience a pressure at the base on order of tens of MPa. This nicely falls within the regime of parameters used in this study. […] it might be worth mentioning a few key papers.”

Action (Q3): We thank the reviewer for pointing out this connection. We completely agree and will add the following to the newly added Limitations and Range of Applicability section:

“The triaxial behavior has been compared only at relatively high mean stress levels, on the order of 10~MPa, due to the limited availability of experimental data. These stress magnitudes are particularly relevant, as they are comparable to basal stress levels typically encountered beneath glaciers and known to influence both sliding dynamics and glacial erosion (Herman et al., 2015, 2021; Cook et al., 2020).”

Reviewer 1 (Q4): “It would be nice to include some discussion about what would be the most useful future studies- especially experimental- that could be done to improve our understanding of ice mechanics in light of the findings in this paper."

Action (Q4): Again, we agree with the reviewer and will add the following passage.

“Dedicated experiments are needed to determine which rheological framework (Maxwell, Kelvin, Burgers, etc.) most closely represents ice behavior. Moreover, although the viscous-driven healing process proposed here is physically appealing, it remains unverified by experiments. Targeted experiments, particularly load–hold–reload and low-frequency cyclic triaxial tests conducted under controlled confinement and temperature, are essential to confirm the potential correlation between strain rate and stiffness recovery.”

We sincerely thank the reviewer for the evaluation and valuable feedback on our manuscript. We have carefully addressed every point raised and made significant changes and additions to the manuscript. 

Reviewer 2 (Q1): “Provide a more detailed discussion of how temperature influences the mechanical properties and overall model response. As an example the authors considered constant elastic modulus and Poisson's ratio. However, in frozen soil, the mechanical properties of the medium vary as a function of temperature or ice saturation."

Response (Q1): It seems from this (and other comments) that the reviewer is comparing the results of our study with frozen soils, most probably because frozen soil literature often involves ice-soil mixtures. However, our model is developed for pure ice and not mixture of ice and soil. Concepts such as ice saturation are only applicable to the latter not the former. The temperature-dependent behavior is already explicitly represented in the model through two mechanisms: (i) the Arrhenius-type viscosity law and (ii) the temperature-dependent damage initiation threshold, see Eqs. 9-12. Regarding the elastic constants, we have clarified explicitly that their potential temperature dependency is ignored in this study, see the discussion following Eq. 6. This is because, as explained in the original manuscript, experimental data show that the variation of elastic parameters with temperature is indeed negligible. Sinha (1989) reported that between 233 K and 273 K, the Young's modulus of ice increases by about 5%, and the Poisson's ratio by only 1%, which is negligible compared with the temperature-dependency of viscosity and damage. Therefore, assuming constant  and  introduces minimal error while simplifying the analysis.

Reviewer 2 (Q2): The incremental form of the constitutive model is solved using an explicit numerical integration scheme. Please specify the computational platform or software (e.g., in-house code, ABAQUS, COMSOL, etc.) used for the implementation. Additionally, discuss whether any stability issues were encountered given the explicit formulation and how these were mitigated.

Response and Action (Q2): The model was implemented in an in-house Python code using an explicit integration scheme. No numerical instabilities were encountered, and sensitivity tests confirmed that further step-size reduction produced no meaningful change in the results. The following sentence will be added to the revised manuscript in response to this comment:

“The system of equations are solved using explicit numerical integration method. Sensitivity tests were performed to confirm that the selected timesteps were sufficiently small, and further step-size reduction produced no meaningful change in the results.”

Reviewer 2 (Q3&4): “Clarify whether the liquid water and ice contents are assumed constant. If the volumetric ratio between the phases changes, how would this affect the governing equations and model predictions? Is the medium assumed to be fully saturated? If so, please elaborate on how the model would perform or need to be modified under unsaturated conditions."

Response (Q3&4): Similar to the first comment, it seems the developed model for pure ice in this study is being conceived as for frozen soil, which has led to some confusion. The present constitutive formulation describes the thermomechanical behavior of single-phase polycrystalline ice. The material is treated as a continuous solid without any soil skeleton, liquid water, or air phases; hence, concepts such as saturation or evolving phase fractions do not apply.

Reviewer 2 (Q5): “A dedicated section discussing the limitations and applicability range of the proposed model would be beneficial.”

Action (Q5): Thank you for this suggestion. We will add the following section to the revised manuscript:

“The proposed model reproduces the main features of the triaxial response within the calibrated window, including initial stiffness, peak strength, and post-peak softening. The formulation nevertheless relies on simplifying assumptions that limit its applicability. It assumes isotropy, while elasticity, creep, and damage in polycrystalline ice depend on crystallographic fabric and on the alignment of microcracks, which can produce direction-dependent stiffness and flow (Alley, 1988; Azuma, 1995; Pralong et al., 2006). Accuracy of prediction is also expected to decline near the pressure melting point because additional thermally activated mechanisms, including grain boundary processes and recrystallization or enhanced mobile dislocation activity, can lead to faster than Arrhenius creep (Morgan, 1991; Cole, 2020). Moreover, microstructural and compositional influences such as salinity and brine content, porosity, and grain size and its evolution are not represented even though they can modify effective stiffness and the operative creep mechanisms in natural ice (Langleben and Pounder, 1963; Goldsby and Kohlstedt, 2001). The proposed viscosity-driven healing mechanism remains yet to be verified and uncalibrated. The triaxial behavior has been compared only at relatively high mean stress levels, on the order of 10 MPa, due to the limited availability of experimental data. These stress magnitudes are particularly relevant as they are comparable to basal overburden pressures typically encountered beneath glaciers (Herman et al., 2015, 2021; Cook et al., 2020).”

Reviewer 2 (Q6): “Experimental data in Figures 3 and 4 are compared visually. Including quantitative error metrics (e.g., RMSE, R²) would significantly strengthen the model validation.”

Response (Q6): We appreciate the reviewer’s suggestion. As explained in the paper, our calibration approach follows a staged, physics-based calibration in which elastic, viscous, and damage parameters are identified from distinct regions of the response rather than from a single global least-squares fit. As a result, aggregate goodness-of-fit metrics such as R² or RMSE are not well aligned with the calibration procedure. Their values depend strongly on how the metric is defined and normalized, and can be biased by the nonuniform distribution of experimental points along the strain axis. For these reasons, such global scores may give a misleading impression of model accuracy. Instead, it is very common in constitutive modeling studies to assess predictive quality primarily through direct visual comparison with experimental curves, and we have chosen to follow the same method in our study.

Reviewer 2 (Q7): Table 1 presents calibrated material parameters. Please discuss any assumptions or uncertainties associated with these values and their impact on model predictions.

Action (Q7): The following passage will be added to the revised manuscript to address this point.

“The parameters in Table 1 are calibrated for isotropic, freshwater polycrystalline ice. Under this assumption the literature indicates a comparatively narrow band for elastic constants, with Sinha (1989) reporting E ≈ 8.93-9.53GPa and ν ≈ 0.308–0.311 over 0 to -50°C. The Maxwell configuration excludes the reversible rate-dependent deformation (primary creep), which is known to lead to lower Young’s modulus values in comparison to commonly reported values. For the Glen–Arrhenius set {A,n,Q}, reported values show bounded variability with n ≈ 2-4, Q ≈ 60-120 kJ⋅mol^-1 while A spans orders of magnitude across studies due to microstructure and differences in calibration windows (Glen, 1955; Barnes et al., 1971; Zeitz et al., 2020). Within this context, our calibrated {A,n,Q}lie within established intervals and are consistent with ranges reported by Durham et al. (1983) for polycrystalline ice.”

Reviewer 2 (Q8): To better illustrate the influence of the damage component, please provide comparative results showing model predictions with and without the damage formulation.

Response (Q8): The damage mechanism is essential for reproducing the post-peak brittle softening observed in the experiments. Without it, the model reduces to a simple viscoelastic Maxwell formulation, which produces only monotonic stress–strain behavior under constant strain rate. For the reviewer’s reference, Figure provided as a supplement to this response illustrates this behavior and highlights the necessity of the damage component. Because the response of such basic viscoelastic models is already well studies in the literature, we propose not to include this comparison in the revised manuscript.