This article provides a detailed review of how resistive stress is calculated for crevasse depth evaluation in various articles in the literature and show that they result in different crevasse depth predictions. As a part of the study, the authors considered both idealized cases and real cases of Antarctic ice shelves, especially those with shear margins. I really appreciated the systematic study in this article, given that my research has focused on exploring the crevasse models in the literature. I believe the article is of interest to the ice sheet modeling community and I will be happy to see it published in this journal. The article is well written and organized (and amazingly I could not find a single typo in text). However, I have several comments listed below, most requiring minor changes or clarifications but I have a couple major comments at the end related to stress evaluation using planar remote sensing data. While I have several articles of mine referenced in my review, I leave it up to the author’s discretion to cite or not cite them in their paper.

-- Ravindra Duddu

Detailed Comments:

Line 51 - I suggest your say zero stress theory instead of Nye crevasse formulation. Originally, Nye (1957) did not include the effect of water pressure, but was later introduced by Jezek (1984).

Line 61 – If I remember correctly, Enderlin and Bartholomaus (2020) used observed surface strain rates to calculate stress in grounded glaciers. If the basal boundary is not free slip, then the stress variation with depth in grounded glaciers is not linear. Please clarify this point that the resistive stress is not a constant if boundary condition is not free slip.

Line 81 – The term “stress differences” could be misunderstood. I think it is better to say, “viscous flow is driven by deviatoric stress, which is the component of Cauchy stress that does not cause volume change during deformation; whereas brittle failure is governed by the Cauchy stress.”

Line 82 – Consider replacing the terms “biaxial” with “equi-biaxial” and “triaxial” with “equi-triaxial” in the paper, as you are referring to the cases when the stress magnitudes are equal. The terms biaxial and triaxial do not imply the stresses are the same in various directions.

Line 90 – I do not agree with the statement “in areas with simple stress states, such as on an ice shelf or near an ice cliff …” Due to free surface effects near an ice cliff of a grounded glacier or floating ice shelf the stress state is not simple. Only in the far field do the stress becomes independent of the horizontal direction and one can use force balance calculation to determine resistive stress. Not sure what I am missing. Please clarify this point.

Line 94 and 98 – Change in ice rigidity due to firn layer are important, as you noted. Sorry for the self-promotion, but I would encourage you to read our recent papers (Gao et al., 2023; Clayton et al., 2024), which examine the influence of firn layer on crevasse propagation in glaciers and ice shelves. In Gao et al. (2023) we show that in ice shelves considering depth-varying density due to the firn layer changes the buoyancy depth and leads to deeper penetration. In Clayton et al., (2024) we consider both depth-varying Young’s modulus and density and derive analytical solutions and conduct analytical LEFM studies. We found that the inclusion of depth-dependent density influences the resistive stress and can thwart or promote deeper crevasse propagation depending on the glacier and ocean water heights, which is more nuanced than the description by van der Veen (1998a).

Gao, Y., Ghosh, G., Jiménez, S., & Duddu, R. (2023). A finite-element-based cohesive zone model of water-filled surface crevasse propagation in floating ice tongues. Computing in Science & Engineering, vol. 25, no. 3, pp. 8-16, May-June 2023, doi: 10.1109/MCSE.2023.3315661

Clayton, T., Duddu, R., Hageman, T., & Martínez-Pañeda, E. (2024). The influence of firn layer material properties on surface crevasse propagation in glaciers and ice shelves. The Cryosphere, 18(12), 5573-5593.https://doi.org/10.5194/tc-18-5573-2024, 2024.

Line 124 – You write height above buoyancy twice. Instead say: “where rho_pw is the density of the proglacial water (lake or ocean) and H_ab is the height above buoyancy defined as” and get rid of repetition.

Line 130 – Perhaps you should mention that Eq. (4) is only valid in 1D.

Line 136 – Better to say “… but considers the stress singularity at the crevasse tip …” Stress concentration is bounded and occurs around holes and U-shaped notches. At the sharp crack tip in LEFM theoretically there is a stress singularity.

Line 175 – In Eq. (8), (9) and others, whenever the subscripts are not indices but rather descriptors like “eff” or “eef, planar” you must use \text{} or \mathrm{}. Only indices that are symbols taking numerical values are italicized.

Line 197 – You can mention that tau_1 and tau_2 in Eq. (13) are invariants with coordinate transformation, whereas tau_flow dir is not in Eq. (12).

Line 266 to 268 – Calculating surface and basal depths using different rigidities based on temperature differences is a bit ad hoc. To see how reasonable it is, a full Stokes FEM simulation could be conducted to obtain the stress field and then the depth where tensile stress becomes zero can be taken as the crevasse depth.

Line 303 to 304 – In this discussion about inverting for damage, it can be noted that damage in the form of crevasses introduces anisotropy. In Huth et al. (2021) we show that when this anisotropy is considered we get more realistic rift propagation. A comment can be added to state that inversion for isotropic viscosity has the limitation attention that it applies the effect of crevasse damage equally in all directions.

Huth, A., Duddu, R., & Smith, B. (2021). A generalized interpolation material point method for shallow ice shelves. 2: Anisotropic nonlocal damage mechanics and rift propagation. Journal of Advances in Modeling Earth Systems, 13(8), e2020MS002292. https://doi.org/10.1029/2020MS002292

Table 3 – You are calculating these at a point on the glacier by assuming that the strain rate is uniform. Is that right or did I misunderstand, please clarify.

Line 333 – Please clarify the phrase “… but the ratio of depths between calculation will be identical” Does mean ratios taken column wise or row wise.

Line 351 – Would be good to include the function plotted in Figure 4 in the Appendix, for the sake of reproducibility. Or maybe you can share the code used to generate these plots.

Line 388 - Perhaps, my only major concern is that the evaluation of stress from observed strain rates in the region of a crack is physically meaningless. This has not been particularly mentioned in the paper. To elaborate, the observed strain rate in an area where there is a crack will be large and the stress evaluated using the Glen’s law will be large, so you may rightly predict a full depth crevasse. However, the true Cauchy stress there will be zero because the ice rigidity becomes zero in an open crevasse. Therefore, it is important to point out that while one can use this approach, the evaluated stress is a trial stress (borrowing this term from plasticity) and not the true stress.

Line 400 to 410 – The argument comparing Larsen B and Stancomb-Wills shelves are a bit difficult to understand. Also, why did you not include the cross-section thickness and crevasse penetration depth plots in Figure 8, just like in Figure 7. This would perhaps make it easy to follow the differences.

Figure 9 could also be clearer if the crevasse depth plot like in Figure 7E was included. Also, by looking at the REMA, how are you able to tell whether crevasses are full depth or not. Please clarify this for the general reader.

Line 443 – The statement “The assumption that the stress calculation that gives the best modeled velocity …” is fine but going to back to my previous comment this is just a trial stress whereas the true stress in a fully crevassed region is zero as the rigidity goes to zero, whereas the strain rate will be large.  

In Figure 10A, I would recommend plotting the root mean squared of the velocity misfit rather than the average. The average would not be accurate measure. Also, if the nodal velocity misfit is the least with inversion, they why do we need to use calculation F using observed strain rates. Can we not just invert for damage and obtain damage and make estimates of crevasse depth.

Line 462 to 465 – Please clarify what you mean by excess velocity. I am also confused by the comment “This may indicate that the damage in the spreading flow region …”  Velocity is not a measure of strain rate or damage, instead the symmetric gradient of velocity could be used.

Figure 11 – Explain why in the 0 - 400 m/yr range the observed velocity is greater than the modeled velocity and why this happens in both cases A and C.

Section 5.2 – As I am reading, I feel like there is a distinction between studies using planar remote sensing data and ice sheet models that was not clearly identified. While you are right about the usage of resistive stress in the formulas with planar remote sensing data, with modeling studies one can evaluate the full stress and identify where it becomes zero in the domain and determine the crevasse depth directly based on zero stress theory. This is what is done I believe in Todd et al. (2018), where they use a full Stokes model that calculates velocity and pressure from mass and momentum balance. In Clayton et al. (2022), we use the momentum balance to calculate the elastic stress in a Maxwell viscoelastic solid and then determine the zero-stress based crevasse depth directly from the stress distribution.

Clayton, T., Duddu, R., Siegert, M., & Martinez-Paneda, E. (2022). A stress-based poro-damage phase field model for hydrofracturing of creeping glaciers and ice shelves. Engineering Fracture Mechanics, 272, 108693.https://doi.org/10.1016/j.engfracmech.2022.108693

Line 534 – The calculation F is recommended for use by the authors, which is reasonable if dealing with planar remote sensing data. With SSA models once can use resistive stress and the crevasse depth formula based on the approach of Sun et al. (2017) or simply once can calculate the depth varying stress using the pressure formula (Huth et al., 2023). Further, principal stress can be obtained from the eigenvalues of the 3 x 3 full stress matrix.

Huth, A., Duddu, R., Smith, B., & Sergienko, O. (2023). Simulating the processes controlling ice-shelf rift paths using damage mechanics. Journal of Glaciology, 69(278):1915-1928. https://doi.org/10.1017/jog.2023.71

Line 541 – Perhaps, I am repeating this statement. I believe the stress calculations are not valid in Mottram and Benn (2009) or Enderlin and Bartholomaus (2020) as they study grounded glaciers, unless they are free slip at the base. For glaciers frozen to the bed the stress is not linear. In Jimenez and Duddu (2018) we used a cubic function to fit to the stress profile. I would encourage the authors to study grounded glaciers and the effect of boundary conditions and how this changes the Nye depth. This could be a quick study that could be added to this paper.

In the acknowledgements, only funding for the DEMs generation is mentioned, but it is not clear how the authors were funded to conduct this study.