Home Journals Articles Updates About
CN
image
收藏切换
Modeling calving process of glacier with dilated polyhedral discrete element method
收藏切换
PDF
Lu Liu1, Ji Li1, Qizhen Sun2, Chunhua Li2, Sue Cook3, Shunying Ji1, *
Acta Oceanologica Sinica | 2021, 40(7) : 159 - 169
Less
收藏切换
Acta Oceanologica Sinica | 2021, 40(7): 159-169
• Marine Geology •
Modeling calving process of glacier with dilated polyhedral discrete element method
Full
Lu Liu1, Ji Li1, Qizhen Sun2, Chunhua Li2, Sue Cook3, Shunying Ji1, *
Affiliations
  • 1 State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology, Dalian 116024, China
  • 2 National Marine Environmental Forecasting Center, Beijing 100081, China
  • 3 Institute for Marine and Antarctic Studies, University of Tasmania, Hobart 7001, Australia
Published: 2021-07-25 doi: 10.1007/s13131-021-1819-x
Outline
收藏切换

Mass loss caused by glacier calving is one of the direct contributors to global sea level rise. Reliable calving laws are required for accurate modelling of ice sheet mass balance. Both continuous and discontinuous methods have been used for glacial calving simulations. In this study, the discrete element method (DEM) based on dilated polyhedral elements is introduced to simulate the calving process of a tidewater glacier. Dilated polyhedrons can be obtained from the Minkowski sum of a sphere and a core polyhedron. These elements can be utilized to generate a continuum ice material, where the interaction force between adjacent elements is modeled by constructing bonds at the joints of the common faces. A hybrid fracture model considering fracture energy is introduced. The viscous creep behavior of glaciers on long-term scales is not considered. By applying buoyancy and gravity to the modelled glacier, DEM results show that the calving process is caused by cracks which are initialized at the top of the glacier and spread to the bottom. The results demonstrate the feasibility of using the dilated polyhedral DEM method in glacier simulations, additionally allowing the fragment size of the breaking fragments to be counted. The relationship between crack propagation and internal stress in the glacier is analyzed during calving process. Through the analysis of the Mises stress and the normal stress between the elements, it is found that geometric changes caused by the glacier calving lead to the redistribution of the stress. The tensile stress between the elements is the main influencing factor of glacier ice failure. In addition, the element shape, glacier base friction and buoyancy are studied, the results show that the glacier model based on the dilated polyhedral DEM is sensitive to the above conditions.

glacier calving  /  discrete element method  /  dilated polyhedral element  /  bond and fracture model
Lu Liu, Ji Li, Qizhen Sun, Chunhua Li, Sue Cook, Shunying Ji. Modeling calving process of glacier with dilated polyhedral discrete element method[J]. Acta Oceanologica Sinica, 2021 , 40 (7) : 159 -169 . DOI: 10.1007/s13131-021-1819-x
Full Text
Less
1 Introduction
Less
The glacial calving process occurs when the front edge of a glacier or ice shelf (also known as its terminus) breaks away and is released into the ocean. This dynamic ice loss process is considered alongside thermodynamic ice melting when using ice sheet models to obtain an accurate prediction of sea level rise (Dyer et al., 2013; Zemp et al., 2019). In addition, coupled simulations of atmospheric, oceanic, and ice-sheet processes show that the ice-sheet mass loss can affect global ocean currents. The meltwater generated by glacial disintegration from Greenland may thus accelerate the further calving of glaciers in the Antarctic by slowing the Atlantic overturning circulation (Golledge et al., 2019). This positive feedback will have a more profound impact on climate change.
Introducing a glacial calving law into existing ice sheet models is an important method of understanding the role of iceberg calving in ice sheet mass loss (Schlemm and Levermann, 2018). A complete law of calving should include environmental forcing, glacier geometry, material properties, and preexisting flaws (Alley et al., 2008). However, an accurate and sufficiently complex calving law describing the location or rate of calving in all glaciers and ice shelves has yet to be achieved (Benn et al., 2007). One approach to develop a more accurate calving law is to parameterize the calving process from simulations, so the results can be adapted to large spatial and time scale simulations.
With increasing interest in glaciers in the Arctic and the Antarctic region recently, both continuous and discontinuous numerical methods have been applied to explore calving laws. Simple laws have previously been based on field data relating to the calving process, such as the velocity of the glacier along-flow (Alley et al., 2008) or the relationship between water depth and ice thickness (Pelto and Warren, 1991). The tensile and shear failure of a rectangular glacier has also been investigated using the stress field obtained from a 2-D full-Stokes finite element method (Ma et al., 2017). Schlemm and Levermann (2018) derived the stress field of calving glacier from the two-dimensional Stokes equations and proposed a calving law based on yield stress (Schlemm and Levermann, 2018). Many previous studies use continuous models such as FEM for glacial simulations including iceberg calving (Cook et al., 2018; Pralong and Funk, 2005; Todd et al., 2018; van Dongen et al., 2018), but the disadvantage of continuous models is that fractures cannot form explicitly, therefore the interaction of fractures and the surrounding stress field must be based on a parameterization (Benn and Aström, 2018).
A discrete model based on a hybrid-lattice spring network approach has also been applied to investigate the unstable and shear damage phenomena of glaciers (Koehn and Sachau, 2014; Åström et al., 2014; Benn and Åström, 2018). This method combines the brittleness of short-term damage and long-term viscous flow (Koehn and Sachau, 2014). The most commonly used discontinuous numerical method in glacier simulations is the discrete element method (DEM). DEM, proposed by Cundall and Strack in 1979 (Cundall and Strack, 1979), has a wide range of applications in geotechnical engineering, ice engineering, and particle physics (Cleary, 2010; Liu and Ji, 2018). Due to the need to study the fundamental mechanisms of the calving process, this method of providing macroscopic and microscopic perspectives has received extensive attention recently. Bassis and Jacobs (2013) simulated the calving process by the discrete element method, finding that glacier geometry, which controls the stress state within the glacier, is the main factor leading to different calving regimes (Bassis and Jacobs, 2013). The sphere-based discrete element model has also been developed to simulate the flow and calving of glacier geometries, which includes elastic, brittle-fracture, and slow viscous behavior (Åström et al., 2013; Riikilä et al., 2015). Benn et al. (2017) used the Elmer/Ice finite element model and Helsinki DEM (HiDEM) discrete element model to simulate the breaking process of a tidewater glacier. A series of factors such as melt under-cutting, buoyancy and cliff geometry were found to contribute to the calving process (Benn et al., 2017).
Continuous models have the advantage of being able to represent the lateral hydrostatic pressure of the surrounding ocean water on tidewater glaciers. The geometric deformation of the glacier over a long time scale because of the larger timestep, so that a more reliable state of stress within the glacier can be obtained. But their disadvantage is that it is difficult to represent the effect of crevasses on the surrounding stress field, and therefore it is difficult to predict their propagation. Discontinuous models such as the discrete element method accurately simulate the crack propagation process during disintegration. Another advantage of the discontinuous method is that it can consider the buttressing force on ice cliffs from the ice mélange (a mixture of ice fragments and sea ice often found near calving glaciers), which is proved important in calving (Burton et al., 2018; Robel, 2017). A simulation of glacier fracture based on the discrete element model also shows that fracture and transport of detached icebergs both play an important role in the disintegration of glaciers (Bassis and Jacobs, 2013). On the other hand, the element shapes used in glacial discrete element models are mainly sphere (circles for 2D) and cube (Åström et al., 2014; Benn et al., 2017; Burton et al., 2018; Riikilä et al., 2015; Robel, 2017). It is believed that the element shape has a significant influence on the dynamic behavior of granular matter (Lu et al., 2015). But the influence of element shape in the simulation of glacier is unclear. From the aspect of computational efficiency, the continuous models can simulate longer timescale cases comparing with the discontinuous models because the continuous models can use a much larger time step in the motion integral.
The DEM with dilated polyhedral elements can be applied to simulate the calving process of tidewater glaciers, providing the potential for irregular convex-shaped elements. The interaction between ice and water is represented by assuming the drag force acting on each element. The elements are bonded together to form a seamless glacier. Factors affecting the disintegration of glaciers are still unclear, but due to their inherent instability, the tidewater glaciers usually experience calving regularly and the disintegration of the calving front can be simulated on short-term scales. In this study, a dilated polyhedral DEM (DPDEM) combined with a bond and fracture model is applied to simulate tidewater glacier calving. The internal stress changes can be obtained at the same time, and the size distribution of icebergs produced during calving process is counted. Based on the simulation of different element shapes, bottom friction and buoyancy conditions, we proved the sensitivity of this model to these factors. Simulation results show that the DPDEM is a potential and effective method for glacial calving simulation.
2 Dilated polyhedron discrete element method
Less
2.1 Dilated polyhedral element and its contact model
A dilated polyhedral element is constructed by the Minkowski sum of a sphere and a polyhedron. Using elements of this shape can prevent several complex and inefficient geometric contact detections on each time step. The Minkowski sum between two geometric shapes A and B can be defined as
${\boldsymbol{A}} \oplus {\boldsymbol{B}} = \left\{ {{\boldsymbol{x}} + {\boldsymbol{y}}|{\boldsymbol{x}} \in {\boldsymbol{A}},{\boldsymbol{y}} \in {\boldsymbol{B}}} \right\},$
where $ {\boldsymbol{x}} $ and $ {\boldsymbol{y}} $ are vectors of sphere and polyhedral surface points, $ \oplus $ is the Minkowski sum operator which means sweeping one set around the profile of the other. Figure 1 shows the process of polyhedral dilation by the Minkowski sum. The polyhedron becomes smooth at the corners while convexity does not change. There are smooth spheres, cylinders and flat faces on the surface of dilated polyhedron, which leads to six types of contact patterns including face-face, cylinder-cylinder, sphere-sphere, sphere-cylinder, sphere-face, and cylinder-face. Different contact force models are established for different contact patterns, which overcomes the limitations of the linear contact force model and avoids the empirical determination of the stiffness parameter (Liu and Ji, 2018).
The normal force of different contact patterns is modeled by the normal stiffness parameter $ {k}_{\rm{n}} $ given in the following formula, which is defined as the normal contact stiffness:
${{\boldsymbol{F}}_{\rm{n}}} ={k_{\rm{n}}}{{\boldsymbol{\text{δ}}}^{\kappa}_{\rm{n}}}-{C_{\rm{n}}}{\left|{{{\boldsymbol{\text{δ}}}_{\rm{n}}}} \right|^{\kappa-1}}{{\dot{\boldsymbol{\text{δ}}}}_{\rm{n}}},$
where $ \kappa $ is usually 1.5 in the Hertzian model, which does not have particular physical significance; $ {\boldsymbol{\text{δ}}}_{n} $ is the normal overlap; and ${\dot{\boldsymbol{\text{δ}}}}_{\rm{n}}$ is the normal relative velocity. ${\boldsymbol{\text{δ}}}_{\rm{n}}$ and ${\dot{\boldsymbol{\text{δ}}}}_{\rm{n}}$ are projections of related variables in the vector of contact normal; $ {k}_{\rm{n}}=4{E}^{*}\sqrt{{R}^{*}}/3 $, where $ {E}^{*} $ and $ {R}^{*} $ are equivalent elastic modulus and equivalent radius respectively; $ {C}_{\rm{n}} $ is the normal damping coefficient, the variable is assigned by the following formula:
${C}_{\rm{n}}={\zeta }_{\rm{n}}\sqrt{{m}_{{\rm{AB}}}{k}_{\rm{n}}},$
where $ {m}_{{\rm{AB}}} $ is the equivalent mass of element A and element B, calculated by $ {m_{{{\rm{AB}}}}} = {m_{\rm{A}}}{m_{\rm{B}}}/({m_{\rm{A}}} + {m_{\rm{B}}}) $, and $ {\zeta }_{\rm{n}} $ is dimensionless damping coefficient given by:
${\zeta _{\rm{n}}} = \frac{{ - {\rm{ln}}e}}{{\sqrt {{\pi ^2} + {\rm{l}}{{\rm{n}}^2}e} }},$
where e is the restitution coefficient, which equals 0.1 here. The tangential force $ {{\boldsymbol{F}}_{\rm{s}}} $ of each contact candidate is defined by the tangential elastic force $ {\boldsymbol{F}}_{\rm{s}}^{\rm{*}} $ without a tangential viscous force, and then related with the friction force:
${\boldsymbol{F}}_{\rm{s}}^{*}={k}_{\rm{s}}{\left|{\boldsymbol{\text{δ}}}_{\rm{s}}\right|}^{\kappa -1}{{\text{δ}}}_{\rm{s}},$
${{\boldsymbol{F}}_{\rm{s}}} = {\rm{min}}\left({\left| {{\boldsymbol{F}}_{\rm{s}}^{\rm{*}}} \right|,\mu \left| {{{\boldsymbol{F}}_{\rm{n}}}} \right|} \right){{\text{δ}}_{\rm{s}}}/\left| {{{\text{δ}}_{\rm{s}}}} \right|,$
where $ {k}_{\rm{s}} $ is the tangential stiffness proportional to the normal stiffness as $ {k_{\rm{s}}} = {r_{{\rm{sn}}}}{k_{\rm{n}}} $; ${\boldsymbol{\text{δ}}}_{\rm{s}}$ is the vector of tangential elastic deformation; $ \mu $ is the sliding friction coefficient between elements. According to the relationship between the elastic modulus and shear modulus of the isotropic material, $ {r_{{\rm{sn}}}} $ is defined as $ 1/2(1+\nu) $, where $ \nu $ is Poisson’s ratio.
2.2 Bond and fracture model of dilated polyhedron element
A brief introduction of bond–failure model in the dilated polyhedral DEM is given in Liu and Ji (2019). The weak form of the governing equation can be derived for each bond point between two rigid elements based on elastic theory:
$\int _S^{}{\text{δ}} {{\boldsymbol{u}}^{\rm{T}}}\left({\text{σ}{\boldsymbol{n}}} \right){\rm{d}}S + \int _V^{}\rho {\text{δ}} {{\boldsymbol{u}}^{\rm{T}}}{\boldsymbol{f}}{\rm{d}}V = \int _V^{}{\text{δ}} {{\boldsymbol{u}}^{\rm{T}}}\rho {\ddot{\boldsymbol u}}{\rm{d}}V,$
where $\text{δ}$ is the variational symbol in the functional analysis; u= {u v w}T is the displacement of the mass center in x, y, z directions, $\ddot u $ is the acceleration of the mass center, n is normal of face S, $ \rho $ is the mass density and the first part on the left indicates the cohesion on the bonded interface, the second part represents the body force f and the right term is force related to acceleration. To avoid this complex calculation of surface integral, an evaluation of the stress and strain is developed by considering a similar approach in the generalized contact model (GCM) (Azevedo and Lemos, 2005). According to the dynamic relaxation method, the continuum can be discretized in space, i.e., the stiffness matrix can be regarded as discretized bond springs. Figure 2 shows the construction of our bond model for the glacier structure, and the contact force model between bonded elements.
Bond springs connect the corresponding vertices on the adjacent faces that are initially in a parallel state. Due to the relative displacement of the bonded vertices, the normal strain $ {\varepsilon }_{\rm{n}} $ and shear strain $ {\varepsilon }_{\rm{t}} $ are obtained from Eqs (8) and (9):
${\varepsilon _{\rm{n}}} = \frac{{{\boldsymbol{d}} \cdot {\boldsymbol{n}}}}{{{C_{ij}}}},$
${\varepsilon _{\rm{t}}} = \frac{{|{\boldsymbol{d}} - ({\boldsymbol{d}} \cdot {\boldsymbol{n}}){\boldsymbol{n}}|}}{{{C_{ij}}}},$
where d is the vector linking two bond points, ${\boldsymbol{n}} = \left({{{\boldsymbol{n}}_i} - {{\boldsymbol{n}}_j}} \right)/ {\rm{norm}} \left({{{\boldsymbol{n}}_i} - {{\boldsymbol{n}}_j}} \right)$ is the unit normal vector of bonded face i and j, norm$(\boldsymbol X) $ means the module of vector X. The characteristic length $ {C}_{ij} $ is the sum of the distances from the center of mass of each element to the bonding surface. The bond stress is determined based on the elastic law of anisotropic material, and both the elastic component $ {\text{σ}_{\rm{e}}} = {\{ {\sigma _{\rm{e}}},{\tau _{\rm{e}}}\} ^{\rm{T}}} $ and viscous component $ {\text{σ}_{\rm{v}}} = {\{ {\sigma _{\rm{v}}},{\tau _{\rm{v}}}\} ^{\rm{T}}} $ are calculated in directions perpendicular and parallel to the bonded faces, $ \sigma $ and $ \tau $ are the normal and tangential stresses respectively. ${\sigma _{\rm{e}}} = {{{{{k}}}}_{\rm{n}}}{\varepsilon _{\rm{n}}}$, $ {\tau _{\rm{e}}} = {k_{\rm{s}}}{\varepsilon _{\rm{t}}} $, $ {\sigma _{\rm{v}}} = \beta {k_{\rm{n}}}{\dot \varepsilon _{\rm{n}}} $, and $ {\tau _{\rm{v}}} = \beta {k_{\rm{s}}}{\dot \varepsilon _{\rm{t}}} $, $ \beta $ is constant scalar taken as 2.0×10−6 here to introduce damping related to the stiffness. The subscripts e and v above represent the elastic and viscous components. The bond force of each bond spring is then considered to be:
${{\boldsymbol{F}}_{\rm{b}}} = \left({{\text{σ}_{\rm{e}}} + {\text{σ}_{\rm{v}}}} \right) \cdot A/n,$
where A is the area of the bonded face and n is the number of bonded points on the bonded face.
For each pair of bonded points, the normal strength $ \bar{\sigma } $ is preset and shear strength $ \bar{\tau } $ is chosen by Mohr-Coulomb criterion, and in our hybrid strength criterion, damage happens as:
${\left({\frac{{\left\langle \sigma \right\rangle }}{{\bar \sigma }}} \right)^2} + {\left({\dfrac{\tau }{{\bar \tau }}} \right)^2} \geqslant 1,$
where $ \left\langle \cdot \right\rangle $ is the Macaulay bracket, and $ \left\langle x \right\rangle=x $ for $ x\geqslant 0 $, $ \left\langle x \right\rangle=0 $ for $ x < 0 $. The variable describing the degree of material damage is $ D $, defined as:
$D = \frac{{\text{δ} _{\rm{m}}^{\rm{f}}(\text{δ} _{\rm{m}}^{{\rm{max}}} - \text{δ} _{\rm{m}}^0)}}{{\text{δ} _{\rm{m}}^{{\rm{max}}}(\text{δ} _{\rm{m}}^{\rm{f}} - \text{δ} _{\rm{m}}^0)}},$
where the equivalent deformation ${\text{δ} _{\rm{m}}} = \sqrt {{{\left\langle {{\text{δ} _{\rm{n}}}} \right\rangle }^2} + \text{δ} _{\rm{s}}^2}$. The superscripts 0 and f denote the initial moment of damage and the moment of complete fracture, respectively. Here $\text{δ} _{\rm{m}}^{\rm{f}} = 2{G^{\rm{c}}}/\sigma _{\rm{m}}^0$, the equivalent stress $ {\sigma _{\rm{m}}} = \sqrt {{{\left\langle \sigma \right\rangle }^2} + {\tau ^2}} $. $\text{δ} _{\rm{m}}^{{\rm{max}}}$ is the maximum equivalent deformation during the loading process, as shown in Fig. 3. The critical energy $ {G^{\rm{c}}} $ is related to the fracture energy in tensile mode, $ G_{\rm{I}}^{\rm{c}} $, and the fracture energy in the shear mode, $ G_{{\rm{II}}}^{\rm{c}} $, which is defined as:
${G^{\rm{c}}} = G_{\rm{I}}^{\rm{c}} + (G_{{\rm{II}}}^{\rm{c}} - G_{\rm{I}}^{\rm{c}}){\left({\frac{{G_{{\rm{II}}}^{\rm{c}}}}{{G_{\rm{I}}^{\rm{c}} + G_{{\rm{II}}}^{\rm{c}}}}} \right)^{\eta} },$
where $ \eta $ is a scalar constant, normally $ \eta =1.75 $ for brittle materials (Camanho et al., 2003).
2.3 Element motion
Each dilated polyhedral element may be subject to gravity, inter-element contact forces, bond forces, buoyancy, and the motion-induced drag force of water (Liu and Ji, 2018). $ {\boldsymbol{F}}_{t} $ is the summation of all forces experienced by the element at time t. The Verlet algorithm was used to carry out the translation in the global coordinate system, through which the position and velocity of the element are calculated at each DEM time step. We use quaternions to store the local position of each element point relative to the initial position after rotation, since the Euler angles may lead to a singularity in special circumstances (Lu et al., 2015). The Verlet formula is defined as follows:
${{\boldsymbol{V}}_{t + \Delta t}} = {{\boldsymbol{V}}_t} + \frac{1}{2}\left({{{\boldsymbol{a}}_t} + {{\boldsymbol{a}}_{t + \Delta t}}} \right)\Delta t,$
${{\boldsymbol{x}}_{t + \Delta t}} = {{\boldsymbol{x}}_t} + {{\boldsymbol{V}}_{t + \Delta t}} \cdot \Delta t + \frac{1}{2} \cdot {{\boldsymbol{a}}_t} \cdot \Delta {t^2},$
${{\boldsymbol{a}}_{t + \Delta t}} = {{\boldsymbol{F}}_{t + \Delta t}}/m,$
where $ m $ is the element mass. The quaternion $ {\boldsymbol{q}} = \left\{ {{q_1},{q_2},{q_3},{q_4}} \right\} $ is employed to calculate the rotational motion of each element. The correspondence between the torque and the rotational acceleration $ \dot{\omega } $ in local coordinate is defined as:
$\left\{\begin{aligned}& \dot \omega _x^{\rm{b}} = \frac{{M_x^{\rm{b}}}}{{{I_{xx}}}} + \left({\frac{{{I_{yy}} - {I_{zz}}}}{{{I_{zz}}}}} \right)\omega _y^{\rm{b}}\omega _z^{\rm{b}},\\& \dot \omega _y^{\rm{b}} = \frac{{M_y^{\rm{b}}}}{{{I_{yy}}}} + \left({\frac{{{I_{zz}} - {I_{xx}}}}{{{I_{yy}}}}} \right)\omega _z^{\rm{b}}\omega _x^{\rm{b}},\\& \dot \omega _z^{\rm{b}} = \frac{{M_z^{\rm{b}}}}{{{I_{zz}}}} + \left({\frac{{{I_{xx}} - {I_{yy}}}}{{{I_{zz}}}}} \right)\omega _x^{\rm{b}}\omega _y^{\rm{b}},\end{aligned} \right.$
where $ {\boldsymbol{\omega }} $ represents the rotational velocity; $ {I}_{\alpha \alpha } $ is the principal moment of inertia around the $ \alpha $ axis including x, y and z; b represents the body-fixed (or local) coordinate system. The relationship between the local coordinator and the global coordinate of each point $ {\boldsymbol{e}}^{\rm{b}} $ can be transformed by matrix ${\boldsymbol{A}}$, G represents the global coordinate system:
${{\boldsymbol{e}}^{\rm{b}}} = {\boldsymbol{A}} \cdot {{\boldsymbol{e}}^{\rm{G}}},$
where the transition matrix $ {\boldsymbol{A}} $ is defined as
${\boldsymbol{A}} = \left( {\begin{array}{*{20}{c}}{q_0^2 + q_1^2 - q_2^2 - q_3^2}&{2({q_1}{q_2} + {q_0}{q_3})}&{2({q_1}{q_3} - {q_0}{q_2})}\\{2({q_1}{q_2} - {q_0}{q_3})}&{q_0^2 - q_1^2 + q_2^2 - q_3^2}&{2({q_2}{q_3} + {q_0}{q_1})}\\{2({q_1}{q_3} + {q_0}{q_2})}&{2({q_2}{q_3} - {q_0}{q_1})}&{q_0^2 - q_1^2 - q_2^2 + q_3^2}\\\end{array}} \right).$
The moment in the global coordinate system can be converted to the local coordinate system with Eq. (16). The quaternion varies with the rotational velocity in the local coordinate system according to:
$\left({\begin{array}{*{20}{c}}{{{\dot q}_0}}\\{{{\dot q}_1}}\\{{{\dot q}_2}}\\{{{\dot q}_3}}\end{array}} \right) = \frac{1}{2}{\boldsymbol{W}}\left({\begin{array}{*{20}{c}}0\\{\omega _x^{\rm{b}}}\\{\omega _y^{\rm{b}}}\\{\omega _z^{\rm{b}}}\end{array}} \right),$
where ${{\text{ω}}^{\rm{b}}} = {\left\{ {\omega _x^{\rm{b}}\;\omega _x^{\rm{b}}\;\omega _x^{\rm{b}}} \right\}^{\rm{T}}}$ is the rotational velocity in the local coordinate system; $ {\boldsymbol{W}} $ is defined as:
${\boldsymbol{W}} = \left({\begin{array}{*{20}{c}}{{q_0}}&{ - {q_1}}&{ - {q_2}}&{ - {q_3}}\\{{q_1}}&{{q_0}}&{ - {q_3}}&{{q_2}}\\{{q_2}}&{{q_3}}&{{q_0}}&{ - {q_1}}\\{{q_3}}&{ - {q_2}}&{{q_1}}&{{q_0}}\end{array}} \right).$
The quaternion can be obtained by step-by-step integration. Accordingly, the transition matrix, rotational velocity, and the new position of each vertex on the element can be updated.
3 Calving process simulation with DPDEM
Less
The contact model and bond–failure model of dilated polyhedral DEM have been used and validated for sea ice simulation in our study (Liu and Ji, 2018, 2019). These numerical models are also used for rock breakage. Thus this method is proved to be an effective approach in the simulation of brittle materials. The glacier calving process over a short time due to inherent instability is simulated by DPDEM. The bond strength is initialized with Weibull distribution on the glacier.
3.1 Glacier model
The tidewater glacier is initialized by extruding triangular elements meshed on a preset geometry. Since the third dimension is not considered in the extruding direction, the following simulation is performed on the plane. The calving of an ocean-terminating glacier is severely affected by the geometry of its terminus (Mercenier et al., 2018). The approximate natural shape of the glacier front is set to gradually thin in the along-flow direction, with under-cutting at the terminus caused by basal melting (Bartholomaus et al., 2013), as shown in Fig. 4a. It is assumed that the part below the sea level on the right-hand side of the glacier will melt due to warm ocean currents, and the length of the melt undercut, LU, increases linearly with the distance from the sea. Several parameters that control the geometry of the glacial front are defined: glacier height H and length L, undercut length LU and height HU, water depth DW, and flotation depth DF. The model is grounded on a flat bottom face with no friction, and we have a fixed vertical boundary on the left side to prevent it from disintegrating backward. Figure 4b illustrates the dilated polyhedral elements in this case.
Geometric parameters and physical parameters in the simulation used in this study are shown in Table 1. The research about the fracture energy of the glacier is rather rare so that similar researches on sea ice are considered. The size-independent fracture energy is 15 N/m according to the experimental and analytical research of Dempsey et al. (1999). But the fracture energy used in numerical simulations has different values, including 12 N/m by Paavilainen et al. (2011) and Paavilainen and Tuhkuri (2012), 15 N/m by Lu et al. (2014), 20 N/m, 28 N/m and 50 N/m by Gürtner et al. (2010). On the other hand, considering the hybrid model of the cohesive element for the breakage simulation of the rock, which is also the brittle material that is similar to ice, the mode-II fracture energy can be 2–5 times the mode-I fracture energy (Ma et al., 2016). Hence, the mode-I and mode-II fracture energy are 20 N/m and 40 N/m respectively. According to the research of Camanho et al. (2003), the power coefficient η in Eq. (11) for brittle material can be 1.75. Thus η = 1.75 is used in this study.
The bond strength represents the tensile and shear resistance of bonds, which is an important property of the material. Although the glacier is typically a nearly continuous structure, towards the terminus, where the glacier flows rapidly, the ice on the leading edge wears under natural conditions, resulting in a difference in the strength of the ice. Thus, the bond strength between each bonded element follows the Weibull distribution with a shape parameter of 2.0 and a size parameter of 4.0. The normal bond strength is between 1.0 MPa and 1.3 MPa in this case. The statistical probability density distributions of normal tensile bond strengths used in the simulation are displayed in Fig. 5, where each square in the histogram has a width of 0.01 MPa. The cohesion of each pair of bonded points is three times the tensile bond strength.
3.2 Calving simulations
This case is performed on an HP Z8 G4 workstation, which is equipped with Intel Xeon Silver 4114 processor and NVIDIA Quadro GP100. We simulated the calving process of 9.389 s, which took an average of 2.527 h. Figures 6 and 7 show the results of the calving process simulation up to 9.389 s. In Fig. 6, green rods represent the unbroken bonds in the glacier model, and red color rods represent bond pairs that have broken at this time step. After the bonds are broken, the further the relative movement distance is, the darker color shown in each figure. The time slice of the calving process and velocity distribution of polyhedron elements is shown in Fig. 7. It can be seen from these two figures that the calving of glacier, in this case, begins with the generation of a top crack until the crack spreads to the bottom of the ice, and then a shear band is formed throughout the fracture zone. There are crisscross shear bands in the glacier which produce a similar grinding effect.
The size distribution of the broken iceberg fragments after disintegration can be important for understanding their future behavior (Savage et al., 2000). Since the model geometry is two-dimensional, the size is measured as the area on the 2D vertical plane. The size distribution of the ice fragments produced in the disintegration can be fitted by a power function, $ N\left(A\right)\propto {A}^{-D} $, with the exponent $ D\approx 1.44 $, as shown in Fig. 8. According to the Arctic ice floe measurements of Rothrock and Thorndike (1984), the exponents are in a broad range, about 1.7 < D < 2.5. The difference between numerical and field results can be complicated. The size in the field data is calculated as $ A=\sqrt[3]{m} $ where m is the mass of the ice fragment, while the size in the numerical data is calculated as the maximum length of each ice fragment. Besides, the 3D simulation can provide the size on the other dimension, but the simulation in this study is 2D.
In this study, the von Mises stress ($ {\sigma }_{s}) $ can be calculated by:
${\sigma }_{s}=\sqrt{{\sigma }^{2}+{3\tau }^{2}},$
where $ \sigma $ is the normal stress; $ \tau $ is the shear stress; and $ {\sigma }_{s} $ is calculated at each bonded vertex. We can obtain the stress results in the computing domain with an interpolation method. Figure 9 shows the von Mises stress at different times in the glacier (0–0.65 s), and bonds that have been broken are indicated by red dots. The results indicate that the von Mises stress at the right bottom of the glacier gradually increases with the separation of the left side of the glacier when fractures form along the boundary. When the ice surface at the glacier terminus begins to break due to excessive stress, the changed glacier geometry causes the internal stress to redistribute.
However, the von Mises stress cannot determine whether the fracture between the bonded elements is caused by stretching or compression. Figure 10 shows the time slice of the normal stress $ \sigma $ in the glacier. It indicates that the generation and propagation of glacial cracks in this example are mainly caused by the tensile stress exceeding the ice strength.
Based on the glacier geometry and physical parameters completely consistent with the above experiment, we use the cubic element to simulate the calving process under the same buoyance conditions. As shown in the following Fig. 11 , the glacier geomtry is composed of cubic elements, the diagonal connection between the square units is only for display and does not mean that the element is split into two. The front end of the glacier is composed of 9 369 elements.
Figure 12 shows the comparison of the glacier velocity distribution at the end of the simulation, i.e., t = 9.389 s using cubic elements and triangular elements. The retreat distances at the bottom of the glacier in Figs 12a and b are 54 m and 53 m, respectively, while the retreat distances at the top are 49 m and 57 m, respectively. The retreat distance is defined as the distance between the front end of the glacier before calving and the front end after calving.
The overall retreat of the glacier is not much different in the two elements, and the retreat distance of the top of the glacier is smaller than that of the bottom of the glacier, but a more detailed analysis will find that due to the limitation of the element shape, the direction of the grinding of the fragments generated by the glacier collapse is affected by the element shape. As shown in the red and black boxes in Figs 12a and b, the crack direction is affected by the shape of the element. The cubic element causes the fragments to follow the vertical sea surface direction, and the use of triangular elements avoids the possibility of vertical crack propagation due to element mesh to a certain extent, and provides a more random propagation direction for crack propagation, but the cracks still grow in the same vertical direction as square elements in general.
3.3 Parameter sensitivity
The bottom friction and buoyancy are important factors that affect the process of glacier calving. In the following research, we conducted a series of simulations on bottom friction and buoyancy. If the parameters used in cases are not listed separately, they are consistent with the above experiment.
3.3.1 Bottom friction
In order to study the influence of bottom friction on the glacier disintegration process, we place a glacier on a slope with an inclination of 0.05 radians, and set a fixed boundary condition on the left side, as shown in Fig. 13. The whole glacier is divided into 6454 triangular elements, the influence of seawater buoyancy is not taken into consideration. The physical parameters in the calculation are the same as those in the above research, but the friction between the bottom of the glacier and the base is set to 0.2, 0.4, and 0.6. The bond strength between elements follows the Weibull distribution with a shape parameter of 2.0 and a size parameter of 4.0 as before, and the normal bond strength is between 1.0 MPa and 1.6 MPa in this case.
With different friction coefficients on the base, different forms of crack propagation and disintegration occur in glaciers. When the friction coefficient is small, such as 0.2 and 0.4, the glacier cracks expand, but the glacier stabilizes finally, as shown in Figs 14a and b. When the friction coefficient is taken as 0.6, the glacier will continue to retreat, as shown in Fig. 14c. By comparing the results of glacier disintegration under different bottom friction coefficients, we know that glaciers tend to disintegrate under larger friction coefficients, and smaller friction coefficients are more likely to cause cracks to growl while the glaciers remain stable. This is because the force provided by the base to the glacier through friction plays a certain role in maintaining the stability of the glacier.
The smaller friction force makes the stress inside the glacier along the advancing direction of the glacier more evenly distribute, and the glacier is more difficult to disintegrate. That is to say, when the friction is smaller, the bottom of the glacier front moves more distance under the action of gravity, and the tensile stress on the upper part of the glacier decreases. The larger coefficient of friction enables the bottom to provide greater friction. The advancement distance of the bottom of the glacier front becomes smaller, the upper part of the glacier is subject to increased tensile stress, and the bottom of the glacier front crush under greater friction and gravity, making it easier to disintegrate.
The propagation of cracks also validates this view. In the case of 0.2 (friction coefficient), cracks occur in the center of the glacier with a lower density, while the crack propagation in Fig. 14b gradually extends to the front edge of the glacier. When the crack propagates further, the part on the right side of the glacier in contact with the basement will not be able to withstand, causing the glacier to collapse and retreat.
3.3.2 Buoyancy
Based on the glacier geometry used in the above-mentioned study of bottom friction, we conducted a study on the influence of sea level position, that is, buoyancy on the calving process. The base friction coefficients used in this example are 0.6. The results are shown in Fig. 15. The sea level is set to 10 m, 20 m, and 30 m from the bottom of the glacier front respectively.
It can be seen from the results that the higher the water level, the less likely the glacier will collapse. In the case of Fig. 15a, the sea level is 10 m from the bottom of the glacier front, and the glacier collapses. In the case of Fig. 15b, the sea level is 20 m from the bottom of the glacier front, the glacier does not collapse, but a longitudinal crack that penetrates the glacier is produced in this case. In the case of Fig. 15c, the sea level is 30 m from the bottom of the glacier front, the glacier does not disintegrate, and the cracks cannot penetrate the entire glacier in the thickness direction.
4 Conclusions
Less
A practical DEM with dilated polyhedral elements for a glacial calving simulation has been developed in this study. The dilated polyhedron element has the potential for creating irregular, convex element shapes. It is also suitable for establishing seamless glacier structures, and the fracture process of glacier ice is simulated by using a hybrid fracture model based on critical energy. The von Mises stress and the normal stress between the elements are analyzed so that the relationship between the fracture of the glacier and the internal stress can be effectively combined to facilitate the observation. In our research, the glacier ice has broken due to excessive tensile stress, and it broke off from the top and spread to the bottom, then the entire glacier became unstable and disintegrated. During the calving process, the grinding between the ice blocks near the shear zone causes the glacier to break down into smaller pieces, while the foremost part of the glacier is broken into larger blocks, which is consistent with some situation of field observation. It indicates that the model is able to simulate a realistic tidewater glacier calving process in a short time, confirming the feasibility of the proposed model for calving simulations. By studying the glacier disintegration process of different element shapes, it is found that cubic elements are more likely to cause fractures along with the element division (vertical to the sea level), while the direction of fragmentation of triangular elements is more random, but in general, the element shape has little effect on the retreat distance of the glacier. The friction at the bottom of the glacier will affect the propagation of cracks and calving in the glacier, and in the case of this article, a large friction coefficient is more likely to cause the glacier failure. Buoyancy conditions have a significant impact on whether glaciers disintegrate, because greater buoyancy will make it more difficult for cracks to penetrate more glaciers. In further research, different glacier geometries and 3D conditions will be considered to explore the performance of calving laws used in ice sheet models.
Funding
Less
  • The National Key R&D Program of China under contract Nos 2018YFA0605902, 2016YFC1402705, 2016YFC1402706 and 2016YFC1401505; the National Natural Science Foundation of China under contract Nos 41576179 and 51639004; the fund of Australian Research Council’s Special Research Initiative for Antarctic Gateway Partnership under contract No. SR140300001; the China Postdoctoral Science Foundation under contract No. 2020M670746.
References
Less
Alley R B, Horgan H J, Joughin I, et al. 2008. A simple law for ice-shelf calving. Science, 322(5906): 1344, doi: 10.1126/science.1162543
Åström J A, Riikilä T I, Tallinen T, et al. 2013. A particle based simulation model for glacier dynamics. The Cryosphere, 7(5): 1591–1602, doi: 10.5194/tc-7-1591-2013
Åström J A, Vallot D, Schäfer M, et al. 2014. Termini of calving glaciers as self-organized critical systems. Nature Geoscience, 7(12): 874–878, doi: 10.1038/ngeo2290
Azevedo N M, Lemos J V. 2005. A generalized rigid particle contact model for fracture analysis. International Journal for Numerical and Analytical Methods in Geomechanics, 29(3): 269–285, doi: 10.1002/nag.414
Bartholomaus T C, Larsen C F, O’ Neel S. 2013. Does calving matter? Evidence for significant submarine melt. Earth and Planetary Science Letters, 380: 21–30, doi: 10.1016/j.jpgl.2013.08.014
Bassis J N, Jacobs S. 2013. Diverse calving patterns linked to glacier geometry. Nature Geoscience, 6(10): 833–836, doi: 10.1038/ngeo1887
Benn D I, Åström J A. 2018. Calving glaciers and ice shelves. Advances in Physics: X, 3(1): 1513819, doi: 10.1080/23746149.2018.1513819
Benn D I, Åström J, Zwinger T, et al. 2017. Melt-under-cutting and buoyancy-driven calving from tidewater glaciers: new insights from discrete element and continuum model simulations. Journal of Glaciology, 63(240): 691–702, doi: 10.1017/jog.2017.41
Benn D I, Warren C R, Mottram R H. 2007. Calving processes and the dynamics of calving glaciers. Earth-Science Reviews, 82(3–4): 143–179, doi: 10.1016/j.earscirev.2007.02.002
Burton J C, Amundson J M, Cassotto R, et al. 2018. Quantifying flow and stress in ice mélange, the world's largest granular material. Proceedings of the National Academy of Sciences of the United States of America, 115(20): 5105–5110, doi: 10.1073/pnas.1715136115
Camanho P P, Davila C G, de Moura M F. 2003. Numerical simulation of mixed-mode progressive delamination in composite materials. Journal of Composite Materials, 37(16): 1415–1438, doi: 10.1177/0021998303034505
Cleary P W. 2010. DEM prediction of industrial and geophysical particle flows. Particuology, 8(2): 106–118, doi: 10.1016/j.partic.2009.05.006
Cook S, Åström J, Zwinger T, et al. 2018. Modelled fracture and calving on the Totten Ice Shelf. The Cryosphere, 12(7): 2401–2411, doi: 10.5194/tc-12-2401-2018
Cundall P A, Strack O D L. 1979. A discrete numerical model for granular assemblies. Géotechnique, 29(1): 47–65
Dempsey J P, Adamson R M, Mulmule S V. 1999. Scale effects on the in-situ tensile strength and fracture of ice. Part II: First-year sea ice at Resolute, NWT. International Journal of Fracture, 95(1): 347–366, doi: 10.1023/A:1018650303385
Dyer M S, Collins C, Hodgeman D, et al. 2013. Computationally assisted identification of functional inorganic materials. Science, 340(6134): 847–852, doi: 10.1126/science.1226558
Golledge N R, Keller E D, Gomez N, et al. 2019. Global environmental consequences of twenty-first-century ice-sheet melt. Nature, 566(7742): 65–72, doi: 10.1038/s41586-019-0889-9
Gürtner A, Bjerkås M, Forsberg J, et al. 2010. Numerical modelling of a full scale ice event. In: 20th IAHR International Symposium on Ice. Lahti, Finland: IAHR
Koehn D, Sachau T. 2014. Two-dimensional numerical modeling of fracturing and shear band development in glacier fronts. Journal of Structural Geology, 61: 133–142, doi: 10.1016/j.jsg.2012.11.002
Liu Lu, Ji Shunying. 2018. Ice load on floating structure simulated with dilated polyhedral discrete element method in broken ice field. Applied Ocean Research, 75: 53–65, doi: 10.1016/j.apor.2018.02.022
Liu Lu, Ji Shunying. 2019. Bond and fracture model in dilated polyhedral DEM and its application to simulate breakage of brittle materials. Granular Matter, 21(3): 41, doi: 10.1007/s10035-019-0896-4
Lu W, Lubbad R, Løset S. 2014. Simulating ice-sloping structure interactions with the cohesive element method. Journal of Offshore Mechanics & Arctic Engineering, 136(3): 16, doi: 10.1115/1.4026959
Lu G, Third J R, Müller C R. 2015. Discrete element models for non-spherical particle systems: from theoretical developments to applications. Chemical Engineering Science, 127: 425–465, doi: 10.1016/j.ces.2014.11.050
Ma Gang, Zhou Wei, Chang Xiaolin, et al. 2016. A hybrid approach for modeling of breakable granular materials using combined finite-discrete element method. Granular Matter, 18: 7, doi: 10.1007/s10035-016-0615-3
Ma Yue, Tripathy C S, Bassis J N. 2017. Bounds on the calving cliff height of marine terminating glaciers. Geophysical Research Letters, 44(3): 1369–1375, doi: 10.1002/2016GL071560
Mercenier R, Lüthi M P, Vieli A. 2018. Calving relation for tidewater glaciers based on detailed stress field analysis. The Cryosphere, 12(2): 721–739, doi: 10.5194/tc-12-721-2018
Paavilainen J, Tuhkuri J. 2012. Parameter effects on simulated ice rubbling forces on a wide sloping structure. Cold Regions Science and Technology, 81: 1–10, doi: 10.1016/j.coldregions.2012.04.005
Paavilainen J, Tuhkuri J, Polojärvi A. 2011. 2D numerical simulations of ice rubble formation process against an inclined structure. Cold Regions Science and Technology, 68(1−2): 20–34, doi: 10.1016/j.coldregions.2011.05.003
Pelto M S, Warren C R. 1991. Relationship between tidewater glacier calving velocity and water depth at the calving front. Annals of Glaciology, 15: 115–118, doi: 10.3189/S0260305500009617
Pralong A, Funk M. 2005. Dynamic damage model of crevasse opening and application to glacier calving. Journal of Geophysical Research: Solid Earth, 110(B1): B01309, doi: 10.1029/2004JB003104
Riikilä T I, Tallinen T, Åström J, et al. 2015. A discrete-element model for viscoelastic deformation and fracture of glacial ice. Computer Physics Communications, 195: 14–22, doi: 10.1016/j.cpc.2015.04.009
Robel A A. 2017. Thinning sea ice weakens buttressing force of iceberg mélange and promotes calving. Nature Communications, 8(1): 14596, doi: 10.1038/ncomms14596
Rothrock D A, Thorndike A S. 1984. Measuring the sea ice floe size distribution. Journal of Geophysical Research: Oceans, 89(C4): 6477–6486, doi: 10.1029/JC089iC04p06477
Savage S B, Crocker G B, Sayed M, et al. 2000. Size distributions of small ice pieces calved from icebergs. Cold Regions Science and Technology, 31(2): 163–172, doi: 10.1016/S0165-232X(00)00010-0
Schlemm T, Levermann A. 2018. A simple stress-based cliff-calving law. The Cryosphere, 13: 2475–2488, doi: 10.5194/tc-13-2475-2019
Todd J, Christoffersen P, Zwinger T, et al. 2018. A full-stokes 3-D calving model applied to a large greenlandic glacier. Journal of Geophysical Research: Earth Surface, 123(3): 410–432, doi: 10.1002/2017JF004349
van Dongen E C H, Kirchner N, van Gijzen M B, et al. 2018. Dynamically coupling full Stokes and shallow shelf approximation for marine ice sheet flow using Elmer/Ice (v8.3). Geoscientific Model Development, 11(11): 4563–4576, doi: 10.5194/gmd-11-4563-2018
Zemp M, Huss M, Thibert E, et al. 2019. Global glacier mass changes and their contributions to sea-level rise from 1961 to 2016. Nature, 568(7752): 382–386, doi: 10.1038/s41586-019-1071-0
Appendix
Less
Year 2021 volume 40 Issue 7
PDF
39
22
Cite this Article
BibTeX
Article Info
doi: 10.1007/s13131-021-1819-x
  • Receive Date：2020-12-19
  • Online Date：2026-03-03
  • Published：2021-07-25
Article Data
Affiliations
History
  • Received：2020-12-19
  • Accepted：2021-01-28
Funding
The National Key R&D Program of China under contract Nos 2018YFA0605902, 2016YFC1402705, 2016YFC1402706 and 2016YFC1401505; the National Natural Science Foundation of China under contract Nos 41576179 and 51639004; the fund of Australian Research Council’s Special Research Initiative for Antarctic Gateway Partnership under contract No. SR140300001; the China Postdoctoral Science Foundation under contract No. 2020M670746.
Affiliations
    1 State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology, Dalian 116024, China
    2 National Marine Environmental Forecasting Center, Beijing 100081, China
    3 Institute for Marine and Antarctic Studies, University of Tasmania, Hobart 7001, Australia

Corresponding:

References
Share
https://castjournals.cast.org.cn/joweb/aos/EN/10.1007/s13131-021-1819-x
Share to
QR

Scan QR to access full text

Cite this article
BibTeX
Citations
表12种不同金属材料的力学参数
科
Family		 属数
Number of
genus		 种数
Number of
species		 占总种数比例
Percentage of
total species (%)		 属
Genus		 种数
Number of
species		 占总种数比例
Percentage of total
species (%)
鹅膏菌科Amanitaceae 2 11 5.26 鹅膏菌属 Amanita 10 4.78
小菇科 Mycenaceae 2 12 5.74 丝盖伞属 Inocybe 5 2.39
多孔菌科 Polyporaceae 8 14 6.70 蜡蘑属 Laccaria 5 2.39
红菇科 Russulaceae 3 23 11.00 小皮伞属 Marasmius 6 2.87
小菇属 Mycena 11 5.26
光柄菇属 Pluteus 5 2.39
红菇属 Russula 17 8.13
栓菌属 Trametes 5 2.39
关闭全屏
  • BibTeX
  • EndNote
  • RefWorks
  • TxT
Links
Articles
Latest Articles
Most Read
Collections
Updates
Events
News
Multimedia
About
About Us
Contact
No. 86 Xueyuan South Road, Haidian District, Beijing
100081
010-62199257
qkjq@cast.org.cn

Copyright © 2025 China Association for Science and Technology. All rights reserved. For all open access content, the relevant licensing terms apply.
Sponsored by the Office of the Leading Group for Cybersecurity and Informatization of CAST, and supported by Science and Technology Review Publishing House