Freeze fracturing of elastic porous media: a mathematical model

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.


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 3 rspa.royalsocietypublishing.org Proc. R. Soc

R(t)
z r 2b(r, t) ice water porous medium l(t) Figure 1. A sketch of a penny-shaped crack propagating in a porous rock, where the tip is at r = R(t) and the ice extends to r = λ(t)R(t). The thickness of the cavity is 2b(r, t).
parameters influencing frost fracturing and interpret experimental observations. It is worth noting that several experimental studies [15][16][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.

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 ∞ < T m , where T m is the melting temperature of the ice measured at the reference pressure p m = p ∞ .
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 ( 2.2) and the functions 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 Linear elasticity predicts that this stress causes the cavity to fracture when it reaches a critical value 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 which is a condition on the shape of the tip when the crack propagates [21].

(b) The water pressure distribution
The flow of water through the pores of the rock is described by Darcy's equation where μ is the dynamic viscosity of water and Π 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 everywhere outside the crack. There is a water flux 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 The liquid pressure field in the thin premelted film can therefore be written as 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 where ρ s is the density of ice, L the latent heat per unit mass and k l the thermal conductivity of water. We define the undercooling θ (r, z, t) = T(r, z, t) − T ∞ (t) > 0 which decays as r → ∞ and satisfies the same equations as 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 r → ∞, we see that the interface undercooling θ I can be expressed in a similar way to p l as an integral of ∂b/∂t We note that in the special case λ = 1, the temperature field is simply proportional to the liquid pressure field 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 (2.14) Here, we have ignored the effect of pressure melting (proportional to 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 where γ is the surface energy of water per unit area.

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 T m − T ∞ = T = Ht −1/4 then there is a similarity solution in which the cavity propagates as for some constant 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 and 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 ( 3.8) In this new formulation, quasi-steady propagation occurs for 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 stressintensity 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 ofK,H andΠ . We investigate the existence of propagating solutions by plotting phase planes ofK versusH for a given value of the permeabilityΠ. We concentrate on the stable (fast) solution, although similar phase planes    are found for the unstable (slow) solution. Figure 4 shows two such phase planes, forΠ = 0.02 andΠ = 200. In general, we find four different solution regimes. Two of these describe fracturing cavities, while the other two represent situations when no propagation occurs, but the reason behind this is different in each case. These four regions are also summarized in the sketch in figure 5.
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 ifK decreases or if H increases. A lower value ofK corresponds to either a weaker medium, which requires less pressure to fracture, or a less compliant medium, in which the pressure build-up takes less time as the cavity does not deform much and hence less ice is required to maintain the pressure. The process takes, the smaller the disjoining pressure. An increase of the undercoolingH means more ice build-up inside the cavity, which can overcome the toughness of the medium.
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 (H,K) for which we can find sets of solutions (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 byK where β 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 3/2 . Substituting these relations in the pressure balance we find This relation gives the propagation rate k in terms of the dimensionless parameters of the problem. The left-hand side of the expression has a minimum at k = 1/ 3(1 +Π ). Hence, if the value ofHΠ 1/4 /K is smaller than that minimum, no solution for k exists. This threshold 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.  For largeΠ , we see that β(Π) → 1, which shows that the permeability does not have a strong effect on the propagation boundary. For small permeabilities, we see that the effect ofΠ on β becomes important (β ∼Π 1/4 ) and hence the fracturing potential of cavities is limited.

(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 asK increases, as shown in figure 6a. There is a similar explanation for why λ increases withH, since a colder cavity causes more water to freeze, resulting in a wider cavity and the ice being able to extend further towards the tip of the crack. 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 6b where, for a givenK, the propagation rate is seen to increase withH. The propagation rate also increases as the fracture toughnessK decreases: for a more brittle rock, less pressure build-up is required inside the cavity for the tip condition to be met, resulting in faster propagation.
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.

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 This shows how the maximum pressure exerted by the ice on the rock is limited by the undercooling. As the critical value of the stress at the tip is reached and the crack extends, more water must freeze to fill the now larger cavity, and build up the stress again. At lower  temperatures, the solidification is faster, meaning that a faster cavity propagation rate can be sustained. 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 as the liquid pressure 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 We can intuitively guess that rocks with smaller pre-existing faults will be less prone to fracturing at a certain undercooling. There is no propagation for K I < K so, for a given undercooling T, we need for the initial cavity to start fracturing. Expressing the stress at the tip as the integral of the pressure over the crack shows that, if the pressure has a maximum value, the stress intensity factor will not reach the critical value when the radius of the crack is too small. As mentioned before, this estimate corresponds to completely ice-filled cracks, since it requires the disjoining pressure p T to be applied along the whole crack surface. If the ice extends to only λR < R, then the condition becomes Of course, the factor λ is a complicated function of the undercooling and the rock properties, and this criterion becomes less straightforward. We can view R min as a lower bound for the minimum radius for propagation but remember that not all faults of initial radius greater than R min are guaranteed to fracture. For the example presented in figure 9, which is a limestone with K = 0.87 MPa m 1/2 for undercooling of T = 15 K, we find that R min ≈ 0.21 cm. This analysis agrees with our conclusions from the similarity solution and demonstrates again the linear relationship between 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 > R min . In figure 9, we have plotted growth curves for three different values of 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 R min conclusions above; the maximum pressure in a cavity is reached faster when the cavity is smaller as less freezing is required. If R 0 < R min , the tip stress caused by this maximum pressure is not enough to fracture the cavity. For R 0 > R min , the propagation will occur as soon as 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 2R 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) = 2R 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Π in figure 10. As expected, the permeability has a big effect on both the timescale of the initial phase as well as the growth rate of the cavity. This effect is larger for very small permeabilities: in theΠ = 0.01 case (solid curve), the build-up time is about five times longer than in theΠ = 0.1 case (dashed curve) and the growth rate about five times slower. By contrast, theΠ = 1 (dotted curve) build up time is just twice longer than thẽ Π = 10 one (dashed-dotted curve) and the growth rate only slower by a factor of 1.5. From this, we deduce that for media with very low permeability, such as granite, the effect ofΠ is especially strong and dominates the timescales of the problem. As the medium becomes more permeable, the resistance to the flow of water is less strong and we see that forΠ greater than about 10, it becomes negligible.

Discussion
We now apply our model to some typical rocks to obtain an understanding of the physical timescales 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    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 timescale of propagation is heavily influenced by parameters of the problem such as the permeability, the elastic modulus and the initial cavity size according to which results in timescales of the order of seconds or minutes for rocks, or up to hours for clays.   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 [4,14] 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 temperaturedependent effective permeability Π T . We write an empirical rule for this of the form which is based on the fact that the permeability of the medium is proportional to the square of the average pore diameter [29]. The parameter 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).

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 flowpath 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 undercoolingdependent 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 r → ∞. Different geometries or finite boundaries could be studied by using the appropriate solutions to the liquid and solid pressure fields. Our study provides a framework for modelling the important physical interactions of frost fracturing, and a base on which to incorporate further effects that are found to be important.