Introduction
Quarrying, the removal of bedrock blocks from the bed of a glacier, is a dominant mechanism of subglacial erosion (Bennett and Glasser, Reference Bennett and Glasser2009; Alley and others, Reference Alley, Zoet and Cuffey2019). At present, however, limited understanding of quarrying mechanics impedes efforts to analyze and predict rates of glacial erosion, development of glacial landforms, production of glacial sediments, and maintenance of the topographic roughness that restrains glacier sliding.
Quarrying proceeds by a sequence of processes, starting with fracturing of bedrock due to tectonic and unloading stresses, additional fracturing in the subglacial environment, loosening of distinct blocks, and finally entrainment and removal. The relative importance and difficulty of these steps varies from situation to situation. There is, none the less, agreement that pre-existing fractures in bedrock facilitate quarrying (Matthes, Reference Matthes1930; Krabbendam and Glasser, Reference Krabbendam and Glasser2011; Hooyer and others, Reference Hooyer, Cohen and Iverson2012; Woodard and others, Reference Woodard, Zoet, Iverson and Helanow2019; Theiss and others, Reference Theiss, Cuffey and Zhao2024). Without them, quarrying must be sparse or absent; massive bedrock tends to be abraded, not quarried. New fracturing in the subglacial environment might, in principle, be driven by glacier loading, by freezing, by hydraulic pressure gradients (sometimes accompanied by wedging of sediment into cracks), or by tectonic and topographic stresses during glaciation. These drivers could interact with one another, and also act as agents of loosening.
The present study focuses on subglacial fracturing due to glacier loading. A sliding glacier presses unevenly on its substrate, especially if cavities form in the lee of bedrock protrusions. Prior work attributes quarrying to fracture growth driven by differential stresses caused by glacial loading. The broad hypothesis is that the weight of the glacier—properly, the effective pressure, meaning the fraction of weight unsupported by subglacial water pressure—concentrates normal stress on the rock lip immediately upstream of a cavity and interacts with pre-existing discontinuities to grow fractures in the bedrock (Iverson, Reference Iverson1991; Hallet, Reference Hallet1996; Woodard and others, Reference Woodard, Zoet, Iverson and Helanow2019; Li and others, Reference Li, Huang and Su2023; Theiss and others, Reference Theiss, Cuffey and Zhao2024). The focus on load-induced subglacial fracturing in the present study and these earlier works is motivated by observations of recently exposed forelands of glaciers eroding into bedrock, such as those described by Hooyer and others, (Reference Hooyer, Cohen and Iverson2012). The bedrock commonly displays zones of rough and angular form, but in-place fragments are absent or rare. Intersecting fracture systems are often prominent and important, but can be irregular, partly discontinuous, and/or widely spaced compared to the scale of observed angular facets. Bedding planes and foliation can be fractured but often remain unbroken. In rock masses with discontinuous fracture networks, distinct blocks are not always available for lifting away, and making them so requires a good deal of hammering and chiseling to create new fractures. The same conditions are found in the channels of actively incising bedrock rivers (Whipple and others, Reference Whipple, Hancock and Anderson2000), and on human-excavated surfaces in unweathered rock (Elmo and others, Reference Elmo, Stead, Yang, Marcato and Borgatti2022).
We note that the foregoing description is by no means universal. Glacial forelands of very different character exist. Some are dominated by smooth abraded surfaces, suggesting that abrasion was the primary erosion mechanism. Others are largely buried in abundant fragmentary material dropped by ablating ice, especially in rugged mountain landscapes with abundant landsliding. In yet other places, the bedrock is intensely fractured in situ, in particular in rock masses with well-developed pre-existing fracture networks and ample evidence of glacial quarrying. Many of the best examples of these are not recently exposed forelands but Pleistocene surfaces. Pleistocene surfaces reflect post-glacial weathering to some extent, and the initial condition of rock prior to glaciation is highly uncertain except where erosion excavated deeply.
The focus in this work is rock masses that do not possess a dense, well-connected fracture network, termed here ‘partially intact rock’. The emphasis of our work on new fracture growth means that a key element in our analysis is sections of intact rock between pre-existing fractures, sometimes referred to as ‘rock bridges’. Unconnected fractures are bordered at least on one side by rock bridges that might control failure. Rock bridges are recognized to be an important factor for rock mass stability and permeability in field settings (Terzaghi, Reference Terzaghi1962; Elmo and others, Reference Elmo, Donati and Stead2018; Nilsson, Reference Nilsson2020), and for this reason are also a focus of experimental investigations (e.g. Wan and others, Reference Wan, Liu, Zhao and Fan2019; Qin and others, Reference Qin, Wang and Wang2025; Wang and others, Reference Wang2025) and numerical simulations (Strauhal and Zangerl, Reference Strauhal and Zangerl2021). Hooyer and others, (Reference Hooyer, Cohen and Iverson2012) recognized their importance in the context of subglacial quarrying. In our modelling, rock bridges can also be regarded as a proxy for continuous fractures in a well-connected network that are partially sealed by vein minerals.
Concerning driving processes
Of all possible driving forces, we focus here on glacier loading because it arises, perforce, from the widespread phenomenon of cavitation when glaciers slide and the requirements of force balance to support the glacier’s weight. Other driving forces deserve to be the topic of separate studies. Here, we note a few difficulties with them to explain our different focus.
Freezing is unlikely to fracture rock at temperate glacier beds, for reasons noted by Hallet (Reference Hallet1996). Given the pressure-dependence of ice melting point, elevated pressures caused by ice growth within pre-existing fractures would melt ice, drawing latent heat and reducing the temperature. Lower temperature, in turn, would draw heat from the surroundings and prevent further ice growth. (This argument does not apply near a glacier’s margin, however, where the thermal regime responds to fluctuations of the outside environment. Here, fracturing could occur by freeze-thaw processes characteristic of periglacial settings.) At the glacier’s sole, freeze-on might be an important process for pulling already loosened blocks from the bed (Rothlisberger and Iken, Reference Rӧthlisberger and Iken1981), but the low strength of ice relative to rock implies that new rock fracturing would be unlikely.
Rock fracturing and loosening by transient gradients of water pressure may occur in some situations with pervasive pre-existing fractures (Krabbendam and others, Reference Krabbendam, Palamakumbura, Arnhardt and Hall2021, Reference Krabbendam, Dioguardi, Arnhardt, Roberson and Hall2023), and one component of our own investigation examines how such loadings work together with glacier loading in the context of partially intact rock (described below). This combination was investigated by Iverson (Reference Iverson1991), who pointed out that a reduction of water pressure in a cavity would simultaneously increase glacier loading (by shifting the weight of overburden onto rock adjacent to the cavity) and create an elevated pressure in voids within the rock compared to the cavity, in the case of low-permeability rock lacking transmissive fractures. Iverson’s study did not, however, consider the pattern of loading arising from sequences of cavities, or address the pervasive role of pre-existing macro-fractures.
Rapid subglacial water pressure variations associated with diurnal source fluctuations and lake drainage events have frequently been observed (Cuffey and Paterson, Reference Cuffey and Paterson2010, Chapter 6; Das and others, Reference Das, Joughin, Behn, Howat, King, Lizarralde and Bhatia2008; Andrews and others, Reference Andrews2014; Wright and others, Reference Wright, Harper, Humphrey and Meierbachtol2016). In the presence of ‘hydraulically open’, transmissive macro-fractures, water pressure loadings are limited as a general mechanism by communication of water between the rock fractures and the nearby cavities and water films along the ice-rock interface. Only in the case of rapid transients and low hydraulic transmissivity would significant water pressure differences arise. The transients would be driven by rapid decreases of water pressure in the cavity system, which also immediately increase the normal loading on rock adjacent to the cavities. Because such loading concentrates on the immediate lip of the cavity, resulting deformation of the rock would widen vertically oriented fractures upglacier, favoring drainage and dissipation of pressure gradients (Theiss and others, Reference Theiss, Cuffey and Zhao2024). Increased normal loading also directly suppresses fracturing in some regions of the rock via increased compression while enhancing it in others. If the fracture network is well-connected, transient pressure gradients can lift rock as a loosening mechanism, or break isolated blocks into smaller ones (Krabbendam and others, Reference Krabbendam, Palamakumbura, Arnhardt and Hall2021; Krabbendam and others, Reference Krabbendam, Hall, Palamakumbura and Finlayson2022; Krabbendam and others, Reference Krabbendam, Dioguardi, Arnhardt, Roberson and Hall2023). The very low compressibility of water, however, prevents this from being an effective mechanism for displacing rock unless new water can flow into the fracture. But if the fracture system is permeable enough to allow refilling, why wouldn’t it also be connected enough to the cavity system to prevent excess pressures from developing? None of these concerns are dispositive, and there are no doubt situations where hydraulic jacking and fracturing are favored. Observations in natural bedrock subject to fluctuating water pressures at its surface show that some zones within the same rock mass rapidly communicate hydraulically with the surface whereas others are unconnected (Neupane and others, Reference Neupane, Panthi and Vereide2020). To model this system properly deserves explicit computation of the evolving water system, which itself depends on unknown and highly variable factors. It is not part of the present study, but deserves a future focus.
Also not the subject of the present paper is the possible role of hydraulic wedging in loosening subglacial rock. A contributing factor for hydraulic damage of rock might be transport of fine particles into fractures, preventing their re-closure after periods of forced opening. High water pressures combined with this wedging action have produced clastic dikes in subglacial sediments. The occurrences of fine sediments filling bedrock fracture networks, observed in some formerly subglacial settings, are likely analogous feature (Krabbendam and others, Reference Krabbendam, Palamakumbura, Arnhardt and Hall2021).
A final reason for our focus on subglacial fracturing by glacier loading is the considerable attention given to the prior literature that uses this framework (Iverson, Reference Iverson1991; Hallet, Reference Hallet1996). Part of the motivation for our study is to follow this widely cited work to its conclusion and ask under what conditions quarrying is encouraged or suppressed.
Foci of the present study
Returning to the prior work on load-driven fracturing near cavities, some of the unresolved aspects of this framework render quarrying difficult. First, to complete quarrying, the downward-growing fractures must connect to the cavity via horizontal fractures. Vertical normal loading imposed on the overlying rock lip, however, strongly suppresses formation or propagation of horizontal fractures. Second, fractures grown by normal loading tend to propagate not only downward but also upglacier, away from the cavity, following the direction of maximum tension (Theiss and others, Reference Theiss, Cuffey and Zhao2024). Third, if a down-going fracture does angle down-glacier to connect with a cavity, as seen in one case and inferred in another (Cohen and others, Reference Cohen, Iverson, Thomason, Jackson and Hooyer2006; fig. 3 of Alley and others, Reference Alley, Zoet and Cuffey2019), the residual geometry buttresses the remaining rock and significantly hinders subsequent fracturing upstream (Theiss and others, Reference Theiss, Cuffey and Zhao2024).
Prior work assumes that the first difficulty would be avoided by the existence of pre-existing, continuous horizontal fractures or broken bedding planes (e.g. Hallet, Reference Hallet1996), and such discontinuities are indeed abundant in some settings. In addition, prior work posits that the second difficulty might be avoided by positioning the normal loading downglacier of pre-existing fractures, close to the cavity lip (Theiss and others, Reference Theiss, Cuffey and Zhao2024). Although a single subglacial block-removal event has been induced experimentally (Cohen and others, Reference Cohen, Iverson, Thomason, Jackson and Hooyer2006), such an event has not been modeled numerically. Block separation that produces a favorable geometry for further quarrying events has never been achieved either experimentally or computationally. At present, it is not clear if persistent quarrying can happen at all unless the prior fracture network is dense and well-connected.
In the present study, we explore what factors might allow persistent quarrying to occur in partially intact rock and give examples in which complete separation of blocks is achieved, albeit dependent on the existence of significant prior fractures. Our arguments are based primarily on calculations with a state-of-the-art numerical model that simultaneously evolves fracture growth and stress state within rock, a capability that has not previously been applied to the quarrying problem.
After introducing the modeling method, we will show a set of calculations representing essentially massive rock. The goal is partly to assess what sort of fracturing is possible in this least favorable situation, partly to illustrate how the model behaves in a simple case as a control, and partly to examine the scenario invoked by seminal papers (Iverson, Reference Iverson1991; Hallet, Reference Hallet1996). These pioneering studies envisioned the bed as a step in a rock mass with no nearby macro-fractures explicitly represented (although their existence can be implied, as fractures continuing into the rock beneath the level of the step will not influence stresses and fracturing in the step itself). The step represents the location of a block previously removed, and these studies asked whether the concentrated glacier load engendered by the step can allow further quarrying and upstream propagation of the step face. Our massive-rock case mimics the Iverson and Hallet scenarios, but we include a single, short, vertical pre-existing fracture, which can result from a previous episode of glacier loading (Hallet, Reference Hallet1996) or can be a pre-existing structural flaw (as in Woodard and others, Reference Woodard, Zoet, Iverson and Helanow2019). All rock bodies contain flaws of various scales. Our scenario simply assumes that one intersects the surface and hence localizes the start of new fracturing.
The Hallet and Iverson models envision that quarrying faces develop perpendicular to ice flow and propagate upglacier. Detailed examination of recently exposed surfaces indicates, however, that quarried faces align with pre-existing fracture sets in the rock rather than with ice flow (Hooyer and others, Reference Hooyer, Cohen and Iverson2012). This is one strand of evidence for a dominant role of pre-existing fractures in controlling the quarrying process, beyond the common-sense notion that more pre-existing fractures implies easier rock removal by the glacier, regardless of quarrying mechanism (e.g. Matthes, Reference Matthes1930). Other evidence comes from the presence of quarried surfaces on well-fractured rock units but not on adjacent massive rock units (Krabbendam and Glasser, Reference Krabbendam and Glasser2011), from the interaction of crescentic fractures with pre-existing joints (Krabbendam and others, Reference Krabbendam, Bradwell, Everest and Eyles2017), from the focusing of erosion on zones of fracture clusters rather than individual discontinuities (Skyttä and others, Reference Skyttä2023), and from the general alignment of eroded zones with rock structures (Addison, Reference Addison1981; Krabbendam and Bradwell, Reference Krabbendam and Glasser2011). Yet more evidence comes from the spatial association between fracture density in bedrock, on the one hand, and the occurrence of abraded vs. quarried surfaces and slow vs. rapid erosion at geological timescales, on the other (Dühnforth and others, Reference Dühnforth, Anderson, Ward and Stock2010).
Consequently, the rest of our calculations are designed to simulate fracture growth in bedrock that already contains pre-existing discontinuities, but does not form a well-connected fracture network. In one scenario, we consider cases for which the pre-existing horizontal discontinuities are nearly connected, and in another scenario we consider what conditions would allow fracturing through larger zones of intact rock, especially in the horizontal direction. In a final scenario, we introduce elevated water pressures in isolated fractures, as an additional loading that acts together with glacier loading. Our target scenarios are highly idealized and should be interpreted as a means to demonstrate tendencies and possibilities. These scenarios were motivated by numerous observations we have made in rock deeply eroded by glaciers and rivers, which match the cases documented by Hooyer and others (Reference Hooyer, Cohen and Iverson2012). New fractures are guided by pre-existing ones, which are extensive but also discontinuous and partly coherent across rock bridges. We acknowledge that there are locations where pre-existing fracture networks are so complete as to render the quarrying problem one of loosening and entrainment, processes not modeled here; the importance of new fracture growth is not supposed to be universal.
Methods
The stress and fracture model
We use commercially available software, Irazu 2D (Geomechanica Inc, 2023), to calculate stress and deformation domains representing bedrock loaded by a glacier. The model simultaneously calculates linear-elastic continuum behavior using a finite-element mesh, fracturing along mesh boundaries, and discrete-element interactions between broken sections (essentially, friction and repulsion). Thus, stresses at fracture boundaries are calculated, not imposed as boundary conditions, except when an internal water pressure is added. Fracture growth calculations account for a zone of interlocking microfractures at a growing crack tip (Munjiza and others, Reference Munjiza, Owen and Bicanic1995). Starting with specified initial configurations for the model domain, a load is imposed and the model steps forward in time, iteratively calculating stresses, elastic deformations, and fracturing. The following paragraphs provide more detail about these essential statements.
The model domain takes the form of a single step on an otherwise rectangular, two-dimensional region, large enough that boundaries do not affect results in the region of interest (as described in Theiss and others, Reference Theiss, Cuffey and Zhao2024, Fig. S1). Normal and shear loadings are applied to a region on the top of the step. Neighboring steps are not represented, which means they are assumed to be spaced apart by approximately twice the step height or more (Theiss and others, Reference Theiss, Cuffey and Zhao2024, Fig. S6). The upper surface of the domain is free to displace in any direction. The left and right boundaries are free to displace vertically but not horizontally. The bottom boundary is free to displace horizontally but not vertically.
The model mesh is two-dimensional, and consists of triangular finite elements with nodes at the three vertices. Along each linear edge between two adjacent elements, fractures can be specified initially or, if the stress level is sufficiently high, can newly form. To simulate fractures, the two nodes at each end of the linear edge are duplicated to form a deformable element, known as a cohesive crack element (CCE). CCEs can be specified as already broken or weakened, while the ones initially representing intact rock are fully cohesive. The CCEs first deform elastically under loading. As the deformation increases, the CCEs can yield, reducing their tensile and shear strengths. With further straining, they can break entirely, reducing tensile strength to zero but still allowing tangential friction. Fig. S1 provides a visual representation of the mesh categories.
For intact parts of the model domain, calculated stresses and strains obey an isotropic, linear-elastic constitutive relationship, solved by iteratively calculating uniform strain of triangular elements and resulting forces at nodes (Munjiza, Reference Munjiza2004; Lisjak and others, Reference Lisjak, Tatone, Grasselli and Vietor2014). For crack elements, boundary-normal convergence generates compression that acts to repel neighboring nodes and edges, while boundary-normal divergence generates tension, and boundary-parallel slip displacement generates shear stress. The formulae governing these behaviors for intact, yielded, and broken crack elements are provided in Supplement S1. In some simulations, a water pressure is specified for a broken crack element, and this is directly added as a force acting on the nodes on both sides of the crack.
Key to calculating fracture development is that the local tensile and shear stresses on the crack element rise with displacement until critical deformation values are reached, analogous to the macroscopic yield strengths of materials. Shear strength obeys a Mohr-Coulomb relationship, hence dependent on a cohesion and friction angle, whereas tensile strength is a single parameter. These critical deformation values are defined by energy-based criteria as shown in Fig. S2.
After critical values are reached, yielding begins and crack-boundary displacements grow. The boundary stresses correspondingly decrease, representing the effects of microstructural damage. Finally, breakage occurs when critical values for boundary separation and slip are attained (Formulae and a diagrammatic representation of this evolution are given in Supplement S1). Subsequently, the fully broken crack element has zero tensile strength, but compression still occurs as boundary-normal repulsion and Coulomb-frictional tangential forces likewise occur, but only if specified as a model setting.
Importantly, the criteria governing strength and deformation behaviors of the model crack elements use parameters that must be calibrated against experiments on rock samples, by matching model output with observed bulk behaviors. Calibration is not yet available for a wide variety of rocks. In this study, we used parameters for two contrasting rock types.
Calibration and load magnitudes
We use two rock types, a strong granite and a weak sandstone, the former previously calibrated by other researchers, and the latter calibrated against laboratory experiments in Hagengruber and others (Reference Hagengruber, Reda Taha, Knight, Rougier and Stormont2021) in the present study (Supplement S2). Subcritical fracture growth occurs at stresses as small as one-third of critical values (Atkinson, Reference Atkinson1987). Thus, the results we report for specified loading magnitudes apply to subcritical fracturing at smaller loads, ranging from about 30% up to 100% of the applied loads.
In zones of ice-bedrock contact, normal stresses are limited to about 10 MPa by the crushing strength of ice (Hallet, Reference Hallet1996). Shear stresses between bedrock and melting ice are zero, but rocks entrained in basal ice and dragged across the bed produce rock-rock friction (Iverson and others, Reference Iverson2003; Cohen and others, Reference Cohen, Moore, Hooyer, Fischer, Jackson and Iverson2005; Hansen and Zoet, Reference Hansen and Zoet2019), with observed values as high as 0.5 MPa. Average basal drag on glaciers is typically only 0.1 MPa, but rock-rock friction can be highly localized in space and time, allowing for higher values during isolated events.
Our model calculations use a variety of load magnitudes, but we allow them to go as high as 10 MPa for normal stresses and 0.5 MPa for shear stress. These upper limits are not well constrained. The 10 MPa limit derives from the measured tensile strength of ice being approximately 1 MPa, for small samples at 0 oC (Haynes, Reference Haynes1978; Petrovic, Reference Petrovic2003). Idealizing the ice-rock contact zone on a subglacial step as an elastic line-load, the normal stress can rise to about 10 MPa before 1 MPa of tension is generated in the adjacent ice (Hallet, Reference Hallet1996). Direct laboratory measurements of the compressive strength of ice depend on temperature, loading rate, and other variables, but extrapolation of experimental results to 0 oC also suggests maximum normal stresses of order 10 MPa. On the other hand, large samples of ice are considerably weakened by the presence of pre-existing flaws (Petrovic, Reference Petrovic2003), and at the scale of entire crevasse fields, the apparent tensile strength of glacial ice is only of order 0.1 MPa (Vaughan, Reference Vaughan1993). In the subglacial setting, perhaps some patches of ice are capable of exerting 10 MPa of normal load while others cannot.
The 0.5 MPa upper value for shear stress is based on one observation in a subglacial laboratory that remains unexplained. There are high-quality measurements of 0.2 MPa shear in this setting (Iverson and others, Reference Iverson2003). Both theoretical considerations and laboratory measurements indicate that widely dispersed rocks embedded in basal ice press lightly on the bed (Hallet, Reference Hallet1979; Thompson and others, Reference Thompson, Iverson and Zoet2020). Thus, large shear stress values likely arise only in special cases, such as large boulders dragging across the bed, basal ice layers with high concentrations of rock debris, or thin layers of till.
The foregoing uncertainties should be kept in mind when interpreting our results. In general, a rock fracture sequence calculated for a specified imposed load also represents the rock fracturing that would occur at a smaller load acting on a weaker rock. In that sense, our results are qualitatively general. Quantitatively, our results apply to subcritical fracture growth at loads that are as small as approximately one-third the size of the loads actually applied (Walder and Hallet, Reference Walder and Hallet1985; Atkinson, Reference Atkinson1987). Given the probable importance of subcritical behavior in the subglacial environment, we have noted, on each relevant figure, the range of load magnitudes corresponding to equivalent subcritical loadings (i.e. 30–100% of the nominal imposed loads). Further, given the importance of subcritical behavior, in running the simulations we were careful to impose loads slowly enough to avoid inertial effects (Supplement S3). Further details of how the model was applied are provided in the same Supplement.
The role of water
Rapid decreases of subglacial water pressure are probably the principal driver of large glacier loads on bedrock steps because they shift the weight of overburden from water-filled cavities to adjacent ice-rock contact zones (Iverson, Reference Iverson1991; Hallet, Reference Hallet1996). The loads enhanced this way include effective normal stresses imposed by the ice and frictional shear stresses imposed by debris-rich basal layers. Loads imposed by large boulders or other isolated rocks are not increased in this situation, however (Hallet, Reference Hallet1979). As loads are simply imposed in our calculations, the results apply to both steady and non-steady conditions as snapshots of ongoing fracture development.
Our model does not calculate water pressures or fluxes. In our calculations, normal stress imposed on modeled bedrock steps should be regarded as effective normal stress, the difference between true normal stress on the rock and water pressure in the subglacial system of cavities and interconnected films. We assume that water fills the cavities and the rock fractures. In the first three sets of calculations (as noted in the Introduction), we further assume that water moves freely enough in the rock fractures to maintain an essentially uniform value, with no pressure difference from the cavity. In the fourth set, we impose an elevated water pressure in the isolated part of the pre-existing fracture network close to the bedrock step and cavity. This represents, specifically, an excess of water pressure relative to that in the cavity and films, nominally induced by a rapid decrease in cavity pressure. The model does not evolve the enhanced water pressure but keeps its magnitude constant. In our standard case, the location of the enhanced pressure also stays constant, but as a sensitivity test, we also perform contrasting simulations in which the enhanced pressure propagates with growing cracks. In nature, rapid fracture growth in a single event would reduce pressure in the fracture, while, in contrast, fracture growth by numerous separate episodes separated by time for re-filling would imply an effective spatial spreading of the elevated pressure condition.
Because reductions of cavity pressure simultaneously enhance normal loading and create elevated pressures within low-transmissive rock fractures, the magnitudes of these two loadings are linked. For our calculations, we consider a simple scenario wherein the transient water pressure drop occurs rapidly enough that the configuration of cavities, and hence the distribution of ice-rock contact areas, does not change at the time-scale considered. Denote the spatial fraction of the bed where ice contacts rock as fr . Further, call the water pressure Pw , the ice overburden pressure Pi (largely controlled by ice thickness), and the initial value of the former a fraction η of the latter (initial Pw = η Pi). Balancing the weight of overburden implies the effective normal stress (i.e. true stress minus water pressure) localized on ice-rock contacts is initially Lo = (1 − η) Pi /fr . Then, if cavity water pressure decreases by a fraction γ of Pi , it simultaneously imposes an internal fracture pressure of that magnitude (γ Pi) and increases the effective normal load to Lo + γ Pi /fr . We use these simple relations to specify loading scenarios involving enhanced fracture pressures.
Scenarios in mostly intact rock
Can quarrying occur in rock with minimal pre-existing fractures? Figure 1a describes a scenario including a bedrock step and a single, short, vertical prior fracture. We refer to the distance of the fracture from the step as the ‘backset’. Normal loading of the entire surface between the fracture and the step lip (LL/λ = 1; the basic scenario of Hallet, Reference Hallet1996) produces trivial fracture growth or growth upglacier (Theiss and others, Reference Theiss, Cuffey and Zhao2024), but shifting the loading zone closer to the step lip induces elastic bending of the rock and downward fracture growth (Fig. 1b, 1c). For a tough granite and large but realistic load magnitudes, however, the new fracture does not reach the bottom of the step and does not trend downglacier toward the cavity (LL/λ = 0.7 and 0.4 in Fig. 1b). To maximize the potential for quarrying, we then switched rock properties to a weak sandstone and allowed the loaded zone to progressively shrink from LL/λ = 0.5 to 0.05. In this case (Fig. 1c) more fracturing occurs, and at a lower stress, and it even turns toward the downglacier cavity. It cannot complete the quarrying process, however, because the fracture cannot propagate across a zone of vertical compression caused by the normal loading at the step surface.
(a) Model setup and geometrical variables. Along the ice-bedrock contact, σ is the effective normal stress and τ the shear stress. (b) Calculated fracture growth in granite for four scenarios, three with fixed L L/λ and one in which LL/λ progressively reduces. Shear loading is included in the third and fourth scenarios (turquoise and magenta). (c) Calculated fracture growth for a sandstone, with LL/λ decreasing in steps. (d) Repeat of (c) but including shear loading. Coloration in (d) and arrows in both indicate displacement (maximum = 0.8 mm in c).
To complete the quarrying event in this situation, we have found that it is necessary to add a shear load to the step surface. In the granite (Fig. 1b, magenta), a combination of normal and shear loads progressively concentrated toward the step lip successfully connects the growing fracture to the lee cavity. In the weak sandstone (Fig. 1d), the same combination succeeds at a shear stress magnitude that is readily achievable subglacially. This combined loading scenario succeeds because it induces a clockwise rotation of the whole step corner (Fig. 1d).
In these ‘successful’ scenarios, if the glacier removes the fractured block, the new bedrock surface is not rectilinear but instead curves downglacier, a geometry that buttresses the newly chiseled bedrock step and suppresses further fracturing (Theiss and others, Reference Theiss, Cuffey and Zhao2024), reducing or precluding further quarrying. Such concave failure planes are observed in the lee of some roche moutonnée (fig. 3 of Alley and others, Reference Alley, Zoet and Cuffey2019), were achieved experimentally (figs. 8 and 9 of Cohen and others, Reference Cohen, Iverson, Thomason, Jackson and Hooyer2006), and might be explained by the scenario outlined here. We suggest, however, that ongoing quarrying would be greatly facilitated by a more rectilinear geometry, and that this generally arises from intersecting prior fractures.
Scenarios with intersecting prior fractures
In some rock masses, intersecting planes of weakness (including fractures, bedding planes, and foliation) can define the angular forms of eroded surfaces and thus provide for subsequent quarrying events. Figure 2 illustrates several calculations for a weak sandstone with intersecting pre-existing fractures, separated by rock bridges (Fig. 2a). Without accompanying shear load, a spatially fixed normal load drives fracturing whose direction depends on how strongly the load focuses on the lip of the bedrock step (Fig. 2b). The most focused case (LL/λ = 0.2) successfully quarries a block, but none of the cases break the horizontal rock bridge (defined in Fig. 2a) to activate the prior intersecting network for quarrying. Progressively narrowing the loading zone toward the step lip, however, not only breaks the horizontal bridge, creating a square quarried block, but also extensively fractures and separates pieces of unweakened rock from a narrow zone along the step riser (Fig. 2c). For a granite, identical conditions as in Fig. 2c failed to sever the horizontal rock bridge even with the normal load increased to σ = 20 MPa (not shown).
(a) Multiple-fracture initial condition and definition of terms. As in Figure 1, σ and τ are imposed effective normal stress and shear stress, respectively. (b) Calculated fracture growth for three fixed values of LL/λ. In (c) and (d), LL/λ decreases progressively. (d) includes a shear stress, with a value 0.05 σ. Arrows indicate displacement (maximum of 1.8 mm). Coloration of new cracks in (d) indicate LL/λ at failure. Darker blue values than at LL/λ = 0.675 indicate vertical fracture growth that occurred before maximum load was achieved.
The addition of shear load greatly facilitates the quarrying process (Fig. 2d). In weak sandstone, a combination of moderate load magnitudes (as low as σ = 2 MPa and τ = 0.1 MPa for subcritical fracture growth) and a progressively narrowing load zone liberates a rectangular clast and leaves behind an angular bedrock morphology conducive to further quarrying events. A clockwise torque accounts for the ability of this case to break the horizontal rock bridge in tension (cf. fig. 6 of Theiss and others, Reference Theiss, Cuffey and Zhao2024). Raising normal and shear stresses to maximal plausible values permits granite to break and quarry in the same fashion (not shown), but only if the pre-existing fracture network is at least as extensive as the one used in Fig. 2.
Promoting horizontal fractures
The calculations behind Figs. 1 and 2 indicate that the most difficult aspect of fracture growth that can lead to quarrying is the completion of horizontal segments. This difficulty, unsurprisingly, reflects the dominance of the vertically directed glacier weight as a loading force, if concentrated on contact zones adjacent to cavities (fig. S4A of Theiss and others, Reference Theiss, Cuffey and Zhao2024). Of course, this difficulty does not exist where the bed already possesses continuous exfoliation or sheeting joints, or already-separated bedding planes or foliation.
The horizontal rock bridges in Fig. 2 are quite small, so the quarrying relies on extensive pre-existing horizontal fractures. One situation that relaxes this difficulty is a geometrical asymmetry wherein vertical fractures are more narrowly spaced than horizontal ones, creating step heights considerably greater than the backset of fractures. Figure 3a illustrates an example for granite with the largest horizontal fractures achieved by our model: the step heights are the same, but the initial horizontal rock bridge is roughly twice as wide as in Fig. 2, while the spacing of vertical fractures relative to the step height is half of that in Fig. 2.
As in prior figures, σ and τ are effective normal and shear stresses, λ is the crack backset distance, and LL is the length of the ice-bedrock contact zone. (a) and (b) Calculated fracture growth given a sequentially decreasing LL/λ, colored according to LL/λ at failure. (a) step is 1 m high, λ = 0.5 m, and the horizontal bridge = 0.38 m. All magenta segments indicate fracture growth at LL/λ = 0.64. (b) λ = 1 m and the horizontal rock bridge = 0.6 m. However, all new fractures occurred along pre-specified planes of weakness, half as strong as the ambient granite, linking the initial fractures. (c) Step height = 1 m, and initial cracks occur as multiple sets of dimensions 0.4 by 0.2 m. A horizontal load compresses the domain from the left, while 0.25 MPa of effective normal stress acts on its upper surfaces. The colour of new fractures indicates the horizontal load at time of failure.
Another way to reduce dependence on prior horizontal breaks is for new fracturing to exploit other structural weaknesses in the rock; bedding and foliation planes that are not yet broken at macro-scale but that lie in the plane of existing horizontal fracture segments, or fractures that have been sealed by chemical precipitate (Fig. 3b). This separation scenario does not require as extensive a fracture network as the previous unweakened cases and produces truly planar surfaces, avoiding the irregular surfaces produced by the co-evolution of stresses and fracture growth in homogeneous rock (Figs. 2b–d and 3a). Unfortunately, at present, we cannot interpret these scenarios quantitatively in terms of required stress magnitudes because the models have not yet been calibrated against experimental rock properties for such structural weaknesses.
As another possibility, horizontal fractures might grow from horizontal loading created by tectonic, topographic, or unloading stresses. It has long been recognized that sheeting joints might facilitate glacial quarrying (ch. 8 of Sugden and John, Reference Sugden and John1976), and that they are abundant features in the shallow zones of otherwise homogeneous and massive bedrock (Martel, Reference Martel2017). Can such fractures continue to develop as the glacier erodes down, as in Leith and others (Reference Leith, Moore, Amann and Loew2014)? Growth of such fractures can only occur without significant vertical compressive loading. While easiest to achieve during interglacial periods, such conditions might prevail subglacially during periods of minimal variability and cavity development but high water pressure, as occurs during wintertime beneath some glaciers (Mathews, Reference Mathews1964; Wright and others, Reference Wright, Harper, Humphrey and Meierbachtol2016). It also requires that the rock be permeable enough for pressurized water to fill widely distributed voids and produce an effective vertical stress much smaller than the overburden load. Figure 3c illustrates a case, near a bedrock step, of tensile fractures growing from an array of small pre-existing fractures due to horizontal loading, opposed in the vertical by 0.25 MPa of effective stress. New quasi-horizontal fractures develop far from a bedrock step (similar to fracture growth calculated from random arrays of flaws, as in Olson and Pollard, Reference Olson and Pollard1989), and these may be exploited subsequently by quarrying as step faces migrate up-glacier through repeated block removal. Near the step itself, however, the stress field is oriented obliquely such that new fractures do not grow horizontally.
Enhancement by hydrofracturing
With sufficient magnitude, elevated water pressures within fractures may drive new fracture growth in various directions favored by the co-evolving stress field. Such hydrofracturing provides another potential way to connect fractures and loosen the rock mass. Fluid pressure reduces the effective stress in the horizontal fractures and hence reduces vertical compression. Hydrofracturing could also help to link vertical fracture segments, because stress concentration at the fracture tips promotes their propagation. High pressure in fractures of other orientations and positions might also contribute to torque on the step corner. In the present scenario, such effects are favored because the bedrock step can deform toward the open cavity.
Figure 4 illustrates calculations for a representative case corresponding to an initial overburden pressure of Pi = 5 MPa (i.e. glacier thickness of about 500 m), an initial water pressure of 90% of overburden, and a rock-ice contact fraction of 1/3 (see Methods). As cavity water pressure drops, normal loading at the step surface increases simultaneously, as specified in Methods. We specify that water pressure load builds in the pre-existing fracture set closest to the step (Fig. 4a), without any dissipation by drainage or crack deformation. We stop the calculation when cavity water pressure has dropped by a magnitude equal to 50% of overburden (55% of initial water pressure).
Results of experiments that include forcing by enhanced water pressure. (a) The model set-up, which includes pressurized water in the two intersecting fractures as indicated. (b) Simulated fracture growth as loading increases, up until the moment that the vertical fractures link. Colour bar indicates load magnitudes. (c) Continuation of fracture growth as in (b), until the horizontal fracture links. (d) A separate case without the initial vertical fracture at the step tread. In (a), (b), and (c), white arrows indicate displacements. In case (d), allowing the water pressure enhancement to propagate with the growing fracture reduces the load required for completion by about 20% (case not shown).
For sandstone, vertical fractures grow readily and link to the upper surface at stress values corresponding to a cavity water pressure reduction of approximately 11% of initial value (Fig. 4b). With this connection made, fracture water pressure would likely dissipate and prevent further hydrofracturing. If we continue the simulation with water pressure anyway, a horizontal fracture connection to the cavity subsequently occurs at stresses equivalent to cavity water pressure reduced by approximately 28% (Fig. 4c). A vertical fracture connection also occurs in the case of granite, but requires a larger cavity water pressure reduction, equivalent to 47%. No horizontal connection occurs in granite for the maximum water pressure loading. In sandstone, vertical fractures connect to the surface before horizontal ones connect to the cavity, even if no pre-existing vertical fracture exists on the upper step tread (Fig. 4d), although in this case a significantly larger (a multiple of ∼1.8) reduction in cavity water pressure is required. Favoring of vertical fracturing compared to horizontal is a consistent behavior across all simulations we have done, independent of parameters describing the initial glaciological conditions and fracture configuration. This reflects the persistent effect of vertical compression induced by overburden, which suppresses horizontal fracturing, combined with the bedrock step’s freedom to deform toward the cavity, which favors tension on vertical planes. An interesting additional aspect of the vertical fracture linkages simulated here is that they occur primarily by downward growth of the upper fracture originating at the step surface, rather than by upward growth of the pressurized, vertically oriented internal fracture.
How significant is the role of hydrofracturing in situations such as in Fig. 4? For the large water pressure loadings imposed here, it is considerable. Averaged across a variety of initial fracture spacings, it reduces the normal load required to link vertical fractures by nearly a factor of two, compared to the case of a well-drained system with no water pressure loading (Fig. 5a). The hydrofracturing influence is much smaller if we consider the vertical extent of fracture growth rather than actual linking. The key distinction here is that, in many cases without hydraulic loading, new fractures grow down from the surface but deviate from vertical, such that they do not link with the tips of underlying, pre-existing, fractures. The new fractures reach the same elevation (and hence the same depth-below-surface) as the underlying tips, but at a different horizontal coordinate.
(a) The loading required to link vertical fractures in the same scenario as shown in Figures 4a and b, except using a variety of sizes of initial intersecting fractures at depth, and hence different initial sizes of the vertical rock bridge. Blue points use the pressure, simultaneous normal load, and pressure loading as described in the text and depicted in Figure 4b. Red points are a sensitivity test demonstrating what happens if the pressure loading fixed in space is replaced by one that propagates along with growing fractures. Green points illustrate the same calculation as for the blue and red ones, but without any pressure loading. (b) Same as panel a, but instead of the load required to link fractures, the load required for the downgoing vertical fracture to reach the same vertical coordinate as the upgrowing one (‘equivalent depth’).
Figure 5b plots the normal load required for the downward growing fracture to reach this depth, which we denote ‘equivalent depth’ in the figure. There is a significant influence only if co-planar rock bridges between initial vertical fractures are large. The greater influence of hydrofracturing for linking (Fig. 5a compared to 5b) arises because the lower pressurized fracture primarily acts by making the co-planar stress field more tensile horizontally, which guides the lower tip of the upper fracture as it approaches. Without this hydraulically induced modification to the stress field, localized compressive stress concentrations in the residual rock bridge divert the down-growing crack and slow its propagation before linkage. There is a random element to the latter process (probably related to the interaction of stress concentration with mesh topology) which produces variability as seen in Fig. 5a, and only general features of these results should be interpreted.
Concerning horizontal fracturing (Fig. 4c), hydraulic loading is essential in these scenarios, as no horizontal connection is achieved without it for either sandstone or granite. On the other hand, it is not clear that the horizontal fracturing simulated here is possible, given the likelihood of pressure dissipation when vertical fractures connect to the surface first. Alternatively, if a complete horizontal fracture already exists, drainage along it would prevent the build-up of any effective loading and hydrofracturing at all, unless the fracture’s transmissivity is low.
Discussion
The following discussion applies to the context of this study, which is to say situations wherein glaciers erode into bedrock that is partially intact and must be fractured for quarrying to occur.
All of the ‘quarrying’ events that appear in our calculations require a small contact area between ice and subjacent bedrock, partly to enhance normal load magnitude and partly to induce elastic bending (Figs. 1–4). Small contact areas, in turn, correspond to large subglacial cavities and hence glaciological conditions of rapid sliding motion and basal water pressures approaching the weight of overburden. Hallet (Reference Hallet1996) and Iverson (Reference Iverson2012) both reached the same conclusion, but our own work (here and Theiss and others, Reference Theiss, Cuffey and Zhao2024) has identified additional fundamental difficulties for quarrying of partially intact rock. Studies relying on scalings derived from these earlier analyses should be viewed with due skepticism (e.g. Ugelvig and others, Reference Ugelvig, Egholm and Iverson2016; Ugelvig and Egholm, Reference Ugelvig and Egholm2018).
Rapid declines of subglacial water pressure likely enhance quarrying by shifting more overburden load onto ice-bedrock contact zones (Hooke, Reference Hooke1991; Iverson, Reference Iverson1991; Hallet, Reference Hallet1996). Our calculations suggest, in addition, that episodes of increasing water volume at the ice-bedrock interface (Harper and others, Reference Harper, Humphrey, Pfeffer and Lazar2007) are particularly favorable for driving quarrying. In such times, cavities expand and normal loading spatially focuses more on the lips of bedrock steps (Figs. 1, 2, 3a, 3b). Our calculations further suggest that very different periods, those with minimal cavity extent and sliding, could favor horizontal fracture growth in settings with significant tectonic, unloading, or topographic stresses (Fig. 3c; fig. 26 of Holzhausen, Reference Holzhausen1989). Alternation of such conditions with periods of rapid sliding and large cavities might be ideal for quarrying, if horizontal fractures are not already continuous in the rock mass.
Rapid declines of subglacial water pressure could enable quarrying by inducing hydraulic fracturing, but only if the permeability of the rock mass and transmissivity of fractures are low (Iverson, Reference Iverson1991; Krabbendam and others, Reference Krabbendam, Palamakumbura, Arnhardt and Hall2021). Our calculations suggest that hydraulic loading reduces the loading necessary to complete the connection of vertical fractures by roughly a factor of two (Figs. 4 and 5). If pre-existing horizontal fractures are continuous, such vertical fracturing can sever layers and complete the separation of a block from its parent rock mass. In the contrary case of discontinuous prior horizontal fractures, hydraulic fracturing can connect them together even in situations where the largest plausible loads fail to do so in the absence of hydraulic loading (Fig. 4c). In our calculations, however, such horizontal connections always occur after the vertical connections are complete. Dissipation of water pressure might therefore stall the process before achieving horizontal connection.
Our calculations indicate that shear loading might be a key driver of quarrying (Figs. 1, 2c, 2d). Shear traction at an ice-bedrock interface must arise from rock-rock friction via particles entrained in the basal ice or forming thin granular layers. Such a role for rock fragments implies that—analogous to abrasion—rates of quarrying depend on the supply of rock material from upglacier. In alpine landscapes, rockfall and landslides from cirque headwalls and lateral cliffs may thus facilitate quarrying down-glacier. Subglacially, a zone of bedrock that is particularly susceptible to quarrying could add debris that fosters erosion in a zone downglacier. In the case of abrasion, unless abrading tools consist of harder rock than the substrate they act on, basal debris expends itself at the same rate as erosion of the subjacent bedrock. In contrast, basal debris driving quarrying through shear traction would not suffer this fate, in general. In this sense, particles in basal ice are more-efficient agents of quarrying than of abrasion.
A bedrock mass with vertical fractures more closely spaced than horizontal fractures is particularly susceptible to quarrying, according to our calculations depicted in Fig. 3b. This implies that zones within mountain belts of steeply dipping layers of metamorphic rock or thin-bedded sedimentary rock might be particularly susceptible to glacial erosion.
Our calculations show how weak rock tends to break into small fragments instead of the large chunks expected for a tough rock (Fig. 2c). Though perhaps unsurprising, this is useful because it can be readily observed if true, and we illustrate one possible example in Supplement S4. This fracturing occurs even in the absence of a pre-existing microfracture to nucleate it. The FDEM method includes the capability to create completely new fractures, as explained in the Methods; it only requires that tensile or deviatoric stresses between two elements be sufficiently large. In the case of Fig. 2c, such large stresses occur because of the combination of spatially concentrated loading and its position atop the free surface of the step riser.
Key limitations of our study include lack of empirical tests for inferred behaviors and lack of theoretical three-dimensional constraints. The latter problem exceeds current computational abilities, but the former could be targets for state-of-the-art glaciological laboratories (e.g. Iverson and Petersen, Reference Iverson and Petersen2011). Another key limitation is the suite of uncertainties precluding confident calculation of transients of water pressure in subglacial fractures: especially the transmissivity of rock fracture networks, how transmissivity changes with rock deformation, and how water pressures evolve as cracks grow. Without the ability to confidently calculate the behavior of water pressure transients, important processes of hydraulic fracturing cannot be modeled well. A related limitation of the present work is that we do not provide simulations for different bedrock configurations, such as isolated roches moutonées, for which vertical glacier loadings are smaller but horizontal loads possibly important (Hättestrand and Stroeven, Reference Hättestrand and Stroeven2002; Krabbendam and Hall, Reference Krabbendam and Hall2019).
We have made no attempt to link our calculations to a model for the rate of erosion by subglacial quarrying. That would require a lesser degree of idealization in our model configurations, a greater accounting for the diversity of subglacial conditions, and the incorporation of a broader set of processes, such as loosening and entrainment that complete the quarrying process and likely rate-limit it in some situations. Nonetheless, our work should be regarded as a step forward in the much longer-term effort to develop a process-based understanding of observed magnitudes and patterns of glacial erosion (Cook and others, Reference Cook, Swift, Kirkbride, Knight and Waller2020; Norris and others, Reference Norris2025).
Conclusion
Thirty years ago, theoretical analyses identified concentrated normal loading on subglacial bedrock steps as the primary driver of quarrying of intact rock, amplified by variability of water pressures in cavities. Geological analyses before and since have shown that pre-existing fracture networks greatly facilitate the process. Can glacial loading drive new fracture growth to allow quarrying, as envisioned in foundational theoretical work, or must such prior networks be continuous? We investigated the two-dimensional form of this problem using a model that rigorously calculates the co-evolution of stresses, elastic deformations, and fracture growth in loaded subglacial steps.
Vertical fracture growth can occur even if pre-existing fractures are minor, but requires normal loading to focus very closely on the corners of steps. Including shear loading can significantly facilitate vertical growth and, in special cases, even allow new fractures to angle downglacier and reach the adjacent cavity. If the glacier removes this block, however, the resulting geometry buttresses the new step, rendering subsequent quarrying even more difficult. In our calculations, intersections of pre-existing planar fractures are necessary for producing geometries favorable to repetitive quarrying.
Horizontal fracture growth is difficult to achieve by glacial loading, and generally requires at least partial pre-existing fractures of this orientation. Rock bridges between horizontal fractures can break if, again, normal loading focuses closely on step corners and if shear loading contributes. Pre-existing horizontal fractures needed for quarrying might be pre-glacial sheeting joints or broken bedding planes, but can also continue to generate subglacially, upglacier of bedrock steps, if tectonic and topographic stresses include horizontal compression.
Additional special cases that facilitate quarrying of partially intact rock include bedrock with vertical fractures spaced more closely than horizontal ones (as sometimes occurs in vertically dipping sedimentary strata), and bedrock that is partially intact because, although they are densely fractured, some of the prior fractures are chemically sealed.
Hydraulic loading (here, in particular, the excess of water pressures in fractures compared to nearby cavities) can significantly reduce the magnitudes of glacial loading needed to connect fractures to quarry. In our calculations, however, this effect primarily arises not because pressurized fractures grow easily themselves but because they make stresses in the surrounding rock more tensile and thereby allows fractures growing downward from the bedrock surface to connect with internal, pressurized fractures.
The key role of pre-existing fracture networks should not obscure the importance of additional bedrock characteristics. In our calculations, the inherent brittleness of the bedrock influences the propensity for quarrying greatly, with a comparatively weak sandstone requiring significantly smaller loads than a tough granite. Hydraulic properties of the rock and its fractures also matter greatly, as Iverson emphasized 30 years ago, because permeability and transmissivity need to be low for hydraulic loadings to develop.
All of our results are consistent with prevailing ideas that subglacial quarrying of partially intact rock is favored by rapid sliding and low effective pressures (which make cavities large and localize glacier-bedrock contact onto step corners) and by variability of the subglacial water system (rapid draining of cavities increases glacial loading and generates hydraulic loading). Our work indicates, however, that earlier efforts to derive scaling relationships between quarrying rate and these variables need to be viewed with skepticism. In addition, we propose that new emphasis should be placed on the possible role of basal debris and the shear loadings it allows.
In the context of our study, the most important questions needing to be addressed are three. What magnitudes of shear loadings arise subglacially? How does the third dimension (transverse to glacier flow) constrain quarrying? In what situations can we expect fracture networks to possess low enough transmissivity to allow strong hydraulic loadings?
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2026.10141.
Data availability statement
The dataset(s) used in this paper can be found on https://doi.org/10.5281/zenodo.18521104.
Acknowledgements
We are grateful to Neal Iverson and Maarten Krabbendam for their critical reviews of our manuscript, which helped turn the ugly caterpillar of our original submission into the beautiful butterfly you have just read. A gift from The Martin Family Foundation to UC Berkeley funded this work. Q. Zhao and C. Xia are supported by RGC GRF grant No. PolyU 15227222.
Competing interests
The author(s) declare none.
Open access
