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.
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 .
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%  or 80% . 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 . 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  simply assumed that the water in the thin films at solid–liquid interfaces experiences an attractive force from the solid boundary. Vignes-Adler  suggested that water in thin films has different properties to water in bulk, but experimental work by Raviv & Klein  does not support this suggestion. A model was developed by Rempel et al.  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  and Murton et al.  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  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  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  but formalize the discussion of physical processes using the premelting theory presented in . 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 , where Tm 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  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.1 where m=G/(1−ν), G is the shear modulus and ν=λ/[2(λ+G)] is the Poisson ratio of the rock. The kernel M is given by 2.2 and the functions Kell and Eell 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 KI given by 2.3 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 KI=K, equation (2.3) together with the expression for the pressure field (2.1) are equivalent to 2.4 which is a condition on the shape of the tip when the crack propagates .
(b) The water pressure distribution
The flow of water through the pores of the rock is described by Darcy's equation 2.5 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 2.6 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 , hence we can express the flux q(r,t) simply in terms of the rate of opening of the cavity 2.7 The liquid pressure field in the thin premelted film can therefore be written as 2.8
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 2.9 Latent heat is released at the solidification front, which is described by the Stefan condition 2.10 where ρs is the density of ice, the latent heat per unit mass and kl the thermal conductivity of water. We define the undercooling which decays as and satisfies the same equations as T, (2.9) and (2.10). By comparing equations (2.7) and (2.10), and noting that both pl(r,t)=p(r,z≈0,t) and θI(r,t)=θ(r,z≈0,t) satisfy Laplace's equation and decay as , we see that the interface undercooling θI can be expressed in a similar way to pl as an integral of ∂b/∂t 2.11 We note that in the special case λ=1, the temperature field is simply proportional to the liquid pressure field 2.12
(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 pT is called disjoining pressure. The liquid pressure pl inside the thin film, together with the disjoining pressure pT acting across the film, must balance the solid pressure ps in the ice  2.13
Owing to the disjoining pressure pT, the liquid and solid pressures across the freezing interface are not equal. This pressure difference affects the freezing temperature TI as described by the Gibbs–Thompson relation 2.14 Here, we have ignored the effect of pressure melting (proportional to pl−pm) 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 2.15 where γ 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 then there is a similarity solution in which the cavity propagates as 3.1 for some constant k to be determined. In terms of the similarity variable η=r/kt1/2, the tip of the cavity is at η=1, the ice extends to η=λ and the cavity thickness is B(η)=b(r,t)/kt1/4. Then, the similarity functions ps(η)=t1/4ps(r,t), pl(η)=t1/4pl(r,t) and ΘI(η)=t1/4θI(r,t) satisfy the following system of equations 3.2 3.3 3.4 3.5 3.6 3.7 where
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 k1 and k2, with k1<k2, 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, ps(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 pl 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 k1<k<k2 the stress intensity factor is higher than its critical value, and for k>k2 or k<k1 it is lower. Therefore, if the cavity propagates just faster than the slow solution, at k=k1+ϵ (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 k1, we find K*<1 therefore the propagation slows even further. Conversely, for a cavity propagating at k=k2+ϵ, 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 (k1,λ1) and (k2,λ2) cannot be found for all values of , and . We investigate the existence of propagating solutions by plotting phase planes of versus 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 and . 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 if decreases or if increases. A lower value of 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 timescale of pressure build-up is important as the undercooling reduces with time. The longer the process takes, the smaller the disjoining pressure. An increase of the undercooling 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. λ=11), 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 for which we can find sets of solutions (k1,λ) and (k2,λ) 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 3.9 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 k3/2 and the temperature field as . Substituting these relations in the pressure balance we find 3.10 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 . Hence, if the value of is smaller than that minimum, no solution for k exists. This threshold is the linear boundary described by equation (3.9) with a slope β. Substituting and in equation (3.10), we find that 3.11 For large , we see that , 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 () 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 ps(r,t) on the inner walls of the crack can be expressed as an integral over ps, 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 increases, as shown in figure 6a. There is a similar explanation for why λ increases with , 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 given , the propagation rate is seen to increase with . The propagation rate also increases as the fracture toughness 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.
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 ps=pl 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 KI 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 KI 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.
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 4.1 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 KI=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 4.2 as the liquid pressure pl 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 ps(r,t)=pT pressure distribution, such that 4.3 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 KI<K so, for a given undercooling ΔT, we need 4.4 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 pT to be applied along the whole crack surface. If the ice extends to only λR<R, then the condition becomes 4.5 Of course, the factor λ is a complicated function of the undercooling and the rock properties, and this criterion becomes less straightforward. We can view as a lower bound for the minimum radius for propagation but remember that not all faults of initial radius greater than are guaranteed to fracture. For the example presented in figure 9, which is a limestone with K=0.87 MPa m1/2 for undercooling of ΔT=15 K, we find that . 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 . In figure 9, we have plotted growth curves for three different values of R0. The curves for R0=1 cm and R0=2 cm are shifted in time to match the R0=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 conclusions above; the maximum pressure in a cavity is reached faster when the cavity is smaller as less freezing is required. If , the tip stress caused by this maximum pressure is not enough to fracture the cavity. For , the propagation will occur as soon as KI=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 R0-sized fault and one that started from a 2R0-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)=2R0. 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 case (solid curve), the build-up time is about five times longer than in the case (dashed curve) and the growth rate about five times slower. By contrast, the (dotted curve) build up time is just twice longer than the 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.
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 , while porosity by 0.2–0.3  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 R0=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 5.1 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 [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 . 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 5.2 which is based on the fact that the permeability of the medium is proportional to the square of the average pore diameter . 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 cm2. We take c=0.33 K, corresponding to an average permeability of ΠΔT=10−16 cm2 at an undercooling of ΔT=1 K, which is in agreement with the values used by Walder & Hallet .
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 , 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 .
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).
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 cm2 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 R0≈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 .
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 . 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.
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.  and, as such, are fully reproducible.
I.V. conducted this research towards a PhD, supervised by M.G.W. Both authors gave final approval for publication.
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  to be proportional to the square of the average pore diameter A 1 where 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 §2d. 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 A 2 will be almost completely ice-filled, reducing the overall average pore diameter of the medium to . We expect this new average diameter to be a function of both the initial average pore diameter d as well as the freezing pore diameter dΔT. When ΔT=0, hence , no ice grows inside the pores and therefore the average diameter remains unchanged, . As the temperature drops, 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 A 3 where α, to be determined, relates to the distribution of pore sizes of the medium.
The effective permeability ΠΔT at an undercooling ΔT is proportional to , according to equation (A 1). We can then express the permeability of the medium ΠΔT in terms of the undercooling ΔT and the unfrozen permeability of the medium (figure 13) A 4 The parameters c1 and c2 have units of K2 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 ΔT2, as in . One of the parameters c1 and c2 can be easily fixed by the condition that ΠΔT=Π for ΔT=0, from which we deduce that . To fix the second parameter, we need to determine the effective permeability at a given undercooling. For example, if 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.
↵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 . 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.