Viscous and viscoelastic stress states at the calving front of Antarctic ice shelves

Abstract

Calving mechanisms are still poorly understood and stress states in the vicinity of ice-shelf fronts are insufficiently known for the development of physically motivated calving laws that match observations. A calving model requires the knowledge of maximum tensile stresses. These stresses depend on different simulation approaches and material models. Therefore, this study compares results of a two-dimensional (2-D) continuum approach using finite elements with results of a one-dimensional (1-D) beam model elaborated in Reeh (1968). A purely viscous model, as well as a viscoelastic Maxwell model, is applied for the 2-D case. The maximum tensile stress usually appears at the top surface of an ice shelf. Its location and magnitude are predominantly influenced by the thickness of the ice shelf and the height of the freeboard, the traction-free part at the ice front. More precisely, doubling the thickness leads to twice the stress maximum, while doubling the freeboard, based on changes of the ice density, results in an increase of the stress maximum by 61%. Poisson's ratio controls the evolution of the maximum stress with time. The viscosity and Young's modulus define the characteristic time of the Maxwell model and thus the time to reach the maximum principal stress.

INTRODUCTION

The calving of icebergs along ice-shelf fronts is an ongoing topic of research that is still not well understood. Consequently, physically motivated calving laws are established neither for calving of small icebergs with sizes smaller than 1 km², nor for giant icebergs. Potential calving laws might use stress criteria for the determination of the position where an iceberg breaks off. Thus, it is of major importance to have a proper understanding of stress states in the vicinity of the calving front. The knowledge of the temporal evolution of the stress state is essential to assess the time between two calving events. This time is needed to establish models for calving processes. A remarkable study in this field was presented by Reeh (1968), who transferred common methods of elastic beam theory to a linear viscous material and utilized this to determine the stress states at the surface of ice shelves. However, Reeh's approach uses major simplifications. First of all, the 1-D theory can only calculate the normal stress in the flow direction and the vertical shear stress. Therefore, the normal stress in the vertical direction has to be approximated independent of the beam theory by a linearly decreasing stress based on the hydrostatic pressure and surface conditions. Furthermore, the spreading in flow direction, as well as the stress-free part of the boundary condition acting on the ice front, is not realizable by the assumptions of the beam theory. Additionally, Reeh (1968) modeled the ice rheology as an incompressible and purely linear viscous material, neglecting elastic deformations. A similar approach was pursued by Fastook and Schmidt (1982), who computed the stress at the top surface by a 2-D finite element simulation using an incompressible linear viscous material. However, the simulation time was restricted by the numerical capabilities at their time and hence not sufficient to detect time-dependent changes in the stress distribution. Iken (1977) and Fastook and Schmidt (1982) argued that the most likely position for the next calving event in grounded tidewater glaciers and floating ice shelves corresponds to the location of the maximum tensile stress.

Some of the ice shelves in Antarctica exhibit a homogeneous, crevasse-free surface structure in the vicinity of the calving front. Calving events also occur at these ice shelves, although no surface crevasses and basal crevasses exist. Those are the types of ice shelves for which we aim to investigate the stress situation. Other situations lead to the formation of medium-size or giant tabular icebergs, where cracks turn into rifts, which propagate episodically and detach an iceberg after years to decades. In order to understand the different types of calving, we start with a very basic situation and analyze calving due to bending stresses as proposed by Reeh (1968) to obtain a general knowledge of the stress situation in the vicinity of the calving front.

Observations of calving are sparse and often only large calving events on Antarctic ice shelves are monitored. Therefore, no time series of calving events of the bending type are available. However, Wesche and others (2013) provided a classification of the surface types of calving fronts in Antarctica. They found that 7.4% of the ice-shelf fronts have no surface structures as for example, crevasses. For these ice shelves, calving presumably happens due to bending stresses based on the boundary disturbances or the formation of giant tabular icebergs. Larour and others (2004) discussed calving events near Hemmen Ice Rise on the Ronne Ice Shelf, Antarctica. Some of these events could also be explained by bending due to the boundary disturbance at the ice front. This bending leads to longitudinal surface stresses and results in the fracture of ice and the calving of icebergs.

At sufficiently long timescales, ice is commonly assumed to behave like a viscous material. This is based on a considerable number of investigations comprising laboratory tests as well as field studies on ice shelves, ice sheets and glaciers. Nevertheless, recent observations indicate an elastic response of ice on shorter timescales. For example, Gudmundsson (2011) showed that the phase shift between the load and reaction of ice streams due to tidal forces cannot be represented by a purely viscous material. Therefore, a more complex viscoelastic material is needed to satisfy the short-term elastic as well as the long-term viscous ice behavior.

The investigation of stress states in ice shelves is a prerequisite for obtaining physically based calving laws and computing calving rates. The aim of this work is to study the influence of different material models on the maximum tensile stress at the top surface. Therefore, the geometry introduced by Reeh (1968) is used. In a first step, the comparison of the stresses after 1 d reveals the importance of the choice of boundary conditions along the vertical cliff face. Then, the temporal evolution of the stresses resulting from the 2-D finite element model with linear viscous material behavior is compared with the results derived by Reeh (1968). The viscous study is concluded by an analysis of the surface strains. In a next step, the time-dependent finite element model is extended to a viscoelastic material. Also a linear elastic simulation is considered to complete the study. The surface stresses of the three rheological models are compared for incompressible (ν ≈ 0.5) as well as compressible materials to point out important differences. The viscoelastic model is used to analyze the influence of the ice-shelf thickness and the material parameters (Young's modulus, Poisson's ratio, viscosity and ice density) on the position and the magnitude of the maximum surface stress. A comprehensive discussion reveals the consequences using the different modeling approaches presented in this study. All finite element models considered in this investigation discard the past history of the ice, which can also affect the formation of cracks. Hence, the results are applicable for ice shelves with rather homogeneous surfaces (e.g. the Ekstroem Ice Shelf studied by Lohse, 2012). However, the following results illustrate strong dependencies of the stress states on quantities influencing the ice-shelf geometry. In the case of uneven surfaces (see Humbert and others (2015) concerning a viscoelastic model of the Jelbart Ice Shelf), it is important to investigate the differences between the application of a measured geometry instead of an idealized geometry.

ICE-SHELF MODELS

Reeh (1968) analyzed the stress state and the vertical deformation of an ice shelf near the ice front. He suggested a viscous material law assuming that stress variations appear very slowly in the analyzed case and thus, elastic responses of the ice are negligibly small. Therefore, the well-known elastic beam model is adapted to the viscous case. Hooke's law implies the stress–strain relationship for the 1-D linear elastic case and reads σ = Eε with the Young's modulus E. The corresponding viscous equation $\sigma = 4\eta \dot \varepsilon $ , with the viscosity η, leads to the stress–strain rate relation for a 1-D linear viscous material. It is derived according to considerations similar to the linear elastic case. The differential equation of an infinitely wide viscous beam on an elastic foundation is given by

(1) $$4\eta I{\dot w^{IV}} + {\rho _{{\rm sw}}}gw = 0,$$

with the moment of inertia per unit width I = H ³/12, the ice thickness H, the sea water density ρ _sw and the deflection w. The superposed dot indicates the first derivative with respect to time, while ( · ) ^IV is the fourth differentiation with respect to space. In order to transfer the results to different glaciers, Reeh (1968) considered all variables dimensionless, which results in the dimensionless differential equation for a viscous beam, see Eqn (17) of Reeh (1968). The viscous beam equation is solved for the vertical deflection w for an idealized rectangular ice shelf, as shown in Figure 1, and suitable boundary conditions (see Eqns (18), (19), (21) and (22) of Reeh, 1968). The stress resultants, normal force N, shear force V and bending moment M (see Fig. 1b), are computed from the deflection w. The normal stress σ _xx and the shear stress σ _xz (Reeh, Eqns (26) and (28)) directly result according to well-known formulas from beam theory. The normal stress σ _zz (Reeh, Eqn (27)) has to be approximated independent of beam theory by the application of a linear decreasing stress from zero at the top surface to the negative water pressure at the bottom surface. Further details on the implementation of the viscous beam theory and the resulting deflections and stresses can be found in Reeh's original work.

Fig. 1. (a) Vertical cross section of the idealized geometry according to Reeh (1968). (b) Stress resultants.

In this work, we are mainly interested in the stresses near the ice front. The stress distribution in this part of the ice shelf is predominantly influenced by the boundary conditions at the front, which cannot be fully captured by the beam model. To detect differences, the conclusions of the viscous beam theory are compared with results using a 2-D finite element model. To establish a direct comparison, the same geometrical configuration of an idealized ice shelf as presented in Reeh (1968) is analyzed, see Figure 1a. The displacements and hence the strain and stress states are computed by solving the quasi-static momentum balance

(2) $${\rm div}{\bi \sigma} + {\bi f} = {\bf 0},$$

where σ denotes the Cauchy stress tensor and f represents volume forces, such as gravity. In order to describe the processes in an infinitely wide ice shelf it is sufficient to consider a 2-D problem in the xz-plane under the assumption of plane strain conditions in y-direction. Thus, the non-trivial equations of the momentum balance reduce to

(3) $$\eqalign{& \displaystyle{{\partial {\sigma _{xx}}} \over {\partial x}} + \displaystyle{{\partial {\sigma _{xz}}} \over {\partial z}} = 0, \cr & \displaystyle{{\partial {\sigma _{xz}}} \over {\partial x}} + \displaystyle{{\partial {\sigma _{zz}}} \over {\partial z}} - {\rho _{{\rm ice}}}\;g = 0}, $$

with the volume force of the ice weight ${f_z} = - {\rho _{{\rm ice}}}g$ , the ice density ρ _ice and the acceleration g due to gravity. Assuming small strains, the kinematic relation is linearized and the three-dimensional strain tensor ε is given by

(4) $${\bi \varepsilon} = \displaystyle{1 \over 2}(\nabla {\bf u} + {\nabla ^T}{\bf u}),$$

with the vector of the unknown displacements u = (u _x , u _y , u _z ) ^T . A suitable material law completes the system of differential equations. The stress and strain tensors are split into volumetric and deviatoric parts

(5) $${\bi \sigma} = \displaystyle{1 \over 3}{\rm tr}({\bi \sigma} ){\bi I} + {{\bi \sigma}^D}\quad {\rm and}\quad {\bi \varepsilon} = \displaystyle{1 \over 3}{\rm tr}({\bi \varepsilon} ){\bi I} + {{\bi \varepsilon}^D},$$

where I is the second order identity tensor and ( · ) ^D denotes the deviator. The trace of an arbitrary second order tensor A is given by tr( A ) = A _xx  + A _yy  + A _zz . In the considered case, the trace of the strain tensor tr( ε ) only consists of the sum of ε _xx and ε _zz , as ε _yy is zero due to the plane strain assumption. In order to implement an incompressible, viscous 2-D continuum, the fluid pressure is often solved for as an additional unknown instead of using a constitutive relation. In this work an approximation using an elastic isometric stress instead of the thermodynamic pressure is introduced, see Altenbach (2012). To achieve incompressibility a Poisson's ratio of ν ≈ 0.5 is used. This allows for a comparison with the incompressible, viscous beam. The corresponding material law reads

(6) $${\bi \sigma} = \left( {\lambda + \displaystyle{2 \over 3}\mu} \right){\rm tr}({\bi \varepsilon} ){\bi I} + {{\bi \sigma}^D}\quad {\rm with}\quad {{\bi \sigma}^D} = 2\eta {\bi {\dot \varepsilon}}^D.$$

Herein, the elastic parameters λ = Eν/[(1 + ν)(1 − 2ν)] and μ = E/[2(1 + ν)] denote the Lamé constants for an isotropic, homogeneous material. The viscous and the later discussed viscoelastic material are only characterized by the stress deviator. At first, the Lamé constants are calculated using a Young's modulus of $E = 1$  GPa, a common lower bound given for ice in literature (see Rist and others, 2002).

For the viscoelastic Maxwell material introduced in the following equation, the short-term material behavior is linear elastic while the long-term behavior resembles a viscous fluid. Consequently, this material provides a more accurate description of the ice behavior. Therefore, the constitutive relation is modified to

(7) $${{\bi \sigma} ^D} = 2\eta {{\bi{\dot\varepsilon}}_{\rm v}^D} = 2\mu {\bi \varepsilon} _{{\rm el}}^D \quad {\rm with}\quad {{\bi \varepsilon}^D} = {\bi \varepsilon}_{\rm el}^D + {\bi \varepsilon}_{\rm v}^D, $$

for the deviatoric stress tensor in Eqn (6). In this case, the deviatoric stresses in the elastic and viscous elements are equal to the total deviatoric stress, and the elastic ( ${\bi \varepsilon} _{{\rm el}}^D $ ) and viscous ( ${\bi \varepsilon} _{\rm v}^D $ ) deviatoric strains combine to result in the total deviatoric strain. For more details on continuum mechanics and rheological models characterizing the material behavior, see for example, Altenbach (2012) or any other textbook on continuum mechanics.

The system of coupled differential equations consists of the balance of linear momentum, the kinematic relation and the material equations. It is solved for discrete nodal displacements using the commercial finite element software COMSOL (http://www.comsol.com). The idealized ice-shelf domain has a length of $L = 5000 $  m and is discretized using triangular elements with Lagrange shape functions of quadratic polynomial order. The maximum element edge length is 10 m. Additionally, the mesh is refined to a maximum element size of 1 m at the ice front and at the ice-shelf surface within the first 1000 m from the front. This results in 215 942 degrees of freedom. The independence of computed results from the mesh has been verified. A coarsening of 5 m instead of 1 m leads to a maximum difference during the simulation time of 0.9 % of the maximum stress value at the upper surface. Furthermore, it was examined that the boundary condition at the inflow has no effect on the stress distribution near the ice front, if the length L is larger than 3500 m. Time stepping is auto-controlled by the time-dependent solver of COMSOL. The following section points out the boundary conditions, which are different from the beam theory. Especially, the crucial influence of the traction-free part at the ice front is considered.

BOUNDARY CONDITIONS AT GROUNDING LINE AND FREEBOARD

In order to obtain a unique solution, boundary conditions are necessary. For the 1-D beam theory in Reeh (1968) the boundary conditions are illustrated in Figure 2a, where the deflection w and the slope of the deflection ∂w/∂x are set to zero at the left boundary (grounding line, see also Fig. 1a). The remaining boundary conditions at the ice front consist of a vanishing shear force V(L) and a bending moment M(L), which is computed using the off-centre resultant force R(L) due to the depth dependent water pressure. As a consequence, the corresponding normal stress σ _xx is linearly varying along the vertical ice front and the shear stress σ _xz is zero. The springs at the bottom of the beam illustrate an elastic foundation (Winkler foundation) modeling the buoyancy forces of the water. The purple line in Figure 3d depicts the dimensionless tensile stress at the upper surface for these boundary conditions. This result is obtained by assuming a constant thickness of $H = 100$  m, see Figure 1a, a gravity acceleration of $g = 9.81\;{\rm m}\;{{\rm s}^{ - 2}}$ , and a density of sea water of ρ _sw = 1028 kg m⁻³. The axes are chosen according to Reeh (1968): the horizontal axis depicts (x − L)/H, the dimensionless distance to the ice front while the vertical axis reports (σ _xx  − σ _zz )/ρ _sw gH, the dimensionless stress difference after 1 d. The time of 1 d is arbitrarily chosen, but is identical for all methods and boundary conditions. The stress component σ _zz is negligibly small at the upper surface. However, Reeh (1968) included this quantity since he also considered the stress difference at the bottom of the ice shelf, where the component σ _zz corresponds to the water pressure.

Fig. 2. (a) Boundary conditions for the beam theory used in Reeh (1968). (b) Boundary conditions for the 2-D continuum approach using finite element simulation.

Fig. 3. (a–c) Separate illustration of different boundary conditions at the ice front for (a) Reeh (1968); (b and c) the finite element simulation. (d) Viscous stress states at the upper surface for the different boundary conditions after the simulation time of 1 d.

The boundary conditions for the 2-D finite element model are illustrated in Figure 2b. In order to represent the beam boundary conditions at the grounding line, the displacement u in the flow direction is set to zero at the inflow boundary. However, the subsequently shown stress distributions are independent of other constant values of u applied along the inflow boundary. The upper surface of the ice shelf is traction-free. The traction σ _n that balances the weight of the ice acts in the normal direction at the bottom of the ice shelf and is given by a Robin-type boundary condition

$${\sigma _{\rm n}} = {\rho_{{\rm sw}}}g({d_{\rm i}}H - w),$$

with the density ratio d _i = ρ _ice/ρ _sw (see Reeh, 1968) such that d _i H represents the depth below sea level. Thus, as the only volume force, gravity is compensated by buoyancy forces. Figures 3b, c illustrate the different boundary conditions at the ice front for the 2-D finite element simulations. The resulting stress states at the upper surface are depicted in the respective colors in Figure 3d. Significant differences arise from this traction-free part, the freeboard. The boundary condition in Figure 3b accurately models the stress state at the ice front of an ice shelf, where water pressure increases with depth and the upper vertical surface (highlighted by the red circle) is traction-free. The accurate modeling of the freeboard forces the stress to be zero at x = L due to the boundary condition at the ice front (see green curve in Fig. 3d). The boundary condition in Figure 3c is a 2-D representation of the boundary condition of Reeh (1968) and has the same stress resultant R(L) as the other two boundary conditions. As a result, the traction-free point is no longer located at sea level but is shifted upwards by Δh = ((1 − d _i )²/2) H. This can be derived by the comparison of stress resultants. Consequently, the water pressure at the bottom of the ice front for the condition in Figure 3c is given by p _o = ρ _sw gH(d _i + (1 − d _i)²/2) for the initial geometry, while the pressure for the case in Figure 3b is p _g = ρ _sw gHd _i. Traction boundary conditions yielding identical stress resultants R(L) at the ice front are considered for the following reason: the stresses, computed in the 2-D finite element simulation, are in good agreement with the value $N/({\rho _{{\rm sw}}}g{H^2}) = d_{\rm i}^{\rm 2}/2 $ obtained by the beam theory some distance away from the ice front boundary. For the following analysis of the stress distribution in a viscous and viscoelastic material, the freeboard boundary condition Figure 3b is applied, as this boundary condition is suitable for accurately describing the conditions at the calving front.

VISCOUS BEAM THEORY VERSUS VISCOUS PLANE STRAIN MODEL

Independent of the applied modeling approach and material law, the tensile stresses are located in the upper part of the ice shelf. This is caused by the boundary disturbance of the freeboard and the increasing water pressure with depth that induce a bending moment. Specifically, the largest tensile stress is located on the upper surface. In this section the deformations and the corresponding stresses and strains are computed using an incompressible (ν ≈ 0.5) viscous material model. Figure 4 depicts the surface stresses at different simulation times t for a density ratio of d _i = 0.9. The considered points in time are identical to the ones analyzed in Reeh (1968) to illustrate similarities and differences. The corresponding timescales are depending on the time factor f = η/(ρ _sw gH) to obtain results independent of the viscosity η. For these computations, the viscosity is set to $\eta { = 10^{14}} $  Pa s (see Greve and Blatter, 2009). Thus, assuming $H = 100$  m, $g = 9.81\;{\rm m}\;{{\rm s}^{ - 2}}$ and ρ _sw = 1028 kg m⁻³, corresponding to the values given in the previous section, the final time step 12.5 f corresponds to ~39 a. The dashed lines in Figure 4 depict the dimensionless stresses at the surface of a viscous beam according to the results shown in Figure 9 in Reeh (1968). The effect of using a 2-D model instead of the beam theory is significant, as illustrated by solid and dashed lines in Figure 4. In contrast to the results by Reeh (1968), the solid curves continuously decrease with time, and the location of the maximum tensile stress shifts towards the ice front. Normalized stresses similar to the results of the beam theory only appear for the first time step (t = 0) and in a certain distance from the front. Closer to the front, considerable differences due to the implemented boundary conditions can be seen for all time steps.

Fig. 4. Comparison of stress differences at the upper surface; solid lines indicate the results of the 2-D viscous material model and dashed lines indicate the results of the 1-D viscous beam.

The solid lines in Figure 5 depict the stress states for earlier points in time ( $t \le 0.31$ f) while the red dashed line repeats the result of the beam theory for the final time $t = 12.5f$ from Figure 4. It is found that the stress evolution of the 2-D model for early points in time is well described by the almost time-independent beam results: 10 d correspond to $t = 8.5 \cdot {10^{ - 3}} f$ (brown line), 30 d to $t = 2.5 \cdot {10^{ - 2}} f$ (gray line), and 1 a to $t = 0.31 f$ (black line). Figure 6 provides a possible explanation for this behavior. The strain component in the flow direction ε _xx monotonically increases with time and develops a maximum in the vicinity of the ice front. This strain component is related to the elongation of the geometry in the flow direction, which is not considered in the beam theory. Since the maximum value of ε _xx increases with time, the discrepancies become larger.

Fig. 5. Comparison of stress differences at the upper surface; solid lines indicate the 2-D viscous material model and red dashed line indicate the 1-D viscous beam for $t = 12.5f$ .

Fig. 6. Time evolution of the strain component ε _xx at the upper surface for a 2-D viscous material model.

VISCOELASTIC PLANE STRAIN MODEL

In the following, the same geometrical setup, with the freeboard boundary condition (Fig. 3b), is used to analyze the temporal evolution of the stress difference for the viscoelastic Maxwell material. The evolution of the stress difference at the upper surface is shown in Figure 7 for an incompressible material (ν ≈ 0.5), concentrating on three different instances in time during 1 a. As for the viscous material, the stress decreases with time. This behavior is illustrated in Figure 8, where the maximum of the dimensionless stress difference is plotted versus time. For comparability, two additional curves are added in this figure: the maximum stress differences of the purely viscous material from the previous section and of a purely elastic material. The stress state at the beginning of the viscoelastic simulation corresponds to the purely elastic response with σ ^D  = 2μ ε ^D (see Eqn (6)), which is rate independent and therefore constant. The viscous and viscoelastic curves monotonically decrease with time and virtually coincide after a certain time period.

Fig. 7. Viscoelastic stress states for an incompressible material for different times. The direction of the arrow indicates an increase in time.

Fig. 8. Maximum stress differences at the surface for incompressible linear elastic, viscous and viscoelastic material versus time.

Figures 9, 10 illustrate the stress states for the same setup as in Figures 7, 8 but for a lower Poisson's ratio of ν = 0.325, a common value given for elastic compressibility of ice in literature (Greve and Blatter, 2009). An initial increase with a subsequent monotonic decrease of the stress response is observed for the viscoelastic material, as indicated by the direction of the arrow in Figure 9. This behavior is also confirmed by the temporal evolution of the maximum stresses in Figure 10 for purely viscous and viscoelastic materials. The short-term behavior for the viscoelastic material is close to the elastic response. For longer time periods, the viscoelastic response converges towards the viscous response. For both materials, incompressible and compressible, the dimensionless distance of the stress maximum to the ice front is almost the same.

Fig. 9. Viscoelastic stress states for ν = 0.325 for different times. The direction of the arrow indicates an increase in time.

Fig. 10. Comparison of the maximum stress differences at the surface for compressible linear elastic, viscous and viscoelastic material versus time for ν = 0.325.

Bassis and others (2008) argued that the calving of ice shelves is controlled by the first (most tensile) principal stress. A final study therefore analyzes the influence of the most important geometric and material parameters on the maximum first principal stress in space and time $\sigma _1^{\mathop {\max} \nolimits_{x,t}} $ . Reasonable limits of the parameter values for typical Antarctic ice shelves are listed in Table 1. As the stress σ _zz in the vertical direction and the shear stress σ _xz are negligibly small at the top surface, the first principal stress

(8) $${\sigma _1} = \displaystyle{{{\sigma _{xx}} + {\sigma _{zz}}} \over 2} + \sqrt {{{\left( {\displaystyle{{{\sigma _{xx}} - {\sigma _{zz}}} \over 2}} \right)}^2} + \sigma _{xz}^2}, $$

is equal to the stress σ _xx in the flow direction. A graphical representation of $\sigma _1^{\mathop {\max} \nolimits_{x,t}} $ is already given in Figure 10. Results in Table 1 indicate that $\sigma _1^{\mathop {\max} \nolimits_{x,t}} $ increases with growing elastic material properties, Young's modulus E and Poisson's ratio ν. Furthermore, the maximum principal stress depends linearly on the geometric parameter H. The density ratio d _i, which controls the extent of the freeboard, influences $\sigma _1^{\mathop {\max} \nolimits_{x,t}} $ such that the maximum principal stress becomes larger if d _i decreases (increasing freeboard). The distance of the maximum principal stress to the ice front increases linearly with the elastic parameters and is inversely correlated to H and d _i. Increasing d _i due to the assumption of surrounding freshwater hence leads to a decreasing of both, the maximum stress and the distance of the respective location from the calving front. This especially applies to calving from floating glacier tongues into freshwater lakes (Trüssel and others, 2013). The value and the position of $\sigma _1^{\mathop {\max} \nolimits_{x,t}} $ , are independent of the viscosity, which is the only material parameter that affects the viscous flow. However, the viscosity has an influence on the time period that is needed to reach the maximum stress. This time period is also influenced by the elastic properties (Young's modulus E, Poisson's ratio ν) while the impact of the geometric quantities is negligible.

Table 1. Magnitude, position and evaluation time of the maximal principal stress dependent on relevant parameter for a viscoelastic material

Bold numbers mark deviations in the material parameters from the standard values used in this study.

In order to show that differences in the stress distribution obtained by the beam theory or the 2-D continuum approach, are independent of the uncertainty related to material parameters, Figures 11a, b illustrate the changes due to the variation of reasonable elastic material parameters for ice. For a better comparison, the geometric parameters, the thickness $H = 100 $  m, and the density ratio d _i = 0.9, are identical to the reference parameters used before. At later points in time, there is virtually no difference between the viscous and viscoelastic response in Figures 8, 10. Therefore, the solid and dashed lines in Figure 11 correspond to those with the respective colors in Figure 4; in fact, the blue ones fit $t = 1.25 f$ and the orange ones fit $t = 4.5 f$ . In Figures 11a, b, the shaded domain is limited by the upper bound with $E = 10 $  GPa and ν = 0.4, and the lower bound with $E = 0.1 $  GPa and ν = 0.1. All other stress distributions, related to reasonable variations of the elastic material parameters according to Table 1, are located within the shaded area.

Fig. 11. Shading indicates the stress variation due to elastic material parameter variations comparable with Table 1 for (a) $t = 1.25 f$ and (b) $t = 4.5 f$ of the viscoelastic 2-D approach; solid and dashed lines correspond to the lines in Figure 4 (dashed: beam theory, solid: viscoelastic 2-D continuum model).

DISCUSSION

It is repeatedly suggested in literature to assume ice as incompressible and to use a stress criterion for crack formation or calving of small icebergs due to bending. The applicability of such a criterion is checked by calculating the stress distribution at the upper surface and its temporal evolution with a 2-D finite element model. The stresses are then evaluated for a variety of different material parameters.

One assumption in this investigation is that temperature and viscosity in the ice shelf are constant. A mean annual surface temperature of − 5 to − 20°C indicates an ice viscosity between 10¹⁴ and 10¹⁵ Pa s (Cuffey and Paterson, 2010). Within this viscosity interval, Table 1 demonstrates the vanishing influence of the viscosity on the value of the maximum tensile stress and its distance to the ice front. The viscosity primarily influences the time interval needed to yield the maximum stress. Further parameter studies, not shown here, indicate that a nonlinear Glen-type viscosity also results in decreasing stresses. The amount and the spatial distribution of the stresses are independent of changes in viscosity. Therefore it is sufficient to assume a constant viscosity in the presented context. This is also described by Vaughan and others (2012), who stated that the temperature variation influencing the ice viscosity appears to have little impact on the magnitude of stresses at the top surface near the calving front of an ice shelf.

The largest tensile stress is located on the upper surface for all applied modeling approaches and material laws. The range of the computed maximum stress values overlaps with suggested stresses for crack nucleation (90–320 kPa), see Vaughan (1993). Resulting stresses of the beam theory approach show only small changes with respect to time (Fig. 4) and the maximum stress values only differ from the initial maximum by at most 9%. Starting with the same initial stress value and considering the same time period, the stresses monotonically decrease with time by more than 56% for the incompressible viscous or viscoelastic material computed by the 2-D continuum approach using finite elements. A closer comparison of surface stresses in Figure 5 indicates that the beam theory is only able to sufficiently represent the stresses at early points in time. Statements about stresses at later times require an approach including the spreading in flow direction. Figure 6 provides a possible explanation for the stress decrease using a 2-D simulation. Since the curves of the strain component ε _xx converges within some distance to the ice front (Fig. 6), the strain rates, as well as the stresses in the flow direction, decrease as a result of the viscous material law. If ε _xx reaches more than 1%, the stress response of the 2-D continuum model increasingly differs from the one computed with the beam theory. This is caused by the longitudinal flow, which cannot be modeled using beam theory assumptions. The stress relaxation for longer time periods was not noticed in the first 2-D finite element analysis of ice shelves by Fastook and Schmidt (1982). Their computational restrictions only allowed for limited time spans (26 d at most), simple shape functions (Q1P0 elements), and a relatively coarse grid. Hence, for this short-time interval the results presented in Figure 5 confirm the misleading observation of only small changes in the stress distribution, which led Fastook and Schmidt (1982) to the assumption that the time evolution does not need to be considered to study crack formation and propagation. At later times, either an other criterion must be applied or mechanisms must be included that lead to an increasing stress. An example of such a mechanism may be changes in the geometry at the ice front, as suggested by O'Leary and Christoffersen (2013).

In the following, the values that influence the stress maximum and its position are discussed. The focus is on the quantities, which lead to increasing stresses and thus provide an explanation for calving. The freeboard is the crucial parameter that considerably influences the formation of a stress maximum at the surface due to the boundary disturbance. This maximum is located within the range of 0.5–1.0 times the thickness upstream of the ice front, an interval commonly given in literature (for example, see Fastook and Schmidt, 1982). Thereby, the stress maximum directly depends on the size of the freeboard and increases linearly with the thickness of the ice shelf. Twice the thickness leads to a doubling of the maximum stress value, while the distance of the associated position to the ice front is influenced by circa 7%. As shown in Table 1, the ratio of the densities d _i, which also controls the extent of the freeboard, has the second largest influence on this stress. Consequently, the study with the thickness $H = 200 $  m and d _i = 0.9 as well as the study with $H = 100 $  m and d _i = 0.8 leads to similar stress values, as the freeboard is 20 m in both cases. The stress value for d _i = 0.8 increases by 61% in comparison with d _i = 0.9, given in the first line of Table 1 and its position to the ice front increases by 25%. The viscous beam theory accounts for the influence of freeboard through the magnitude of the applied bending moment, but it cannot realistically reproduce the stress-free boundary condition, particularly in the upper corner of the ice front.

Many of the parameters of the 2-D viscoelastic model are poorly constrained, and we therefore explore the sensitivity of our results to variations in material properties. Young's modulus has a limited impact of <13% on the magnitude of the maximum principal stress, as a softening only implies a slight decrease of the stresses. A more significant influence of Young's modulus is observed regarding the position of the maximum stress by up to 27%, see Table 1. This indicates that Young's modulus affects the size of a possible iceberg. In comparison with other parameters varied in this study, Poisson's ratio has only little influence on the magnitude of the maximum tensile surface stress (<4%) and the position to the ice front (<7%). However, the impact on the temporal evolution of $\sigma _1^{\mathop {\max} \nolimits_{x,t}} $ is remarkable. A Poisson's ratio of ν ≈ 0.5 (incompressible) leads to monotonically decreasing surface stresses with time, which means that the maximum stress is reached immediately at the beginning of the simulation (see Fig. 8). A calving law based on in compressible ice would therefore lead to essentially continuous calving. Only compressible ice leads to stress evolutions that result in maximum stresses after some finite period of time. The time needed to reach the maximum first principal stress in space and time is mostly influenced by the viscosity and Young's modulus, the two parameters defining the characteristic time of a Maxwell model. More precisely, increasing the viscosity leads to a longer evolution time of the maximum stress while the correlation between Young's modulus and the evolution time is inversely proportional. In the specific case of ν = 0.325, the maximum tensile stress is reached after ~20 d, a relatively short period of time. Ice-shelf calving based on a critical stress criterion seems unlikely during decreasing maximum principal stresses. Therefore, the result ${t_{{\sigma _{\max}}}} = 20 $ d can be seen as an indicator for the most probable timespan of new crack initiation at the surface. During this time period only small deformations occur in the flow direction. The beam theory, which neglects horizontal spreading, leads to a similar maximum stress for a very long period of time compared with the overall stress maximum of the finite element simulation with compressible viscous or viscoelastic material behavior (see Figs 6, 10).

Figure 11 illustrates the effect of uncertainties in elastic parameters on the stress distribution. Thereby, the shaded areas for $t = 1.25 f$ (Fig. 11a) and $t = 4.5 f$ (Fig. 11b) are at smaller stress values than the stress curves for the beam theory. The decrease of the stresses for longer time periods is thus not influenced by the uncertainties, since the decrease appears for all applied material parameters. The stress curves differ only slightly for different Poisson's ratios in the case of $E = 10 $  GPa. For smaller Young's moduli at early times, the influence of Poisson's ratio increases. Further simulations show that the stresses for smaller density ratios or thicker ice shelves are shifted to larger values (compare Table 1). However, the stress decrease with time is also recognizable for changes in thicknesses and density ratios. Moreover, for these results, the viscosity has no influence on the stress states at the top surface for the same considered points in time.

In conclusion, the application of a stress criterion under the assumption of incompressibility implies that the nucleus for the next calving process forms directly during the previous calving event. Then, an ice shelf without surface cracks should not be possible. However, if observations show that some time passes between a calving event and the next crack nucleation, the following cases are likely. The short-term behavior is assumed to be compressible for ice considered as a brittle material during crack formation and propagation, see for example, discussions of Rist and others (1999) or Schulson and Duval (2009). If ice is assumed to be incompressible, only a time-dependent change in geometry or boundary conditions at the ice front will lead to increasing stresses and therefore to crack nucleation.

CONCLUSIONS

In a fundamental work on ice-shelf calving Reeh (1968) expanded the linear elastic beam theory to a viscous material setting. However, the beam theory can neither consider the extension in flow direction nor accurately model the freeboard at the ice front. Therefore, the extension of the ice-shelf model to a 2-D continuum is necessary to model the spreading in flow direction and the consequent stress decrease with time. The influence of different boundary conditions indicates that the formation of the maximum tensile stress within a certain distance to the ice front only appears in the case of a stress-free upper part for the 2-D approach using finite elements. The application of boundary conditions associated with the beam theory approach leads to a maximum stress directly at the ice front of the 2-D continuum model. The boundary disturbance at the ice front hence dominates the stress state at the surface in the vicinity of the ice front. Thus, the thickness and the density ratio are especially important for the maximum principal stress and its position. This leads to the presumption that larger stress values can be observed, if the thickness or freeboard changes due to melting or freezing, which was not part of this study but should be investigated in future works and has already been studied for tidewater glaciers, see, for example, O'Leary and Christoffersen (2013).

Viscoelasticity is important to model crack formation or calving happening on short timescales. The instantaneous elastic response is crucial to decide whether a crack is stable or instable. On the other hand if a crack is stable and therefore does not expand due to the instantaneously existing elastic strain, the viscous stresses might lead to crack propagation at later points in time. The development of the maximum stress is regulated by Poisson's ratio while the time evolution of the stress is influenced by the viscosity and Young's modulus controlling the characteristic time of the Maxwell model. In case of longer time periods than in the presented work, the total strain may exceed 10% and thus violates the assumption of small deformations and a large deformation model must be considered. However, in the present context, it is sufficient to analyze the stresses for small strains, as the maximum stress response occurs in a rather short time. Therefore, a restriction of the maximum time ${t_{\max}} = 12.5 f$ was reasonable.

The application of stress criteria is not only for crack initiation but rather for calving implies several problems. On the one hand the stresses increase with the thickness, i.e. thicker ice shelves are more prone to calving by bending than thinner ice shelves. On the other hand the time evolution of the stress in an incompressible viscous or viscoelastic fluid does not yield a stress maximum, as the tensile stresses monotonically decrease over time. Hence, if a critical stress criterion was applied, the next calving event would occur directly at the beginning of the considered timespan (here t = 0). This might result in the collapse of the whole ice shelf. We therefore conclude that a stress criterion for discrete small calving is only applicable for a compressible ice behavior or additional, uninvestigated, time-dependent geometric evolution. Possibly the best criterion is a self similarity criterion, which we are still working on. If it is known that some stress distribution was critical in the past and is reached again, the next calving event will take place.