A thermomechanical model for frost heave and subglacial frozen fringe

Abstract Ice-infiltrated sediment, or frozen fringe, is responsible for phenomena such as frost heave, ice lenses and metres of debris-rich ice under glaciers. Understanding frozen fringes is important as frost heave is responsible for damaging infrastructure at high latitudes and sediment freeze-on at the base of glaciers can modulate subglacial friction, influencing the rate of global sea level rise. Here we describe the thermomechanics of liquid water flow and freezing in ice-saturated sediments, focusing on the conditions relevant for subglacial environments. The force balance that governs the frozen fringe thickness depends on the weight of the overlying material, the thermomolecular force between ice and sediments across liquid premelted films and the water pressure required by Darcy flow. We combine this mechanical model with an enthalpy method that conserves energy across phase change interfaces on a fixed computational grid. The force balance and enthalpy model together determine the evolution of the frozen fringe thickness and our simulations predict frost heave rates and ice lens spacing. Our model accounts for premelting at ice–sediment contacts, partial ice saturation of the pore space, water flow through the fringe, the thermodynamics of the ice–water–sediment interface and vertical force balance. We explicitly account for the formation of ice lenses, regions of pure ice that cleave the fringe at the depth where the interparticle force vanishes. Our model results allow us to predict the thickness of a frozen fringe and the spacing of ice lenses in subaerial and subglacial sediments.


Introduction
Freezing of interstitial water in sediments commonly occurs in subaerial and subglacial environments, contributing to effects such as frost heave, needle ice, and the transport of subglacial debris (Hemming, 2004;Dash et al., 2006;Wettlaufer and Worster, 2006).Through multiple cycles of freeze and thaw in high latitude environments, patterns can develop such as the arctic stone circles (Kessler and Werner, 2003).In this paper, we consider the thermodynamical and fluid dynamical processes that occur as water freezes in a porous medium.We describe the melting and freezing processes using an enthalpy formulation, which facilitates our numerical method as we avoid tracking phase change interfaces and elucidates the role of the frozen fringe as a mushy zone between ice-and water-saturated sediments.Our treatment is general enough to apply in a variety of industrial and environmental contexts where interstitial freezing occurs, yet here we primarily focus on geophysical applications.
Consider frost heave, the common freeze/thaw phenomenon that takes place throughout high latitudes.As water held within sediments freezes, ice lenses cleave the sediment and expand, causing vertical displacement of the ground surface.Such surface displacement causes significant damage to infrastructure at high latitudes and is due to the growth of distinct ice lenses within the soil rather than the water density change on freezing.Taber (1930) demonstrated this key fact by freezing a sediment pack saturated with benzene, which contracts on freezing; the benzene produced significant heave through the expansion of discrete ice lenses.
Early models for frost heave relied on surface tension to draw water to the lowest ice lens, i.e. the so-called "primary model for frost heave".This model suffers from several deficiencies, most importantly that surface tension acts tangential to the ice surface and cannot provide the upward force to drive heave.In addition, there is no mechanism to form distinct ice lenses in primary heave, which led O' Neill and Miller (1985) to derive the "secondary model of frost heave," wherein a zone of partially ice saturated sediment extends below the lowest ice lens.Fowler and Krantz (1994) clarified the mathematical model for secondary frost heave and analysed an asymptotically reduced form of the model.Rempel et al. (2004) highlighted the role of premelting at the interface between ice and sediment grains.The disjoining pressure across the liquid between ice and sediment grains relates the local melting temperature to the vertical force balance.Fowler and Krantz (1994), on the other hand, choose the liquid pressure to be given as a function of soil water content as set by surface tension, reminiscent of a primary frost heave model.In what follows, we build on the Rempel et al. (2004) formulation, highlighting an alternate derivation, writing out the equations, and systematically reducing the equations asymptotically, similar to Fowler and Krantz (1994).
Field observations show several meters of frozen sediment are commonly attached to the base of glaciers, motivating efforts to understand glaciohydraulic supercooling (e.g.Röthlisberger and Lang, 1987;Lawson et al., 1998;Creyts et al., 2013) and frost heave.In the fastest flowing reaches of glaciers and ice sheets, sliding dominates glacier motion and the rate of sliding is tied to the temperature and water pressure at the glacier base -the same things that control the growth of a frozen fringe.Rempel (2008) treats the overlying glacier as a large lowest ice lens and predicts meters-scale frozen fringes below glaciers for typical parameters, in line with observations.Anderson andWorster (2012, 2014) analyse freezing colloid suspensions using a directional solidification experiments with aqueous suspensions of alumina particles.Anderson and Worster (2014) developed a model built on the Rempel et al. (2004) framework that include compaction and cohesion of the colloid suspension.The Anderson and Worster (2014) model assumes a steady state linear temperature profile throughout the experiment and boils down to a system of ordinary differential equations for the location of the compaction front and frozen fringe extent, reminiscent of Fowler and Noon (1993).Based on their model, Anderson and Worster (2014) describe a regime diagram showing the three primary freezing regimes observed in their experiments: periodic ice lenses, disordered ice lenses, and periodic ice banding.
Frozen fringes and freeze-thaw cycles are inherently problems of phase change and partial melting.In these types of problems, it can be valuable to solve the energy conservation equation in an enthalpy form rather than for the temperature to avoid explicitly tracking phase T f < l a t e x i t s h a 1 _ b a s e 6 4 = " n P J T b 1 M y T J U F 1 U 3 o o x 3 E e G B 9 m W M = " > A A A B 7 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L L a C p 5 L 0 o s e C I B 4 r 9 A v a U D b b T b t 0 s w m 7 E 6 G E / g g v H h T x 6 u / x 5 r 9 x 2 + a g r Q 8 G H u / N M D M v S K Q w 6 L r f T m F r e 2 d 3 r 7 h f O j g 8 O j 4 p n 5 5 1 T J x q x t s s l r H u B d R w K R R v o 0 D J e 4 n m N A o k 7 w b T u 4 X f f e L a i F i 1 c J Z w P 6 J j J U L B K F q p W 2 0 N s 3 B e H Z Y r b s 1 d g m w S L y c V y N E c l r 8 G o 5 i l E V f I J D W m 7 7 k J + h n V K J j k 8 9 I g N T y h b E r H v G + p o h E 3 f r Y 8 d 0 6 u r D I i Y a x t K S R L 9 f d E R i N j Z l F g O y O K E 7 P u L c T / v H 6 K 4 a 2 f C Z W k y B V b L Q p T S T A m i 9 / J S G j O U M 4 s o U w L e y t h E 6 o p Q 5 t Q y Y b g r b + 8 S T r 1 m u f W v M d 6 p X G f x 1 G E C 7 i E a / D g B h r w A E 1 o A 4 M p P M M r v D m J 8 + K 8 O x + r 1 o K T z 5 z D H z i f P 6 Q 1 j x s = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " n P J T b 1 M y T J U F 1 U 3 o o x 3 E e G B 9 m W M = " > A A A B 7 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L L a C p 5 L 0 o s e C I B 4 r 9 A v a U D b b T b t 0 s w m 7 E 6 G E / g g v H h T x 6 u / x 5 r 9 x 2 + a g r Q 8 G H u / N M D M v S K Q w 6 L r f T m F r e 2 d 3 r 7 h f O j g 8 O j 4 p n 5 5 1 T J x q x t s s l r H u B d R w K R R v o 0 D J e 4 n m N A o k 7 w b T u 4 X f f e L a i F i 1 c J Z w P 6 J j J U L B K F q p W 2 0 N s 3 B e H Z Y r b s 1 d g m w S L y c V y N E c l r 8 G o 5 i l E V f I J D W m 7 7 k J + h n V K J j k 8 9 I g N T y h b E r H v G + p o h E 3 f r Y 8 d 0 6 u r D I i Y a x t K S R L 9 f d E R i N j Z l F g O y O K E 7 P u L c T / v H 6 K 4 a 2 f C Z W k y B V b L Q p T S T A m i 9 / J S G j O U M 4 s o U w L e y t h E 6 o p Q 5 t Q y Y b g r b + 8 S T r 1 m u f W v M d 6 p X G f x 1 G E C 7 i E a / D g B h r w A E 1 o A 4 M p P M M r v D m J 8 + K 8 O x + r 1 o K T z 5 z D H z i f P 6 Q 1 j x s = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " n P J T b 1 M y T J U F 1 U 3 o o x 3 E e G B 9 m W M = " > A A A B 7 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L L a C p 5 L 0 o s e C I B 4 r 9 A v a U D b b T b t 0 s w m 7 E 6 G E / g g v H h T x 6 u / x 5 r 9 x 2 + a g r Q 8 G H u / N M D M v S K Q w 6 L r f T m F r e 2 d 3 r 7 h f O j g 8 O j 4 p n 5 5 1 T J x q x t s s l r H u B d R w K R R v o 0 D J e 4 n m N A o k 7 w b T u 4 X f f e L a i F i 1 c J Z w P 6 J j J U L B K F q p W 2 0 N s 3 B e H Z Y r b s 1 d g m w S L y c V y N E c l r 8 G o 5 i l E V f I J D W m 7 7 k J + h n V K J j k 8 9 I g N T y h H e n Y 9 l a 8 H J Z 0 7 h D 5 z P H 5 6 K j V Y = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " w Z t e t r H o e R F F V N w R P H e n Y 9 l a 8 H J Z 0 7 h D 5 z P H 5 6 K j V Y = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " w Z t e t r H o e R F F V N w R P H e n Y 9 l a 8 H J Z 0 7 h D 5 z P H 5 6 K j V Y = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " w Z t e t r H o e R F F V N w R P e X H e n Y 9 F a 8 H J Z 4 7 h D 5 z P H 5 0 F j V U = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " t X e X H e n Y 9 F a 8 H J Z 4 7 h D 5 z P H 5 0 F j V U = < / l a t e x i t > z < l a t e x i t s h a 1 _ b a s e 6 4 = " s a F F g 9 T e X H e n Y 9 F a 8 H J Z 4 7 h D 5 z P H 6 A P j V c = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " s a F F g 9 T e X H e n Y 9 F a 8 H J Z 4 7 h D 5 z P H 6 A P j V c = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " s a F F g 9 T e X H e n Y 9 F a 8 H J Z 4 7 h D 5 z P H 6 A P j V c = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " s a F F g 9 T e X H e n Y 9 F a 8 H J Z 4 7 h D 5 z P H 6 A P j V c = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 7 Y 7 R i g W j F j s 9 M q l water and sediments Figure 1: Schematic of the a frozen fringe system: (a) components including pure ice (lowermost ice lens or the bottom of a glacier), frozen fringe, and the unfrozen porous mixture of water-saturated sediments (after Rempel et al., 2004;Anderson and Worster, 2014).(b) surface integral path Γ and inward-pointing normal vector dΓ.(c) volume integral domain Ω with the outward-pointing normal vector dΩ.change interfaces.The enthalpy (i.e.sum of the sensible and latent heat) accommodates the phase change, which facilitates numerical solutions.From sea ice (Katz and Worster, 2008), permafrost (Clow, 2018) and meltwater percolation through snow (Meyer and Hewitt, 2017) to industrial processes (Voller and Prakash, 1987), the enthalpy approach to phase change problems is useful for many applications.Enthalpy methods have been used extensively for polythermal glaciers, where part of the glacier is below the melting point and the rest of the glacier is at the melting point, i.e. temperate ice (Aschwanden et al., 2012;Schoof and Hewitt, 2016).Here we use the enthalpy method to solve for energy conservation within a frozen fringe.
In this paper, we focus on the conditions relevant for subaerial frost heave and subglacial environments.Although numerous treatments of frost heave exist in the literature (e.g.O' Neill and Miller, 1985;Fowler and Krantz, 1994;Rempel et al., 2004), here we derive our model from scratch for completeness and clarity.In section 2, we start by writing down mass, momentum, and energy conservation equations for a frozen fringe.Then, we nondimensionalise and systematically reduce the equations by exploiting the small density difference between ice and water as well as the large latent heat of fusion upon freezing water (i.e. a large Stefan number).We solve our reduced model using an enthalpy method, where phase-change boundaries are determined implicitly on a fixed grid.In section 3, we demonstrate the results of our enthalpy model.We analyse a steady state frozen fringe thickness in melting and balanced thermodynamic conditions in both a semi-analytical model and an enthalpy framework.Then, we examine the local effective pressure for melting and freezing conditions, highlighting ice lens formation.Lastly, we show the formation of periodic ice lenses and map out the different behaviour in a regime diagram for the heave rate and effective pressure.Finally, we offer conclusions and discuss future directions in section 4.

Model
Inside the frozen fringe, which is shown schematically in figure 1, we define a coordinate system that is fixed with respect to the immobile, water-saturated sediment below, with z vertical, x the lateral coordinate, and y pointing into the page.We label the deepest extent of the fringe as z = z f and the top of the fringe as z = z , or equivalently z = z f + h, where h is the fringe thickness.

Mass conservation
The frozen fringe is partitioned into three components: ice, water, and sediment.The porosity φ denotes the volume of voids (i.e.ice and water) within a representative control volume.The fraction of the voids that is taken up by ice is the ice saturation S. Mass conservation for sediment, ice, and water implies that where m is the rate at which ice (density ρ i ) is melted, i.e. converted into liquid water (density ρ w ), V is the speed at which the ice moves through the fringe due to heaving at the ice lens above the fringe, and V s is the sediment (density ρ s ) velocity.The water flux through the fringe is U , which is given by Darcy's law as where the permeability k depends on the ice saturation as well as the porosity and other properties of the sediment matrix.Here g is the acceleration due to gravity.Values for all of the parameters are given in table 1.
In this paper, we assume that the porosity is constant throughout the fringe.Both below the frozen fringe and in its interior, there may be significant compaction due to the reduced water pressure as the lens pulls in water to freeze (Fowler and Noon, 1999;Anderson andWorster, 2012, 2014).We neglect such complications for now and take the entire sediment pack to maintain a constant porosity φ that jumps to φ = 1 at ice lenses.A new ice lens forms when the force between sediment grains reaches zero, as described in the next section.

Force balance
The force balance within the fringe is composed of three components: the weight of the material above and within the fringe, the water pressure within and below the fringe, as well as the thermomolecular force between sediment grains and interstitial ice, which acts across a thin film of premelted water (Dash et al., 2006).Integrating these components over the surface area of the fringe gives in which p i is the isotropic ice pressure (i.e.local normal stress) at the ice-water interface within the fringe, p w is the water pressure at the boundary, and n is the inward-pointing unit normal to the boundary Γ (cf.figure 1).The difference between the ice and water pressure, i.e. the first term on the right hand side of equation ( 5), is accommodated by a thermomolecular force (Rempel et al., 2001(Rempel et al., , 2004;;Wettlaufer and Worster, 2006).
We convert these surface integrals over Γ to volume integrals over the unfrozen component of the fringe Ω.We construct a closed surface by adding surface integrals with flat surfaces at z f and z, as shown schematically in figure 1.That is, for a generic pressure field p, we have the integrals where n is the outward-pointing normal for the volume Ω.In other words, n = k on the upper cap at z, the lowest ice lens, and n = −k on the lower cap at the bottom of the fringe z f , where k is the unit vector in the z-direction.
We take the surface Γ z to be the ice boundary at some height z, which we can write |Γ z | = (1 − φS)A where A is the cross-sectional area.This surface has the two limits of |Γ z | = 0 at the bottom boundary of the lowest active ice lens (φ = 1, S = 1) and Γ z f = A at the bottom of the fringe (S = 0).For this reason, no upper cap is necessary when integrating across the entire fringe, and equation ( 6) reduces to We assume that water flow through ice-saturated porous fringe is governed by Darcy's law and that the water pressure varies only on a lengthscale set by the fringe and not on the scale of individual grains.This assumption allows us to define the water pressure throughout the volume Ω, even though part of the domain is filled with ice and sediment.Crucially, we follow Rempel et al. (2004) and assume that the microscale pressure is the homogenised Darcy pressure.In other words, we treat the water pressure in the thin films between sediment and ice as well as the water pressure in the pore throats between sediment grains as determined by Darcy's law, which is a key difference between Fowler and Krantz (1994) and Rempel et al. (2004).
At this stage, we restrict our focus to a one-dimensional water pressure that only depends on the vertical coordinate z and assume that there are no transverse pressure gradients.Therefore, inserting the water pressure p w into equation ( 6), we find that where the porosity φ and saturation S can also depend on the vertical coordinate z.Now integrating across the entire fringe from z = z f to z = z gives The analogous equation for p i − p w represents the thermomolecular contribution to the force balance (Rempel and Worster, 1999;Rempel et al., 2004) and is given by When integrated across the entire fringe, the water pressure and thermomolecular force balance the total normal stress σ n at z f resulting from the weight of the overlying material (e.g.equation ( 5) and figure 1).With this in mind, we write where we have included the weight of fringe material as the integral over each constituent.
We now combine all of these equations including the effects of gravity and arrive at We define the effective pressure at the base of the fringe N as the total normal stress σ n at z f supported by the fringe less the water pressure at the base of the fringe p w (z f ) so that We recognise N as the load supported by contacts between sediment grains at z f .Additionally, we define the local effective pressure N loc (z) as the portion of the overlying load that is supported at a height z by sediment grain contacts.The rest of the overlying load is supported by thermomolecular forces or water pressure acting on the ice fringe below the height z as well as water pressure at height z.The thermomolecular and water pressure contributions from below z are given as We assume that sediment grains have infinitesimal contacts and, therefore, the water pressure at the height z supports the force A (1 − φS) p w (z)k, which excludes the areas occupied by ice.The total force supported by grain contacts at a height z is the overburden σ n minus both equation ( 14) and the water pressure at z. Thus, we can write the effective pressure N loc (z) as A new ice lens initiates at the height z n where the local effective pressure is zero, i.e.N loc (z n ) = 0, as there is no longer any force on the sediment grains (O'Neill and Miller, 1985;Rempel et al., 2004;Anderson and Worster, 2014).We treat the effective pressure at the bottom of the fringe N as an input to the model that is determined by groundwater hydrology or subglacial drainage (e.g.Schoof, 2010).

Generalized Clausius-Clapeyron and Gibbs-Thomson
The pressure difference between ice and water is related to temperature through the generalised Clausius-Clapeyron equation, which in its linearised form is given by where the bulk melting temperature at the reference pressure p m is T m , the specific latent heat of fusion for ice is L , and the densities of ice and water are given as ρ i and ρ w , respectively (Worster and Wettlaufer, 1999;Worster, 2000;Clarke, 2005).We choose the reference pressure to be the overburden σ n , so that at the bottom of the fringe z = z f , we have and in the interior of the fringe, we have Note that in this formulation changes to σ n affect the value of N and T m .At the bottom of the fringe, the ice is in contact with water in between pore throats, as shown in the schematic in figure 1.The curvature induced by the space between sediment grains, leads to a difference in the pressure in ice and water phases due to the Gibbs-Thomson effect (Worster, 2000), which is given by where γ is the ice-water surface energy and κ is the curvature of the ice-water interface.Moreover, the curvature κ at the bottom of the fringe is related to the radius of curvature for sediment pore throats r p as κ = 2/r p .The critical effective pressure required to overcome the pore throat curvature is given by the Gibbs-Thomson effect and defined by the right-hand side of equation ( 19), i.e.
Combining the generalised Claussius-Clapeyron equation ( 18) and the Gibbs-Thomson effect (19) at the bottom of the fringe, we can relate the curvature induced by pore throats to the temperature at the interface, which is given by The contribution from surface energy on the right hand side typically dominates the term proportional to the density difference, since ice and water densities differ by less than 10%.
Therefore, it is useful to define the undercooling temperature T f that supports this balance as which allows us to write equation ( 21) as Rempel ( 2008) drops the second term on the right arguing that it is small, which is consistent with our dominant balance above.We, however, keep all terms for now and reduce the model systematically in §2.7.2.
At this stage, we can now insert the generalised Clausius-Clapeyron equation ( 18) and the Gibbs-Thomson effect (19) into the effective pressure integrals ( 13) and ( 15).The effective pressure at the bottom of the fringe N is then In the interior of the fringe, the local effective pressure is given by which connects the local pressure and temperature.If there is not a fringe, the development of the force balance at the base of the fringe still holds, except that z f = z and the integrals vanish.Also, the curvature at the bottom of the ice is no longer κ = 2/r p and so the effective pressure is given by while the temperature at bottom of the ice is given by which is modulated by the effective pressure.

Energy conservation
Mass exchange between the liquid and solid phases within the fringe leads to changes in ice saturation.Conservation of energy determines the temperature and phase change within the fringe and can be expressed in terms of specific enthalpy, h α .For the constituents α that make up the fringe, sediment (s), ice (i), and water (w), conservation of energy in enthalpy form is given as where K e is the effective thermal conductivity (Appendix B of Rempel, 2008), which can be represented as for the thermal conductivities of the fringe constituents, sediment K s , ice K i , and water K w (Clauser and Huenges, 1995).Using the divergence theorem, we can write equation ( 28) as We define the difference between the water flux U and the heave rate V as u, i.e.
which will typically be small as it is the flow of water that allows for heave.Water flow is given by Darcy's law (4) as Thus, we combine the mass conservation equations ( 2) and (3) as which allows us to define the total ice and water W (e.g.Meyer and Hewitt, 2017) that is given by so that in the rigid-ice limit with ∇ • V = 0 where the heave rate is spatially independent, mass conservation (33) can be written succinctly as Similarly, for conservation of energy we can define the total enthalpy H as which is specified up to a constant, reference enthalpy H 0 , which we choose to make H = 0 at T = T f Thus, equation (30) reduces to Expanding out the total enthalpy (36) and using the fact that the specific latent heat L is equal to the difference between the liquid and ice specific enthalpies, i.e.L = h w − h i , we find that Figure 2: Schematic for the enthalpy below an ice lens: (a) Temperatures below T occur above the lowest ice lens; within the fringe, the temperature is tied to the ice saturation curve; and below the fringe, the water and sediment temperatures rise above T f .(b) Nondimensional version of (a) showing that the enthalpy is zero at θ = 0 and negative nondimensional temperatures correspond to unfrozen water with positive enthalpy.
which is a sum of sensible and latent heat contributions to energy.We leave equation ( 37) in enthalpy form to facilitate the description of the numerical method.For an incompressible medium, the specific enthalpy h α is equivalent to a change in temperature, i.e. dh α = c α dT , where c α is the specific heat capacity with α representing ice, liquid, or sediments.Thus, the enthalpy can be related to temperature as T f ≤ T (water and sediments) which we show schematically in figure 2. As we describe in the next section, the ice saturation in the fringe is a function of the temperature.

Constitutive relations for saturation and permeability
The ice saturation S and, therefore, the permeability k of the frozen fringe depend on the local thermodynamics with the pressure difference between the phases controlled by the Gibbs-Thomson effect and interfacial premelting (Andersland and Ladanyi, 2004;Hansen-Goos and Wettlaufer, 2010;Rempel, 2012).Here we use the generalized Clausius-Claperyon relation to specify the ice saturation S i and permeability k as functions of the local difference between ice pressure and water pressure.That is, we generalize the relationships used by Rempel (2007Rempel ( , 2008) ) and define the function Ξ as the ratio of the critical effective stress to the local pressure difference p i − p w , i.e.
This reduces to the expression used by Rempel (2007Rempel ( , 2008) ) if we ignore the density difference between ice and water.For now, we proceed with the definition in equation ( 40) and write the ice saturation and permeability as where the empirical exponents are β > 0 and α > β (typically α > 1; Rempel, 2008).

Boundary conditions
Now that we have specified the governing equations for enthalpy H and total water W , we describe the boundary conditions.At the top of the lowest fringe, i.e. z = z , a finite jump in saturation can occur as the ice lens is fully occupied by ice (φ = 1 and S = 1) whereas ice only partially saturates the interstices in the underlying fringe (φ < 1 and S < 1).
Integrating mass and energy conservation across the jump at the lowest ice lens boundary gives the following conditions Taking the heat flux into the ice lens and material above as q, we can simplify these conditions to where mass conservation implies that at the top of the fringe, all of the liquid water must freeze onto the lowest ice lens, and conservation of energy implies that all the heat that enters the fringe at the bottom must leave through the top.At the bottom of the fringe (or at any point below the fringe), we impose a conductive heat flux q which includes contributions from geothermal heat and friction from sliding, i.e.
which is the full heat flux as H = 0 at the base of the fringe (i.e.where we have chosen H = H 0 such that H = 0 when T = T f ).The ice saturation transitions smoothly from 0 < S < 1 within the fringe to S = 0 below the fringe and, therefore, no jump condition is required at z = z f .The water pressure p w at the bottom of the fringe is set by the effective pressure N and is given by

Nondimensionalisation
We now scale our model to find the dominant physical balances.For example, we write t = [t]t * for the time t, where [t] is the scale and t * is the nondimensional variable.Proceeding in this way, we write all variables as and choose the scales for the variables based on the expected physical balances.
A scale for the effective pressure within the fringe is the threshold entry pressure, i.e.
[N ] = N c , and a scale for the heat flux into the fringe is the geothermal heat flux, i.e.
[q] = q • n.We choose the temperature scale to be the temperature difference implied by premelting, i.e.
[N ] = A scale for the vertical distance [z] comes from the heat flux scale, i.e.
[q] = K The rate of heaving is determined by water percolation, so we choose and time can be scaled for the solidification as For the permeability, we choose the scale to be the prefactor as We also define the dimensionless variables where δ is the scaled density difference, ν is the ratio of the sediment density to water density, Pe is the Péclet number, Gr is the ratio of gravitational hydrostatic pressure to infiltration pressure, and St is the Stefan number.

Full model
We now write the model in nondimensional variables and for concision, we drop the asterisks.
Rewriting equation ( 24), the force balance across the fringe is if a fringe exists, or the effective pressure is constrained by N < 1 if there is not a fringe.The local effective pressure in the fringe is The total water in the fringe is and the enthalpy is θ ≤ 0 (water and sediments) Nondimensionalising equation ( 37) results in The scaled total water equation is given by which depends on both the flow of water u and the rate of heave V .The constitutive laws for permeability and saturation are written nondimensionally as Finally, the nondimensional boundary conditions are given as

Model reduction
Typical values for the nondimensional variables based on the parameters are given in table 1.The scaled density difference δ is a small value and therefore it is reasonable to neglect terms that are multiplied by δ (Rempel, 2008).Taking this limit, we find that the vertical force balance reduces to unless N < 1, in which case there is not a fringe.In the same way, the local effective pressure N loc (z) is to leading order in 1/St in the fringe.Only sensible heat terms persist in the pure ice and water and sediments and we retain the 1/St dependence to meet flux boundary conditions.The enthalpy variation in the frozen fringe is tied to the temperature θ through the ice saturation S as analogous to the liquidus condition in a mushy zone (Worster, 2000).
Given that frozen fringes are often much wider than thick, we now restrict our attention to one vertical dimension for conservation of mass and energy.Thus, in the same large Stefan number limit, the evolution equation for enthalpy is where we have neglected terms proportional to Pe/St as well.In the limit δ = 0, the total water W reduces to which is a constant, meaning that the pore space is entirely occupied by ice and water, yet there is no distinction in this limit due to the small density difference.Therefore, mass conservation implies The boundary conditions for mass, momentum, and energy conservation reduce to We now integrate equation ( 74) and impose the boundary condition ( 75), which implies that the water pressure gradient is We can now insert this water pressure gradient into the vertical force balance (68) to find that the heave rate V is given by as shown previously by Rempel (2008).This prescription of the heave rate is determined by force balance as well as conservation of mass and requires integrating the temperature field θ and the attendant saturation S. Thus, we can summarise our full model for the transient evolution of a frozen fringe as: enthalpy evolution (72) with heave rate (80) subject to boundary conditions ( 76) and (77).

Enthalpy numerical method
We write equation ( 72) in conservative form, defining the flux as the sum of the advective and diffusive components.We then discretise the conserved fluxes in space using a finite volume method implemented in python.In this numerical method, we divide the domain into cells and each variable is constant within a cell whereas velocities and fluxes are evaluated at cell edges.For advection, we use an upwinding scheme where the advective fluxes on cell edges are given by the cell values 'upwind', which is determined by the sign of the heave rate.We evolve explicitly equation (72) in time using solve ivp and the method of lines in python.With these choices, the finite volume implementation is conservative, meaning that the flux transferred between cells respects conservation of energy and phase change.The code is included in the supplemental information as well as in a github repository (link/doi to be added in proofs).
The two input parameters for the model are the effective pressure N and the heave rate V .Thus, coupling the governing equation ( 72) with the heave rate equation (80), this problem takes an integro-differential equation form, where at each timestep we integrate equation (80).Rather than implementing the top boundary condition (76) as a total flux, we apply a boundary condition to the diffusive part as where V input is the input value.The steady state is diagnosed when the value of V computed through equation ( 80) is equal to V input .We study the steady state problem in more detail in the next section.The outputs for the model are the frozen fringe thickness h = z − z f , the temperature profile through the fringe and the ice saturation profile in the fringe, as well as the initiation, timing, and spacing of ice lenses.

Steady state frozen fringe
To understand the development of a frozen fringe as well as the relationship between a free boundary representation and the enthalpy method, we start by considering a steady state frozen fringe with constant thermal conductivity.The problem is then: for a fixed location of the lowest ice lens z (e.g. the glacier-sediment interface), a constant heave rate V , and a known effective pressure N , what is the steady state temperature profile θ(z) and fringe-front location z f ?Conservation of energy in the fringe and the water-saturated region in front of it is given by where the enthalpy H in the fringe is given by ( 71) as H = −φS.
The boundary conditions are with the internal conditions dθ dz We start by deriving the temperature profile for the region below the fringe.By integrating equation ( 83), we have which satisfies the geothermal heat flux boundary condition (85) and the scaled temperature goes to zero at the bottom of the fringe to satisfy the internal condition (86).Now in the fringe, we integrate equation ( 82) once and apply the boundary condition (84), which results in which is the governing ordinary differential equation for the temperature profile in the fringe, subject to the boundary condition θ = 0 at z = z f .Thus, for a given heave rate V and effective pressure N , the temperature profile is specified by ( 89).The only piece of information that is missing, is the location of the bottom of the fringe z f , which we determine through force balance (80).To solve for θ and z f , we numerically integrate (89) for the temperature, insert the solution into the force balance (80), and find the fringe front z f using a root-finding algorithm (i.e.similar to a shooting method).The result of this numerical procedure with the parameters given in table 1 is shown as the red line in figure 3. The black line shows the solution to the same problem using the enthalpy method, where the domain runs from z = 0 to z = z .Here we set z = 1.We find the steady state through a relaxation method, i.e. we integrate equation ( 72) in time until the difference between the computed heave rate and the target value is less than 10 −3 .The enthalpy in the frozen fringe is negative, taking on its smallest value at the base of the lowest ice lens z and rising monotonically up to zero at the base of the fringe z f .The enthalpy is positive in the water-saturated sediments, yet is very small due to the large Stefan number St.The nondimensional temperature θ is shown in the right panel of figure 3.In line with T = T f − [T ]θ, we see that θ is positive in the fringe and negative below.The temperature profile is close to linear, which makes sense given that the heave rate V is small, and is exactly linear in the thermodynamically balanced case where V = 0 and when the thermal conductivity is constant.At the bottom of the fringe θ = 0 and the location z f is determined through force balance.The nondimensional fringe thickness is given as h = z − z f and we find h = 0.36 for the parameters in figure 3.
In figure 4, we show the dimensional fringe thickness h = [h](z − z f ) in meters as a function of the nondimensional effective pressure N/N c for balanced (V = 0, left) and temperate melting (V = −q/(ρ i L [V ]), right) thermodynamics using 5 different solution methods.The first two techniques are what we just described: 'ode' is the solution to (89) subject to (80) and 'enthalpy' is the conserved finite volume method for solving (72).In the balanced, V = 0 melting, V = −1.1 Figure 4: Compilation of benchmark solutions for dimensional fringe thickness h for two nondimensional heave rates (left) V = 0 and (right) V = −q/(ρ i L [V ]) = −1.1.The 'uega' and 'shooting' are solution methods from Rempel (2008).The 'enthalpy' solution is from the finite volume method, 'ode' is the steady state ode solution, and 'root' is the semi-analytical root-finding method.All methods the same result.thermodynamically balanced case, (80) can be integrated exactly and the fringe thickness can be found with a root-finding algorithm, which we name 'root' in figure 4. The last two methods 'uega' and 'shooting' are from Rempel (2008) and are run here using the parameters in table 1.In the uniform external gradient approximation, i.e. uega, method Rempel (2008) maps the fringe with imposed external heat fluxes to a Stefan problem domain whereas the 'shooting' method searches for a consistent temperature at the base of the lowest ice lens.As expected, all 5 of these solution methods give the same result.The fringe thicknesses are much lower in the temperate melting case because the liquid pressure distribution in the fringe needed to expel meltwater supports a larger portion of the overburden and the fringe melts as the ice infiltrates into the sediments.

Ice lens initiation
In the calculation of the steady state fringe thicknesses, we focused on melting (V < 0) and balanced (V = 0) thermodynamics.Although steady states do exist for relatively small freezing rates and relatively small effective pressures (Meyer et al., 2018), transient behaviour such as ice lens nucleation can occur when there is net freezing (V > 0).In figure 5, we show the local effective pressure N loc as a function of depth z for melting (left) and freezing (right).Here the top of the domain is the lowest ice lens z = 20 and N loc is evaluated in the fringe region above z f ≈ 18 (left) and z f ≈ 6 (right).In the steady melting case, N loc monotonically increases from the ice lens to the effective pressure at the bottom of the fringe N .The low pressure at the base of the ice lens draws in water as the ice lens infiltrates into the sediments through regelation (Gilpin, 1980;Fowler and Krantz, 1994;Rempel and Meyer, 2019).In the transient freezing case, water is drawn into the fringe and freezes onto the base of the lowest ice lens.The ice saturation increases, which lowers the permeability and requires a larger pressure difference to continue freezing.At some point in time, the local effective pressure reaches zero N loc (z n ) = 0 and a new ice lens forms at z n .We determine the time when a new ice lens forms using the 'events' functionality built into 'solve ivp', which flags the location and time of the new ice lens as an event and stops the integration.We then shift our domain up so that the new ice lens is at the top and pad the bottom of the domain with water-saturated sediments following the same incoming heat flux.Then, we restart the integration until the next ice lens forms.Written inside a loop, we generate sequences of ice lenses with an interlens time t .
Using the constant heave rate V and interlens time t , we can reconstruct the porosity structure a posteriori.We treat the porosity as a constant φ in the fringe and watersaturated sediment.In the ice lenses, the porosity is also constant with φ = 1.In one vertical dimension, mass conservation for sediments from equation (1) is for a constant sediment density ρ s .If we say that the sediment advection V s is given by the heave according to we find that ∂φ ∂t above z f since V is constant in this region.This model assures that the ice lenses and interstitial fringe advect vertically as one unit of rigid ice (O'Neill and Miller, 1985).Importantly, this treatment assumes that there is no volume expansion upon freezing and the water freezes in place, requiring that ice and water have the same density.This is true to leading order in the scaled and simplified model we describe in §2.7.2.We solve equation (92) equation analytically using the d'Alembert form, i.e. which we augment with the location of each new ice, thereby increase the region of applicability for equation ( 92) and use the fact that φ = 1 at z = z f .In figure 6, we show the evolution of the porosity with time as 5 new ice lenses sequentially nucleate and grow.With these parameters, the lenses form periodically with equal spacing and interlens times.
As we have seen, the two primary control parameters for the system are the heave rate V and the effective pressure N .So far we have shown steady states for balanced (V = 0) and melting (V < 0) conditions as well as ice lens nucleation for freezing conditions (V > 0).In figure 7, we show a regime diagram for the system behaviour as a function of V and N .Each point on the figure is a simulation and the colour denotes the grouping.For N < N c , no fringe forms and the lens either melts or grows, depending on the sign of V .When N > N c , ice infiltrates into the sediment forming a fringe.In melting or balanced cases where V ≤ 0, a steady state fringe thickness emerges (e.g.figures 3 and 4).For positive heave rates V > 0, steady states for relatively small effective pressures give way to periodic lenses as the heave rate increases.At the boundary between the steady and periodic regimes, there is a zone of hysteresis, where depending on the initial conditions the system either relaxes to a steady state or periodically generates new ice lenses.This same hysteresis was observed by Rempel et al. (2004).Although we observe some small variability in the interlens times, we do not see anything reminiscent of the chaotic regime found by Anderson and Worster (2014).Based on figure 7, if we have estimates for V and N in a geophysical context such as below a glacier, we can predict the system behaviour such as whether ice lenses will form.

Conclusions
In this paper, we derived the thermodynamical and fluid mechanical equations governing frost heave in a geophysical context.We systematically reduced the equations using the fact that ice and water have similar densities (i.e.small δ) as well as the large Stefan number St. We solved the reduced set of equations using an enthalpy method where the interstitial ice saturation acts like a liquidus condition in the frozen fringe and conservation < l a t e x i t s h a 1 _ b a s e 6 4 = " / g g l l 5 6 a h H 3 j u p R J s W N E d x r o / P u c e + y k 0 K w Y 0 N w 9 9 r 6 x v B g 4 e P N h / 3 n j x 9 9 n x r e + f F m V G l Z j B g S i h 9 k V A D g k s Y W G 4 F X B Q a a J 4 I O E + u P / n 9 8 x v Q h i t 5 a q c F j H I 6 l j z j j F o H f d 8 N d 6 + 2 + + F + 2 A y 8 H E S z o I 9 m 4 + R q Z 8 O S V L E y B 2 m Z o M Y M o 7 C w o 4 p q y 5 m A u k d K A w V l 1 3 Q M Q x d K m o M Z V c 1 Z a / z G I S n O l H a v t L h B 5 x U V z Y 2 Z 5 o l j 5 t R O T H f P g 6 v 2 h q X N P o w q L o v S g m S t U V Y K b B X 2 F 8 c p 1 8 C s m L q A M s 3 d W T G b U E 2 Z d e V Z c P G 5 t c m M u 4 m E W 6 b y n M q 0 I q 6 s o h 5 G o 4 p 8 A y I g s 6 T q R 0 T z 8 c S S u s M 9 P a 6 H s W N m L n 9 1 3 I 9 q N 8 V d k k 3 n S K k n p X i Z V c y z S O G r T A X 2 7 L + L F a J V q s u 4 q 8 P / z W H v l + N y S Q p 3 x U 1 b K l 8 n T A S V Y w F O i N t 6 u U + D d A s C 0 i h t W 2 U p U 9 D + t 5 4 P n f W S l 8 5 a Q Z 6 o u 8 q v W t c 9 7 N m t 3 9 s V I r a o Y v e R u c n U F S b t 0 0 i x w 3 A 9 g z p 0 U 9 w q n T Y + i x p v 8 S 9 R Y t q D J U q k v h G U a G 7 d c 3 0 a d b t y O T i L 9 6 N 3 + / H X u H 9 0 O O v Y T f Q K v U Z 7 K E L v 0 R H 6 g k 7 Q A D E 0 R j / Q T / Q r 2 A o O g s P g Y 0 t d X 5 t p X q K F E X z + A 2 c w h 0 M = < / l a t e x i t > 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " X w 6 p P M M 8 f w q y x b / l X R + B E x I x X 6 4 = " > A A A E j X i c j V N d b 9 M w F P W 2 A K N 8 7 I N H X i w 6 p P F S N Z E G P E z T J C b B 4 4 B 1 m 1 R 3 k + P c t N Y c O 7 K d b V W U n 8 A r / D b + D X Z S i T Y t a I 7 i X B + f c 4 9 9 l R v n g h v b 7 / 9 e W 9 8 I H j 1 + s v m 0 8 + z 5 i 5 d b 2 z u 7 5 0 Y V m s G A K a H 0 Z U w N C C 5 h Y L k V c J l r o F k s 4 C K + + e T 3 L 2 5 B G 6 7 k m Z 3 m M M r o W P K U M 2 o d 9 H 0 v 3 L v e 7 v Z 7 / X r g 5 S C c B V 0 0 G 6 f X O x u W J I o V G U j L B D V m G P Z z O y q p t p w J q D q k M J B T d k P H M H S h p B m Y U V m f t c J v H Z L g V G n 3 S o t r d F 5 R 0 s y Y a R Y 7 Z k b t x L T 3 P L h q b 1 j Y 9 O O o 5 D I v L E j W G K W F w F Z h f 3 G c c A 3 M i q k L K N P c n R W z C d W U W V e e B R e f W 5 v U u J t I u G M q y 6 h M S u L K K q p h O C r J N y A C U k v K b k g 0 H 0 8 s q V r c s 5 N q G D l m 6 v K X J 9 2 w c l P U J t l k j p R 4 U o K X W f k 8 i + S + y l R g z / 6 7 W C F a p b q K 2 j r 8 3 x z 2 Y T m u l q R w n 9 8 2 p f J 1 w k R Q O R b g h L i p l / v U S L s g I I 3 S t l E W M g H t f + v 5 0 F k v e e m 0 E W S x u i / 9 q n H d x 5 7 d + L 1 b I W K L K v Y Q m Z t M V W L S P L U U O w x X M 6 h F N / m d 0 k n t s 6 j x F v 8 S x a Y 5 W K x E 4 h t B i f r W H d e n Y b s r l 4 P z q B e + 7 0 V f o + 7 x 4 a x j N 9 W V e f J R e f W 5 v U u J t I u G M q y 6 h M S u L q K q p R O C 7 J d y A C U k v K X k g 0 n 0 w t q V r c s + N q F D l m 6 v K X x 7 2 w c l P U J t l k g Z R 4 U o J X W f k i i + S + y l R g z / 6 7 W C N a p 7 q K 2 j r 8 3 x z 2 Y T m u V q R w n 9 8 2 p f J 1 w k R Q O R H g h L i p l / v U S L s g I I 3 S t l E W M g H t / + v F 0 F m v e O m 0 E W S x u i / 9 q n H d w of enthalpy allows us to determine the fringe interface implicitly.For melting and balanced thermodynamics, we compared our enthalpy method to a steady state cast as an ordinary differential equation for temperature and found excellent agreement.In freezing cases, we found that the local effective pressure can go to zero within the fringe and nucleate a new ice lens.We accommodate this process in our enthalpy model using an 'events' function that stops the integration when a new ice lens forms and restarts the integration with a domain shifted below the new ice lens.Based on our solutions for the time between lenses, we can reconstruct the porosity profile showing the sequence of ice lenses.Finally, we compiled a regime diagram of our simulation results, showing the onset of periodic lensing and behaviour including hysteresis.These results will inform the regime of geophysical systems and future investigations into the role of compaction as well as comparisons with laboratory experiments.

A Steady heave near thermal balance
For an effective pressure N greater than the infiltration threshold N c , a steady fringe will form for melting (V < 0), balanced (V = 0), and weakly freezing (0 < V < V max (N )) thermodynamics, where V max depends on the effective pressure N (viz.figures 4 and 7).Alternatively, the maximum freezing heave rate can be thought of in the reverse: for a given heave rate (V > 0), there is a maximum effective pressure N max (V ), which is the largest load that can be supported by a steady fringe.For larger effective pressures, periodic lenses form.Here we will show how to calculate V max (N ) and N max (V ).
We start by rearranging (80) for N , which is To find the maximum effective pressure, we treat z f as fixed and set the derivative of N with respect to z equal to zero, i.e.
Treating the thermal conductivity variation as negligible, we insert the heat flux from equation (89) to find which can be solved using a root-finding algorithm for the lens temperature θ * (V ), i.e. the lowest lens temperature that can be supported for a given heave rate V .We now change the integration variable from z to θ in the force balance equation (94) using equation ( 89) and use the minimum lens temperature θ * in the limits of integration, i.e.
For a given value of N , we use another root-finding algorithm to find the heave rate V that satisfies equation ( 97), which is the maximum heave rate with a steady fringe.Repeated application of this root-finding algorithm results in the black curve V max (N ) shown on figure 7 and partitions the periodic lens regime from the hysteresis/steady lens regime.

Figure 3 :
Figure 3: Comparison between the finite volume enthalpy method and the steady state ode solution (nondimensional parameters V = −0.055and N = 1.5).(left) enthalpy H with height z below the lowest ice lens at z = 1.The frozen fringe extends from z f = 0.64 to z = 1, with water and sediments below.The nondimensional fringe thickness is h = z − z f = 0.36.The enthalpy in the water and sediments portion is close to zero throughout the depth.(right) temperature θ as a function of height z as calculated from the enthalpy.

Figure 5 :
Figure 5: Local effective pressure N loc (z) with height z scaled with the effective pressure N at the base of the fringe for: (left) a melting steady state (V = −0.01 and N = 2.9) where N loc > 0 throughout the fringe and (right) a transient freezing simulation (V = 0.20 and N = 2.9) where the local effective pressure goes to zero N loc (z n ) = 0, at point in the fringe z n , initiating a new ice lens.

Figure 6 :
Figure6: Porosity structure of a freezing sediment pack: (left) initial porosity profile with z = 25.The region above z is the lowest ice lens and the region below is fringe as well as water-saturated sediments.(right) a posteriori d'Alembert advection solutions for porosity after nucleation and growth of 4 periodic ice lenses.The nondimensional parameters V = 0.5 and N = 1.5 and the nondimensional interlens time is t = 9.6.
F r 9 A b t o x B 9 Q M f o C z p F A 8 T Q G P 1 A P 9 G v Y C s 4 C A 6 D o 4 a 6 v j b T v E I L I / j 8 B 2 t Y h 0 Q = < / l a t e x i t > 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " R o 5 j y k I A M / u i B 4 e + x S H Q 3 5 L + G d E = " > A A A E j n i c j V N d b 9 M w F P W 2 A K N 8 d f D I i 0 W H N B 6 o m j w M H j Y x i T 3 s c a B 1 m 1 R 3 k + P c t N Y c O 7 K d b V W U v 8 A r / D X + D X Z S i T Y t a I 7 i X B + f c 4 9 9 l R v n g h s 7 G P z e 2 N w K H j 1 + s v 2 0 8 + z 5 i 5 e v u j u v z 4 0 q N I M h U 0 L p y 5 g a E F z C 0 H I r 4 D L X Q L N Y w E V 8 8 9 X v X 9 y C N l z J M z v L Y Z z R i e Q p Z 9 R 6 6 O N u u H v d 7 Q 3 6 g 3 r g 1 S C c B z 0 0 H 6 f X O 1 u W J I o V G U j L B D V m F A 5 y O y 6 p t p w J q D q k M J B T d k M n M H K h p B m Y c V k f t s L v H Z L g V G n 3 S o t r d F F R 0 s y Y W R Y 7 Z k b t 1 L T 3 P L h u b 1 T Y 9 P O 4 5 D I v L E j W G K W F w F Z h f 3 O c c A 3 M i p k L K N P c n R W z K d W U

Figure 7 :
Figure 7: Regime diagram showing the system behaviour as a function of the heave rate V and effective pressure N .The maximum steady fringe thickness is shown as a solid black line and theory is described in the appendix.

Table 1 :
Table of parameters for frozen fringe.Now since the Stefan number St is large, the sensible heat contributions to the enthalpy H within the fringe can be ignored.Thus, we have that