Determining the fate of surface and basal cracks on ice shelves is an important yet unsolved problem in the cryosphere. In this paper, following a prior work focused on the development of the boundary element model and LEFM by the same authors, a novel dynamical system approach is employed to generalize the predicted outcomes of surface and basal cracks for a range of relevant parameters. The final result is subject to a few assumptions, such as the periodic boundary condition, and the neglect of the buoyancy restoring force related to the ice-ocean boundary condition. However, the dynamical system analysis and interpretation of the phase space are generally applicable pending improvements to the LEFM model itself. My comments primarily pertain to the clarity of the model explanations. With these minor revisions, I recommend this paper for publication.



First, as this is the first time a phase plane approach for crack propagation is presented, the authors can include some examples where clear theoretical expectations exist to facilitate the explanation of the phase plane. For example, the authors mention that their previous analytical result that \tau = 0.039 would cause calving by basal crevasse. Theoretically without surface water we would simply see zero surface crack growth and basal crack growth. I think Fig 3 a4 attempted to illustrate this, but I’d think for eta = 0.1, any surface crack depth <0.1 corresponds to dry surface crack. Thus dt < 0.1 would all correspond to arrows pointing towards the right. Could the authors explain the curved trajectories near [dt, db] = [0, 0]?

Would you have a simpler phase plane if water height is rescaled by the crack depth rather than ice thickness?



On the other hand, when the surface crack is fully filled with water, and no resistive stress tau = 0, it makes sense that surface crack always leads to full-depth penetration (figure 3e1). However, why by increasing the resistive stress to tau = 0.04 (figure 3a1) the surface crack dt < 0.15 would not propagate? This seems to have an important implication, that when the ice surface is fully filled with meltwater, as long as the resistive stress is large enough and the surface crack under certain depth, calving could be induced by basal crack rather than surface crack.



Second, in the propagation rate model described in equation 4, the crack reaches steady state when the stress intensity factor reaches the fracture toughness. As demonstrated in figure 6, generally when the stress intensity factor curve bends back down at larger crack depth, there exists two steady state crack depths. Are the K = Kc steady states at the smaller dt, db simply treated as unstable steady states  in the dynamical system approach? Can the author specify this?



Finally, please note that the left hand side of equation 4 is still dimensional; the units on the left-hand side don't seem to align with the unit on the right-hand side. What other system parameters are missing to ensure a matching unit and determine the time scale of crack propagation?



Line by line suggestions/comments:



Eqn 1: Briefly explain the resistive stress Rxx and give a reference



Line 115: Add “Poisson’s ratio” before \nu

Line 117: “written in the form” -> “written in the dimensionless form”



Line 122: As surface crack propagates the dimensionless water level generally won’t stay constant. Please add a justification or acknowledge the model limitation caused by this simplification.



Line 123: “are constant during crack propagation” -> “are assumed to be constant during crack propagation.”



Line 125: “as t increases” -> “as time t increases”



Line 134: Have the authors checked that when W further increases the result doesn’t change for this aligned surface-basal crack case?



Line 141: “maximum ensures not only ensures that cracks cannot shrink” The first “ensure” appear to be a typo



Line 142: What variable is non-differentiable against what variable? K against db, dt?



Line 144: “intensity factors is equal to the fracture toughness.” -> “intensity factors is equal to the fracture toughness (green and yellow lines in figure 2).”



Figure 7: Is the critical stress to drive calving sensitive to the resolution of your numerical model? If yes. Have the authors checked that the result had already converged with higher resolution?



Line 385: “We anticipate that incorporating the feedback between displacement and fluid pressure at the boundary will lead to additional torques generated by vertical displacements in the far field, suppressing crack growth for very large crack spacings”. The effect of vertical elastic deformation on buoyancy at the ice-ocean interface was included in Buck and Lai (2021), although their result corresponds to zero fracture toughness. The ice-ocean restoring buoyancy force can increase the critical stress \tau_{crit} for basal crack to reach the sea level.

Buck, W. Roger, and Ching‐Yao Lai. "Flexural control of basal crevasse opening under ice shelves." Geophysical Research Letters 48, no. 8 (2021): e2021GL093110.

This study is the second (or third) in a trilogy of papers by the same author team that examines propagation of crevasses in freely floating ice shelves using a boundary element model.  The contribution of this manuscript is to study the effect of surface water filled and longitudinal extension on the interaction between surface and bottom crevasses.It has long been assumed that surface and basal crevasses intersect to form rifts, but the dynamics of surface and basal crevasse propagation and interaction have rarely been studied.  Hence, this is a welcome study into an old, but important problem. 

The problem of the interaction between adjacent cracks also has a long history in the fracture mechanics literature.  This is a surprisingly challenging problem because, even under pure model I loading, the crack tip stress field from interaction results in mixed-mode loading.  As a consequence, experiments and theory indicate that en echelon cracks coalesce in a “kink” rather than in a straight intersection.  The analytic, numerical and experimental results that I am familiar with for this problem are typically done under idealized pure mode I propagation so it isn’t clear that this necessarily applies to the ice shelf problem. Nonetheless, it would be reassuring if the authors can use their model to reproduce some of the classic results of en en-echelon elastic fracture propagation—or at least touch base—with some of the literature.  I would be surprised if the interaction between surface and bottom crevasses resulted in pure mode-I behavior.  The fact that this problem has a lot of history, in my opinion, deserves a little bit more attention in the introduction.

Overall, I think this is an interesting study that merits publication.  I have a few comments.  However, I also reviewed a separate study and some of my comments repeat previous comments associated with that review as well.  The authors and editors should take care to determine if those comments need to be addressed here or are a merely a restatement of my previously expressed prejudices.

Overarching comments.

1. I mentioned this in my previous review of a different manuscript, but I find the non-dimensionalization counter-intuitive and hard to track.  The two main dimensionless numbers are tau and eta.  The parameter tau is a measure of longitudinal extension and eta is a measure of the water pressure filling crevasses.  A more natural (to me) definition of tau would define the non-dimensional  longitudinal extension stress based on the reduced gravitational acceleration (g’=(1-r)*g) or, equivalently, based on the resistive stress associated with a freely spreading ice shelf.  This would imply that a value near unity corresponds to an ice shelf spreading under its own weight and values larger or smaller would correspond to extensional stresses that are larger or smaller compared to an ice shelf spreading purely under its own weight.   I have a hard time visualizing what a tau of 0.02 means physically without resorting to using my calculator to mess with densities.  I think the eta parameter is even more difficult for me to visualize.  The situation most relevant for most ice tongues is the surface-water free case.  Previous studies have defined water depth in crevasses as a fraction of the crevasse depth, which is a bit more intuitive to visualize (brim-full vs empty).  I would encourage the authors to consider their non-dimensionalization and to connect the values as much to physical situations as possible (i.e., water-free crevasses, extension larger/smaller than the gravitationally induced spreading, etc.) to make it as easy as possible for readers to understand the underlying physical situation the authors envision.

1. One of the novelties of this study is the display of basal and surface crevasse depths in a phase plane.  I think this is an interesting way of displaying the results with a lot of potential.  This method introduces a slightly different perspective than the way we typically think of these problems.  The way we normally think of the system is how deep will a crevasse penetrate given a small “starter crack” of some pre-determined size.  The phase plane encourages us to think about pre-existing crevasses of a variety of sizes, including those that aren’t necessarily “small”.  The question that this introduces is what processes introduce large-ish crevasses that seed the initial conditions?  Is the idea that crevasses advect from a region where the stress was larger?  I see the authors come back to this on page 15.  It might be helpful to foreshadow or mention this earlier.  This is especially relevant because what we typically see is that rifts and crevasses initiate along the margins and propagate from the margins into the interior of the ice shelf.  This requires a more 3D treatment of fracture, but it seems relevant that the starting depth for basal or surface crevasses here might be related to the horizontal propagation of a crevasse or rift with some stress that includes stress concentrations associated with the horizontal fracture.

1. I think this might be addressed in one of the other manuscripts, but when the authors  introduce a crevasse into a freely floating ice shelf, the ice shelf has a flexural response that is not incorporated by  the
   
   "viscous  pre-stress".  The flexural response tends to reduce the stress concentrated ahead of crevasses.  Is this included in the boundary element model?  What effect would neglecting it have on model results?  What does the flexural stress do to the lateral boundary conditions?  I assume this is negligible for domains that are very large compared to the flexural wavelength, but the domain sizes here seem roughly comparable to the flexural wavelength or smaller.  This seems especially relevant to the interaction between crevasses.  I see this is returned to near lines 385.  I think it might be worth introducing this earlier, perhaps in the methods/model section as it seems quite important.  

1. Is it true that vertical propagation is the most optimal orientation for crevasse propagation?  If the direction of propagation is determined by the direction of maximum principal stress, are crevasses expected to kink or turn based on the direction of maximum principal stress?  A relevant physical question is what happens to crevasses that are slightly offset from each other?  It would be surprising if crevasses were exactly aligned, but what if they are mis-aligned by a small fraction of the width?  Would they never intersect?  Is it possible that the phase space is not well resolved if crevasses are allowed to kink or turn?

Minutia:

Line 18-35.  I think the more relevant comparison is between boundary element models and damage mechanics.  Damage mechanics can be used so simulate failure under a wide variety of circumstances.  Judicious choice of the damage production function allows damage mechanics to reproduce LEFM results, creep rupture or any heuristic method of simulating failure.  One of the reasons that damage mechanics is so popular is that it avoid the need to remesh that is the bane of many LEFM simulations.  Damage mechanics has been used to simulate the growth of both isolated surface and basal crevasses and arrays of crevasses.  It would be nice to a more detailed comparison between the results considered here and those previous results.

Line 33.  It is true that discrete element models do have a dependency on the packing orientation, but it has been shown that these models to converge to the continuum elastic limit under some circumstances.  One of the open questions, however, is how to specify the bond strength.  Conventional discrete element models include two fracture parameters and this allows mixed-mode failure.  Mixed-mode failure is something that can also be difficult to simulate within a linear elastic fracture mechanics framework because it requires an additional criterion to allow cracks to kink or bend.  Typically, one assumes that cracks propagate in the direction of the largest principal stress.  It seems like this study, however, assumes single mode loading.  

Line 135:  Vanishing elastic traction implies that elastic strains vanish at the domain boundaries, but least displacements are allowed, right?

Equation (4).  What are the units and numerical value of  K’(0)?  I apologize if I missed this in the manuscript.  If K’(0) is dimensionless, it is unclear how the units of the equation work out.  If it is dimensional, then we need to know the numerical value.

Line 245: Placing cracks a distance of W/4 and 3W/4 depends on the width of the domain.  What about slightly offset crevasses?  There is, in theory, two distances in the problem, right?  The distance between the crevasses and the length of the (periodic) domain.  What happens if the distance between crevasses remains the same, but the length of the domain increases?

Line 87: Punctuation?  Is the semi colon supposed to be there?

I think lines 295 are saying that you need a larger stress to propagate an array of crevasses all the way through compared to isolated crevasses.  This is consistent with previous analytic calculations by Weertman and others.