Subcritical crack propagation in glacial quarrying during subglacial water pressure variation

Abstract

Glacial quarrying stems from the fracturing of subglacial bedrock. Much evidence shows that subcritical crack propagation of bedrock is closely related to subglacial water pressure fluctuations. Here we employ a model that assesses the impact of subglacial water pressure fluctuation on cavity length and subcritical crack propagation, while analyzing the effect of a pre-existing crack location using a phase-field model (PFM). Our results indicate that the cavity length is reduced during diurnal fluctuations in water pressure. There are two patterns of subcritical crack propagation on the corner of the step. The first stems from a rapid drop in water pressure. The second occurs after the water pressure recovers from the fluctuation to the initial steady state. This pattern is a consequence of enhanced stress concentration on the step since the modeled cavity length exceeds its steady value and has higher efficiency in promoting subcritical crack propagation. Additionally, based on the PFM results, we speculate that the subcritical crack initiation and propagation happen on a broader scale, including the ice-bed contact region and its adjacent region. Our findings imply that the duration of subcritical crack propagation is short and typically ceases once the cavity length adjusts to reduced water pressure levels.

1. Introduction

Subglacial erosion causes mountain surfaces to be worn down and produces vast amounts of sediments that are transported into rivers and oceans, changing the biogeochemical balance (Wadham and others, 2019; Herman and others, 2021). Subglacial erosion is strongly linked to landscape evolution (Egholm and others, 2009). Two processes govern subglacial erosion rates: glacial quarrying, which involves the dislodging of large fragments from the subglacial bed, and glacial abrasion, defined as the wear of rock fragments against the subglacial bed. Glacial quarrying is recognized as a dominant process (Hallet, 1996; Hildes and others, 2004; Loso and others, 2004; Riihimaki and others, 2005). Field observations reveal that glacial quarrying was significant with a pre-existing crack or joint set (Rea and Whalley, 1996; Krabbendam and Glasser, 2011; Hooyer and others, 2012; Glasser and others, 2020). Additionally, it is suggested that quarrying rates increase with subglacial water pressure fluctuations (Iverson, 1991; Alley and others, 1999; Cohen and others, 2006).

In the past three decades, a series of studies of glacial quarrying has further supported the view that quarrying stems from fracturing the subglacial bed adjacent to water-filled cavities, which expand or contract, depending on whether the water pressure is high or low (Iverson, 1991, 2012; Hallet, 1996; Rea and Whalley, 1996; Cohen and others, 2006; Anderson, 2014). However, direct observations of glacial quarrying are limited due to the difficulty of digging tunnels into the glacier base to study glacial erosion. Some studies analyze fracture growth on a subglacial bed through scaled drone-based photos and finite element models (Woodard and others, 2019).

In general, several quarrying models describe the relationship between the geometry of subglacial bed, subcritical crack growth, cavity length and erosion rates (Iverson, 1991, 2012; Hallet, 1996; Hildes and others, 2004). These quarrying models estimate the basal stresses on the corner of subglacial bedrock steps. As the water level falls, the basal stresses increase with the water pressure decrease in the cavities, and then the subcritical crack growth is estimated by Linear Elastic Fracture Mechanics (LEFM) (Atkinson, 1984; Iverson, 1991; Hallet, 1996). The crack growth determines the subglacial erosion rates in these quarrying models.

As far as glacial quarrying is concerned, the bedrock geometry is commonly chosen (Iverson, 1991, 2012; Hallet, 1996; Rea and Whalley, 1996; Cohen and others, 2006; Hooyer and others, 2012; Anderson, 2014), as shown in Figure 1. There are two models to describe the cavity roof, Kamb's model and Nye's theory. Kamb's model supposes that the ice rheology is linear and the slope of the cavity roof is small (Kamb, 1987); the cavity height is calculated based on the cavity length. In comparison, Nye's theory is based on the closure rate of the cylindrical cavity in infinite glacier ice (Nye, 1953). According to Kamb's model and Nye's theory, the cavity roof geometry is obtained in a steady state. Nevertheless, the glacier is only sometimes steady and changes every moment. Water level records from boreholes show dynamic water pressure character, fluctuating on sub-daily scales during the melt season (Hooke, 1991; Fudge and others, 2008; Sugiyama and others, 2019). This character inevitably changes cavity geometry and influences subcritical crack propagation. Up to now, there have been few systematic analyses of the intrinsic connection between the diurnal fluctuation of water pressure, cavity geometry and subcritical crack propagation, which is the key to enhancing the knowledge of glacial quarrying (Iverson, 1991).

Fig. 1. The cavity profile is assumed to be a part of a circle of radius R. OA and OB are equal in length to R. The velocities of all the points of curve AB are decomposed into ice flow sliding velocity U _s and cavity closure rate U _c. L(t) is the horizontal displacement of these points, and h(t) is the vertical displacement. Their velocities are L′(t) and h′(t). A pre-existing crack extends 0.1 m into the bedrock. The step height is one meter, and the tread length is 10 m. The melting of the cavity roof is not included in the computational model. In the steady state, points on the cavity roof have a velocity that is tangential to the roof. But in the dynamic state, these points don't move on their tangent due to the deformation of the cavity roof.

This study presents a model that assesses the impact of subglacial water pressure fluctuations on cavity length and subcritical crack propagation. The study is organized as follows. Section 2 gives the technical parameters, calculating method, numerical implementation and phase-field method (PFM) simulation, which involves the cavity length change, the subcritical crack propagation and the effect of a pre-existing crack location. Section 3 shows the results, and section 4 carries on the discussion and model limitations. Finally, the research conclusion is given with recommendations in section 5.

2. Technical parameters and model

The normalized cavity extent and effective stress determine the far-field tensile stress that drives subcritical crack propagation. In order to investigate this process, we approximate cavity geometry under steady state conditions, where the cavity radius is defined by a unique function of cavity length and step height that is not sensitive to other details of the cavity height profile. Next, we examine the impact of dynamic water pressure fluctuations on cavity geometry using a time-varying creep closure rate that captures the effects of changes in effective stress. Based on cavity geometry characteristics, we obtain the ice normal stress acting on the rock step by cavity length and then use LEFM to compute subcritical crack propagation. In a dynamic state, the range of ice flow sliding velocity is from 100 to 500 m yr⁻¹. The amplitude of water level fluctuations is from 50 to 100 m. The computation of steady and dynamic states is implemented by using the FORTRAN program. Finally, we use PFM to study the effect of pre-existing crack location in the COMSOL Multiphysics software.

2.1. Steady state

The calculation of cavity length is grounded in Nye's theory (Nye, 1953). It yields the cavity closure rate,

(1)$$U_{\rm c} = \left({\displaystyle{{P_{\rm e}} \over {nB}}} \right)^nR, \;$$

where U _c is the cavity closure rate; P _e is the effective pressure (ice-overburden pressure P _i minus water pressure P _w); B and n are the ice flow-law parameters in Glen's flow law (Glen, 1955). Finally, R is the radius of the cavity.

The parameter value B is based on the field data collected at the bed of Engabreen in Norway (Cohen, 2000). The value of n is from laboratory and field data. In this study, we choose n = 3 for Eqn (1) according to the previous studies (Cuffey and Paterson, 2010). There are two different ways of estimating the value of the cavity radius R: one way is to consider R as the step height (Walder, 1986; Schoof, 2010), and the other is to consider R as the cavity length (Humphrey, 1987; Iverson and Petersen, 2011). In order to reflect the geometry features of cavities more comprehensively, an alternative way is developed (Petersen, 2012), as shown in Figure 1. At a steady state (Fig. 1a), The cavity profile is assumed to be a part of a circle of radius R. The circle center is point O. A is where ice separates from bedrock, and B is where ice reconnects. OA and OB are equal in length to R.

Theoretically, as glaciers move, the velocities of all the points of curve AB are decomposed into two parts: ice flow sliding velocity U_s in the horizontal direction and the approximately considered cavity closure rate U_c in the direction toward the center of the circle. As a result, we get two equations for the deformation of the cavity in the horizontal and vertical directions:

(2)$$\mathop \smallint \limits_0^t \left({U_{\rm s}-\displaystyle{{L( x ) } \over R}U_{\rm c}} \right){\rm d}x = L( t ) , \;$$

(3)$$\mathop \smallint \limits_0^t \displaystyle{{\sqrt {R^2-L^2( x ) } } \over R}U_{\rm c}{\rm d}x = h( t ) , \;$$

where t is the time interval in which an arbitrary point of the curve AB moves from point A to the current location, L(t) is the horizontal displacement of these points, and h(t) is their vertical displacement. By taking the derivative of Eqns (2) and (3) with respect to the time interval t, we obtain two differential equations, Eqns (4) and (5).

(4)$$U_{\rm s}-\displaystyle{{L( t ) } \over R}U_{\rm c} = {L}^{\prime}( t ) , \;$$

(5)$$\displaystyle{{\sqrt {R^2-L^2( t ) } } \over R}U_{\rm c} = {h}^{\prime}( t ).$$

The Euler method is chosen for the differential equations above, and then we approximate L ^′(t _i) resulting in (L(t _i+1) − L(t _i))/(t _i+1 − t _i) and h ^′(t _i) resulting in (h(t _i+1) − h(t _i))/(t _i+1 − t _i).

(6)$$L( {t_{{\rm i} + 1}} ) = L( {t_{\rm i}} ) + \left[{U_{\rm s}-\displaystyle{{L( {t_{\rm i}} ) } \over R}U_{\rm c}} \right]( {t_{{\rm i} + 1}-t_{\rm i}} ) , \;$$

(7)$$h( {t_{{\rm i} + 1}} ) = h( {t_{\rm i}} ) + \left[{\displaystyle{{\sqrt {R^2-L^2( {t_{\rm i}} ) } } \over R}U_{\rm c}} \right]( {t_{{\rm i} + 1}-t_{\rm i}} ).$$

The value of time interval (t _i+1 − t _i) is the size of every step, yet the radius R is undetermined. In Eqn (7), the value of h(t) increases over time; when h(t _k) = AC in Figure 1, it is equivalent to a particle moving from point A to point B during the time of k × (t _i+1 − t _i), and a corresponding value of L(t _k) can be obtained as cavity length; then we have a formula about R from Figure 1:

(8)$$R = \displaystyle{{h^2( {t_{\rm k}} ) + L^2( {t_{\rm k}} ) } \over {2\;h( {t_{\rm k}} ) }}, \;$$

where h(t _k) is the step height and L(t _k) is the cavity length. Furthermore, it is required that L(t _k) > h(t _k) ; otherwise, the circle's center is on AC.

At last, based on Eqns (6–8), we have a cavity geometry in a steady state. The calculation procedure is presented in Table 1.

Table 1. The basic steps of the algorithm in the steady state (An example in the supplementary material)

2.2. Dynamic state

Previous research finds that diurnal borehole water level fluctuations are remarkable characteristics during the melt season, whereas water level is high and stable for the rest of the year (Hooke, 1991; Fudge and others, 2008; Sugiyama and others, 2019).

This study divides the water level into two parts: dynamic level and stable level. We apply a cosine function model to calculate effective pressure during water level fluctuations, as shown in Eqns (9) and (10):

(9)$$P_{\rm e}^{\rm ^{\prime}} = P_{\rm i}-P_{\rm w}^{\rm ^{\prime}} , \;$$

(10)$$P_{\rm e} = P_{\rm e}^{\rm ^{\prime}} + \left[{0.5-\;0.5\cos \left({\displaystyle{{2\;\pi t{\rm^{\prime}}} \over {3600\;\times \;24}}} \right)} \right]w\rho g, \;$$

where t′ is the time span as the fluctuation of water level takes place, as can be seen in Figure 3b, w is the amplitude of water level fluctuations, ρ is the density of water, g is the gravitational acceleration, P _i is the ice-overburden pressure, $P_{\rm w}^{\rm ^{\prime}}$ is the stable water pressure, $P_{\rm e}^{\rm ^{\prime}}$ is the effective pressure in stable water level, P _e is the effective pressure in dynamic water level.

In Nye's theory, Eqn (1) is derived based on Glen's experiments about the steady state creep of ice (Nye, 1953; Glen, 1955). Therefore it is a reasonable approximation that different effective pressure P _e corresponds to different cavity closure rates U _c. That means the fluctuation of the effective pressure P _e can cause the fluctuation of the cavity closure rate U _c by Eqns (1) and (10). Further, this change of U _c influences the value of the cavity length in Figure 1, thereby changing the value of the cavity radius R by Eqn (8). Combined with the steady state analysis method above, the cavity geometry can be presented in a dynamic state by Eqns (1–10). Figure 1 shows a schematic diagram of cavity geometry in a steady and dynamic state. The calculation procedure is presented in Table 2.

Table 2. The basic steps of the algorithm in the dynamic state (An example in the supplementary material)

2.3. Subcritical crack propagation

Subcritical crack propagation in rocks is directly affected by the surrounding environment. According to fracture mechanics, a crack occurs when the applied force breaks the atomic bonds of rocks; that is, effective pressure P _e induces subcritical crack propagation. In this study, crack propagation velocities are calculated by the quarrying model of Hallet (Hallet, 1996), based on LEFM (Atkinson, 1984). Suppose there is a pre-existing crack in the rock step. Then subcritical crack propagation can be described by the following equations:

(11)$$V = V_{\rm I}[ {{\rm e}^{\gamma ( {( {K_{\rm I}^2 /K_{\rm c}^2 } ) -\;1} ) }-{\rm e}^{\gamma ( {-\;8/9} ) }} ] \;, \;\quad K_{\rm I} > \;\displaystyle{1 \over 3}\;K_{\rm c}, \;$$

$$V = 0, \;\quad K_{\rm I} \le \displaystyle{1 \over 3}K_{\rm c}, \;$$

(12)$$K_{\rm I} = \sigma _{\rm d}\sqrt {\displaystyle{{4l_{\rm c}} \over \pi }} , \;$$

(13)$$\sigma _{\rm d} = \displaystyle{2 \over 3}( {\sigma_{\rm n}-P_{\rm w}} ) = \displaystyle{2 \over 3}P_{\rm e}\left({\displaystyle{1 \over {1-{S}^{\prime}}}} \right), \;\quad \sigma _{\rm n}-P_{\rm w} < \sigma _{\rm n}^\ast , \;$$

where V is the subcritical crack propagation velocity, V _I and γ are the growth-law parameters of rock, K _c is the fracture toughness of rock, K _I is the stress intensity factor of rock, σ _d is the far-field tensile stress, σ _n is the ice normal stress acting on the rock step, P _w is the water pressure, (σ _n − P _w) is the deviatoric stress, $\sigma _{\rm n}^\ast$ is the finite strength of the ice, S′ is the normalized extent of the subglacial cavity which is equal to cavity length divided by tread length in Figure 1, l _c is the crack length, the (1/3) K _c is the stress-corrosion limit which determines whether the crack propagation occurs or not.

The effective pressure P _e and the cavity geometry are consistent when the system is in a steady state. Therefore, according to Eqns (11–13), if K _I > (1/3) K _c, crack propagation occurs and does not stop, whereas crack propagation never happens if K _I ≤ (1/3) K _c.

In dynamic state, water level fluctuations change the effective pressure P _e and the cavity geometry. These changes influence the subcritical crack propagation velocity V and crack propagation increases the crack length l _c. In return, a more extended crack length l _c also induces a higher subcritical crack propagation velocity V. However, if the value of the effective pressure P _e is low enough, crack propagation stops.

2.4. Experiment design

2.4.1. Numerical implementation

Figure 1 illustrates ice sliding over the subglacial bed, with a pre-existing crack extending 0.1 m into the bedrock. We considered two scenarios. In the first scenario, it was assumed that the water level was in a steady state at first and then began to fluctuate with an amplitude w = 100 m, continuing for several months. In the second scenario, it was also assumed that the water level was steady initially and then began to fluctuate. However, this fluctuation continued for just 12 days. After that time, the water level recovered to the initial level in a steady state. In the two scenarios, the initial effective pressure P _e was considered 0.4 MPa, the ice flow sliding velocity was U _s, and the amplitude of water level fluctuations was from 50 to 100 m. The other parameters of computation are given in Table 3. Finally, the solution algorithm's basic steps are presented in Tables 1 and 2.

Table 3. Parameter used for the model

To investigate the effect of water pressure fluctuations on cavity geometry and subcritical crack propagation, we mainly present three types of graphs. The first is about cavity length evolution with the fluctuation of water pressure. The second is about cavity roof shape, which influenced cavity length, and the third is about subcritical crack propagation.

2.4.2. Phase-field model (PFM) simulation

Based on assessing the impact of subglacial water pressure fluctuations on cavity length and subcritical crack propagation, we study the effect of a pre-existing crack location using the PFM.

The PFM is considered a convenient method to predict crack initiation and propagation (Bourding and others, 2011; Miehe and others, 2015). In the PFM, the sharp crack topology is described as a phase-field parameter:

(14)$$\phi ( x ) = \left\{{\matrix{ {1, \;} \hfill & {if\;material\;is\;cracked} \hfill \cr {0, \;} \hfill & {if\;material\;is\;intact} \hfill \cr } } \right..$$

A dimensional phase-field is given approximately by a standard exponential function:

(15)$$\phi ( x ) = {\rm e}^{-\vert x \vert /l_0}, \;$$

where l ₀ is the length scale parameter about the crack, which indirectly reflects the width of the crack. In this study, we implemented the PFM in the COMSOL Multiphysics software. The schematic diagram of the model can be seen in Figure 2. The model's dimensions equal 8 m in width and 2 m in height. In order to simulate the actual condition, the normal stress σ _n is calculated by Eqn (13). Basal shear stress τ _b has a linear dependence on debris concentration (Cohen and others, 2005). The values of debris-bed friction do not typically exceed 0.1 MPa (Cohen and others, 2005; Cuffey and Paterson, 2010). This paper assumes that basal shear stress τ _b is 0.1 MPa.

Fig. 2. Geometry and boundary condition of PFM.

In the COMSOL simulation, we chose six pre-existing cracks along with the bedrock. Granite's material parameters can be seen in Table 3, including elasticity modulus E, Poisson ratio v, length scale parameter l ₀ and critical energy release rate G _c. We set a high-density mesh in the vicinity of the propagating crack in order to properly resolve the phase field. The maximum mesh element is 0.3 mm. In addition, we selected the model's bottom as a fixed constraint. The boundary loads can be seen in Figure 2. There are no constraints and no loads on the other boundaries. To improve the convergence, we set up a segregated solver sequence. The number of iterations is set to be 3 to meet computation accuracy.

With the same stress and critical energy release rate, these cracks' propagation can be compared to each other and give us valuable information about the quarrying process.

3. Results

Figure 3a shows the relationship between the cavity length and the ice flow sliding velocity during subglacial water pressure variations. Figure 3b shows the corresponding change at the water level. We suppose the initial effective pressure P _e is 0.4 MPa. It is approximately equal to an ice thickness of 404 m with a water level of 330 m. At first, the cavity length remains constant at a stable water level. Then, after the 10th day of the simulated calculation, the water level begins to have a diurnal fluctuation with an amplitude of 100 m, which causes a change in the cavity length. As shown in Figure 3a, the cavity length is much less during the water level fluctuation, and higher ice velocities lead to greater cavity lengths. Furthermore, changes in cavity length and water level fluctuations are not synchronized in Figure 5, the former lags behind the latter.

Fig. 3. The relationship between cavity length, ice flow sliding velocity and water level.

Figure 4 shows the cavity roof shape on the 11th, the 12th and the 20th day corresponding to the water levels shown in Figure 3b. The cavity length significantly decreases on the 11th day when the water level begins to decline in Figures 4a–4c. On the 12th day, the cavity length is generally stable in Figures 4d–4f. Combining Figures 3a and 4g–4i, we can see that diurnal changes in cavity length are the same every day after the 20th day.

Fig. 4. Cavity roof shape on the 11th, 12th and 20th day with the water level in Figure 3b (U _s = 500 m yr⁻¹).

Fig. 5. The lag effect between the cavity length change and the water level fluctuation (U _s = 500 m yr⁻¹).

Figure 6a shows cavity length changes with water level going from a steady state to a fluctuating state, then back to a steady state. The water level is in a steady state at first and then begins to fluctuate with an amplitude of 100 m. This fluctuation continues for just 12 days. After that time, the water level recovers to the initial level in a steady state. Comparing different states in Figure 6a, we can find that the cavity length is longer after recovering the water level. Figure 8a shows the cavity closure rate change based on water level fluctuation in Figure 8c, the ice flow sliding velocity is 500 m yr⁻¹, and the change of the cavity radius is in Figure 8b. In Figures 7a–7i, we can see the change in the cavity roof shape; the corresponding water level fluctuation can be seen in Figure 6b. At first, the water level is steady, and the cavity roof shape is an approximate triangle in Figure 7a. When the water level begins to fluctuate, the cavity is compressed, and its roof shape is a wavy pattern because of the diurnal change of the water level, as shown in Figures 7b–7d. After the fluctuation of 12 days, the water level recovers to the initial steady state. However, during the recovery process, with the water pressure rising, the front part of the cavity roof shape remains much plumper than that in the steady state. This significant difference can be seen between Figures 7a and 7f. When the ice flow moves from left to right, the plump front part causes a gradual increase in the cavity length, which exceeds that in the steady state on the 28th day. The fundamental reason for the overshoot in the cavity length is the closure rate in Figure 8a. After the water level recovers, from the 23rd to the 26th day, the cavity closure rate is continuously lower than that in the first ten days. In Figure 1, the vertical velocities of the points of curve AB are components of the cavity closure rate. This indicates that the cavity roof's contact with bedrock is delayed by the reduced vertical velocities. However, the cavity length is too long to maintain this unique cavity roof shape for a long time because the increased cavity length leads to a higher cavity closure rate in Eqn (1). In the end, the cavity length returns to that in the steady state. The evolutionary processes above can be seen in Figures 7a–7i.

Fig. 6. The relationship between cavity length, ice flow sliding velocity and water level.

Fig. 7. Cavity roof shape with the water level in Figure 6b (U _s = 500 m yr⁻¹).

Fig. 8. The cavity closure rate with water level (The dashed line is the average cavity closure rate, U _s = 500 m yr⁻¹).

Using Hallet's model (Hallet, 1996) and LEFM, we calculated the subcritical crack propagation velocity, described in Figure 1. In the first ten days, there is no subcritical crack propagation in the steady state due to the low effective pressure. With the water level fluctuation, two subcritical crack propagation patterns are presented. The first occurs when the water level falls in the first fluctuation cycle on the 11th day in Figure 9a. At that time, the effective pressure is higher than that of the steady state. Figures 9a and 9b show that the first pattern of subcritical fracture propagation coincides with the initial decrease in water level. When the water level drops for the first time on the 11th day, the effective pressure and the ice-bed contact area increase simultaneously. However, the more significant effective pressure prompts subcritical crack propagation, while the greater ice-bed contact area produces opposite effects. Only in the first drop in water level, the effective pressure increases more. Meanwhile, the ice-bed contact area increases a little. After the 12th day, the effective pressure increases periodically. However the shorter cavity length leads to greater ice-bed contact area in Figure 3a. This is why the first pattern of subcritical crack propagation happens in the first drop of water level on the 11th day.

Fig. 9. Subcritical crack propagation with water level fluctuation. The first pattern of crack propagation occurs when the water level falls in the first fluctuation cycle.

The second pattern of subcritical crack propagation occurs after the water level recovers from the fluctuation to the initial steady state on the 28th day in Figure 10a. In this pattern, the effective pressure is the same as that of the steady state, while the contact area between ice flow and bedrock is smaller. The reason is due to the excessive length of the cavity, which exceeds that in the steady state, as shown in Figures 7a and 7g. The excessive length derives from the lower cavity closure in Figure 8a, which is described above. Because of the smaller contact area, the subsequent growing stress intensity factor triggers subcritical crack propagation. As illustrated in Figures 6, 7g and 10a, the water level fluctuation ceases on the 22nd day, but the cavity length continues to grow until the 28th day, delaying the spread of the subcritical crack.

Fig. 10. Subcritical crack propagation with water level fluctuation.

In Figure 10, the water level amplitudes 60, 70 and 80 m have their maximum crack extension length. Nevertheless, when the water level amplitudes are 90 and 100 m, according to Eqns (11–13), higher ice normal stress σ _n leads to a larger intensity factor K _I. Subsequently, the larger factor K _I increases the subcritical crack propagation velocity V and the crack length l _c. In turn, the crack length l _c promotes the larger intensity factor K _I. At last, the positive feedback leads to unstable crack propagation. We note that there is only one crack propagation pattern in Figure 9a, while two patterns are shown in Figure 10a. This is because, with the continuation of diurnal fluctuations, the cavity length never returns to recover to the pre-fluctuation level in Figure 3a, and the subcritical crack propagation of the first patterns stops growing in Figures 9a and 10a.

After finishing the calculation of the cavity length change with water level fluctuation, we chose the geometry of Figure 7g as our simulation object in the PFM because the pre-existing crack most likely extends in that environment. For the geometry shown in Figure 7g, the normal stress σ _n is 6.38 MPa, the basal shear stress τ _b is 0.1 MPa, and the length of the ice-bed contact region is 1.3 m. Through the COMSOL Multiphysics software, we obtained the evolution paths of six pre-existing cracks. Figure 11 shows the schematic diagram of PFM. Figure 12 shows the difference in the evolution paths of the six pre-existing cracks. When the pre-existing crack is in the ice-bed contact region, as seen in Figures 12e and 12f, crack propagation does not occur. However, if the pre-existing crack is on the cavity side of the ice-bed contact region or its edge in Figures 12b–12d, the crack grows, while the extension length is negatively correlated with the distance of the pre-existing crack from the ice-bed contact region.

Fig. 11. Crack paths in PFM (The figure is a portion of Fig. 12d).

Fig. 12. Crack paths in PFM. From A to F, The distance between the pre-existing cracks and the rock step edge is 3, 2.5, 1.9, 1.3, 0.8 and 0.3 m.

4. Discussion and model limitations

4.1 Model limitations

In this study, we consider one meter of step height and 10 m of tread length. The outcomes of crack propagation will differ depending on how the subglacial bedrock step is shaped. For example, a shorter tread length leads the cavity to span to the next step, and effective pressure would change correspondingly. Therefore, further studies should focus on the different kinds of subglacial bedrock steps and their crack propagation. Moreover, this study is based on an approximation that different effective pressure P _e corresponds to different cavity closure rates U _c. If there is a time delay between P _e and U _c, this asynchrony would change the results. Therefore, we recommend conducting additional research to test the relationship of effective pressure P _e and cavity closure rate U _c with water level fluctuation, which is essential to evaluate crack propagation.

In the model of Nye, cavity roof is arc shaped when water pressure is constant. Our study shows a wave-like shape because of subglacial water pressure fluctuation in the cavity roof. Although a change in water pressure can affect the shape of the cavity roof, there is currently no evidence to support this wave-like pattern. It is necessary to investigate the subglacial cavity during the melt season, in which water pressure changes daily. In addition, the cavity radius R is calculated based on step height and cavity length in Eqn (8). However, the cavity roof is not a symmetrical circular borehole in a dynamic state, as shown in Figures 4 and 7. This means the cavity radius R is just an approximation of the cavity geometry. Another uncertainty is the cavity closure rate of the points near A and B in Figure 1. We assume that all points of the cavity roof have the same closure rate according to Eqn (1). As the cavity roof evolves in the model, there is no consideration of mass or mass balance. The cavity roof can freely ‘sink’ beneath the bedrock level: at this stage, the points on the cavity roof beneath the bedrock are disregarded and not used in further computations. Therefore, more accurate models should be proposed in the future.

The PFM in the study gives us a valuable reference on the effect of a pre-existing crack location. We chose the cavity length and stress in which the pre-existing crack most likely expands after the water level recovers from the fluctuation to the initial steady state. The normal stress σ _n is 6.38 MPa, and the basal shear stress τ _b is 0.1 MPa at that time. Nevertheless, these stresses are insufficient to cause crack initiation or propagation in the PFM if the critical energy release rate G _c is too high. Hence we chose a relatively low value. Figures 12b–12d indicate that the crack will probably spread in this area. Currently, no studies have been done on determining the reasonable value of the critical energy release rate G _c in subcritical crack propagation.

4.2 Cavity length

In this study, we chose Nye's theory. Kamb's model also applies to the dynamic state with water level fluctuation (Kamb, 1987). However, in Kamb's model, cavity roof shape is determined by cavity length. It means that the same cavity length leads to the same cavity roof shape at any time. Yet, in Nye's theory, when the water level fluctuates, the cavity roof shape changes with time, even with the same cavity length. This suggests that Nye's theory is more suited to describe the dynamic state.

Through numerical modeling, the study presents the details of cavity length change and subcritical crack propagation during subglacial water pressure variations. This result is consistent with field observations. Research has also found that cavity length changes lag behind the water level fluctuation in Figure 5. The lag effect stems from the change in the cavity closure rate. In Eqn (1), two parameters determine the cavity closure rate, the effective pressure P _e and the cavity radius R. When the water level falls, or rather, the effective pressure P _e begins to rise, it should trigger the cavity to close. Nevertheless, the cavity length or radius R is still small, so the cavity closure rate is not large enough to compress this cavity, according to Eqn (1). In other words, cavity length keeps increasing until there is a reversal, as shown in Figure 5. When the water level rises, the same lag effect can be seen between the change in cavity length and the water level fluctuation.

4.3 Subcritical crack propagation

Based on a dynamic simulation of water level fluctuation, this research identified two crack propagation patterns. As shown in Figure 10, the second one has higher efficiency in promoting subcritical crack propagation. The mechanism of two crack propagation patterns means glacial quarrying happens within a short time. Most times in a year, it stops. This study assumes only a short melt season with diurnal water level fluctuation. However, if there are several short melt seasons at intervals of a few days, crack propagation will occur much more often. This implies that glacial hydrology has a more complex impact on quarrying rates.

Previous studies mention the effect of water pressure within cracks (Iverson, 1991; Hooyer and others, 2012). Field observation indicates that the water pressure within the crack is likely important (Hooyer and others, 2012). Water pressure in the cavity and the crack are different from each other during water level fluctuation. When water pressure in the cavity drops abruptly, the internal water pressure of the crack does not reduce immediately, and the pressure difference promotes subcritical crack propagation. In our LEFM, the effect of water pressure within the crack most probably happens during the first crack pattern. If the difference of water pressure between the cavity and the crack significantly impacts crack extension, we might underestimate the first crack pattern's length. However, when water pressure in the cavity rises rapidly, the water pressure change does not affect the second crack pattern because this difference in the value of water pressure happens at that moment of water pressure rising in the cavity. The delayed subcritical crack propagation indicates that there is sufficient time for the water pressure to equalize inside and outside of the crack.

To examine the effect of pre-existing crack location on subcritical crack propagation, we used PFM in COMSOL Multiphysics software. According to Figures 12b–12d, pre-existing crack propagation happens on the cavity side of the ice-bed contact region or its edge. However, this location of crack propagation is beyond the scope of LEFM of glacial quarrying. The LEFM is based on tensile stress close to the corner of the ice-bed contact region (Hallet, 1996), while PFM predicts crack initiation and propagation using a scalar phase field. Although the two methods produce different predictions for crack evolution, the validity of LEFM and PFM has been confirmed by experiments (Cohen and others, 2006; Miehe and others, 2015). We emphasize that the LEFM and PFM models are approximate theories about crack propagation. A reasonable speculation is that the subcritical crack initiation and propagation happen on a broader scale, including the ice-bed contact region and its adjacent region.

Glacial quarrying is a complicated process. Some earlier research uses the Weibull probability distribution of bedrock strength to calculate the erosion rate of quarrying, which assumes that the strength depends on its weakest component (Iverson, 2012; Ugelvig and others, 2018). In these theories, the mean volume eroded is half the ice-bed contact region. By comparison, our model analyzes subcritical crack initiation and propagation based on sliding velocity, bedrock characteristics and water pressure fluctuations. Although this crack propagation model cannot be applied directly to the calculation of quarrying rates, it does provide a description of crack initiation and subcritical crack length. Moreover, the PFM finding suggests the location of the pre-existing crack propagation beyond the ice-bed contact region. It implies that the location of crack propagation on the surface of bedrock is subject to a probability distribution. We now know relatively little about the ice-bed contact region and its adjacent region. Further research might contribute to calculating the mean volume eroded of bedrock based on statistical probabilities.

4.4 Evidence from a glacier surge

Currently, fast glacier motion, such as glacier surges, provides clues to the crack propagation patterns mentioned above. In some glacier surges, high velocities are correlated with a high water level. Field survey results also indicate a high sediment output during glacier surges (Hallet and others, 1996; Benn and Evans, 2010). An essential feature of the surge of Variegated Glacier is that the water pressure rose rapidly to a level near or greater than the ice-overburden pressure and then began a slow decline after reaching the peak. The process continued for a day or two, with a sliding velocity of 100–300 cm d⁻¹ (Kamb and Engelhardt, 1987; Raymond, 1987). For example, during the 5th mini-surge of Variegated Glacier in 1980 (Raymond, 1987), water levels quickly rose almost 150 m and maintained at 100 m for two days; this is exactly the characteristic of water level fluctuation concerned with the second crack propagation pattern in Figure 10. It could be speculated that if the subcritical crack propagation happened during the surge, the second pattern might account for a considerable proportion of damage. Additionally, compared to the first pattern, the second pattern's water level changes are more effective at promoting crack propagation, which should, in theory, result in a high quarrying rate. Despite a high sediment output in the surge of Variegated Glacier (Humphrey and Raymond, 1994), more research is needed to analyze the distribution of sediment's grain sizes (e.g. Crompton and Flowers, 2016) and identify how much of it was produced by quarrying and abrasion (Riihimaki and others, 2005).

5. Conclusion

We have presented an assessment model to analyze cavity length and subcritical crack propagation in glacial quarrying with water level fluctuation. In this model, the Euler method was used to solve the differential equations that describe cavity size, and subcritical crack propagation was estimated by LEFM. The results have shown that the cavity length is much less in periods of water level fluctuation because of the higher cavity closure rate. The water level fluctuation and the cavity length variation are not synchronized. Moreover, suppose the water level recovers from the fluctuation to the initial steady state. In that case, for a short time the cavity is longer than it was during the initial steady state. With the water level fluctuation and recovery, there are two patterns of subcritical crack propagation. The first one stems from the rapidly increasing effective pressure when the water level falls during the first fluctuation cycle. The second one happens once the water level returns to the initial steady state following the fluctuation. By comparing the two patterns, we have shown that the second is more effective at promoting subcritical crack propagation. The findings suggest that if crack propagation relies on water level change, this propagation will happen in the short term and cease once cavity length is reduced sufficiently for the stress intensity factor to fall below the stress corrosion limit.

Based on assessing the impact of water pressure fluctuations on cavity length and subcritical crack propagation, the effect of a pre-existing crack location was analyzed. The result of PFM implies that the subcritical crack initiation and propagation happen on a broader scale, including the ice-bed contact region and its adjacent region.

Overall, this study's results provide new insights into subcritical crack propagation during subglacial water pressure variations. Further research is required to understand better the relationship between the two patterns of subcritical crack propagation, the subglacial water pressure variations and the interior features of the ice-bed contact region.