Evaluation of the criticality of cracks in ice shelves using finite element simulations

The ongoing disintegration of large ice shelf parts in Antarctica raise the need for a better understanding of the physical processes that trigger critical crack growth in ice shelves. Finite elements in combination with configurational forces facilitate the analysis of single surface fractures in ice under various boundary conditions and material parameters. The principles of linear elastic fracture mechanics are applied to show the strong influence of different depth dependent functions for the density and the Young’s modulus on the stress intensity factor KI at the crack tip. Ice, for this purpose, is treated as an elastically compressible solid and the consequences of this choice in comparison to the predominant incompressible approaches are discussed. The computed stress intensity factors KI for dry and water filled cracks are compared to critical values KIc from measurements that can be found in literature.


Introduction
Eight of twelve ice shelves in the Antarctic Peninsula have retreated or disintegrated in the past decades (Cook and Vaughan, 2010;Braun et al., 2009).The processes that lead to break-up events at ice shelves are all linked to fracturing of weakened polycrystalline ice.Causes for the weakening are surface cracks forming due to bending stresses at the surface and crevasses, formed by tensile stresses, originating at shear margins or along the ice front.These cracks might initially be stable, but additional loads may let them become critical.Our analysis of cracks, based on well established fracture mechanical concepts, is focused on simplified scenarios which we derived from the break-up events that happened at the Wilkins Ice Shelf in 2008/2009(Braun et al., 2009)).
Fracture mechanical concepts investigate the criticality of cracks by determining the stress intensity factor K I at the crack tip and comparing it with critical values K Ic , obtained by experiments.Rist et al. (2002) performed three-pointbending and short-rod fracture tests on samples taken from a core of the Ronne Ice Shelf that contained meteoric as well as marine ice.The measured K Ic show a strong dependence on density or porosity, respectively.Further experimental values for K Ic were assembled by Schulson and Duval (2009).
The vertical propagation of crevasses has been investigated for more than 50 yr.Nye (1955) argued that crevasses propagate vertically to the point where the normal stress σ xx changes sign.The first one to assume elastic material behaviour for the analysis of vertical crevasses in glaciers was Weertman (1973).He used dislocation distribution functions to evaluate the characteristics of dry and water filled single crevasses.Smith (1976) was the first who applied methods of linear elastic fracture mechanics for the evaluation of stress intensity factors of dry and water filled surface cracks in ice shelves.He simplified the crack geometry and boundary conditions (BCs) to facilitate the use of tabulated values gained from semianalytical methods by Tada et al. (1973) and Sih (1973).The method of Smith (1976) was adapted and extended by Van der Veen (1997), who discussed the importance of depth dependent density profiles and Rist et al. (2002) who additionally analysed depth dependent tensile stresses.The model by Rist et al. (2002) has been repeatedly used for the analysis of bottom and surface crevasses, e.g. by Nath and Vaughan (2003) and Luckman et al. (2011), Published by Copernicus Publications on behalf of the European Geosciences Union.
and will be used here as a benchmark for our finite element simulations.
A different approach was introduced by Weiss (2004) who argued that critical crack growth can not explain slow crevasse propagation.He therefore analysed subcritical crack growth for very simplified geometries, boundary conditions and material parameters.
Stress intensity factors for the analysis of horizontal crack propagation have been used by Hammann and Sandhäger (2003) and Hulbe et al. (2010).Hammann and Sandhäger (2003) used a numerical model of ice shelf dynamics to compute stresses within the ice shelf that are then used as input for the evaluation of stress intensity factors.Hulbe et al. (2010) treated the ice as linear elastic solid and used the boundary element method for the analysis of crack criticality and growth.
The qualification of a linear elastic model to describe the ice rheology for fracture mechanical purposes is discussed controversially within the glaciological community.Therefore, different approaches exist to investigate the flow induced strain rate and/or stress fields in ice shelves and glaciers and to formulate yield stress or strain rate criteria for the nucleation and propagation of cracks.Larour et al. (2004) and Albrecht and Levermann (2012) analyse the evolution of crevasse fields or single crevasses based on ice dynamical simulations.Fracture is then understood as a softening parameter or an enhancement factor (Humbert et al., 2009) leading to higher velocities within the ice shelf or glacier.These methods require less knowledge about the material parameters of the ice and are a valid approximation for fields of closely spaced crevasses, where the stress concentration at the crack tip is reduced.However, these approaches do not provide a physical examination of fracture processes which is the purpose of our analysis.
This study investigates the effect of different BCs, loads, density profiles, Young's moduli and Poisson's ratios on the criticality of a single surface crevasse.For this purpose, the adequate model is presented in Sect. 2. The plane strain model with the corresponding equations is introduced in Sect.2.1.Sections 2.2 and 2.3 explain the finite element discretisation and the resulting discrete configurational forces.Valid BCs and a satisfying numerical model are identified in Sects.2.4 and 2.5 followed by the validation using the well known model of Rist et al. (2002) in Sect.2.6.The results of the numerical simulations are presented and discussed in Sect. 3 and 4. In Sect.3, we study the influence of depth varying material properties on the stress intensity factor of dry cracks.For this purpose, we first show the effect of different BCs in Sect.3.1.The results for the different varied material parameters are presented in Sect.3.2 to 3.4.Two scenarios for water filled cracks are discussed in Sect. 4 followed be the summary in Sect. 5.

Model
In order to analyse the dependence of the criticality at the crack tip on different depth dependent material parameters and various BCs, finite element (FE) simulations are used.This provides the basis for further simulations with advanced 2-D and 3-D geometries.The advantage of the FE method in comparison to semi-analytical methods or other numerical methods like finite differences lies primarily in its flexibility and in its selectable accuracy via mesh refinement.

Basic equations
Creep dominates the long-term behaviour of ice.Therefore, in most ice dynamical analyses, ice is modelled as viscous fluid.On the short time scale, ice can be considered as linear elastic solid.A combination of both rheologies results in a visco-elastic model.A first and very simplified viscoelastic representation for a fluid-like visco-elastic material is the Maxwell element consisting of a serial coupling of a spring and a dashpot.The equilibrium condition yields that the stress in both components is the same.If a Maxwell element is subjected to a stress, the elastic component instantaneously responds by a time independent elastic strain while the viscous strain is a function of time.Therefore, on large time scales, the viscous answer predominates, while on the short time scale, the elastic response is important.The present analysis concerns the influence of different parameters on brittle fracture, which is a short time process.For this reason, only the elastic response of ice, e.g.fractures due to elastic strain, is considered.
Our analysis of the criticality of certain crack scenarios is based on the evaluation of the crack driving force at the tip of a sharp Griffith crack (Lawn and Wilshaw, 1977), where the maximum distance between the crack faces is much smaller than the crack depth.
Thus, we consider a static, linear elastic plane strain model of an edge crack as depicted in Fig. 1a.The equation for the solution of the boundary value problem for a linear elastic solid in equilibrium is with the volume forces f and the Cauchy stress σ .The Cauchy stress is obtained by the constitutive equation where C is the stiffness tensor and ε is the symmetric part of the displacement gradient For the isotropic case, the stiffness tensor C depends on only two independent constants, the Young's modulus E and the Poisson's ratio ν.Gross and Seelig (2011).
The solution of Eq. ( 1) with Eqs. ( 2) and (3) in conjunction with proper BCs yields the displacements u and the stress field σ , which are required for the subsequent calculation of the crack driving force and the stress intensity factor.For the evaluation of the crack driving force, configurational forces are used.Configurational forces can be interpreted as a negative energy release rate.The benefit of configurational forces in comparison with the J-integral is that these forces can be evaluated at every nodal point within the FE mesh.They provide a measure for the integrity of the material structure and allow the consideration of inclusions, cracks or inhomogeneous changes of material properties.The evaluation of the configurational forces follows the method presented in Müller et al. (2002).Here, the authors introduce the Eshelby stress tensor as a function of the strain energy, U =1 2 ε : (Cε) 1 , the identity tensor I, the transposed displacement gradient (∇u) T and the Cauchy stress tensor σ .With the definition of a configurational volume force the configurational balance equation can be written as The configurational volume force g vol considers the contribution of the physical volume force f to the configurational force balance while the contributions of inhomogeneous material properties as, e.g.due to a spatial dependence of C = C(x, y), are given by Insertion of Eqs. ( 4), ( 7) and (8) in Eq. ( 6) leads to a form of the configurational balance equation that holds within the boundaries of a continuous material.

Finite element discretisation
The solution of Eq. ( 1) is obtained by using the FE method.For this purpose, Eq. ( 1) is transformed into the weak form by multiplication with a test function η and integration over the body B, Integration by parts with application of the Gauss' theorem leads to with the applied traction vector t * = σ •n on the stress boundaries ∂B t .FE discretisation of the test function and the displacement vector u yields where N I and N J are the standard shape functions for the applied elements.Insertion of Eq. ( 12) in Eq. ( 11) results in The integral on the left hand side represents the vector of internal forces (F int I ) and the integrals on the right hand side the vector of the external (F ext I ) and volume forces (F vol I ), respectively.The solution of the residual equation provides the nodal displacements u J .The discretisation of Eq. ( 9) follows analogous.Using the function φ with φ = 0 on ∂B as test function, the weak form of Eq. ( 9) takes the form The FE discretisation of the test functions leads to

C. Plate et al.: Evaluation of the criticality of cracks in ice shelves
As residual equation, Eq. ( 16) can be written as

Interpretation of discrete configurational forces
With the application of FE, the continuous Eq. ( 9) is transformed into a discrete form, where G I (u J ) = 0 is a measure for the discontinuity at every node of the FE mesh.The measure of the discontinuity at the crack tip (index "ct") can be interpreted as the crack driving force G = G ct (u J ).
From the predominant vertical component of the crack driving force at the crack tip, G y = G • e y , the stress intensity factor K I is calculated using the interrelation Further information on configurational forces, stress intensity factors and their relation can be found in Müller et al. (2002), Gross and Seelig (2011), Maugin (1993), Steinmann and Maugin (2010), Gurtin (1999) and Kienzler and Herrmann (2000).

Geometry, load and boundary conditions
In reality, the ice shelf is subjected to gravity, as well as to different boundary tractions (tension, pressure and shear).The modelled ice shelf consists of a vertical cut through a finite section of an "infinite" ice shelf, which is replaced by a sufficiently long rectangular domain (l = 2000 m, b = 250 m) under plane strain conditions, see Fig. 1a.The model only considers the gravity induced pressure and horizontal strain due to the ice flow.In-plane-shear is neglected as fractures tend to align perpendicular to the first principal stress direction, which is shear free.In a sufficient distance from the grounding line, the horizontal velocities and displacements in an ice shelf are depth-independent (Greve and Blatter, 2009).This constraint can not be fulfilled using traction BCs at the vertical boundaries, as those would allow tilting of the vertical boundaries.Therefore, unless stated differently, the vertical boundaries are loaded with prescribed vertically constant displacements u.Using Hooke's law for an uncracked homogeneous body under uniaxial tension, σ = E ε, with E = E ice in the case of plane stress and E = E ice /(1 − ν2 ) for plane strain, the magnitude of the boundary displacement u on one side of the model ice shelf is related to the horizontal stress σ at the ice shelf surface by The stress field at the ice shelf surface is evaluated from the flow velocity in the ice shelf using Glen's flow law (Glen, 1958).The bottom boundary is loaded by the water pressure at the respective depth of the undeformed body.Further traction BCs are eventually applied on the crack faces to consider water filled crevasses.

Numerical model
The stresses and displacements in the rectangular domain are determined by solving Eqs. ( 1)-( 3), using the commercial FE program COMSOL 2 .The crack driving force and the resulting stress intensity factors are evaluated in postprocessing routines in MATLAB3 .As the stress intensity and therefore K I at the crack tip highly depends on the element size in the vicinity of the crack tip, as well as on the distribution of elements within the geometry, the appropriate mesh has to be chosen carefully.Figure 1c shows the difference between the simulated K I at the crack tip and a semianalytical solution (Gross and Seelig, 2011, p. 79) for an edge crack under linear loading for different discretisations.The black curve with an element edge length of 0.0125 m at the crack tip shows satisfying results for reasonable computation time.Computation time is saved by cutting the model geometry along the crack and using half the geometry with symmetry BCs. Figure 1b illustrates the discretisation assigned to the black curve.All simulations are conducted using 6-node triangular elements with quadratic shape functions.The resulting total number of elements for the meshed geometry is about 8400 for a crack tip located in the domain centre.

Benchmark
The FE model is validated using the geometry and material parameters presented in Rist et al. (2002).The authors introduce a semianalytical approach for the evaluation of stress intensity factors for cracks in ice shelves, taking a 422 m thick part of the Ronne Ice Shelf as an example.They assume that the stress intensity factor at the crack tip depends on the total stress acting on the flaw.Using a usual power law for the ice flow (Glen, 1958) and balance equations, following Weertman (1957), a distribution of the normal stress σ xx is derived as a function of the vertical position z in the ice shelf: Here ρ(z) is the depth-dependent density of the ice, which is parametrized by  based on measurements of ice cores, ρ sw = 1028 kg m −3 is the density of salt water and B is the temperature dependent and therefore depth-dependent rate factor.
Further information on the applied temperature and density profile can be found in Rist et al. (2002).Rist et al. (2002) use the weight function method presented in Bueckner (1970) to evaluate the stress intensity factor K I based on σ xx .The normal stress σ xx and the resulting stress intensity factor K I are shown in Fig. 2a.The diagram shows that the critical stress intensity factor K I c , which ranges between (100-400) kPa √ m (Rist et al., 2002), is reached at larger depth than the depth where in the uncracked body the normal stress σ xx changes sign.The assumption that cracks will only propagate to a depth where the stresses change sign, as presented by Nye (1955), turns out to be an underestimation.Figure 2b shows a comparison of the results of Rist et al. (2002) and two FE simulations: (I) a crack which has solely been loaded on its faces by tractions given by Eq. ( 20) and (II) the same geometry loaded by gravity, a vertically constant displacement BC, equivalent to the non-cryostatic part of the stress in Eq. ( 20) and the water pressure as stress BC at the ice shelf bottom.As the stress function of Eq. ( 20) assumes incompressibility, the FE simulation in (II) is conducted with a Poisson's ratio of ν → 0.5.The results are in a good agreement, taking numerical inaccuracies in the semianalytical results and the FE simulation into account.

Dry cracks
Dry cracks are simulated to validate the model and to analyse the influence of different material parameters and loading scenarios on the stress intensity factor at the crack tip.The studies first concentrate on the evaluation of the appropriate type of BCs for further simulations.Then the influence of the applied load, Poisson's ratios, density profiles and Young's moduli is analysed.

Study A: dependence on type of boundary conditions
Rist et al. (2002), Van der Veen (1997) and Weertman (1973) use stress BCs at the vertical boundaries of the ice shelf.Therefore, the constant or depth varying tensile normal stress is superposed by the cryostatic pressure of the ice.This approach requires a hydrostatic stress state within the ice shelf that is only valid if ice can be understood as incompressible (ν = 0.5).This is a good approximation for the long term behaviour of ice.In contrast, a fracture event in a brittle medium occurs on a rather short time scale.The measurements in Rist et al. (2002) indicate brittle material properties for ice.Therefore, it seems more reasonable to take material properties from short time measurements into account, as can be found in Greve and Blatter (2009) and Schulson and Duval (2009), where ν ranges between 0.2-0.4.Poisson's ratios of ν = 0.5 have been used in elastic analyses of ice by Rist et al. (1999), Hulbe et al. (2010) and Konovalov (2011).The influence of the Poisson's ratio on the choice of BCs will be presented in the following chapter, while the influence on the K I is analysed in Sect.3.2.

Boundary conditions and scale effects
Figure 3a shows two uncracked example geometries of a homogeneous isotropic body, i.e. ρ = const., with different horizontal BCs (upper geometry: constrained, lower geometry: unconstrained).Additionally, the vertical displacement at the basal boundary is constrained in both configurations, u z (z = 0) = 0.The geometries are solely loaded by gravity.
Figure 3b shows the resulting the normal stress σ xx for different Poisson's ratios and BCs.The stress component σ zz that can also be referred to as the ice overburden pressure is identical for all simulated Poisson's ratios and BCs.It equals the horizontal stress for the constrained boundaries and ν = 0.499 (dot-and-dashed line).The horizontal stress component σ xx is affected by both the changes in the Poisson's ratio and the BCs.A Poisson's ratio of ν = 0.499, approximating an incompressible material, leads to a cryostatic stress state for the horizontally constrained body, meaning that for every material point in the body, the normal stress components are equal and shear stresses vanish.The value ν = 0.499 is a good approximation of the incompressible case (ν = 0.5), which numerically can not be treated with the applied constitutive law as it yields singularities in the stiffness matrix.For ν = 0.3, the body is compressible.This leads, for the horizontally constrained body, to stresses σ xx which are less than half of the stress component σ zz .For a horizontally unconstrained body, the horizontal stress component σ xx is identical zero.These results show that the assumption of incompressibility overestimates the crack closing pressure due to the weight of the ice by approximately a factor two.
In a next step, we apply different BCs on a cracked geometry to analyse the effect of the type of boundary condition on the stress concentration at the crack tip as well as possible dependencies on the length of the model domain.To ensure comparability to semi-analytical results, the simulations are performed without additional gravity or bottom pressure.Figure 4a shows that simulations with pure stress boundary conditions on a domain of 2 km length and 250 m depth (green drawn through line) reproduce the semi-analytical results (black circles), which assume infinitely long domains (Gross and Seelig, 2011).Further simulations with different depth/length ratios and stress BCs indicated that a ratio larger than 1 2 is sufficient to yield a crack tip field that is independent of the disturbance introduced by the boundaries.Different results can be expected for constrained boundaries (dashed and dotted lines) where the dependence on the domain length affects the stress intensity at the crack tip in two ways.To identify the factors that induce a length dependence, the cracked domain is on the one hand loaded by displacement BCs u and on the other hand loaded by stress BC σ with the additional constraint ∂u ∂z = 0, which prevents tilting and bending of the vertical boundaries.The relation between the applied u and σ is evaluated via Eq.( 19) using the uncracked geometry.The displacements u are kept constant for all crack depths.The dashed and the dotted lines in Fig. 4a show the resulting stress intensity factors K I (d) considering different domain length l for displacement and stress BCs, respectively.For large crack depths d, a length dependent difference between the constrained (dotted line) and the unconstrained stress BC represented by the drawn through line becomes obvious.The difference is caused by the crack closing bending moment induced by constrained boundaries.For l reaching infinity, the influence of the constraints on the deflection in proximity of the crack becomes negligibly small and the K I converge to the solution of the unconstrained problem with stress BCs.
The difference between the dashed and the dotted lines for geometries of equal length l results in the compliance effect of deep cracks loaded by constant displacement BCs.If we imagine the ice shelf as a spring, the compliance due to the crack growth can be understood as a reduction of the spring stiffness.A reduction in stiffness yields a reduction in the resulting stresses and therefore decreasing stress intensity factors.The compliance effect becomes less obvious for values of l approaching infinity and lim l→∞ (K I (d)) yields the solution of the unconstrained problem with stress BCs.
The results show that constrained boundaries imply a dependence on the length of the cracked domain.For infinitely long domains the influence on the type of BC vanishes and The Cryosphere, 6, 973-984, 2012  the solution of the pure stress BC is reached.The crack closing bending moment as well as the compliance due to the application of displacement BCs are reasonable qualities of the chosen model, considering an ice shelf with a given vertically constant velocity field.
Cracks in ice shelves are small-scale defects -crack depth and crack opening are small in comparison to the lateral extend of the ice shelf.As we are interested in the local effect of the crack we have to choose a model size that is large in comparison to the crack length but small in comparison to the characteristic length of the ice shelf.In addition, the stress field in the shelf can vary over length scales of less than 10 km for which reason it is unrealistic to model cracked domains with lengths of more than 10 km.To this end, we chose a finite domain length of l = 2 km for all further simulations, taking the length dependence due to displacement BCs into account.
Figure 4b illustrates the deviation in the resulting K I for different BCs and realistic loading due to additional volume forces (ρ = const.)and water pressure at the bottom of the shelf.
The applied BCs consist of pure tension (a), the equivalent displacement boundary condition given by Eq. ( 19) (b) and the superposition of tension and horizontal pressure (c).As for ν = 0.3, the horizontal pressure is not equivalent to the ice overburden pressure (see Fig. 3b), it has to be evaluated from the horizontal reaction forces of the horizontally constrained body.For a small load of 100 kPa and resulting shallow cracks, the difference between the BCs (b) and (c) is marginal.Case (a), with pure tension represents a totally different loading case.Even though the body is loaded by volume forces, they do not influence the horizontal stresses.The stress intensity factors for higher loads are on the other hand more sensitive to the choice of the BC as (d) and (e) demonstrate.Figure 4c shows a qualitative plot of the horizontal normal stress σ xx and the deformed shape (exaggerated presentation, scaled by a factor 100) for stress BCs and equivalent displacement BCs (no application of gravity or water pressure at the bottom boundary).As displacement BCs prevent the ice from bending, a bending moment that works against the crack opening is induced and the stress intensity at the crack tip, especially for deeper cracks, is smaller.

Load
Next, different loads using displacement BCs are applied on the pre-cracked ice shelf.Here, the density is chosen as constant over depth (ρ = 910 kg m −3 ) and ν = 0.3.For the evaluation of plausible load cases, the principal stresses in a part of the Wilkins Ice Shelf are calculated using the velocity field of Braun et al. (2009).The resulting first principal stress ranges from about −400 kPa to 400 kPa.The purpose of this work is to analyse the criticality of fractures due to tension in the open ice shelf.Therefore, only positive stresses to a maximum of 300 kPa seem relevant.Figure 5a shows the stress intensity factors for displacement BCs equivalent to 0, 100, 200 and 300 kPa.It is obvious that for zero boundary displacement only the ice overburden pressure is acting on the crack, leading to negative values for K I , which can be interpreted as crack closure.Non-zero load leads to positive K I , varying with the crack depth.For very small cracks, the stress intensity factors are low as there is enough unbroken area to absorb the load.For deeper initial cracks, the stress intensity factors grow until a maximum value is reached.Then K I decreases as the influence of the ice overburden pressure starts to compensate the tensile stress arising from the BCs.The dashed red lines represent the range for measured values of the critical stress intensity factor K I c , see Rist et al. (2002).Values of K I beyond the critical value K Ic are interpreted as crack growth, while values lower than K I c imply stable cracks.It appears that none of the simulated load cases leads to penetration of initial cracks through the entire depth, as K I becomes negative before the bottom of the ice shelf www.the-cryosphere.net/6/973/2012/The Cryosphere, 6, 973-984, 2012 is reached.This result is in good agreement with previous findings by Rist et al. (2002).

Study B: influence of Poisson's ratio
Unlike the results for various constant Young's moduli (Sect.3.4), there is a big difference in the stress intensity factors for different vertically constant Poisson's ratios as can be seen in Fig. 5b for a constant density (ρ = 910 kg m −3 ) and a displacement BC equivalent to 100 kP.This is due to the coupling between the stresses in x-and z-direction which is influenced by ν.In other words: for a Poisson's ratio of ν = 0, the normal stress σ xx does not experience stress contributions of σ zz induced by the body loads.For ν → 0.5, the stress state is hydrostatic for constrained vertical boundaries.
There is hardly any reliable data on the depth dependence of the Poisson's ratio.Therefore, only constant distributions were simulated to obtain a general understanding of the relation between K I and ν.

Study C: influence of different density profiles
Previous studies by Rist et al. (2002), Van der Veen (1997), Scambos et al. (2000) and Scambos et al. (2009) motivate the necessity to take depth-dependent density profiles into account.The density of the ice is estimated from the densification model of Herron and Langway (1980).There are different mechanism of densification, which contribute to the depth ranges under consideration here.The densification in the upper regime, driven by grain growth and sintering, down to a density of 550 kg m −3 , depends only on temperature.Below that the grains form bonds, allowing recrystallisation and deformation to become dominant and the density depends on temperature and accumulation rate.As we are lacking in-situ measurements of the mean annual surface temperature and the accumulation rate of the Wilkins Ice Shelf, we choose upper and lower bounds for both variables.The mean annual surface temperature of the Wilkins Ice Shelf was proposed by Morris and Vaughan (2003) to be −8 • C.However, the surface of the Wilkins Ice Shelf melts every summer and most likely refreezes in winter, so that the latent heat increases the firn temperatures.Swithinbank (1988) reports a temperature measured in a drill hole at 5.5 m of −2.5 • C. We thus assume −2.5 • C as a maximum mean annual temperature, leading to an almost isothermal ice shelf.
The accumulation rate was given by Vaughan et al. (1993) to be 0.5 m a −1 WE, based on a stake measurement over a short time period.Thus, we choose for this study additionally 1.0 m a −1 WE, to have an upper estimate.Exponential fits of the estimated density profiles are presented in Fig. 6a.Additionally, a constant density profile is presented for comparison reasons.Figure 6b shows the corresponding stress intensity factors for a displacement load equivalent to 100 kPa and ν = 0.3.It can be seen that small differences in the density profiles lead to only marginal differences in the K I .Nevertheless, we observe that density profiles with larger values at lower depth (270 K, a s = 0.5 m a −1 WE) lead to lower K I than more moderate profiles (264 K, a s = 1 m a −1 WE).This effect becomes more obvious when applying the constant profile, which leads to considerably lower K I .We conclude that higher densities lead to a higher ice overburden pressure at the crack tip and therefore to less tensile stresses which are known to be responsible for larger K I .Rist et al. (1999) and Scambos et al. (2009) motivate a density dependent examination of the critical stress intensity factor K Ic ranging from K Ic = 50 kPa √ m, for low density firn to K Ic = 150 kPa √ m for meteoric ice.This, in our model, results in a change in the critical crack depth of less than 10 m.The variance due to different density profiles or elastic material parameters is larger, hence the depth dependence of K Ic will not be considered.

Study D: influence of Young's modulus variation
The stress field and the consequential K I resulting from stress boundary value problems in linear elastic solid mechanics are invariant to the choice of the material parameters E and ν.In contrast, for displacement BCs as chosen in the present studies, the value of E and ν has to be considered.Different constant Young's moduli do not change the outcome of the stress intensity factors, considering displacement BCs which are calculated by Eq. ( 19).Different results can be expected from depth dependent Young's moduli.Figure 7a shows the different dependencies of the Young's modulus E used in the simulations in which a constant, a linear and an exponential shape are considered.Rist et al. (2002) motivate a density related and therefore exponential dependency of the Young's modulus on the depth, which gives reason for a coupled analysis of the influence of both parameters.However, in order to compare the simulation results with former analyses applying depth dependent density profiles and to separate effects and mechanisms, we decided to look at E and ρ independently.To this end, the Poisson's ratio and the density profile are kept constant with ν = 0.3 and ρ = 910 kg m −3 , respectively.Equation ( 18) shows that the Young's modulus is included in the relation between stress intensity factors and the calculated configurational forces.The simulations presented here use the Young's modulus at the crack tip for the evaluation of K I .As stated before, the Young's modulus contributes to the displacement BC, Eq. ( 19).Therefore, an equivalent method for the evaluation of u has to be found for varying E. It seems reasonable to look at two cases.At first, a displacement BC, for which the surface stress remains the same as in the previous analysis, is simulated.The results are presented in Fig. 7b.It shows that for a small surface load of 100 kPa (solid line) only the exponential function for the Young's modulus leads to considerably different results.For a higher load of 300 kPa (dashed line), the crack reaches into deeper zones where the Young's modulus, and therefore the tensile stress is larger.This results in a strong variation of the stress intensity factors for all three profiles.Secondly, the displacement BC are adjusted by the average of the stress component σ xx over the depth in an uncracked geometry without volume forces.The average stress is equal to the previously chosen surface stress.Figure 7c shows that this choice of BCs leads to smaller differences in the resulting K I for both, a load of 100 kPa (solid line) and 300 kPa (dashed line).

Dry cracks: conclusion
Our results achieved by FE simulations on dry cracks support the general findings from previous studies by Nye (1955), Weertman (1973), Smith (1976), Van der Veen (1997) and Rist et al. (2002): dry surface cracks under reasonable tensile loading won't reach the base of the ice shelf.The importance to consider depth dependent density profiles was affirmed (Van der Veen, 1997;Rist et al., 2002).Our results showed an increase of the crack depth by ≈ 50 % in comparison to the constant profile.Our analyses differ from previous findings in the applied BCs and the associated choice of material parameters.Here the influence of the chosen Poisson's ratio of ν = 0.3 has to be emphasised.The higher stress concentration at the crack tip could only be marginally reduced by the crack stabilizing effect of displacement BCs.This, in conclusion, leads to larger crack depth under comparable geometry and loading conditions.Beyond that, for the first time, the influence of depth varying Young's moduli and different vertically constant Poisson's ratios for cracks in ice was investigated.It shows that the influence of different Young's modulus profiles is significant and strongly depends on the choice of BCs.The Poisson's ratio proves to be the most important parameter for the analysis of crack criticality.Values of ν ranging between 0.2 and 0.5 show results varying from no crack growth to cracks only stabilizing at 70 m depth for equal loading.

Wet cracks
The results above show that even though changing material parameters and loadings have a strong influence on the stress intensity factors at the crack tip, hardly any of the simulated configurations would cause a break through of an initial crevasse.It seems reasonable to consider some additional loading as for example due to water of different origin inside the crevasses.Previous studies by Weertman (1973), Van der Veen (1997), Rist et al. (2002) and Scambos et al. (2000) show that water pressure on the crack faces can lead to a crevasse breaking the ice shelf.However, as these studies have been conducted using a Poisson's ratios of ν = 0.5, the influence of the ice overburden pressure σ zz as a crack closing factor has been overestimated.This leads to the conclusion that for reasonable tensile loadings, a deep crack had to be almost entirely filled to break through.The following simulations show that for different material properties, even less water leads to critical situations.The simulations are conducted using the depth-dependent density profile of Fig. 6 (T s = 264 K, a s = 0.5 m a −1 WE), a constant Young's modulus and a Poisson's ratio of ν = 0.3 unless stated differently.

Study A: surface melt water
Figure 8a visualises the model for rising melt water within a crack.For the depth dependent density profile as shown in Fig. 6 (T s = 264 K, a s = 0.5 m a −1 WE) pore closure is reached at about 50 m depth.We assume that melt water can rise inside a crack as long as the water-level h w is below the pore closure.The simulation is conducted for three different displacement BCs equivalent to 100, 200 and 300 kPa.For these simulations, the crack depth is kept constant at the depth for which the unfilled crack was about to reach crack closure, which is indicated by negative values for K I .This leads to a crack depth of 66 m for 100 kPa, 122 m for 200 kPa and 172 m for 300 kPa. Figure 8b shows the resulting stress intensity factors as a function of the water-level for the applied loading.We find that the critical stress intensity factors are reached for only 10 to 20 m high water columns inside the crack.The crack corresponding to a 100 kPa load can be filled up to h w = 16 m before the water reaches permeable ice.As for this water level the stress intensity factors are critical, it is interesting to evaluate how deep this crack would penetrate.Figure 8c visualises the stress intensity factor for various crack depths, considering an equivalent to 100 kPa tensile loading and 16 m water filling, starting at a crack depth of 66 m.K I increases for few more meters before the influence of the ice overburden pressure starts to dominate the stress state at the crack tip and the stress intensity factor decreases.Crack closure is reached before the crack can break through.Sufficient additional melt water supply at deeper crack depths will, however, lead to crevasse penetration.

Study B: brine infiltration
Cracks in a closer proximity to the calving front can be exposed to brine infiltration through porous firn.These cracks are therefore always filled up to sea level as visualised in Fig. 8d.Nevertheless, the pressure of brine inside the ice shelf, as well as on the crack faces, for z between the sea level and pore closure depth compensates and does not increase the stress intensity at the crack tip.The load on the crack faces therefore rises linearly from the constant load of p c = (h w − h PC )ρ sw g at pore closure depth to the maximum load of p c = h w ρ sw g at the crack tip. Figure 8e shows the resulting stress intensity factors starting at d = 0 m.The stress intensity factors for crack depths less than d = 50 m are equivalent to those for the unfilled cracks with exponentially fitted density profiles.Cracks below pore closure are exposed to water pressure and therefore show increasing stress intensity factors.By trend, all these cracks would break through.This leads to the question: How much water is required for a crevasse to break through?Table 1 shows the water level required to reach critical stress intensity factors for 249 m deep cracks.
We find that increasing the load within the applicable range of 300 kPa only leads to a decrease in the required water level for penetration of less than 30 m for both simulated magnitudes of Poisson's ratio.On the other hand, a decreasing Poisson's ratio from ν = 0.5 to ν = 0.3 leads to a decrease in the required water level of ≈ 90 m representing a decrease by ≈ 50 % for all simulated loads.

Wet cracks: conclusion
The studies on wet cracks could confirm previous findings on water filled cracks: water pressure on the crack faces profoundly increases the criticality of cracks and can lead to crack penetration where unfilled crevasses are stable.However, it has to be mentioned that continuous water supply is needed as crevasses in an ice shelf of 250 m thickness need to be filled up to 91 m-119 m, depending on the load, to reach penetration.The study repeatedly showed that the choice of the Poisson's ratio is more crucial to the evaluation of crack criticality than the applied load, a finding that has not been discussed in previous studies.

Summary
Finite elements in combination with configurational forces proved to be a comfortable numerical tool for the evaluation of crack criticality under different setups.The choice of a valid rheological model for cracks in ice in combination with appropriate BCs and material parameters turned out to be crucial for a physical understanding the of fracture mechanism that lead to disintegration of ice shelves.We showed that especially the Poisson's ratio and the associated compressible or incompressible treatment of ice during fracture plays an important role that has hardly been discussed so far.Schulson and Duval (2009) and Rist et al. (2002) showed a density and therefore depth dependence of the Young's modulus.The use of finite elements allowed us to evaluate the influence of the depth dependence on the criticality of cracks.It showed that for the prescribed surface stress, the depth dependent and therefore exponential function for the Young's modulus doubled the critical crack depth.However, depth dependent Young's moduli raised the questions which choice of BCs is appropriate for cracks in floating ice shelves and how the interaction between ice dynamics and linear elastic fracture mechanics should be formulated.Despite the different choice of BCs and material parameters, the studies showed that for the applied loading, dry surface cracks will not penetrate an ice shelf.Also the general findings of previous studies on water filled surface cracks could be affirmed.Nevertheless our choice of the appropriate Poisson's ratio lead to ≈ 50 % less water required for crevasse penetration.
Fig. 1.(a) Model geometry for dry cracks.(b) Discretization in entire geometry with focus at crack tip.(c) Difference between numerical simulation and analytical result for different mesh sizes.( x = element edge length at crack tip).
Fig. 2. (a) Resulting horizontal stress and stress intensity factors.(b) Comparison of the numerical model to results of Rist et al. (2002).
Fig. 3. (a) Uncracked model for the evaluation of proper BCs, with constrained (upper geometry) and unconstrained (lower geometry) vertical boundaries.(b) Normal stress σ xx for the uncracked model.

Fig. 8 .
Fig. 8. (a) Stress intensity factors for different BC.(b) Qualitative contour plot of the normal stress σxx for displacement and equivalent stress BC.
Fig. 5. (a) Stress intensity factors for different loads.(b) Stress intensity factors resulting from different constant Poisson's ratios.
Fig. 6.(a) Applied density profiles estimated from the densification model of Herron and Langway (1980) and an additional constant profile with ρ = 917 kg m −3 .(b) Stress intensity factors resulting from different density profiles.
Fig. 7. (a) Simulated depth-dependent Young's modulus functions.(b) Stress intensity factors resulting from different Young's modulus profiles with surface stresses equal to previous simulations (drawn through line, u = 100 kPa surface stress, dashed line, u = 300 kPa surface stress).(c) Stress intensity factors resulting from different Young's modulus profiles with resulting depth integrated tensile stresses equal to previous simulations (drawn through line, u = 100 kPa depth integrated tensile stress, dashed line, u = 300 kPa depth integrated tensile stress).

Table 1 .
Critical water level (WL) at 249 m crack depth in a 250 m thick ice shelf for different loads and Poisson's ratios.