the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
A leadingorder viscoelastic model for crevasse propagation and calving in ice shelves
Abstract. We use a leadingorder viscoelastic model for crevasse evolution, in which a purely viscous model for the deformation of the domain couples with linear elastic fracture mechanics models through a viscous prestress. The fracture mechanics model conversely couples with the viscous model by inserting cracks into the domain, which viscous flow subsequently pulls apart. By contrast with prior work, we solve the fracture mechanics problem on the actual domain geometry using a boundary element method, coupled with a finite element solution of the Stokes equations describing the viscous flow. We study a periodic array of surface and basal crevasses on an ice shelf being stretched at a prescribed rate. We find that calving can either occur instantly for large enough stretching rates or sufficiently high surface water levels or through feedbacks between partial fracture propagation and subsequent viscous deformation and adjustment of the viscous prestress acting on crack faces. Our results show that purely stressbased calving laws cannot robustly describe the process of calving, since they cannot account for the gradual evolution of local crevasse and surface geometry, which can be understood at the large scale as being more akin to the evolution of a damage variable. Future work will need to coarsegrain the type of process model we describe here in order to make it applicable to ice sheet simulations.
 Preprint
(6216 KB)  Metadata XML
 BibTeX
 EndNote
Status: final response (author comments only)

RC1: 'Comment on egusphere2023807', Jeremy Bassis, 08 Sep 2023
Review of TCD: Interaction between surface and basal crevasses
This manuscript is the second (or third?) paper in a trilogy of papers by the same authors focused on a boundary element model that is used to simulate the vertical propagation of surface and bottom crevasses in ice shelves and ice tongues. This manuscript focuses on a viscoelastic implementation of the model in which the (elastic) boundary element model is coupled with a full Stokes “flowlline” model. In this approach, the viscous stress is used to drive fracture propagation assuming that fractures propagate on a time scale that is rapid compared to that of viscous flow. This separation of scales allows the large scale viscous deformation of the ice to alter the stress near preexisting crevasses providing a means to simulate the viscous widening and vertical propagation of crevasses.
The authors use this model to map out conditions where a combination of surface water filling surface crevasses and longitudinal stretching eventually permit crevasses to penetrate the entire ice thickness. This ductile deformation driven models result in a version of what is often called a “necking” instability evident in ductile failure of metals and other materials that exhibit creep on a time scale comparable to fracture.
There is an impressive amount of work that has gone into developing the boundary element model and then coupling it into a viscous Stokes flow model, which, as the authors clearly show, is not a trivial exercise and involves both nuance and significant numerical hurdles. Overall, I think this is an interesting manuscript that makes novel contributions to the field. The amount of details and technical exposition is, however, challenging. I had to read it several times and I’m still not entirely sure I understand some of the details. Below I give several suggestions that might ease the exposition for a broad readership. My overarching comment is that, as much as possible, it would be helpful for readers to be provided with a physical picture associated with model results. I also think some of the choices for nondimensionalization are opaque making it harder to interpret the results. For example, the authors provide thorough descriptions of conditions when surface crevasses penetrate the entire ice thickness, but it isn’t clear *why* these conditions lead to full thickness penetration. Is that once the crevasse penetrates deep enough, there is a critical height of water that allows the crevasse to rapidly penetrate the rest of the ice shelf?
Before I delve into some of my more detailed comments, I do want to point out that the mechanism for calving that the authors envision is one in which “starter cracks” in the middle of an ice shelf or ice tongue gradually propagate through the ice thickness is not the way calving occurs for most ice shelves and ice tongues. What we see most often is that fractures initiate along the margins and pinning points where stress is concentrated and then propagate horizontally across the ice shelf. This is not to say that the mechanism that the authors simulate doesn’t happen (or isn’t related), but it might be useful to signpost the fact that the twodimensional model may not be adequate for many interesting situations that aren't well represented by flow lines.
Major comments:
 Nondimensional numbers. There are two main nondimensional numbers in the manuscript tau_star and eta_star and these two numbers play a crucial role in the exposition. I had to constantly flip backandforth to remind myself what the numerical values correspond to. As far as I understand it, both of these are defined with respect to some mean ice thickness. This is perhaps my own issue, but I am especially interested in certain limits that correspond to physical situations. Examining tau_star first, my suspicion is that the relevant quantity is whether the deviator stress (or resistive stress) is larger or smaller than it would be for a freelyfloating shelf/tongue that is spreading solely under its own weight. Intuitively, it would be a little bit surprising if fractures could penetrate the entire thickness of a freely floating ice tongue spreading solely under its own weight without *some* additional forcing. At the very least, it seems as though this should be marginally stable since we see freely floating ice exist stably in many environments. I think it would be more intuitive to define a nondimensional number that measures the ratio of the imposed longitudinal extension stresses relative to that of a freely spreading ice shelf such that a value of unity (or something!) corresponds to a freely spreading ice shelf and we can more easily assess whether it is being pulled apart with more or less vigor. I apologize if I am misinterpreting or imposing my own framework on the authors. One aspect of the nondimensionalization used by the authors that isn’t intuitive to me is that because of the way tau_star is defined, it seems like assuming a reference density of ice that is 800 kg/m^3 (some firn/marine ice) vs 900 kg/m^2 (mostly ice) will result in distinct values of tau_star even for a freely floating ice tongue that is solely spreading under its own weight. Is that intentional? Or am I hung up on own mental rigidity. Similarly, it is hard to ferret out the physical significance of the eta_star parameter. What I would really like to know for the eta_star parameter is when are surface crevasses water free? Does this correspond to eta*=0 (was this not considered)? Are all of the simulations done for surface crevasses with some amount of water (provided they propagate deep enough)? I think eta*=0.05 corresponds to a water table at the mean surface height (but this, again depends on the ratio of the density of ice to water, right?). I think in the past the water table was defined relative to crevasse depth so that some fracture of the crevasse was water filled. I may just be very inflexible on this. Another possibility would be to annotate the graphs more to show where on graphs crevasses are water free and where crevasses are “brimful” and also to show where the ice is being pulled apart with a stress larger than or smaller than the stress associated with a freely floating ice shelf.
 Physical picture/cartoon of the situation considered with water table and other variables defined. I think it would be helpful to readers to provide more details in the cartoon of the setup envisioned. I had a hard time understanding the boundary conditions imposed along the horizontal edges of the domain: are stresses or velocities imposed along the boundaries? This is a bit embarrassing, but I’m still not sure I understand what boundary conditions where imposed, which is my next point.
 Horizontal boundary conditions. I’m still a little bit confused about the horizontal boundary conditions imposed on the full Stokes model. In fracture mechanics we tend to distinguish between displacement and stress controlled loading. In displacement controlled loading, you take a sample and pull it apart with a constant velocity or stepwise displacement. In stress controlled loading a stress is applied to the edges of the sample. Stress controlled loading tends to result in catastrophic failure whereas displacement controlled loading tends to result in quasistable fracture propagation. I’m still not entirely sure which situation the authors are considering here. I think (?) the authors have imposed a horizontal velocity to the edges of the domain, but is a bit confusing because they talk about applying stretching rates. This would seem to be a velocity or displacement controlled boundary condition. For ice shelves and ice tongues, stress controlled seems the more appropriate boundary condition and I think that a more physically obvious boundary condition would be to impose a strain rate (velocity gradient) at the boundary. (I apologize if this is exactly what the authors did!). A stress controlled boundary condition would imply (I think) that the elastic strains at the boundary vanish rather than the displacements. Boundary conditions can be really important and I would encourage the authors introduce a separate section to describe the boundary condition, whether they are displacement or stress controlled. I’m not insisting that the authors change the boundary conditions. Just explain what they are. There is text to explain this, it was just scattered in a couple of different places making it harder to follow and the notation wasn’t always clear.
 Reference to past and current work. This is a pretty minor point, but I found I had to constantly refer back to the authors prior and current work and this made it a little bit more difficult to follow the manuscript. One example of this is that the the authors refer to the nondimensionalization done in a previous manuscripts to justify a result *before* the nondimensionalization is introduced in this manuscript. I think it would be helpful to streamline the references to previous/existing work in the introduction to tell summarize what has been done in those manuscripts and how this manuscript differs and then keep the crossreferencing to a minimum in the methods/results sections. The discussion would be an excellent time to come back to any crossreferencing between results.
 Model domain size. This a bit of numerical issues, but typically for fracture studies we take a domain length that is 10 to 20 times the ice thickness to avoid contamination from edge effects. This is a bit of rule of thumb. I think the authors used a much smaller domain in this case (perhaps out of numerical necessity). I would like to encourage the authors to demonstrate that their results are insensitive to domain size. This doesn’t need to be included in the manuscript, but a single statement that their results were insensitive to model domain size would be helpful in convincing readers that results aren’t sensitive to arbitrary domain choices.
 Finite element model resolution. Similarly, I would encourage the authors to conduct tests to demonstrate that their results are insensitive to the mesh size of the finite element Stokes solver. My understanding is that the authors did this for the boundary element size, but not the finite element sizes. I apologize if I misunderstood. Again, this doesn’t necessarily need to take significant space in the manuscript, but explain any sensitivity to mesh size. Similarly, ’m not sure that the authors told us the number of elements or range of element sizes used for the FEM model domain. This would be helpful information.
 The interplay between viscous and elastic stresses presented here is interesting. Given that one of the requirements for vertical crevasse propagation is that the viscous prestress must be tensile, is it possible to determine the vertical crevasse depth solely using the viscous model by assessing if the stress at the tip of crevasses is tensile? A relevant question is whether the elastic part of the problem is even needed to determine if (when?) crevasses will penetrate the entire thickness. It seems entirely plausible that the elastic part of the problem is entirely determined by the viscous prestress and it may be possible to ignore the elasticity and viscoelasticity completely. I think that is a result that the wider community would find very interesting. I was left wondering if the boundary element model was even necessary and if a simpler full Stokes model could suffice when paired with some remeshing and a Jintegral or some other estimate of energy fluxes.
 Section 2.1 This section provides a detailed description of the contact conditions imposed and is fairly elaborate. These are fairly standard conditions, but as the authors found, these conditions can be numerically challenging to impose. Later in the manuscript we learn that the contact conditions are *not* actually imposed in the Stokes solver. I will be honest that I was little bit miffed to have gone through and followed all of the details only to find out that these conditions were eventually dismissed because of the numerical difficulty imposing them. Fair enough. But there isn’t any detail about how these equations are treated numerically and what the authors end up solving for the viscous part of the problem is a standard Stokes flow problem. I think the manuscript can be simplified by ditching large portions of this section or relegating it to supplementary material. I would be a lot more interested in reading about it in a manuscript that also provided a plausible numerical method/discretization that was implemented.
 Setting the critical stress intensity of zero (near Line 255). The argument around setting the critical stress intensity factor to zero isn’t obvious to me. Or at least it requires more details. To start, the nondimensionalization used here isn’t introduced until later in the manuscript making it a harder to follow the argument. As far as I understand it, the authors are introducing a scale for the stresses that is proportional to ice thickness and this results in a scale for the stress intensity factor that is proportional to ice thickness to the 3/2. But the rho*g*H term provides the magnitude of the hydrostatic pressure, which is compressive. In my experience, the argument for setting the critical stress intensity to zero relies on an argument centered on starter crack lengths. For instance, the critical stress necessary to propagate a crack roughly scales like the critical stress intensity factor divided by the square root of the length of the starter crack size. For large starter crack sizes, this quantity becomes small. I would have thought that the authors require an additional length scale here associated with starter crack or boundary element sizes. Are starter cracks assumed to be the same size irrespective of ice shelf thickness? Or are they assumed to scale with ice thickness? Either way, small cracks are less likely to propagate and large cracks are more likely to propagate. The assumption here seems to be that large cracks are always present and thus the boundary element size is sufficiently large? I think it would be helpful to sketch out the argument in more detail. It actually seems like by setting the critical stress intensity of zero, the situation is analogous to a preexisting vertical fracture that already exists with zero separation between the surface of the cracks. The question the authors ask is when will that preexisting crack open.
Small Details.
Line 15: It is true that fulldepth propagation is a key requirement for calving. However, even when these fractures are present, fractures need to also propagate horizontally in most cases.
Lines 3035: There is a big chunk of literature missing here from methods associated with damage mechanics (and the associated phase field methods). Damage mechanics can be used to simulate the interaction between multiple crevasses in a continuum model. I see that the authors come back to this in the discussion. I think it might be useful to add damage mechanics and phase field methods into the introduction as another method to contrast with the boundary element. The reality is that damage mechanics has become extremely popular in the fracture community because it avoids or sidesteps many of the issues that they authors have discovered with remeshing and mesh tangling.
Lines 30: It has been demonstrated that for sufficiently small particles, discrete element models do converge to elastic models (without fractures) and can reproduce key features of wingcracks and other continuum fracture observations given an appropriate parameterization of the failure strength. There is nonuniqueness inherent in discrete element models description of fracture. But this nonuniqueness is related to failure under mixedmode loading, which is something that also requires additional assumptions to incorporate into LEFM. To me, the bigger issue with discrete element models it that they are fully dynamic and typically limited to time scales of minutes to hours. These models cannot be easily integrated into larger scale and longer time scale simulations.
Title: Leading order in what? Is it leading order in the Maxwell viscoelastic time scale? Or the fracture propagation time scale? I don’t think this is ever stated. One potential subtlety is that the time scale of fracture propagation has to be short enough to be fully elastic, but also slow enough to avoid dynamic effects.
Line 85, just as a historical point the field of dislocation mechanics has long recognized the nonuniqueness associated with displacement and the compatibility equations in the presence of dislocations.
Line 151: missing word?? “Occurs over a time scale??””
Line 216, equation 21. I think I’m confused about the boundary conditions used. I thought that a traction boundary condition was imposed at the edges of the domain, which seems to imply that the gradient in displacements normal to the boundary must vanish and not necessarily the displacement. I think this comes down to whether the authors are considering a displacement controlled or stress controlled boundary condition. The natural choice for this boundary conditions seems to be a traction/stress boundary condition, but I think (?) the authors are using a velocity boundary condition.
Section 2.4. Why restrict crack growth to the vertical direction. These cracks will generally experience mixedmode loading and, for sufficiently slow propagation, it is commonly assumed that cracks will propagate in the direction of maximum tensile stress. This should result in a slight curvature of the crack surface. Is this negligible enough to be neglected? Or this just done for convenience?
Line 235: What is used to initialize cracks here and how is it different from the other submitted study? I got lost in this sentence. I think this might just be a matter of rewriting the clauses in the sentence or separating into two sentences.
Lines 302. In Berg and Bassis (2020), we actually found that the sea spring resulted in a time step size dependent solution. This was generally small for most situations, but where it become apparent was when we remeshed or the domain abruptly changed size. I suspect that those effects are more muted here because of the periodic boundary conditions and the relatively flat domain. The sea spring had the most artifacts when we included the horizontal ocean pressure condition at sloped calving calving front. I think doing convergence studies associated with the FEM element size should demonstrate that this isn’t an issue here.
Line 305. Contact conditions. There is a lot of text on contact conditions and imposing the right boundary conditions. But, as far as I can tell, these are largely abandoned and instead a standard Stokes flow problem is solved. This results in mesh tangling, which is a long standing problem in large deformation problems and one of the reasons that damage mechanics has become so much more popular than LEFM for many problems. Can the authors assess if there is any impact on the solutions? The mesh tangling can have subtle effects on the solution because of the unphysical nature of even small node crossing.
Line 355. The s is defined here based on Archimedes principle using the mean ice thickness, right?
Line 359. When portions of the lower boundary are above sea level, this means that crevasse propagate above the water line? Is there any other situations when this occurs? Can this be explained in a simpler fashion?
Line 375. Note that although the nondimensional equations are independent of length scale, there is still a scale effect associated with starter crack size. The assumption that critical stress intensity factor remains small implies critical crack lengths also scale.
Line 397. What is r again? Can you refer back to equations where these are introduced because I found myself flipping and back forth frequently trying to find definitions.
Line 410. This might be better described using an illustration. I think this is consistent with the way most define the normalized crack length, but I got lost in the technical description. Maybe this should be moved to the methods section.
Figure 2. I suspect that these dimensionless numbers will be meaningless for most readers without a strong feeling. Can they be put into a more intuitive scaling or provide physical numbers for a reference ice thickness? I’m having a hard time understanding why the edges of the boundary are not straight given that I thought that a constant velocity was imposed along the edges? Panels a and b didn’t initially look like the same shape to me, but I see that different scales are used. Can the authors use the same horizontal scale for both sub panels to help readers translate between the two. Also, can the authors use a different color in panel to illustrate which of these is being shown in panel b?
Another minor question that I is whether the episodic propagation nature is inherent to the problem or has that been built in. The authors have assumed that fracture propagation is rapid compared to the viscous time scale and slow and steady propagation would be inconsistent with the approximation. Would omitted terms in the viscoelastic relationship result in more continuous propagation? At the vary least, it is encouraging to see that the results are selfconsistent with the approximation.
Is it possible to benchmark the model against any of the fully viscoelastic models that have been used for bottom or surface crevasses? This might be interesting and build confidence in the approach used.
Line 441. The overshoot is pretty typically of the elastic problem where the stress concentration allows fractures to penetrate slightly deeper than it would otherwise. Actually, you get the exact same overshoot if you use a linear viscous rheology instead of the powerlaw viscous.
How is the longitudinal extension boundary condition determined? I think I get the impression that Rxx decreases as the ice thickness decreases? Is that true?
Line 465: Is there a physical explanation that can be provided? My suspicion is that this is all surface water pressure mediated and then once the water pressure reaches a critical value, the surface crevasse propagates the entire depth, but I think readers would like a little bit more physical description of what is happening this situation to balance the technical description.
Figure 7. There are two blue lines. Which is the water level?
Line 480: *boundary* element size, right? This isn’t done with respect to finite element size, right?
Figure 9: Can authors use the same axis for all of the subpanels? Also, can you use a different color than red and green for the lines as this combination can be difficult for color blind folks to distinguish.
What does it physically mean when the surface water table is at or beneath sea level? Also, what condition would correspond to water free surface crevasses (the most reasonable condition for most ice tongues in the Antarctic)?
Line 640650. I think the interpretation of damage mechanics is not entirely accurate. Damage mechanics is a *method* that can be used to simulate failure in a continuum model that alleviates the need to remesh necessary in fracture mechanics approaches. The damage production function can be chosen to reproduce results from LEFM, micromechanics or be chosen based on a heuristic calibrations.
Line 725. A big difference between the work of Crawford et al., and this study is that the DEM in Crawford et al., is initialized with no memory of its previous state. Hence, the model has no memory of basal/surface crevasses that had propagated over part of the domain when reinitialized.
Citation: https://doi.org/10.5194/egusphere2023807RC1 
AC1: 'Reply on RC1', Maryam Zarrinderakht, 29 Feb 2024
The comment was uploaded in the form of a supplement: https://egusphere.copernicus.org/preprints/2023/egusphere2023807/egusphere2023807AC1supplement.pdf

RC2: 'Comment on egusphere2023807', Anonymous Referee #2, 15 Oct 2023

AC2: 'Reply on RC2', Maryam Zarrinderakht, 29 Feb 2024
The comment was uploaded in the form of a supplement: https://egusphere.copernicus.org/preprints/2023/egusphere2023807/egusphere2023807AC2supplement.pdf

AC2: 'Reply on RC2', Maryam Zarrinderakht, 29 Feb 2024
Viewed
HTML  XML  Total  BibTeX  EndNote  

321  135  22  478  18  18 
 HTML: 321
 PDF: 135
 XML: 22
 Total: 478
 BibTeX: 18
 EndNote: 18
Viewed (geographical distribution)
Country  #  Views  % 

Total:  0 
HTML:  0 
PDF:  0 
XML:  0 
 1