## Abstract

We present a mathematical model of the fracturing of water-saturated rocks and other porous materials in cold climates. Ice growing inside porous rocks causes large pressures to develop that can significantly damage the rock. We study the growth of ice inside a penny-shaped cavity in a water-saturated porous rock and the consequent fracturing of the medium. Premelting of the ice against the rock, which results in thin films of unfrozen water forming between the ice and the rock, is one of the dominant processes of rock fracturing. We find that the fracture toughness of the rock, the size of pre-existing faults and the undercooling of the environment are the main parameters determining the susceptibility of a medium to fracturing. We also explore the dependence of the growth rates on the permeability and elasticity of the medium. Thin and fast-fracturing cracks are found for many types of rocks. We consider how the growth rate can be limited by the existence of pore ice, which decreases the permeability of a medium, and propose an expression for the effective ‘frozen’ permeability.

## 1. Introduction

Large pressures can develop inside water-saturated porous media at sub-zero temperatures. These pressures occur owing to the solidification of water inside the pores and can cause fracturing of pre-existing faults, leading to the degradation of the rock. A similar process, called frost heave, occurs in soils, where segregated ice lenses form devoid of soil particles, causing upward movement of the ground above. This frost-induced deformation of material can destroy building foundations, damage roads and statues, and deform cooled gas pipelines laid in submarine sediments. Phenomena such as patterned ground, potholes, needle ice and pingos play an important role in the development of landscapes [1].

The problem of frost fracturing of rocks is closely related to that of frost heave in soils, which describes the upwards displacement of the ground surface caused by the formation of particle-free lenses within a soil. Frost heave has been studied extensively, starting with Taber [2,3], who performed experiments freezing benzene in a column of soil to demonstrate that the expansion of the freezing liquid is not an important contributor to the heaving, as benzene contracts upon freezing.

The volumetric-expansion model is dependent on the idea that water is forced away from the solidification front and raises the pressure inside the medium, as the space that can be occupied by water is limited. This requires the medium to be completely (or, at least, considerably) saturated with water. If a large proportion of the pores were empty, the water would simply fill them and the pressure would be relaxed. In particular, as freezing water expands by 9%, a minimum saturation level of about 91% is required for fracturing to occur by this mechanism. However, experimental studies have shown otherwise, with fracturing being reported for saturation levels as low as 65% [4] or 80% [5]. In addition to this, the volumetric-expansion model predicts fracturing in bursts, and only during the freezing cycle, which isn't supported by experimental evidence [4]. A scenario in which volumetric expansion would have an important role is that of a rock freezing rapidly from all sides. However, the applications of such a model are limited.

Taber's studies suggested that thin films of melt exist between the ice and the solid particles, with liquid from the surrounding medium being drawn towards them. A lot of subsequent research was focused on exploring the nature of these films. Several studies [6–9] considered a surface-tension driven flow towards the freezing front, while Gilpin [10] simply assumed that the water in the thin films at solid–liquid interfaces experiences an attractive force from the solid boundary. Vignes-Adler [11] suggested that water in thin films has different properties to water in bulk, but experimental work by Raviv & Klein [12] does not support this suggestion. A model was developed by Rempel *et al.* [13] based on the idea that repelling forces between the ice and the soil particles create these premelted thin films, while pushing the two substrates apart. The pressure in the water is lowered as a result of the forces tending to widen the gap, and more water gets sucked in.

The idea of the importance of this *premelting* regime versus the expansion mechanism can also be used in the modelling of frost fracturing, as there are strong similarities between the two processes. As in soils, disjoining (intermolecular) forces between the ice and the rock lower the pressure in unfrozen water films adjacent to the ice surface, which draws water in from the surrounding saturated medium. These disjoining forces cause the ice-filled cracks to grow. The difference between the process in soils and rocks lies in the way each medium deforms under the forces exerted by the ice: soil particles can be rejected from the solidification front if the freezing is slow enough; the cohesion of the rock means that the same is not possible. Instead, the pressure exerted on the rock deforms it elastically and can cause fracturing of the cavity if the stress at its tip is above a critical value. While in soils the pressure exerted by the ice is balanced by the overburden pressure, in rocks it is instead balanced by the elastic pressure of the medium. It is the relative magnitude of these two pressures that determines whether the cavity expands and the ice continues to grow.

Walder & Hallet [14] and Murton *et al.* [4] used the concept of ice-filled cavities in rocks acting as ice lenses, and developed models which incorporate the concept of premelted films and flow towards freezing regions. More specifically, Walder & Hallet [14] discussed how these films exert an ‘attractive force’ on the pore water (hence the flow towards the ice front) and a disjoining pressure that pushes the ice and the rock apart. They showed that the fastest growth rate occurs at temperatures in the range −4 to −15°C as, in colder systems, the transport of water is limited by the large amount of pore ice reducing the permeability of the rock.

The aim of this paper is to develop a complete model for rock fracturing during freezing, coupling linear elasticity and fracture mechanics with the fluid dynamics governing flow through the porous medium. A mathematical model is required to quantify the important parameters influencing frost fracturing and interpret experimental observations. It is worth noting that several experimental studies [15–17] focused on fast freezing rates and water content and ignored the structure of the rock, as they explained fracturing within the framework of volumetric expansion.

Vlahou & Worster [18] showed that the pressure contribution from expansion is in most cases negligible compared with that resulting from premelting. Therefore, in this paper, we concentrate on a premelting regime in which the cavity is filled with ice. We follow ideas developed in [14] but formalize the discussion of physical processes using the premelting theory presented in [19]. We develop a similarity solution for a specific time-dependent undercooling, which serves as a mathematical tool to understand the physical interactions in a simple framework. The time-dependent problem is then solved to confirm and extend the conclusions found from the similarity solution.

We show that it is the maximum value of undercooling rather than the rate of freezing which determines the maximum pressures applied on the medium. Material properties such as permeability and porosity also contribute to the flow of water, hence they need to be taken into account when studying the susceptibility of materials to frost damage. Finally, the ability of the medium to withstand pressures without failing is described by fracture mechanics through the fracture toughness parameter, which we show to be just as important as the undercooling in determining the susceptibility of a material.

## 2. Governing equations

### (a) The elastic pressure

We consider an axisymmetric, penny-shaped crack of radius *R*(*t*) and small half-width *b*(*r*,*t*)<*R*(*t*) inside an elastic porous rock of infinite extent. The rock is saturated with water that has been supercooled to a temperature *T*_{m} is the melting temperature of the ice measured at the reference pressure

Because the cavity is thin, curvature melting can become important towards the tip of the cavity. To take this into account, we assume that the ice in the cavity only extends up to *r*=λ*R*, where λ∈[0,1] is a parameter to be determined, with the rest of the cavity filled with water. A sketch of the cavity is shown in figure 1.

We use linear elasticity to describe the pressure field within the solid. For the configuration considered here, we follow [20] and assume that the crack has zero internal shear stress. By symmetry, the problem for the disc-shaped cavity is equivalent to that for a half space *z*>0 with prescribed pressure field on *z*=0 within the circle *r*<*R* and no normal displacement for *r*>*R*. The pressure field in the surrounding elastic medium can then be expressed in terms of the cavity shape as
*m*=*G*/(1−*ν*), *G* is the shear modulus and *ν*=λ/[2(λ+*G*)] is the Poisson ratio of the rock. The kernel *M* is given by
*K*_{ell} and *E*_{ell} are complete elliptic integrals of the first and second kind, respectively. The pressure applied on the boundary of the cavity causes a stress at the tip *K*_{I} given by
*K*, the fracture toughness of the material. In the critical state *K*_{I}=*K*, equation (2.3) together with the expression for the pressure field (2.1) are equivalent to

### (b) The water pressure distribution

The flow of water through the pores of the rock is described by Darcy's equation
*Π* is the permeability of the medium. We combine (2.5) with the mass continuity equation ∇ . ** u**=0 to find that the water pressure

*p*satisfies Laplace's equation

*q*(

*r*,

*t*) towards the cavity, and the crack, being long and thin, can be approximated as a planar sink. The flow along the film can be shown to be negligible compared with the flow into the film [22], hence we can express the flux

*q*(

*r*,

*t*) simply in terms of the rate of opening of the cavity

Complications arise when λ<1, as the region towards the tip of the crack is not occupied by ice but is instead filled with water. If this gap is large enough to allow easy flow in the radial direction, we can approximate the liquid pressure there as uniform. However, close to the tip, it is likely that the width of the crack will be small, and we shall therefore assume that expression (2.8) applies throughout 0<*r*<*R*.

### (c) The temperature field

We assume that the temperature field is quasi-steady and hence described by Laplace's equation
*ρ*_{s} is the density of ice, *k*_{l} the thermal conductivity of water. We define the undercooling *T*, (2.9) and (2.10). By comparing equations (2.7) and (2.10), and noting that both *p*_{l}(*r*,*t*)=*p*(*r*,*z*≈0,*t*) and *θ*^{I}(*r*,*t*)=*θ*(*r*,*z*≈0,*t*) satisfy Laplace's equation and decay as *θ*^{I} can be expressed in a similar way to *p*_{l} as an integral of ∂*b*/∂*t*

### (d) The dynamics of the premelted film

When a solid is separated from a foreign substrate by a thin film of its melt, repelling intermolecular forces can develop between the two substances, pushing them apart. The resulting pressure *p*_{T} is called *disjoining pressure*. The liquid pressure *p*_{l} inside the thin film, together with the disjoining pressure *p*_{T} acting across the film, must balance the solid pressure *p*_{s} in the ice [19]

Owing to the disjoining pressure *p*_{T}, the liquid and solid pressures across the freezing interface are not equal. This pressure difference affects the freezing temperature *T*^{I} as described by the Gibbs–Thompson relation
*p*_{l}−*p*_{m}) as we take *ρ*_{s}=*ρ*_{l} approximately. We also assume the contribution of curvature melting is negligible, except at the edge of the ice lens, where it determines λ.

### (e) The ice extent

Unlike the rest of the ice surface, the tip is highly curved and this can be important. Here, we assume that the main contribution to a solid–liquid pressure difference at the tip is because of its curvature, i.e. we ignore pressure melting and disjoining forces. We can then use the Gibbs–Thompson equation above, evaluated at *r*=λ*R*, to find an expression for the maximum possible tip curvature for a given undercooling. The ice tip is approximated as circular and hence the curvature can be written as *κ*=1/*b*(λ*R*,*t*). The relation determining the ice extent λ is therefore
*γ* is the surface energy of water per unit area.

## 3. Similarity solution

Before considering the full time-dependent problem from a specific initial condition, we explore the following special solution of the governing equations, from which significant insight can be gained. If the far-field undercooling is described by *k* to be determined. In terms of the similarity variable *η*=*r*/*kt*^{1/2}, the tip of the cavity is at *η*=1, the ice extends to *η*=λ and the cavity thickness is *B*(*η*)=*b*(*r*,*t*)/*kt*^{1/4}. Then, the similarity functions *p*_{s}(*η*)=*t*^{1/4}*p*_{s}(*r*,*t*), *p*_{l}(*η*)=*t*^{1/4}*p*_{l}(*r*,*t*) and *Θ*^{I}(*η*)=*t*^{1/4}*θ*^{I}(*r*,*t*) satisfy the following system of equations

The last three parameters can be thought of as dimensionless permeability, fracture toughness and undercooling, respectively. We find that two solutions exist for each set of parameters. We denote their propagation rates as *k*_{1} and *k*_{2}, with *k*_{1}<*k*_{2}, and refer to them as ‘slow’ and ‘fast’ solutions, respectively. Figure 2 shows the cavity width and pressure distribution for a typical fast solution. The difference in cavity thickness between the part of the cavity occupied by ice and the tip is clear in the *B*(*η*) plot. The solid pressure features a step drop at *η*=λ, and a negative spike at *η*=1. The latter agrees with the prediction of infinite solid pressure at the tip of the crack, caused by the assumption of a parabolic tip as predicted by linear elasticity. In practice, *p*_{s}(1) remains finite due to the limitations of the discretized numerical scheme used, but its value is unimportant as it is not used to determine any characteristics of the cavity. Finally, a rise in the liquid pressure *p*_{l} is seen towards the tip of the crack. This is most likely due to the fact that the cavity growth is driven by the ice-occupied region, where water continues to freeze, and hence the flux of water is larger there.

### (a) Stability analysis

Before continuing with the analysis of the different solution regimes, it is important to determine whether the solutions are stable or unstable. We define a modified stress-intensity factor *K** which corresponds to a propagation rate *k* of the tip
*K**=1, as follows from equation (3.7). Figure 3 shows the values of the modified stress intensity *K** plotted against the corresponding propagation rate, with the points of quasi-steady propagation denoted by the squares.

We can now understand how the two solutions react to small perturbations in propagation rates, therefore determining which solution is stable. We note that for *k*_{1}<*k*<*k*_{2} the stress intensity factor is higher than its critical value, and for *k*>*k*_{2} or *k*<*k*_{1} it is lower. Therefore, if the cavity propagates just faster than the slow solution, at *k*=*k*_{1}+*ϵ* (where *ϵ*>0), the pressure builds up higher than the critical value. This forces the cavity to propagate even faster, making the solution unstable. For a propagation rate just slower than *k*_{1}, we find *K**<1 therefore the propagation slows even further. Conversely, for a cavity propagating at *k*=*k*_{2}+*ϵ*, the stress-intensity factor is lower than the critical value, causing the propagation to slow down and hence making the fast solution stable.

### (b) Phase planes

Solutions for (*k*_{1},λ_{1}) and (*k*_{2},λ_{2}) cannot be found for all values of

The main region with no propagation is the one represented by the crosses. For these cases, the stress at the tip of the crack is less than the critical value *K*, even for completely ice-filled cavities. This is the case where the material is tough enough to withstand the pressure caused by the build up of ice for a given undercooling Δ*T*. Fracturing will occur if

The smaller region shown by the circles also describes situations for which propagating similarity solutions cannot be found. The reason we distinguish this region from the crosses is because in these cases, though the material is too weak to withstand the pressure from a completely ice-filled cavity, the curvature-melting effect does not allow the ice to extend all the way to the tip of the crack. We can see from figure 4 that this occurs for small undercoolings Δ*T*, indicating that the maximum curvature determined by the Gibbs–Thompson relation is small. In addition to this, the dimensionless fracture toughness is generally small, representing media that do not deform much under pressure, resulting in thin cavities. An increase in the undercooling Δ*T* results in ice being able to extend further towards the tip, with the resulting stress causing fracturing.

The non-propagating regimes cannot be properly described within the similarity solution framework since *k*=0 implies *R*=0. In reality, we expect an equilibrium situation were the disjoining pressure balances the elastic pressure of the cavity. As the environment warms up with time, the maximum curvature imposed by the Gibbs–Thompson relation decreases, resulting in the melting of ice towards the tip. We therefore expect the pressure to relax slowly and part of the ice to melt as the environment is warming, which is not captured by the characteristics of the similarity solution.

The propagating solutions are again split into two regions, the ones with ice-filled cavities (i.e. λ=1^{1}), denoted by the diamonds in figure 4, and the ones denoted by squares with λ<1, where the Gibbs–Thompson relation dictates how far into the tip of the cavity the ice extends. The squares describe the sets of *k*_{1},λ) and (*k*_{2},λ) satisfying both the tip and ice extent conditions. The diamonds represent sets of parameters where solutions for *k* exist but the curvature of the ice tip is always less than the maximum value defined by the Gibbs–Thompson relation, i.e. the ice can extend to the tip of the cavity and the curvature does not affect the extent of the ice.

We see from figure 4 that the boundary between propagating and non-propagating solutions seems to be linear and hence given by
*β* is a parameter that depends on the value of the permeability

To understand this relationship, we need to examine what changes across the boundary. Equation (2.14) describes the balance of pressures across the premelted film for a specific propagation rate *k*. If no *k* exists for which this balance is achieved, the cavity does not propagate, at least under the similarity framework. The tip condition gives *B*∼*k*^{−1/2}, hence we find that the solid pressure scales as *k*^{−1/2} while the liquid pressure scales as *k*^{3/2} and the temperature field as *k* in terms of the dimensionless parameters of the problem. The left-hand side of the expression has a minimum at *k* exists. This threshold is the linear boundary described by equation (3.9) with a slope *β*. Substituting *β* becomes important (

### (c) Propagation rate and ice extent

The stress-intensity factor at the tip of a crack induced by a loading *p*_{s}(*r*,*t*) on the inner walls of the crack can be expressed as an integral over *p*_{s}, as described by equation (2.3). We can see that the integral is weighted towards the stress contributions closer to the tip. For a small value of λ, the lack of ice (and therefore pressure) close to the tip has to be counteracted by substantial ice growth, enough to raise the stress intensity at the tip to the critical value. This implies that for smaller λ the crack is wider. However, the wider the cavity, the further the ice can extend towards the tip. The balance of these two processes leads to the result that λ increases as *a*. There is a similar explanation for why λ increases with

The propagation rate depends on the freezing rate, as more water needs to freeze to maintain the stress at the tip. The faster the solidification, the quicker the pressure in the cavity builds up and the stress condition at the tip is met. This can be seen in figure 6*b* where, for a given

An interesting feature that arises from the assumption of a warming environment is the importance of the permeability *Π*, which controls the timescale of flow and water supply. Since the undercooling decreases with time, media with low permeabilities can experience no fracturing, even for parameters for which more permeable media would fracture. This is due to the time-dependence of the undercooling and, since the next section is predominantly involved with constant undercoolings, it is an important point to take away from this study.

## 4. Full time-dependent problem

So far, we have analysed similarity solutions which require specific far-field solutions. In order to investigate general boundary and initial conditions, as well as transient behaviour, we solve the full time-dependent equations (2.1), (2.8), (2.11), (2.14) for *r*≤λ*R* and *p*_{s}=*p*_{l} for *r*>λ*R*, together with the tip-extent condition (2.15) and the tip-propagation condition (2.4). This allows us to study not only the long-time behaviour of the crack, where the process has developed and the pressure build-up is causing the cavity to propagate, but also the initial stage of ice growth. This is characterized by the stress at the tip *K*_{I} being below the critical value *K* and no fracturing occurring. During this stage, the cavity thickness increases as a result of the ice formation and subsequent pressure build-up, while the stress at the tip *K*_{I} approaches the critical value *K*. The ice extent λ is a complicated function of the undercooling of the surrounding rock and the thickness of the cavity, and this framework allows for it to vary with time.

We also find that the similarity solution can be reproduced by the time-dependent problem with the appropriate temperature boundary condition, which serves as a useful check for our numerical scheme.

### (a) Initial condition

We first consider the initial stages of ice formation and pressure build-up. A simple way of thinking about how a situation like this can develop in nature is to imagine a pre-existing fault with some ice growing inside it. Initially, the stress at the tip is sub-critical and hence no fracturing occurs. Instead, water keeps freezing inside the non-propagating cavity, causing the cavity to increase in thickness. This in turn increases the pressure applied on the rock, and subsequently the stress at the tip. When the critical value is reached, the fracturing begins.

As we saw in the study of the similarity solution, depending on the properties of the porous medium and the undercooling of the surroundings, cases exist in which the pressure induced by the ice on the rock is not enough to make the cavity propagate. If the cavity is completely ice-filled, the ice growth is limited by the back pressure from the rock. In that case, the maximum pressure has been achieved but the induced stress at the tip is not large enough to cause fracturing. In cases of small undercoolings or very inelastic rocks (hence very thin cavities), where the extent of the ice is dictated by the curvature-melting at the tip, the equilibrium is reached with the ice extent λ<1. An example can be seen in figure 7 where the ice has reached the maximum extent it can achieve and the back pressure from the rock on the ice is preventing any further freezing. The stress at the tip for this example is about half the value of the critical fracture toughness *K* and therefore no propagation occurs in this case. This effect is independent of the initial cavity thickness used and corresponds to the ‘no propagation’ regions discussed in the similarity solution section.

### (b) Results

Faster propagation occurs for greater undercooling of the surroundings, which can be seen in the growth curves *R*(*t*) versus *t* in figure 8. We see that the change in temperature affects both the initial stage of pressure build-up inside the cavity and the propagation rate of the crack. For smaller undercoolings, it takes longer for the stress at the tip to reach the critical value and for the propagation to begin. We note that the minimum value of the undercooling Δ*T* for propagation to occur is just below 8 K for this example. The maximum disjoining pressure acting between the ice and the rock through the premelted film is a linear function of Δ*T* given by

An increase of the critical fracture toughness *K* has a similar effect to a decrease of Δ*T*, slowing propagation as more water must freeze to maintain the critical state *K*_{I}=*K*. There is a linear relationship between the temperature and the disjoining pressure, demonstrated by equation (4.1), as well as the stress intensity at the tip, shown by equation (2.3). This implies that the ratio *K*/Δ*T* of the fracture toughness to the undercooling is important in determining the propagation rate.

The maximum solid pressure on an ice-filled crack is given by the disjoining pressure as
*p*_{l} is negative during ice growth. This means that the maximum possible pressure is dependent on the undercooling. The stress intensity factor at the tip is given as an integral of the liquid pressure over the crack and its maximum value comes from the uniform *p*_{s}(*r*,*t*)=*p*_{T} pressure distribution, such that
*K*_{I}<*K* so, for a given undercooling Δ*T*, we need
*p*_{T} to be applied along the whole crack surface. If the ice extends to only λ*R*<*R*, then the condition becomes
*K*=0.87 MPa m^{1/2} for undercooling of Δ*T*=15 K, we find that *K* and Δ*T*. In this case though, there is no dependence on the permeability since the undercooling is constant.

We are also interested in how the initial size of the cavity affects the propagation rate for initial radii *R*_{0}. The curves for *R*_{0}=1 cm and *R*_{0}=2 cm are shifted in time to match the *R*_{0}=0.5 cm curve at *R*(*t*)=1 cm and *R*(*t*)=2 cm, respectively. We see that the smaller the initial radius the less time it takes for the critical stress condition to be reached and the propagation to start. This does not contradict the *K*_{I}=*K*.

It is interesting that the three curves appear to coincide. This tells us that the propagation rate at any time is only dependent on the cavity radius at that time and not on the previous evolution of the crack. Hence, a crack that started from an *R*_{0}-sized fault and one that started from a 2*R*_{0}-sized fault will propagate in an almost identical way if we ignore the initial time it takes for the smaller cavity to reach *R*(*t*)=2*R*_{0}. This conclusion means that the initial condition we use only has an important effect on the fracturing potential of the crack rather than the rate of fracturing.

As we saw from the similarity solution, the permeability limits the flow of water towards the freezing front and hence the rate of solidification. In a warming environment, this results in a change in the fracturing potential. If it takes a long time for the ice and the pressure in the cavity to build up (i.e. small permeability), then the undercooling is reduced and hence the ability of the ice to fracture the cavity is limited. When the undercooling of the medium is not time-dependent, the permeability will simply affect the propagation rate rather than the fracturing potential.

We plot growth curves for different values of

## 5. Discussion

We now apply our model to some typical rocks to obtain an understanding of the physical time-scales of fracturing in natural settings.

Table 1 lists the values of the physical constants relevant to water and ice, which are used for all the examples we consider. Typical values for the fracture toughness *K*, the shear modulus *G*, Poisson's ratio *ν*, the porosity *ϕ* and the permeability *Π* for different types of rocks and clays are presented in table 2. Considerable variation in these physical parameters exists within each type of medium, as is clear by the range of values reported in the literature. For example, permeability is likely to vary by roughly 1–2 orders of magnitude [29], while porosity by 0.2–0.3 [30] within each rock type.

The values given for the permeability of the different media so far assume that there is no pore ice, or that its effect on the flow of water through the medium is negligible. An expression for an undercooling-dependent permeability is discussed at the end of this section.

Figure 11 shows the growth of a cavity of initial radius *R*_{0}=1 cm in limestone subjected to an undercooling of Δ*T*=10 K. The pressure of the ice is also plotted. During the initial stage of ice growth, the pressure rises rapidly until the critical fracture toughness is reached. At that point, the fault begins to fracture and the pressure drops quickly. This is consistent with equation (2.3) and the fact that the integral of the pressure over the crack surface needs to remain constant for quasi-steady propagation.

We note that the rate of fracturing is fast, of the order of several cm min^{−1}. In general, the time-scale of propagation is heavily influenced by parameters of the problem such as the permeability, the elastic modulus and the initial cavity size according to

The values given in table 2 are representative of unfrozen rocks. So far, we have assumed that little ice exists in the pores and it does not affect the permeability. Clearly, a change in the permeability has a strong effect on the timescale as they are inversely proportional, and hence it is interesting to study the existence of pore ice further. We expect the change in permeability to have an effect on the timescale of pressure build-up and propagation, but not on the fracturing potential of a fault, at least for constant undercoolings.

Previous studies of rock fracturing [14,4] have assumed the existence of pore ice. Physically, when a water-saturated porous medium is subjected to sub-zero temperatures, water inside the bigger pores can freeze, restricting the hydraulic connectivity of the medium. The extent of the freezing depends on the freezing-point depression, which is a result of curvature melting or the disjoining forces between rock and ice. Quantifying this effect is difficult, as pore size varies from material to material, and there can be a large range of pore diameters within one material. For example, while most limestones have pore sizes of the order of 0.1−1 μm, marbles and lithographic limestones can have an average pore size of approximately 0.005 μm [31]. The critical nucleation radius for an undercooling of Δ*T*=1 K is approximately 0.1 μm and, while inversely proportional to Δ*T*, we see that even for low temperatures, curvature melting can stop ice from growing inside porous material.

We can approximate the effect of pore ice on permeability by introducing a temperature-dependent effective permeability *Π*_{ΔT}. We write an empirical rule for this of the form
*c* is dependent on the properties of the specific medium, and could be determined experimentally (see the appendix for more details). For the following example, we consider a granite with unfrozen permeability of *Π*=10^{−15} cm^{2}. We take *c*=0.33 K, corresponding to an average permeability of *Π*_{ΔT}=10^{−16} cm^{2} at an undercooling of Δ*T*=1 K, which is in agreement with the values used by Walder & Hallet [14].

We saw in §3 that the growth rate of a crack increases with the undercooling Δ*T*. However, if we assume the existence of pore ice and approximate the undercooling-dependent permeability of the medium by expression (5.2), the effect of an increased Δ*T* will be more complicated. While a large undercooling speeds up fracturing, it also reduces the permeability of the medium. A smaller permeability *Π* means that flow through the porous medium is more restricted, hence the propagation slows down. These two effects act in opposing ways and the balance of the two determines how an increase in Δ*T* changes the growth curves. This effect can be seen in figure 12, where we have plotted the instantaneous growth rate of a crack both after 2 h of freezing and 20 h of freezing. We discover a similar behaviour to that found by Walder & Hallet [14], with the most effective fracturing occurring at temperatures higher than −5°C and the fracturing rate decreasing for large values of the undercooling, as there is more pore ice restricting the flow of the water. The rates of freezing are of the order of 10^{−4} mm s^{−1} which is in agreement with the 10^{−4} to 10^{−5} mm s^{−1} rates found by Walder & Hallet [14].

It should be pointed out that the effect of the permeability through the undercooling is dependent on the value of *c*. We have used *c*=0.33 K here but a larger value would result in the effect of *Π*_{ΔT} being weaker and we could see the propagation rate increasing with Δ*T* even with the decreasing permeability. Experimental work could provide further useful insight here, in particular in investigating the existence of pore ice, determining *c* and validating expression (5.2).

## 6. Conclusion

In conclusion, this study of a penny-shaped cavity has provided us with a mathematical model capable of explaining the physical processes governing fracturing, and establishing the relevance of different parameters.

Phase planes of the different types of solutions were discussed in §3, and general patterns of dependence of the cavity growth curves on the different parameters of the problem were established. The ratio of fracture toughness to undercooling not only determines the susceptibility of a medium to fracturing but also strongly affects the propagation rate, with a larger *K*/Δ*T* ratio resulting in slower propagation. The growth rate is also proportional to the elastic modulus of the medium, hence very inelastic rocks are expected to fracture quickly.

We also established that when the surrounding environment is undercooled to a constant Δ*T*, the permeability of the medium does not affect the fracturing potential, but only the timescale of cavity growth. In addition, the permeability must be smaller than about 10^{−12} cm^{2} in order to have an important effect on either the fracturing potential or the fracturing rate.

We found a relation between the key parameters of the problem, given by equation (4.4), which needs to be satisfied for fracturing to occur. For example, for limestone undercooled by Δ*T*=5 K, no faults smaller than *R*_{0}≈1.2 cm will fracture. By contrast, a 1 cm-sized fault in granite will require undercoolings of at least 13 K to crack.

In the case where the condition above is not met, ice can still grow inside the cavity and deform it without fracturing it. The process continues, with the thickness of the cavity increasing, until an equilibrium is reached where the disjoining pressure between the ice and the rock is balanced by the elastic response of the medium.

This discussion assumes that several flow paths to the fault exist, a reasonable assumption for initial fault sizes of millimetres which, for most rocks, is the order of magnitude required for fracturing to occur at undercoolings of a few degrees. For small faults in rocks with low flow-path connectivity (which will be reflected in the low permeability of the medium via equation (A 1)), different mechanisms might develop, especially at the initial stages of freezing. If the flow paths to a pre-existing fault are restricted, it could instead be the expansion of the freezing water causing large pressures to develop, as discussed in [18].

The discussion of the conditions for the existence of a frozen fringe and an undercooling-dependent permeability, as well as its effect on the fracturing, shows the variety and complexity of factors that influence the process. The time-dependent model can also incorporate different background conditions, including time-varying far-field undercoolings. Another assumption essential for the formulation of the integral equations is the spherical symmetry of the pressure fields, with the background conditions applied as

## Data accessibility

The results of this paper are derived from a complete set of mathematical equations given in the paper using a custom numerical algorithm described in the thesis by I.V. [22] and, as such, are fully reproducible.

## Authors contributions

I.V. conducted this research towards a PhD, supervised by M.G.W. Both authors gave final approval for publication.

## Funding statement

This research was supported by a PhD studentship from the EPSRC.

## Conflict of interests

We have no competing interests.

## Appendix A. The dependence of permeability on undercooling

The permeability is closely related to the size of the pores, and can be shown [29] to be proportional to the square of the average pore diameter
*C* is a constant relating to the flow paths connecting the pores. If we consider the growth of ice in a small spherical pore, two effects contribute to the depression of the freezing temperature: the curvature-induced melting, and the disjoining forces which become important when the pore is almost ice-filled, as explained in §2*d*. It can be shown that, at first order, the smallest depression of the freezing temperature is determined by curvature melting (dictated by the Gibbs–Thompson relation). In that case, the ice is close enough to the rock for the curvature to be approximately given by *κ*=1/*d*, where *d* is the diameter of the pore, but not close enough for the disjoining forces to be important, hence causing further depression of the melting temperature.

For a given undercooling Δ*T*, pores that are larger than the corresponding diameter
*d* as well as the freezing pore diameter *d*_{ΔT}. When Δ*T*=0, hence *d*_{ΔT} becomes smaller and hence some of the pores fill with ice. This reduces the average pore size. This analysis implies a relationship between *d* and *d*_{ΔT} of the form
*α*, to be determined, relates to the distribution of pore sizes of the medium.

The effective permeability *Π*_{ΔT} at an undercooling Δ*T* is proportional to *Π*_{ΔT} in terms of the undercooling Δ*T* and the unfrozen permeability of the medium (figure 13)
*c*_{1} and *c*_{2} have units of K^{2} and K, respectively, and depend on the medium characteristics, such as the pore size distribution. We note that the frozen permeability is inversely proportional to Δ*T*^{2}, as in [14]. One of the parameters *c*_{1} and *c*_{2} can be easily fixed by the condition that *Π*_{ΔT}=*Π* for Δ*T*=0, from which we deduce that *c*=0.1 K, an undercooling of just 1°C will cause a permeability 100 times smaller than the permeability of the unfrozen medium, while *c*=0.33 K means that, at Δ*T*=1 K, the permeability is 10 times smaller than the unfrozen value.

## Footnotes

↵1 While we have denoted this as λ=1, this would require the curvature of the tip of the ice to be equal to the curvature of the tip of the crack. In the model we have used here, this is represented by a parabolic tip, which results in a curvature several orders of magnitude larger than the approximate circular curvature of the ice close to the tip. This indicates that the cavity can never be completely ice-filled. Indeed, the section we have labelled as such actually represents solutions for which λ>1−1/

*N*where*N*is the number of panels that the cavity has been split into to solve numerically, i.e. the ice extends inside the last panel. This has a negligible effect on the results, especially for large*N*, and the boundary between λ<1 and λ=1 is useful in representing contours of constant λ. The assumption of a parabolic tip is only approximate and breaks down when within a few nanometres of the tip [23]. In reality, the tip is more likely to be sharp, making the curvature at that point infinite. The details of what happens at such small scales close to the tip are beyond the scope of this paper.

- Received September 30, 2014.
- Accepted January 9, 2015.

© 2015 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.