CVPM 1.1: a flexible heat-transfer modeling system for permafrost

The Control Volume Permafrost Model (CVPM) is a modular heat-transfer modeling system designed for scientific and engineering studies in permafrost terrain, and as an educational tool. CVPM implements the nonlinear heat-transfer equations in 1-D, 2-D, and 3-D Cartesian coordinates, as well as in 1-D radial and 2-D cylindrical coordinates. To accommodate a diversity of geologic settings, a variety of materials can be specified within the model domain, including organic-rich materials, sedimentary rocks and soils, igneous and metamorphic rocks, ice bodies, borehole fluids, and other engineering materials. Porous materials are treated as a matrix of mineral and organic particles with pore spaces filled with liquid water, ice, and air. Liquid water concentrations at temperatures below 0 ^∘C due to interfacial, grain-boundary, and curvature effects are found using relationships from condensed matter physics; pressure and pore-water solute effects are included. A radiogenic heat-production term allows simulations to extend into deep permafrost and underlying bedrock. CVPM can be used over a broad range of depth, temperature, porosity, water saturation, and solute conditions on either the Earth or Mars. The model is suitable for applications at spatial scales ranging from centimeters to hundreds of kilometers and at timescales ranging from seconds to thousands of years. CVPM can act as a stand-alone model or the physics package of a geophysical inverse scheme, or serve as a component within a larger Earth modeling system that may include vegetation, surface water, snowpack, atmospheric, or other modules of varying complexity.

1 Introduction

Given the recent surge of interest in the cryosphere and its role in the Earth's climate system, a large number of permafrost models have been developed over the past few decades. An important characteristic of permafrost, especially in its fine-grained form, is that significant amounts of liquid water can coexist with ice within the pore spaces at temperatures well below 0 ^∘C due to a combination of interfacial, grain-boundary, curvature, solute, and pressure effects (Dash et al., 2006; Davis, 2001). Even at standard atmospheric pressure, liquid water has been observed at temperatures as low as −31 ^∘C in silty soils and −40 ^∘C in glass powders (Watanabe and Mizoguchi, 2002). The existence of liquid water at such low temperatures is of interest biologically (Rothschild and Mancinelli, 2001), particularly in the case of Mars, where permafrost could serve as a microbial refuge from high radiation levels if life ever existed there. From a geoscience perspective, it is also the presence of liquid water that makes permafrost so dynamic and interesting. Since the liquid water content of permafrost is highly temperature dependent and the thermal properties of solid and liquid water are so different (Anderson et al., 1973; Holten et al., 2012; Huber et al., 2012; Yen, 1981), the thermal response of permafrost to any temperature change is complicated by nonlinearities and feedbacks. As with the thermal properties, the mechanical properties of permafrost can be highly sensitive to temperature and the unfrozen water content, particularly within a few degrees of the freezing point T_F. As temperatures approach T_F, the material strength generally declines, increasing the likelihood of downslope creep, slope failures, accelerated lakeshore and coastal erosion, and ultimately thaw settlement (thermokarst) if temperatures become warm enough. If enough liquid water is available and the permafrost is sufficiently permeable, migration of liquid water towards colder temperatures can lead to significant frost heave, damaging buildings, roadways, and other facilities. With a warming climate, the dynamic response of permafrost is expected to be amplified, leading to accelerated landscape changes, disruption of vulnerable habitats and ecosystems, and damage to human infrastructure (ACIA, 2005; USARC, 2003).

A wide range of models has been developed to better understand the occurrence of permafrost and its dynamics in a warming world. These models range from simple analytical models to sophisticated numerical codes with integrated vegetation, snow, and atmospheric layers overlying permafrost (Riseborough et al., 2008; Zhang et al., 2003). The vast majority of these are 1-D vertical models. Although useful for simulating conditions beneath a uniform surface, 1-D models ignore important lateral heat-transfer effects occurring near large land-surface contrasts such as at the boundaries between tundra, rivers, lakes, oceans, glaciers, and human infrastructure. In addition, these models almost universally use empirical equations to predict the unfrozen water content at temperatures below 0 ^∘C. A significant limitation with this approach is that the coefficients appearing in the unfrozen water equations must be “calibrated” using field data for every material type, pressure, water saturation, and solute condition. Even with calibration, the empirical equations remain valid only over a limited range of temperatures. With an emphasis on simulating shallow permafrost and active-layer conditions, most permafrost models currently neglect the freezing-point depression due to pressure and the radiogenic heat-source term, both of which are needed to simulate conditions in deep permafrost.

In this paper, we present the new Control Volume Permafrost Model (CVPM v1.1) which is designed to relax several of the limitations imposed by previous models. CVPM implements the nonlinear heat-transfer equations in 1-D, 2-D, and 3-D Cartesian coordinates, as well as in 1-D radial and 2-D cylindrical coordinates. A variety of materials can be specified within the modeling domain, including organic-rich materials, sedimentary rocks and soils, igneous and metamorphic rocks, ice bodies, borehole fluids, and other engineering materials. Numerical implementation is based on the control-volume method (Anderson et al., 1984; Minkowycz et al., 1988; Patankar, 1980), allowing enthalpy fluxes to be exactly balanced at control-volume interfaces (e.g., at the interface between an ice lens and a siltstone). The unfrozen water content at temperatures below 0 ^∘C is found using relationships from condensed matter physics that utilize physical quantities (e.g., particle radii), rather than non-physical empirical coefficients requiring calibration. Pore pressure and solute effects are included in the unfrozen water equations. A radiogenic heat-production term is also included to allow simulations to extend into deep permafrost and underlying bedrock. CVPM is designed for use over a broad range of depth, temperature, rock and soil types, porosity, water saturation, and solute conditions. These conditions include the coldest temperatures experienced on Earth through the ice ages, as well as conditions on Mars, where the upper crust of the planet consists entirely of permafrost (Squyres et al., 1992). The model is suitable for applications at spatial scales ranging from centimeters to hundreds of kilometers and at timescales ranging from seconds to thousands of years. To achieve the greatest flexibility, CVPM does not include heat-transfer processes within a vegetation canopy, snowpack, or atmospheric boundary layer. Rather, CVPM focuses on permafrost and the underlying Earth materials. In this way, CVPM can act as a stand-alone model or the physics package of a geophysical inverse scheme, or serve as a component within a larger Earth modeling system that may include vegetation, surface water, snowpack, atmospheric, or other modules of varying complexity.

2 Governing equations

The basis for the CVPM model is the conservation of mass and enthalpy over time within any finite volume V. In integral form, the conservation equations take the form

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }\underset{V}{\int }{\displaystyle \frac{\partial \mathit{\rho }}{\partial t}}\mathrm{d}V=-\underset{A}{\int }\mathit{\rho }\mathit{v}\cdot \mathrm{d}\mathit{A}\left(\mathrm{mass}\right)\text{(2)}& {\displaystyle }& {\displaystyle }\underset{V}{\int }{\displaystyle \frac{\partial \left(\mathit{\rho }H\right)}{\partial t}}\mathrm{d}V=-\underset{A}{\int }\mathit{J}\cdot \mathrm{d}\mathit{A}+\underset{V}{\int }S\mathrm{d}V\left(\mathrm{enthalpy}\right),\end{array}$$

where ρ is the bulk density, ρv is the mass flux, A is the area-bounding volume V, H is the specific enthalpy, J is the enthalpy flux, S is the enthalpy production rate, and t is time. For the current version of CVPM, the velocity v is assumed to be sufficiently small so that the advective heat flux is negligible compared to the diffusive heat flux. In this case, the enthalpy flux is simply $\mathit{J}=-k\mathrm{\nabla }T$, where k is the bulk thermal conductivity and T is temperature. The medium within the model domain is assumed to consist of organic-rich materials, rocks and soils, ice bodies, and engineering materials. Porous materials are treated as a matrix (m) of mineral and organic particles with pore spaces filled with liquid water (ℓ), ice (i), and air (a). The porosity at any model location is then $\mathit{\varphi }={\mathit{\varphi }}_{\mathrm{\ell }}+{\mathit{\varphi }}_{\mathrm{i}}+{\mathit{\varphi }}_{\mathrm{a}}$, where ϕ_ℓ, ϕ_i, and ϕ_a are the volume fractions of the pore's constituents.

2.1 Heat capacity

Figure 1 (a) Variation of the specific heat with temperature for liquid water, ice, and two common minerals (quartz, albite). The vertical grey bar shows the specific heat range for most minerals at 300 K. Panel (b) shows the specific heat of liquid water and of ice in more detail. The vertical dashed line shows the temperature of the liquid–liquid critical point (T_c).

In permafrost, the enthalpy at temperature T consists of two components: one associated with the vibrational modes of the molecular lattice and the other due to the latent heat associated with the phase change of water:

$$\begin{array}{}\text{(3)}& \mathit{\rho }H\left(T\right)=\mathit{\rho }\underset{\mathrm{0}}{\overset{T}{\int }}{c}_p\left(T^{\prime }\right)\mathrm{d}{T}^{\prime }+{\mathit{\rho }}_{\mathrm{\ell }}\mathrm{\Delta }{H}_{\mathrm{fus}}{\mathit{\varphi }}_{\mathrm{\ell }}\left(T\right).\end{array}$$

Here, c_p is the specific heat of the bulk material, ρ_ℓ is the density of liquid water, and ΔH_fus is the specific enthalpy of fusion for water. Differentiating Eq. (3), the volumetric heat capacity at constant pressure, defined by $C\equiv \mathit{\rho }(\partial H/\partial T{)}_{\mathrm{P}}$, is given by the sum of the lattice-vibration and latent-heat terms:

$$\begin{array}{}\text{(4)}& C=\mathit{\rho }{c}_p+{\mathit{\rho }}_{\mathrm{\ell }}\mathrm{\Delta }{H}_{\mathrm{fus}}{\displaystyle \frac{\partial {\mathit{\varphi }}_{\mathrm{\ell }}}{\partial T}}.\end{array}$$

Since the density of air is much less than that of the other permafrost constituents, the lattice-vibration term is well approximated by the volume-weighted sum of the specific heats of the matrix materials, liquid water, and ice:

$$\begin{array}{}\text{(5)}& \mathit{\rho }{c}_p=(\mathrm{1}-\mathit{\varphi }){\mathit{\rho }}_{\mathrm{m}}{c}_{p\mathrm{m}}+{\mathit{\varphi }}_{\mathrm{\ell }}{\mathit{\rho }}_{\mathrm{\ell }}{c}_{p\mathrm{\ell }}+{\mathit{\varphi }}_{\mathrm{i}}{\mathit{\rho }}_{\mathrm{i}}{c}_{p\mathrm{i}},\end{array}$$

assuming ϕ_a<1.

Below 500 K, the specific heat of most matrix minerals is strongly temperature dependent, primarily due to the energies associated with the acoustic and optical modes of vibration (Kieffer, 1979, 1980). Due to the wide variety of crystalline mineral structures, simple analytic expressions for the temperature dependence of c_pm encompassing all minerals are unavailable. However, the relatively smooth temperature dependence predicted by detailed models indicates the specific heat for a given mineral can be adequately represented by a three-term Taylor expansion over the temperature range of interest for permafrost (Fig. 1). For most minerals, the specific heat falls within the range 630–870 $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ at 300 K and tends towards zero as T→0 K (Kieffer, 1979; Kittel, 1967; Robertson, 1988). Taylor expansions describing the temperature dependence ω(T) of c_pm for most common mineral groups are incorporated into CVPM. The matrix specific heat is then given by $c_{p\mathrm{m}}=c_{p\mathrm{m}}^{\circ }\mathit{\omega }\left(T\right)$, where $c_{p\mathrm{m}}^{\circ }$ is the specific heat of the dominant minerals at a standard temperature of 293.15 K.

Unlike most materials, experimental data for liquid water show an anomalous increase in specific heat (c_pℓ) with decreasing temperature. Holten et al. (2012) explained this and other peculiar behaviors of supercooled water with a thermodynamic model that assumes the existence of a liquid–liquid critical point at low temperatures. Based on their interpretation of available thermodynamic data, the liquid–liquid critical temperature T_c is near 227 K. A least-squares fit to a composite of data reported by Angell et al. (1982) below 273 K and the International Association for the Properties of Water and Steam (IAPWS) 2008 values above 273 K provide the following relationships:

$$\begin{array}{}\text{(6)}& c_{p\mathrm{\ell }}\left(T\right)=\left[\begin{array}{l}{a}_{\mathrm{1}}+a_{\mathrm{2}}{\left({\displaystyle \frac{T}{T_{\mathrm{c}}}}-\mathrm{1}\right)}^{-\mathrm{1}},\mathrm{235}\mathrm{K}<T\le \mathrm{265}\mathrm{K}\\ \sum _{i=\mathrm{1}}^{\mathrm{5}}{b}_{\mathrm{i}}\left({\displaystyle \frac{T}{\mathrm{310}\mathrm{K}}}-\mathrm{1}\right),\mathrm{265}\mathrm{K}<T\le \mathrm{360}\mathrm{K}.\end{array}\right.\end{array}$$

Values for the coefficients a₁, a₂, and b_i are listed in Table 1.

Table 1 Coefficients a_i and b_i in Eqs. (6)–(7) for the specific heat of liquid water and of ice (c_pℓ and c_pi in $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$).

For water ice (I_h), lattice vibrations lead to a simple linear relationship between the specific heat c_pi and temperature for T>150 K (Yen, 1981). Based on Yen's empirical relationship, c_pi is well represented by

$$\begin{array}{ll}{\displaystyle }{c}_{p\mathrm{i}}\left(T\right)=& {\displaystyle }{a}_{\mathrm{1}}+a_{\mathrm{2}}\left({\displaystyle \frac{T}{\mathrm{273.15}\mathrm{K}}}-\mathrm{1}\right),\\ \text{(7)}& {\displaystyle }& {\displaystyle }\mathrm{150}\mathrm{K}<T<\mathrm{273.15}\mathrm{K},\end{array}$$

with a₁=2096.1 and a₂=1943.8 ($\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$).

2.2 Unfrozen water

Studies dating back to the mid-1800s show that a melt layer can stably exist at the interface between ice and a foreign substrate (e.g., a mineral grain), even at temperatures well below the bulk freezing temperature of water T_f (Dash, 2002; Dash et al., 1995). A melt layer can similarly exist at the grain boundaries within polycrystalline ice. For both interfacial and grain-boundary melting, the liquid phase exists because it reduces the system's total free energy. Electrostatic interactions in molecular substances such as ice tend to be dominated by non-retarded van der Waals forces. In this case, the thickness of the liquid layer adjacent to a planar substrate is $L=\mathit{\lambda }\mathrm{\Delta }{T}^{-\mathrm{1}/\mathrm{3}}$, where $\mathrm{\Delta }T=(T_{\mathrm{f}}-T)$ is the temperature below the bulk freezing point (Dash et al., 2006; Wettlaufer and Worster, 1995). A ramification of this behavior is that the melting and freezing of ice in contact with mineral grains occur over a range of temperatures, rather than at a distinct temperature. Depending on the value for the interfacial melting parameter λ, substantial amounts of liquid water can exist at temperatures well below T_f. For a planar substrate, the interfacial melting parameter is given approximately by

$$\begin{array}{}\text{(8)}& \mathit{\lambda }={\left({\displaystyle \frac{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\Delta }\mathit{\gamma }{T}_{\mathrm{f}}}{{\mathit{\rho }}_{\mathrm{i}}\mathrm{\Delta }{H}_{\mathrm{fus}}}}\right)}^{\mathrm{1}/\mathrm{3}},\end{array}$$

where Δγ is the difference in the interfacial free energy with and without the melt layer, and σ is a constant on the order of a molecular diameter (Wettlaufer and Worster, 1995). Imperfections due to internal disorder (polycrystallinity, point defects, dislocations) within the ice and irregularities (pits, scratches, steps) in the substrate's surface can greatly increase the magnitude of Δγ and thereby the effective interfacial melting parameter λ. Due to the irregular nature of mineral surfaces and the likely disorder within interstitial ice, reliable expressions for parameter λ are currently lacking for most Earth materials. Thus, λ is best determined experimentally for frozen ground.

Surface curvature also affects the interfacial free energy and hence the thickness of liquid water films surrounding mineral grains. By considering the detailed effects of curvature along with interfacial and grain-boundary melting, Cahn et al. (1992) found that the volume fraction of liquid water in a porous medium consisting of spheres with radius r can be described by the sum of two temperature-dependent terms:

$$\begin{array}{}\text{(9)}& {\mathit{\varphi }}_{\mathrm{\ell }}=a_{\mathrm{1}}\left({\displaystyle \frac{\mathit{\lambda }}{r\mathrm{\Delta }{T}^{\mathrm{1}/\mathrm{3}}}}\right)+a_{\mathrm{2}}{\left({\displaystyle \frac{\mathit{\xi }}{r\mathrm{\Delta }T}}\right)}^{\mathrm{2}}.\end{array}$$

The first term is due to the combined effects of interfacial melting at the surface of the spherical particles and grain-boundary melting within polycrystalline pore ice. The second term is associated with high-curvature areas on the ice–liquid interface (e.g., near mineral grain-contact points and where ice-grain boundaries approach mineral surfaces). Coefficients a₁ and a₂ depend on how the particles are packed. For simple cubic packing, a₁=1.893 and a₂=3.367, while for cubic close packing, a₁=2.450 and a₂=8.572 (Cahn et al., 1992). Earth materials are likely to have intermediate values for the packing coefficients. In addition to its dependence on particle size and temperature, the second term depends on the interfacial free energy at the ice–water interface (γ_sℓ). The curvature coefficient at this interface, $\mathit{\xi }={\mathit{\gamma }}_{\mathrm{s}\mathrm{\ell }}{T}_{\mathrm{f}}/\left({\mathit{\rho }}_{\mathrm{i}}\mathrm{\Delta }{H}_{\mathrm{fus}}\right)$, numerically evaluates to 0.0259 µm K (Cahn et al., 1992). Given the temperature dependencies, the interfacial/grain-boundary term ($\mathrm{\Delta }{T}^{-\mathrm{1}/\mathrm{3}}$) dominates for strong undercooling (ΔT), while the second term (ΔT⁻²) dominates as temperatures approach the bulk freezing point (T_f). As the particle size r decreases, the transition between the two behaviors shifts to larger ΔT values (colder temperatures). While the first term is inversely proportional to r, the second term has an even stronger dependence on particle size (${\mathit{\varphi }}_{\mathrm{\ell }}\propto r^{-\mathrm{2}}$). Both terms result in greater amounts of unfrozen water in fine-grained materials (Fig. 2). Experimental data with graphitized carbon black and polystyrene powders confirm the form of Eq. (9) for monosized particles (Cahn et al., 1992).

Figure 2 Sensitivity of the volume fraction of liquid water ϕ_ℓ to particle diameter d=2r using Eq. (9) with λ=0.36 µm K^1∕3 and T_f=273.15 K, where $\mathrm{\Delta }T=T_{\mathrm{f}}-T$. The porosity, ϕ=0.54 in this example, sets the upper limit on ϕ_ℓ.

Although the particle and associated pore-size distributions in sandstones, limestones, and other rocks are often unimodal, those in mudrocks and soils typically are not (Kuila and Prasad, 2013). To accommodate a multimodal pore-size distribution, CVPM finds the total liquid water content by summing the contributions from the dominant modes. This is implemented in CVPM using a variant of Eq. (9):

$$\begin{array}{}\text{(10)}& {\mathit{\varphi }}_{\mathrm{\ell }}=\sum _j{\mathrm{\Psi }}_j\left[a_{\mathrm{1}}\left({\displaystyle \frac{\mathit{\lambda }}{r_j\mathrm{\Delta }{T}^{\mathrm{1}/\mathrm{3}}}}\right)+a_{\mathrm{2}}{\left({\displaystyle \frac{\mathit{\xi }}{r_j\mathrm{\Delta }T}}\right)}^{\mathrm{2}}\right],\end{array}$$

where Ψ_j is the relative volume fraction of pores associated with the mode whose mean particle size is r_j. Tests with sample pore-size distributions show the larger pore sizes contribute by far the most to the unfrozen water content. Thus, CVPM currently limits the number of modes contributing to ϕ_ℓ to two (j≤2). In this case,

$$\begin{array}{}\text{(11)}& {\mathrm{\Psi }}_{\mathrm{1}}={\displaystyle \frac{{\left(\frac{r_{\mathrm{1}}}{r_{\mathrm{2}}}\right)}^{\mathrm{3}}}{{\left(\frac{r_{\mathrm{1}}}{r_{\mathrm{2}}}\right)}^{\mathrm{3}}+\left(\frac{n_{\mathrm{2}}}{n_{\mathrm{1}}}\right)}},{\mathrm{\Psi }}_{\mathrm{2}}={\displaystyle \frac{\left(\frac{n_{\mathrm{2}}}{n_{\mathrm{1}}}\right)}{{\left(\frac{r_{\mathrm{1}}}{r_{\mathrm{2}}}\right)}^{\mathrm{3}}+\left(\frac{n_{\mathrm{2}}}{n_{\mathrm{1}}}\right)}},\end{array}$$

where (n₂∕n₁) is the ratio of the number density of pores with radius r₂ to those with radius r₁. As an example, suppose there are 100 times as many small pores (with a mean modal radius of 0.1 µm) as large pores (mean modal radius of 2 µm). Despite the greater number of small pores, the relative volume fraction of large pores is Ψ₁=0.988, while that of the smaller pores is only Ψ₂=0.012.

The undercooling ΔT used to evaluate the interfacial, grain-boundary, and curvature effects is measured relative to the bulk freezing temperature (T_f). For permafrost, pore pressures and dissolved solutes can significantly reduce T_f below the point at which pure water freezes ($T_{\mathrm{f}}^{*}$): 273.16 K at the triple point pressure P_tp=611.66 Pa (Guildner et al., 1976; Kittel, 1969; Nicholas and White, 2001). If θ_P and θ_s are the freezing-point depressions due to pressure and solutes, respectively, CVPM predicts the bulk freezing temperature using

$$\begin{array}{}\text{(12)}& T_{\mathrm{f}}=T_{\mathrm{f}}^{*}-{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{P}}.\end{array}$$

When solutes remain dilute, the freezing-point depression due to impurities can be approximated using simple relationships such as Blagden's law (Delapaz, 2015). However, due to the insolubility of most solutes in ice, impurities become strongly enriched in the liquid pore water as permafrost freezes. As a result, solute–solute interactions become increasingly important, leading to significant deviations from the ideal behavior exhibited by dilute solutions. To account for the non-ideal behavior of aqueous electrolyte solutions at higher solute concentrations, CVPM uses the relationship

$$\begin{array}{}\text{(13)}& -\ln \left(a_{\mathrm{w}}\right)={\displaystyle \frac{\mathrm{\Delta }{H}_{\mathrm{fus}}^{\circ }}{R{T}^{\circ }}}\left({\displaystyle \frac{{\mathit{\theta }}_{\mathrm{s}}}{T^{\circ }}}\right)+\left({\displaystyle \frac{\mathrm{\Delta }{H}_{\mathrm{fus}}^{\circ }}{R{T}^{\circ }}}-{\displaystyle \frac{\mathrm{\Delta }{c}_p^{\circ }}{\mathrm{2}R}}\right){\left({\displaystyle \frac{{\mathit{\theta }}_{\mathrm{s}}}{T^{\circ }}}\right)}^{\mathrm{2}}\end{array}$$

between the water activity a_w and the solute freezing-point depression θ_s (Robinson and Stokes, 1959). Here, $\mathrm{\Delta }{H}_{\mathrm{fus}}^{\circ }$ is the molar enthalpy of fusion at a standard temperature T^∘, $\mathrm{\Delta }{c}_p^{\circ }$ is the difference between the molar heat capacities of liquid water and ice, and R is the gas constant. Using established values for $\mathrm{\Delta }{H}_{\mathrm{fus}}^{\circ }$, $\mathrm{\Delta }{c}_p^{\circ }$, R, and T^∘, Eq. (13) simplifies to

$$\begin{array}{}\text{(14)}& -\ln a_{\mathrm{w}}=(\mathrm{9.687}\times {\mathrm{10}}^{-\mathrm{3}}){\mathit{\theta }}_{\mathrm{s}}+(\mathrm{4.76}\times {\mathrm{10}}^{-\mathrm{6}}){\mathit{\theta }}_{\mathrm{s}}^{\mathrm{2}},\end{array}$$

where θ_s is in Kelvin. The extent to which aqueous electrolyte solutions deviate from ideal behavior varies greatly, depending on the composition of the solute. As a result, the water activity a_w depends on the particular solute and its mole fraction x_s. Several expressions have been proposed for the water activity of non-ideal electrolyte solutions. CVPM uses the following equation proposed by Miyawaki et al. (1997):

$$\begin{array}{}\text{(15)}& a_{\mathrm{w}}=(\mathrm{1}-x_{\mathrm{s}})\exp \left(\mathit{\alpha }{x}_{\mathrm{s}}^{\mathrm{2}}+\mathit{\beta }{x}_{\mathrm{s}}^{\mathrm{3}}\right),\end{array}$$

where the coefficients α and β are solute dependent. For example, α=1.825, 4.754, and 11.859 for NaCl, KCl, and MgCl₂, respectively, while $\mathit{\beta }=-\mathrm{20.78}$, −49.37, and −404.5 (Miyawaki et al., 1997). Since essentially all the solutes are concentrated in the aqueous solution upon freezing, the solute mole fraction at any stage during freezing or thawing is given by $x_{\mathrm{s}}=x_{\mathrm{s}}^{\star }\left[\right({\mathit{\varphi }}_{\mathrm{\ell }}+{\mathit{\varphi }}_{\mathrm{i}})/{\mathit{\varphi }}_{\mathrm{\ell }}]$, where $x_{\mathrm{s}}^{\star }$ is the solute concentration in the fully melted system (ϕ_i=0). Once x_s and a_w have been established, the freezing-point depression θ_s can be found by solving Eq. (14).

Figure 3 shows the sensitivity of the unfrozen water content ϕ_ℓ to solutes predicted by CVPM for a medium-grained silt along with the temperature sensitivity $\partial {\mathit{\varphi }}_{\mathrm{\ell }}/\partial T$ needed to find the latent-heat component of the volumetric heat capacity C. The results show that even small amounts of solute can significantly affect ϕ_ℓ and $\partial {\mathit{\varphi }}_{\mathrm{\ell }}/\partial T$. As noted by Dash et al. (2006), the great sensitivity of ϕ_ℓ to impurities is a likely cause for the considerable disagreement between the results of various unfrozen water experiments. An additional sensitivity can occur at cold temperatures if solute concentrations are sufficiently high. This occurs when the solution reaches its saturation limit, beyond which the solute begins to precipitate upon further cooling. This leads to the spike in $\partial {\mathit{\varphi }}_{\mathrm{\ell }}/\partial T$ values near −21 ^∘C seen in Fig. 3b when aqueous NaCl concentrations $x_{\mathrm{s}}^{\star }$ exceed 0.005. The temperature at which the solute-saturation spike occurs varies, depending on the particular solute. Outside of the saturation limit, the largest $\partial {\mathit{\varphi }}_{\mathrm{\ell }}/\partial T$ values occur as the last bits of pore ice melt upon warming (ϕ_i→0). The volumetric heat capacity mirrors the temperature sensitivity $\partial {\mathit{\varphi }}_{\mathrm{\ell }}/\partial T$ but with minimum values established by the lattice-vibration term in Eq. (4).

Figure 3 Volume fraction of liquid water ϕ_ℓ predicted for a medium-grained silt (r=15 µm) with NaCl solute concentrations $x_{\mathrm{s}}^{\star }$ ranging 0 to 0.03 (a). The depth z≈0 m, ϕ=0.54, and λ=0.36 µm K^1∕3. Green dots show measured ϕ_ℓ values for Fairbanks silt (ϕ=0.54) reported by Yuanlin and Carbee (1987). Stars correspond to the empirical equation fitting unfrozen water measurements in Fairbanks silt given by Anderson et al. (1973). Panel (b) shows the sensitivity of the liquid water content to temperature, $\partial {\mathit{\varphi }}_{\mathrm{\ell }}/\partial T$.

As previously mentioned, the interfacial melting parameter is best determined experimentally for natural Earth materials. Inversion of unfrozen water data (shown in Fig. 3a) from Yuanlin and Carbee (1987) using CVPM yields a value of λ=0.36 µm K^1∕3 for Fairbanks silt. Other parameters determined by the inversion are Ψ₁≃1 and $x_{\mathrm{s}}^{\star }=\mathrm{0.0008}$. Thus, the pore-size distribution for this material is approximately unimodal. In addition, trace amounts of impurities appear to have been present during the unfrozen water experiments despite efforts to eliminate them. The λ value determined for Fairbanks silt is about 100 times that determined for liquid films adjacent to smooth metal wires (Cahn et al., 1992), testifying to the importance of mineral surface irregularities and imperfections on the interfacial free energy in natural Earth materials. While more work needs to be done to quantify λ for the range of materials expected in permafrost, preliminary inversions for sedimentary materials (e.g., Suffield silty clay and kaolinite) yield values within ±10 % of that found for Fairbanks silt. At this point, the interfacial melting parameter λ does not appear to vary substantially amongst natural Earth materials.

Given that permafrost occurs at depths in excess of 1 km in some high-latitude areas and at 3–4 km beneath the polar ice sheets (Davis, 2001; MacGregor et al., 2016), the effect of pressure can be substantial on the bulk freezing temperature T_f and thereby the unfrozen water content ϕ_ℓ. If the interstitial pores are freely connected to the planet's surface, the pore pressure is equal to the hydrostatic pressure (P=ρ_ℓgz), where g is the acceleration of gravity and z is the depth below the surface. However, if the pore water is trapped, the pore pressure can be nearly equal to or exceed the lithostatic pressure (P=ρgz) (Turcotte and Schubert, 1982). In either case, the freezing-point depression due to pore pressure is

$$\begin{array}{}\text{(16)}& {\mathit{\theta }}_{\mathrm{P}}=a(P-P_{\mathrm{tp}}).\end{array}$$

As water in permafrost is likely to be saturated with air, the appropriate value for coefficient a is $\mathrm{9.8}\times {\mathrm{10}}^{-\mathrm{8}}$ K Pa⁻¹ (Cuffey and Paterson, 2010). Since both pressure situations are known to occur in sedimentary basins, both are implemented in CVPM. The lithostatic effect is generally 2–3 times that of the hydrostatic effect. Not only does the pressure effect increase the unfrozen water content with depth, it also increases the temperature sensitivity ($\partial {\mathit{\varphi }}_{\mathrm{\ell }}/\partial T$) and therefore the volumetric heat capacity C (Fig. 4).

Figure 4 Sensitivity of the volume fraction of liquid water ϕ_ℓ to depth below surface z. Solid lines are for hydrostatic pore pressures, while dashed lines are for lithostatic pressures. In this example, the solute (NaCl) concentration is $x_{\mathrm{s}}^{\star }=\mathrm{0.001}$, r=15 µm, ϕ=0.3, and λ=0.36 µm K^1∕3.

2.3 Thermal conductivity

Several mixing models are available for estimating the bulk thermal conductivity of multi-component systems. Of these, CVPM uses the Brailsford and Major (1964) two-phase random mixture (BM2) and three-phase (BM3) models, which are recommended as being the best for use with in situ Earth materials (Roy et al., 1981). Assuming a random mixture of pores and matrix material, the bulk thermal conductivity of permafrost can be described by the BM2 model:

$$\begin{array}{ll}{\displaystyle }k=& {\displaystyle }{k}_{\mathrm{m}}[{\displaystyle \frac{(\mathrm{2}\mathit{\chi }-\mathrm{1})-\mathrm{3}\mathit{\varphi }(\mathit{\chi }-\mathrm{1})}{\mathrm{4}\mathit{\chi }}}\\ \text{(17)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{{\left\{{\left[(\mathrm{2}\mathit{\chi }-\mathrm{1})-\mathrm{3}\mathit{\varphi }(\mathit{\chi }-\mathrm{1})\right]}^{\mathrm{2}}+\mathrm{8}\mathit{\chi }\right\}}^{\frac{\mathrm{1}}{\mathrm{2}}}}{\mathrm{4}\mathit{\chi }}}],\end{array}$$

where k_m is the conductivity of the matrix material, k_p is the conductivity of the pores, and χ is their ratio (k_m∕k_p).

For matrix minerals, the thermal conductivity depends primarily on the temperature and mineral composition. Using thermal conductivity data obtained by Birch and Clark (1940a, b) over the temperature range 273–473 K, Sass et al. (1992) found the temperature dependence could be separated from the compositional dependence using a function of the form

$$\begin{array}{ll}{\displaystyle }{k}_{\mathrm{m}}\left(T\right)=& {\displaystyle }{k}_{\mathrm{m}}^{\circ }{\left[a_{\mathrm{1}}+\left(T-\mathrm{273.15}\mathrm{K}\right)\left(a_{\mathrm{2}}-{\displaystyle \frac{a_{\mathrm{3}}}{k_{\mathrm{m}}^{\circ }}}\right)\right]}^{-\mathrm{1}},\\ \text{(18)}& {\displaystyle }& {\displaystyle }\mathrm{150}\mathrm{K}<T<\mathrm{570}\mathrm{K},\end{array}$$

where $k_{\mathrm{m}}^{\circ }$ is the value of the matrix conductivity at a standard temperature of 273.15 K. As the a_i coefficients are fairly insensitive to rock type, the effects of mineralogy and texture are almost entirely encapsulated in $k_{\mathrm{m}}^{\circ }$. A more recent analysis indicates the a_i coefficients (Table 2) are slightly different for the mineral assemblages that dominate sedimentary rocks from those that occur in magmatic and metamorphic rocks (Vosteen and Schellschmidt, 2003). The upper temperature limit for Eq. (18) is set below the temperatures at which metamorphosis occurs in sedimentary rocks and well below the point where radiative heat transfer within crystal lattices becomes important (Clauser and Huenges, 1995). As very little thermal conductivity data exist for rocks and minerals below 273 K, the validity of Eq. (18) has yet to be tested at lower temperatures. The little data that do exist suggest that the transition from the intermediate-temperature behavior (Eq. 18) to the low-temperature behavior (k_m∝T³) (Parrott and Stuckes, 1975) generally occurs below 100 K. For example, in garnets, the transition occurs at 20–30 K (Slack and Oliver, 1971). We tentatively set the lower limit of validity for Eq. (18) at 150 K.

Table 2 Coefficients a_i and b_i in Eqs. (18), (21)–(22), (24), and (26)–(27) for the thermal conductivity of matrix minerals, liquid water, ice, air, and CO₂ gas (k_x in $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$). For matrix minerals, $a_{\mathrm{i}}^{\mathrm{mm}}$ and $a_{\mathrm{i}}^{\mathrm{sed}}$ refer to the coefficients appropriate for magmatic/metamorphic and sedimentary mineral assemblages, respectively.

To find the thermal conductivity of the pores (k_p), CVPM utilizes the three-phase BM3 model:

$$\begin{array}{}\text{(19)}& k_{\mathrm{c}x}=k_{\mathrm{1}}\left[{\displaystyle \frac{{\mathit{\psi }}_{\mathrm{1}}+\mathrm{3}{\displaystyle }\left({\displaystyle \frac{{\mathit{\psi }}_{\mathrm{2}}}{\mathrm{2}{\mathit{\chi }}_{\mathrm{2}}+\mathrm{1}}}+{\displaystyle \frac{{\mathit{\psi }}_{\mathrm{3}}}{\mathrm{2}{\mathit{\chi }}_{\mathrm{3}}+\mathrm{1}}}\right)}{{\mathit{\psi }}_{\mathrm{1}}+\mathrm{3}{\displaystyle }\left({\displaystyle \frac{{\mathit{\psi }}_{\mathrm{2}}{\mathit{\chi }}_{\mathrm{2}}}{\mathrm{2}{\mathit{\chi }}_{\mathrm{2}}+\mathrm{1}}}+{\displaystyle \frac{{\mathit{\psi }}_{\mathrm{3}}{\mathit{\chi }}_{\mathrm{3}}}{\mathrm{2}{\mathit{\chi }}_{\mathrm{3}}+\mathrm{1}}}\right)}}\right],\end{array}$$

where the three phases (x) are liquid water, ice, and air. Here, ${\mathit{\chi }}_{\mathrm{2}}=(k_{\mathrm{1}}/k_{\mathrm{2}})$, ${\mathit{\chi }}_{\mathrm{3}}=(k_{\mathrm{1}}/k_{\mathrm{3}})$, and the ψ_x are the relative volume fractions of the pore's constituents (${\mathit{\psi }}_x={\mathit{\varphi }}_x/\mathit{\varphi }$). Similar to other three-phase models, BM3 assumes phases 2 and 3 are randomly distributed within a continuous phase 1. If the relative volume fraction of any of the three constituents exceeds a threshold (ψ_x≥α), it is assumed that component is the continuous phase and k_p is calculated directly from Eq. (19). A comparison of the results from the BM3 model in the limit ψ₂→0 or ψ₃→0 with those of the BM2 model suggests a reasonable choice for α is ∼0.75. If none of the relative volume fractions exceed α, Eq. (19) is used to calculate the conductivity of the pore space assuming each of the three components, in turn, is the continuous phase to produce values k_cℓ (continuous liquid water phase), k_ci (continuous ice phase), and k_ca (continuous air phase). The pore conductivity is then found from a simple weighted average:

$$\begin{array}{}\text{(20)}& k_{\mathrm{p}}=w_{\mathrm{\ell }}{k}_{\mathrm{c}\mathrm{\ell }}+w_i{k}_{\mathrm{ci}}+w_{\mathrm{a}}{k}_{\mathrm{ca}},\end{array}$$

where the weights w_x are based on the relative volume fractions ψ_x and the requirement that k_p be continuous across the lines ψ_ℓ=α, ψ_i=α, and ψ_a=α in three-phase space (Fig. 5).

Figure 5 Variation of the pore thermal conductivity k_p on Earth at −10 ^∘C with the relative volume fractions of liquid water (ψ_ℓ), ice (ψ_i), and air (ψ_a) within the pores. Threshold α is 0.75 in this example.

For the thermal conductivity of liquid water k_ℓ, CVPM uses the simplified correlating equation recommended by Huber et al. (2012) for use at 0.1 MPa:

$$\begin{array}{}\text{(21)}& k_{\mathrm{\ell }}\left(T\right)=\sum _{i=\mathrm{1}}^{\mathrm{4}}{a}_{\mathrm{i}}{\left({\displaystyle \frac{T}{\mathrm{300}\mathrm{K}}}\right)}^{b_{\mathrm{i}}},\mathrm{250}\mathrm{K}<T\le \mathrm{383}\mathrm{K},\end{array}$$

with coefficients a_i and b_i (Table 2). Although the formal lower limit for Eq. (21) is 273.15 K, Huber et al. (2012) find that it extrapolates in a physically reasonable manner down to ∼250 K (Fig. 6), producing results very close to those of the new detailed IAPWS formulation for k_ℓ. Thermal conductivity data for supercooled water do not appear to exist below 250 K at this time, preventing the development of accurate correlating equations at lower temperatures. This is a minor limitation for the thermal model as the relative amount of liquid water is expected to be small at colder temperatures.

Figure 6 Variation of the thermal conductivity with temperature for liquid water, ice, air (terrestrial atmosphere), CO₂ gas (Martian atmosphere), and magmatic/metamorphic matrix minerals with $k_{\mathrm{m}}^{\circ }=\mathrm{3.5}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ (feldspar poor) and $k_{\mathrm{m}}^{\circ }=\mathrm{2.0}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ (feldspar rich).

Experimental data for the thermal conductivity of ice k_i exist at temperatures as cold as 60 K. Based on these data, Yen (1981) recommends the function

$$\begin{array}{}\text{(22)}& k_{\mathrm{i}}\left(T\right)=a_{\mathrm{1}}\exp (-a_{\mathrm{2}}T),\mathrm{60}\mathrm{K}<T\le \mathrm{273.15}\mathrm{K}\end{array}$$

for describing the temperature dependence of k_i, with a₁=9.828 $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ and a₂=0.0057 K⁻¹.

For the terrestrial environment, the thermal conductivity of air k_a can be separated into the sum of two terms: a “dilute gas” term k_o that depends solely on temperature and a “residual” term Δk that depends on air density:

$$\begin{array}{}\text{(23)}& k_{\mathrm{a}}({\mathit{\rho }}_{\mathrm{a}},T)=k_{\mathrm{o}}\left(T\right)+\mathrm{\Delta }k\left({\mathit{\rho }}_{\mathrm{a}}\right).\end{array}$$

For the dilute gas term, Stephan and Laesecke (1985) recommend the correlating equation:

$$\begin{array}{ll}{\displaystyle }{k}_{\mathrm{o}}\left(T\right)=& {\displaystyle }\sum _{i=\mathrm{1}}^{\mathrm{9}}{a}_{\mathrm{i}}{\left({\displaystyle \frac{T}{\mathrm{132.52}\mathrm{K}}}\right)}^{(i-\mathrm{4})/\mathrm{3}},\\ \text{(24)}& {\displaystyle }& {\displaystyle }\mathrm{70}\mathrm{K}<T<{\mathrm{10}}^{\mathrm{3}}\mathrm{K},\end{array}$$

with a_i coefficients (Table 2). At typical terrestrial surface pressures (∼0.1 MPa), the residual term Δk is $\mathrm{5.17}\times {\mathrm{10}}^{-\mathrm{5}}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ (Stephan and Laesecke, 1985).

When considering permafrost on Mars, the thermal properties of a different atmospheric gas must be used. The Martian atmosphere is currently 95 % carbon dioxide, a gas that has a thermodynamic critical point at 304.107 K, 7.3721 MPa. At gas densities below 25 kg m⁻³, the effects of the critical region are small enough that the thermal conductivity can again be described by Eq. (23). In the case of CO₂, the dilute gas contribution to k_a is

$$\begin{array}{ll}{\displaystyle }{k}_{\mathrm{o}}\left(T\right)=& {\displaystyle }{\displaystyle \frac{\mathrm{4.75598}\times {\mathrm{10}}^{-\mathrm{4}}\left(\mathrm{1}+{\displaystyle \frac{\mathrm{2}{c}_{\mathrm{int}}}{\mathrm{5}{k}_{\mathrm{B}}}}\right)T^{\mathrm{1}/\mathrm{2}}}{{\mathcal{C}}_{\mathrm{R}}\left(T\right)}},\\ \text{(25)}& {\displaystyle }& {\displaystyle }\mathrm{200}\mathrm{K}<T<{\mathrm{10}}^{\mathrm{3}}\mathrm{K},\end{array}$$

where c_int is the ideal gas heat capacity, k_B is the Boltzmann constant, and 𝒞_R is the reduced effective cross-section (Vesovic et al., 1990). The correlating equations for 𝒞_R and c_int provided by Vesovic et al. (1990) are

$$\begin{array}{}\text{(26)}& {\displaystyle }& {\displaystyle }{\mathcal{C}}_{\mathrm{R}}\left(T\right)=\sum _{i=\mathrm{0}}^{\mathrm{7}}{a}_{\mathrm{i}}{\left({\displaystyle \frac{T}{\mathrm{251.196}\mathrm{K}}}\right)}^{-i}\text{(27)}& {\displaystyle }& {\displaystyle }{c}_{\mathrm{int}}=k_{\mathrm{B}}\left[\mathrm{1}+\exp \left({\displaystyle \frac{-\mathrm{183.5}\mathrm{K}}{T}}\right)\sum _{i=\mathrm{1}}^{\mathrm{5}}{b}_{\mathrm{i}}{\left({\displaystyle \frac{T}{\mathrm{100}\mathrm{K}}}\right)}^{\mathrm{2}-i}\right]\end{array}$$

(coefficients a_i and b_i given in Table 2). At current Martian surface pressures (0.6 kPa), the residual component Δk is $\mathrm{3.53}\times {\mathrm{10}}^{-\mathrm{7}}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$. Even when the Martian atmosphere was denser, the contribution of Δk to the overall thermal conductivity of the CO₂ gas would have been quite small.

2.4 Compaction and heat-production functions

In sedimentary basins, overburden pressure causes the porosity ϕ to decrease with depth due to pressure solution and mechanical-compaction processes (Revil et al., 2002). The former process changes the mineral shapes in response to grain-contact stresses, while the latter results in the slippage and rotation of the grains. With increasing overburden pressure, the porosity ultimately reaches a residual (or critical) porosity ϕ_c that depends on the grain-shape and grain-size distribution. Shales and mudstones are much more easily compacted than sandstones due to the plate-like shape of the mineral grains. Although fairly sophisticated compaction models now exist, CVPM uses the simple frequently used compaction function attributed to Athy (1930),

$$\begin{array}{}\text{(28)}& \mathit{\varphi }\left(z\right)={\mathit{\varphi }}_{\mathrm{0}}\exp \left(-z/h_{\mathrm{c}}\right),\mathit{\varphi }\ge {\mathit{\varphi }}_{\mathrm{c}},\end{array}$$

to account for overburden pressures. This function has been successfully used in a large number of studies (Burns et al., 2005; Fjeldskaar et al., 2004). Here, ϕ₀ is the porosity extrapolated to the surface, while h_c is the compaction length scale. Parameters ϕ₀ and h_c depend on both the lithology and the effective stress history.

CVPM assumes the enthalpy-production rate S is associated with the decay of radionuclides. In this case,

$$\begin{array}{}\text{(29)}& S\left(z\right)=S_{\mathrm{0}}\exp \left(-z/h_{\mathrm{s}}\right),\end{array}$$

where S₀ is the radioactive heat-production rate extrapolated to the surface and h_s is the heat-production length scale (Turcotte and Schubert, 1982). Surface heat-production rates S₀ can vary from 0.002 to 5.5 µW m⁻³, depending on lithology (Rybach, 1988), while h_s is typically on the order of 10 km.

3 Numerical implementation

Figure 7 Schematic showing the nomenclature associated with a control volume centered on grid point P for 2-D Cartesian (XZ) coordinates. The control volume is bounded by interfaces located at x_w, x_e, z_u, and z_d, through which enthalpy fluxes J_w, J_e, J_u, and J_d pass. Grid points W, E, U, and D are located at the center of the neighboring control volumes (CVs). The nomenclature for 2-D cylindrical coordinates is completely analogous with R replacing X. Three-dimensional (3-D) Cartesian coordinates introduce an additional axis (Y) with neighboring grid points S and N.

The CVPM modeling system implements the governing equations in 1-D, 2-D, and 3-D Cartesian coordinates (X, XZ, XYZ), as well as in 1-D radial (R) and 2-D cylindrical (RZ) coordinates. Discretization follows the control-volume approach (Anderson et al., 1984; Minkowycz et al., 1988; Patankar, 1980) in which the problem domain is divided into a set of contiguous control volumes (CVs). Scalars such as temperature T and thermal conductivity k are computed at grid points located in the center of the CVs, while the enthalpy fluxes J are computed at control-volume interfaces (Fig. 7). Development of the CVPM permafrost model begins by integrating the two conservation equations (Eqs. 1–2) over a time step Δt. With velocity v≃0, the conservation equations become

$$\begin{array}{ll}\text{(30)}& {\displaystyle }& {\displaystyle }\underset{V}{\int }\left({\mathit{\rho }}^{n+\mathrm{1}}-{\mathit{\rho }}^n\right)\mathrm{d}V=\mathrm{0}{\displaystyle }& {\displaystyle }\underset{V}{\int }\left[(\mathit{\rho }H{)}^{n+\mathrm{1}}-(\mathit{\rho }H{)}^n\right]\mathrm{d}V\\ \text{(31)}& {\displaystyle }& {\displaystyle }=-\underset{t^n}{\overset{t^{n+\mathrm{1}}}{\int }}\underset{A}{\int }\mathit{J}\cdot \mathrm{d}\mathit{A}\mathrm{d}t+\underset{t^n}{\overset{t^{n+\mathrm{1}}}{\int }}\underset{V}{\int }S\mathrm{d}V\mathrm{d}t,\end{array}$$

where superscript n refers to the time step following standard numerical nomenclature (e.g., $t^{n+\mathrm{1}}=t^n+\mathrm{\Delta }t$). Equation (30) states the bulk density integrated over any control volume is time invariant. The conservation of enthalpy equation (Eq. 31) can be put in a more convenient form by noting

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }(\mathit{\rho }H{)}^{n+\mathrm{1}}-(\mathit{\rho }H{)}^n=\left[\mathit{\rho }{c}_p+{\mathit{\rho }}_{\mathrm{\ell }}\mathrm{\Delta }{H}_{\mathrm{fus}}{\left.{\displaystyle \frac{\partial {\mathit{\varphi }}_{\mathrm{\ell }}}{\partial T}}\right|}^{T^n}\right]\left(T^{n+\mathrm{1}}-T^n\right)\\ \text{(32)}& {\displaystyle }& {\displaystyle }=C^n\left(T^{n+\mathrm{1}}-T^n\right).\end{array}$$

This follows from Eq. (3) and a Taylor series expansion. To provide the flexibility of running the model in either explicit, implicit, or fully implicit modes, the net heat flux into a control volume is approximated by a linear combination of values at either end of a time step:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }-\underset{t^n}{\overset{t^{n+\mathrm{1}}}{\int }}\underset{A}{\int }\mathit{J}\cdot \mathrm{d}\mathit{A}\mathrm{d}t=-\mathrm{\Delta }t\left[f{\left(\underset{A}{\int }\mathit{J}\cdot \mathrm{d}\mathit{A}\right)}^{n+\mathrm{1}}\right.\\ \text{(33)}& {\displaystyle }& {\displaystyle }\left.+(\mathrm{1}-f){\left(\underset{A}{\int }\mathit{J}\cdot \mathrm{d}\mathit{A}\right)}^n\right].\end{array}$$

The explicit/implicit weighting factor f can take any value between 0 and 1. Following Patankar (1980), the heat fluxes across the z_u and z_d interfaces in the vertical direction (see Fig. 7) are

$$\begin{array}{}\text{(34)}& J_{\mathrm{u}}=-{\stackrel{\mathrm{̃}}{k}}_{\mathrm{u}}\left({\displaystyle \frac{T_{\mathrm{P}}-T_{\mathrm{U}}}{z_{\mathrm{P}}-z_{\mathrm{U}}}}\right),J_{\mathrm{d}}=-{\stackrel{\mathrm{̃}}{k}}_{\mathrm{d}}\left({\displaystyle \frac{T_{\mathrm{D}}-T_{\mathrm{P}}}{z_{\mathrm{D}}-z_{\mathrm{P}}}}\right),\end{array}$$

where ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{u}}$ and ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{d}}$ are the “effective” conductivities at the upper and lower interfaces, defined by

$$\begin{array}{}\text{(35)}& {\stackrel{\mathrm{̃}}{k}}_{\mathrm{u}}=\left[{\displaystyle \frac{\mathrm{1}}{{\displaystyle \frac{(\mathrm{1}-{\mathit{\epsilon }}_{\mathrm{u}})}{k_{\mathrm{U}}}}+{\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{u}}}{k_{\mathrm{P}}}}}}\right],{\stackrel{\mathrm{̃}}{k}}_{\mathrm{d}}=\left[{\displaystyle \frac{\mathrm{1}}{{\displaystyle \frac{(\mathrm{1}-{\mathit{\epsilon }}_{\mathrm{d}})}{k_{\mathrm{P}}}}+{\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{d}}}{k_{\mathrm{D}}}}}}\right],\end{array}$$

with fractional distances

$$\begin{array}{}\text{(36)}& {\mathit{\epsilon }}_{\mathrm{u}}=\left({\displaystyle \frac{z_{\mathrm{P}}-z_{\mathrm{u}}}{z_{\mathrm{P}}-z_{\mathrm{U}}}}\right),{\mathit{\epsilon }}_{\mathrm{d}}=\left({\displaystyle \frac{z_{\mathrm{D}}-z_{\mathrm{d}}}{z_{\mathrm{D}}-z_{\mathrm{P}}}}\right).\end{array}$$

Subscripts used here indicate grid point and interface locations. For example, T_P is the temperature at grid point P, while J_u is the heat flux across the interface located at depth z_u. Fluxes across the other interfaces are defined in a completely analogous way. The use of effective conductivities guarantees that the heat fluxes exactly balance at an interface between materials with very different thermal properties (e.g., between a siltstone and an ice lens). The source-term integral in Eq. (31) is left in a very general form:

$$\begin{array}{}\text{(37)}& {\mathcal{S}}_{\mathrm{P}}=\mathrm{\Delta }t\underset{V}{\int }S\mathrm{d}V.\end{array}$$

Substituting Eqs. (32)–(34) and (37) into Eq. (31), the discrete form of the enthalpy balance for a control volume centered on grid point P can be written as

$$\begin{array}{}\text{(38)}& a_{\mathrm{P}}{T}_{\mathrm{P}}^{n+\mathrm{1}}=\sum a_{\mathrm{nb}}{T}_{\mathrm{nb}}^{n+\mathrm{1}}+\sum a_{\mathrm{nb}}^{\prime }{T}_{\mathrm{nb}}^n+b,\end{array}$$

where the sums are taken over the values at the neighboring (nb) grid points (W, E, N, S, U, D). Putting all the geometric information into factors A_x and V_P (Table 3), the discretization coefficients for the internal control volumes are

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{a}_{\mathrm{W}}=f\mathrm{\Delta }t{A}_{\mathrm{w}}{\stackrel{\mathrm{̃}}{k}}_{\mathrm{w}},a_{\mathrm{E}}=f\mathrm{\Delta }t{A}_{\mathrm{e}}{\stackrel{\mathrm{̃}}{k}}_{\mathrm{e}}\\ {\displaystyle }& {\displaystyle }{a}_{\mathrm{S}}=f\mathrm{\Delta }t{A}_{\mathrm{s}}{\stackrel{\mathrm{̃}}{k}}_{\mathrm{s}},a_{\mathrm{N}}=f\mathrm{\Delta }t{A}_{\mathrm{n}}{\stackrel{\mathrm{̃}}{k}}_{\mathrm{n}}\\ \text{(39)}& {\displaystyle }& {\displaystyle }{a}_{\mathrm{U}}=f\mathrm{\Delta }t{A}_{\mathrm{u}}{\stackrel{\mathrm{̃}}{k}}_{\mathrm{u}},a_{\mathrm{D}}=f\mathrm{\Delta }t{A}_{\mathrm{d}}{\stackrel{\mathrm{̃}}{k}}_{\mathrm{d}}{\displaystyle }& {\displaystyle }{a}_{\mathrm{P}}=V_{\mathrm{P}}{C}_{\mathrm{P}}^n+\sum a_{\mathrm{nb}},a_{\mathrm{P}}^{\prime }=V_{\mathrm{P}}{C}_{\mathrm{P}}^n-\sum a_{\mathrm{nb}}^{\prime }\\ {\displaystyle }& {\displaystyle }b={\mathcal{S}}_{\mathrm{P}}.\end{array}$$

The primed counterparts of a_W, a_E, a_S, a_N, a_U, and a_D are identical, except that f is replaced by (1−f).

Table 3 Geometric factors A_x and V_P appearing in the enthalpy discretization (Eq. 39) for Cartesian and cylindrical coordinate systems. Dimensions of the control volume centered on grid point P are $\mathrm{\Delta }x=(x_{\mathrm{e}}-x_{\mathrm{w}})$, $\mathrm{\Delta }y=(y_{\mathrm{n}}-y_{\mathrm{s}})$, and $\mathrm{\Delta }z=(z_{\mathrm{d}}-z_{\mathrm{u}})$, while the distances between grid points in the X direction are $(\mathit{\delta }x{)}_{\mathrm{w}}=(x_{\mathrm{P}}-x_{\mathrm{W}})$ and $(\mathit{\delta }x{)}_{\mathrm{e}}=(x_{\mathrm{E}}-x_{\mathrm{P}})$. Distances between grid points in Y and Z directions are defined similarly. For radial geometries, $\mathrm{\Lambda }=(r_{\mathrm{e}}^{\mathrm{2}}-r_{\mathrm{w}}^{\mathrm{2}})/\mathrm{2}$.

Consideration of the enthalpy balance shows that the discretization coefficients are slightly different for CVs adjacent to the boundaries of the problem domain. CVPM can be “forced” at the boundaries using either a prescribed temperature (Dirichlet) or heat-flux (Neumann) boundary condition (a convective boundary condition will be introduced in a later version). For a control volume adjacent to a boundary with a Dirichlet boundary condition, a factor of (4∕3) is introduced into the discretization coefficients (Eq. 39) associated with the boundary and the opposing interface. When the heat flux is prescribed on a boundary (Neumann BC), the coefficients associated with the boundary are zero and the specified heat flux appears in discretization coefficient b (Table 4). Boundary conditions along the edges of the problem domain are allowed to vary both spatially and temporally in CVPM.

Table 4 Discretization coefficient b for a control volume adjacent to a prescribed heat-flux (Neumann) boundary condition.

To complete the setup of the discretization coefficients, the material properties must be specified at every grid point within the model domain. Parameters controlling these properties include material type, mean density of matrix particles ρ_m, mineral-grain thermal conductivity $k_{\mathrm{m}}^{\circ }$ at 273.15 K, mineral-grain specific heat $c_p^{\circ }$ at 293.15 K, heat-production rate extrapolated to the surface S₀, heat-production length scale h_s, porosity extrapolated to the surface ϕ₀, critical porosity ϕ_c, compaction length scale h_c, degree of pore saturation S_r, dominant solute type, solute mole fraction in the fully melted system $x_{\mathrm{s}}^{\star }$, interfacial melting parameter λ, mean diameter of larger mode pores d₁, mean diameter of smaller mode pores d₂, and the ratio of the number density of small pores to large pores (n₂∕n₁). The material “type” specifies which governing equations to utilize when finding the heat capacity and thermal conductivity (Sect. 2); types include pure ice and a variety of rocks, soils, organic-rich materials, fluids, and metals. Only a subset of these properties needs to be specified for non-porous layers (e.g., an ice lens or a borehole casing). Examples of the required thermophysical parameters are provided in Sect. 4 and in the CVPM version 1.1 modeling system user's guide (Clow, 2018).

Any temperature field can be used to set the initial temperature condition, including a user-supplied field (e.g., a measured temperature field), a CVPM-determined steady-state field consistent with the boundary conditions and material properties, or a field generated by a previous CVPM modeling experiment.

With the initial condition, boundary conditions, and discretization coefficients specified, the enthalpy-balance equation (Eq. 38) is solved recursively at each time step using the tridiagonal matrix algorithm (TDMA) for 1-D models and the line-by-line method with TDMA for 2-D and 3-D models. Given their temperature sensitivities, the thermophysical properties (ϕ_ℓ, ϕ_i, C, k, etc.) are updated at every time step. In order for the numerical scheme to remain unconditionally stable, all of the discretization coefficients must be non-negative. This consideration leads to the numerical stability condition:

$$\begin{array}{}\text{(40)}& \mathrm{\Delta }t<{\displaystyle \frac{V_{\mathrm{P}}{C}_{\mathrm{P}}^n}{{\displaystyle }(\mathrm{1}-f)\sum A_{\mathrm{nb}}{\stackrel{\mathrm{̃}}{k}}_{\mathrm{nb}}}},\end{array}$$

which must be satisfied in all the CVs at each time step for the scheme to remain unconditionally stable.

Model verification was conducted in two phases. In the first phase, the general model structure and numerical implementation were tested by comparing model results with analytic solutions for a series of simple heat-transfer problems without phase change. Test problems included steady-state and transient boundary conditions, homogeneous and composite media with fixed thermal properties, materials whose thermal properties vary linearly with temperature, and materials with and without radiogenic heating. In all cases, maximum model errors ϵ are on the order of 0.1 mK or less under the test conditions (spatial resolution, time step, etc.). For most cases, max(ϵ) ranges 1 µK to 0.01 pK. Since analytic solutions are unavailable for simultaneously testing all of the model physics, the second testing phase consisted of separately testing each physics module to guarantee it properly simulates the appropriate governing equations (Sect. 2).

4 Example simulations with a sedimentary sequence

Figure 8 Vertical sequence of sedimentary rocks used for the example simulations in Sect. 4.1–4.3. The sequence is assumed to consist entirely of shale below 1 km.

To illustrate the capabilities of the CVPM model, several examples are provided in this section based on the response of a thick vertical sequence of sedimentary rocks to changing boundary conditions. The sequence consists of flat-lying mudrock, carbonate, and sandstone units of various thicknesses (Fig. 8). Values of the parameters controlling the thermophysical properties (Sect. 2) are listed in Table 5. Heat-production rates are from Rybach (1988), while the compaction length scale is loosely derived from values found for a partially exhumed basin on the Arctic slope of Alaska (Burns et al., 2005). Pore spaces are assumed to be fully saturated with water throughout the geologic section [$S_r\equiv (\mathrm{1}-{\mathit{\varphi }}_{\mathrm{a}}/\mathit{\varphi })=\mathrm{1}$]. Sodium chloride, present at relatively low levels when the sediments are completely thawed ($x_{\mathrm{s}}^{\star }\simeq \mathrm{0.003}$), is the dominant pore-water solute.

Table 5 Thermophysical parameters for the sedimentary rock units in Fig. 8. Parameters include thermal conductivity of matrix particles at 0 ^∘C, $k_{\mathrm{m}}^{\circ }$ ($\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$); density of matrix particles, ρ_m (kg m⁻³); specific heat of matrix particles at 20 ^∘C, $c_{p\mathrm{m}}^{\circ }$ ($\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$); heat-production rate extrapolated to surface, S₀ (µW m⁻³); heat-production length scale, h_s (km); porosity extrapolated to surface, ϕ₀; critical porosity, ϕ_c; compaction length scale, h_c (km); degree of pore saturation, S_r; solute mole fraction in the fully melted system, $x_{\mathrm{s}}^{\star }$; interfacial melting parameter, λ (µm K^1∕3); effective diameter of larger mode pores, d₁ (µm); effective diameter of smaller mode pores, d₂ (µm); ratio of the number density of small pores to large pores, n₂∕n₁.

4.1 Permafrost response to ice-age cycles

The first simulation explores the response of the sedimentary sequence to surface-temperature changes over the last ice-age cycle. The upper boundary condition is based on the surface-temperature history determined for the Greenland Ice Sheet during the Holocene and Wisconsin glacial period by Cuffey and Clow (1997) and the Eemian interglacial by the NEEM Community Members (2013). To construct the upper BC for the permafrost simulation, the Greenlandic temperature history was rescaled and shifted to yield an 8 K warming between the Last Glacial Maximum (LGM) and the early Holocene and a mean surface temperature of −11 ^∘C during the 1800s (Fig. 9). For the most recent period, upper boundary temperatures warm 2.5 K between the mid-1800s and 1980 and an additional 2.5 K by the present time, similar to the record for the Arctic Coastal Plain of Alaska (Clow, 2017). The temperature history was replicated back several ice-age cycles to allow for model spin-up. At the lower boundary (depth 2 km), the conductive heat-flux q_b was fixed at 60 mW m⁻², slightly above the continental average. To initialize the model, the subsurface-temperature profile was assumed to be in a steady-state condition at 255 ka with the mean surface temperature over an ice-age cycle. For the spatial grid, 2 m thick control volumes were used above the 550 m depth. Below this, the grid spacing Δz increased progressively to 50 m near the lower boundary. To explore the response of permafrost well below the surface, a 20-year computational time step (Δt) was deemed sufficient.

With the above setup, CVPM was run forward in its 1-D vertical mode from 255 ka to the present. During the last ice-age cycle, the base of permafrost P_d (defined by the 0 ^∘C isotherm) is found to vary by 90 m in this example from 435 to 525 m (Fig. 9). Of greater physical significance is the maximum depth where interstitial ice is present in permafrost. Due to the freezing-point depression caused by interfacial, curvature, pressure, and solute effects, the base of ice-bearing permafrost (B-IBPF) is located 20–27 m above the P_d throughout the simulation. During the most recent ice-age cycle, the B-IBPF and P_d both reached their greatest depths at ∼14 ka, a delay of 10 kyr from the Last Glacial Maximum. Since then, these interfaces have been steadily rising. With the conditions of this simulation, the B-IBPF is currently located at 431 m in a silty claystone, about 22 m below the shallowest depth projected to have occurred following the last interglacial, while the base of permafrost P_d is currently at 467 m in the underlying sandstone unit. Both interfaces are predicted to be rising about 1 cm yr⁻¹ at present, a rate that may be detectable with a carefully designed experiment.

Figure 9 Upper boundary condition (T_s) used to explore the response of a sedimentary sequence (Fig. 8, Table 5) to surface-temperature changes over the last ice-age cycle. Also shown are the depths for the base of permafrost (P_d) and the base of ice-bearing permafrost (B-IBPF) predicted by the CVPM model.

As the simulation confirms, the volume fraction of ice ϕ_i depends strongly on lithology (Fig. 10). This leads to zones with relatively high ice contents in coarse-grained sediments as is often observed in electrical resistivity geophysical logs (Hnatiuk and Randall, 1977; Osterkamp and Payne, 1981). An interesting facet of this behavior is that a pocket of ice-rich sandstone can occur below a fine-grained unit that has little or no ice (e.g., the silty claystone above the sandstone at 450 m in Fig. 10). Except near the B-IBPF, the ice content in coarse-grained materials appears to be relatively stable. In contrast, the volume fractions of ice ϕ_i and unfrozen water ϕ_ℓ are much more dynamic in fine-grained materials over ice-age cycles due to their greater temperature sensitivity. In the upper two-thirds of the permafrost zone, up to 15 % of the water within silty claystones and limestones converts between ice and unfrozen water over ice-age cycles; the percentages are much higher in the lower third of the permafrost zone. Due to the phase change of water, a high heat capacity zone occurs just above the B-IBPF (Fig. 11). This zone, which tracks the B-IBPF over time, is roughly 100 m thick. Pockets of elevated heat capacity also occur in fine-grained materials closer to the surface, especially during periods affected by the warm interglacials. Thermal conductivity variations over an ice-age cycle are also primarily driven by the melting and refreezing of pore ice. Hence, conductivity variations are greatest near the B-IBPF, where conductivity changes of ±15 % occur (Fig. 12). In the upper two-thirds of the permafrost zone, thermal conductivity variations are greater in the fine-grained materials, at least on a percentage basis.

Figure 10 Ice content ϕ_i variations over the last ice-age cycle within the example sedimentary sequence (Fig. 8). Unfrozen water ϕ_ℓ variations (not shown) mirror the ice content fluctuations.

Figure 11 Volumetric heat capacity C variations over the last ice-age cycle within the example sedimentary sequence (Fig. 8). The color mapping corresponds to log(C).

Figure 12 Bulk thermal conductivity k variations over the last ice-age cycle within the example sedimentary sequence (Fig. 8).

4.2 Permafrost response to drilling a deep borehole

We now return to the state of the sedimentary sequence simulated in Sect. 4.1 during the year 1980. By this time, temperatures in the upper 100 m of the sedimentary sequence have warmed in response to the 2.5 K surface warming since the termination of the Little Ice Age (∼1850). At greater depths, temperatures still reflect conditions earlier during the Holocene and the Wisconsin glacial period. With this initial condition, we consider the drilling of a 3 km deep borehole through the example sedimentary sequence (Fig. 8, Table 5) over a 60-day period. Drilling fluids pumped into the proposed 30 cm diameter hole at 30 ^∘C interact thermally with the drill pipe and surrounding rock as they circulate to the bottom of the hole and back to the surface. As a result of the drilling processes, temperatures within the borehole warm within 700 m of the surface. The degree of warming depends on both depth and time as the drill bit advances into the warmer rocks below (Clow, 2015). Figure 13 shows the evolution of temperature changes (ΔT_a) at the borehole wall during drilling based on the Szarka and Bobok (2012) wellbore model. This thermal drilling disturbance, when added to the initial 1980 temperature field, establishes the boundary condition at the borehole wall over the 60-day drilling period. To complete the CVPM setup for a cylindrical 2-D simulation of the permafrost surrounding the borehole, the problem domain was extended radially far enough (40 m) from the hole that the radial heat flux at the outer boundary could be set to zero. The heat flux on the lower boundary was again 60 mW m⁻², while the upper boundary remained at its initial 1980 temperature (−8.5 ^∘C). The vertical grid was identical to that used in Sect. 4.1, while the radial grid spacing increased progressively from 2.5 cm near the borehole to 2 m near the outer boundary. To resolve the rapid temperature changes that occur near the advancing drill bit, the computational time step was set to 0.2 days. The initial condition was provided by the values of all the state variables from the previous example (Sect. 4.1) during 1980.

Figure 13 Temperature change ΔT_a at the borehole wall while drilling a 3 km borehole through the example sedimentary sequence (Fig. 8, Table 5) over a 60-day period. This thermal disturbance, when added to the initial temperature field, provides the boundary condition at the borehole wall during drilling. The drill bit advances from the surface on day 0 to 3 km on day 60.

Running CVPM with the described initial and boundary conditions, the drilling disturbance is found to be great enough in this simulation to melt all of the permafrost ice within 1–2 m of the well by the time the borehole is completed on day 60 (Fig. 14). In this case, borehole electrical resistivity measurements used to detect ice in permafrost must be able to penetrate at least 1.5 m of rock to successfully detect ice. The exact location of the melting front is controlled partially by lithology, with it advancing further from the borehole in fine-grained high-conductivity sediments (e.g., limestones) than in coarser grained low-conductivity mudrocks such as siltstones. Sediment texture is a factor because it affects the volumetric ice content of a material in its undisturbed state. Although the thermal disturbance due to drilling is greatest inboard of the melting front, a substantial disturbance also occurs beyond the front, particularly in the upper couple hundred meters of permafrost, where the initial undisturbed temperatures were −7 to −9 ^∘C (Fig. 15). Outboard of the melting front, the drilling disturbance extends further into the higher conductivity limestone and sandstone units than in the low-conductivity mudrocks. Because of the effect of temperature on the thermal conductivity of minerals, ice, and water, and because of the conversion of ice to liquid water, the bulk thermal conductivity drops about 30 % in the sedimentary units near the wellbore during drilling (Fig. 16). Attempts to infer the thermophysical properties of the sedimentary units from borehole temperature measurements using geophysical inverse methods must carefully consider these changes (Nicolsky et al., 2007).

Figure 14 Volumetric ice content ϕ_i in the sedimentary sequence (Fig. 8) penetrated by a newly completed 3 km deep borehole (day 60). Radial distance is measured from the central axis of the borehole.

Figure 15 Temperatures in the sedimentary sequence (Fig. 8) penetrated by a newly drilled 3 km deep borehole (day 60).

Figure 16 Bulk thermal conductivity k in the sedimentary sequence (Fig. 8) penetrated by a newly drilled 3 km deep borehole (day 60).

4.3 Permafrost response to the formation of a lake

Shallow lakes are ubiquitous on the Arctic Coastal Plain. In thermokarst areas, these lakes are constantly in transition, shrinking, enlarging, draining, and filling new depressions in response to changing temperatures and stream flows. The seasonal ice that forms on these lakes is categorized as “bedfast” ice if it freezes solid to the bottom of the lake, “floating” ice if some liquid remains beneath the ice throughout the winter, and “intermittent” if it is bedfast some years and floating during others. Whether a lake is a bedfast-ice lake or a floating-ice lake depends on whether the maximum seasonal ice-cover thickness ($Z_{\mathrm{ice}}^{\max }$) exceeds the depth of the lake. During the 1970s and 1980s, $Z_{\mathrm{ice}}^{\max }$ for lakes along the Beaufort Sea coast of Alaska was 2.0±0.2 m (Arp et al., 2012; Weeks et al., 1981). By 2010, the maximum seasonal ice thickness had decreased to about 1.5 m due to the warming climate in the region over the last few decades (Bieniek et al., 2014; Wang et al., 2017). Thus, 1.5 m deep lakes that would have been bedfast-ice lakes during the 1970s and 1980s would have transitioned to intermittent-ice lakes by 2010. Arp et al. (2012) recently provided lake–bed temperature data for seasonally ice-covered lakes near the Beaufort Sea coast. Based on these data, lakes whose depth is 0.5–1.0 m less than $Z_{\mathrm{ice}}^{\max }$ have a mean-annual temperature ∼4 K warmer than the surrounding tundra, while intermittent-ice lakes are about 9 K warmer.

Here, we briefly explore the permafrost response over a 35-year period to the instantaneous creation of a 200 m wide lake on the Arctic Coastal Plain of Alaska during 1980. The lake is assumed to be 1.5 m deep with a 100 m wide shallow (1.0 m deep) shelf and is assumed to occur in the same geologic terrain as in the previous two examples. Thus, the lake lies on a 50 m thick silty claystone unit which overlies deeper sedimentary rocks (Fig. 8, Table 5). To investigate the permafrost response, CVPM was used to simulate conditions along a 2-D vertical cross-section passing beneath the lake and surrounding tundra. Initial conditions were identical to those used in the previous example (Sect. 4.2). The upper boundary was set at the 1.5 m depth, i.e., at the bottom of the deeper portion of the lake. Temperatures on this boundary, which were used to force the model, varied depending on location. Beneath the tundra, temperatures were assumed to warm from the initial 1980 surface temperature (−8.5 ^∘C) at 0.75 K decade⁻¹, similar to recent trends observed along the Beaufort coast of Alaska (Clow, 2017; Wang et al., 2017). With these temperatures, the seasonal ice cover is projected to have been bedfast across the entire lake when it was first created and to have transitioned to intermittent over the deeper lake section towards the end of the simulation. Thus, lake–bed temperatures beneath the deeper portion of the lake, assumed to have initially been 4 K warmer than the surrounding tundra, were prescribed to have warmed until they were 9 K warmer than the tundra by 2015. Temperatures beneath the shelf were prescribed to be 4 K warmer than the tundra throughout the simulation. Temperatures beneath the tundra, shelf, and deeper portion of the lake provided the upper BC for the simulation. A heat-flux BC was prescribed on the lower boundary (q_b=60 mW m⁻²), taken to occur at the 300 m depth in this example. The problem domain was extended laterally far enough beyond the lake that the heat flux across the outer boundaries could be set to zero. Vertical grid spacing was Δz=0.1 m in the upper silty claystone unit. Beneath this, Δz increased progressively to 10 m near the lower boundary. Horizontal grid spacing was 1 m beneath the lake and 5 m beneath the surrounding tundra. Although the prescribed upper BC does not include seasonal effects, the computational time step Δt was set to 0.1 years to satisfy the numerical stability condition (Eq. 40).

Figure 17 Simulated permafrost temperatures over a 35-year period following the creation of a 200 m wide lake on a silty claystone unit. Depth is measured relative to the bottom of the deeper portion of the lake.

Running CVPM in its 2-D Cartesian mode for 35 years with the described boundary conditions, temperatures beneath the deeper portion of the lake are found to become warm enough to melt all of the pore ice at the lake–bed interface 19 years after the lake is created (Fig. 17). Thereafter, the melting front propagates downward in the upper silty claystone unit at ∼19 cm yr⁻¹, creating a thaw bulb in its wake (Fig. 18). Lateral migration of the melting front is much more modest: ∼4 cm yr⁻¹ where the deep section adjoins the tundra and ∼6 cm yr⁻¹ where it meets the shallow shelf. By the end of the simulation, a thaw bulb has not yet begun to form beneath the shelf. There are two ramifications of the growing thaw bulb beneath the deeper section of the lake. (1) Fine-grained materials such as silty clay lose much of their mechanical strength once they thaw. In this state, the sides of the lake are much more vulnerable to erosion, which may lead to eventual drainage of the lake. (2) Old carbon stocks stored in the previously frozen permafrost are likely to decompose in the anaerobic thaw bulb and contribute greenhouse gases to the atmosphere (Anthony et al., 2016).

5 Summary and conclusions

This paper presents the governing equations and numerical methods underlying the Control Volume Permafrost Model v1.1 which was designed to relax several of the limitations imposed by previous models. CVPM implements the nonlinear heat-transfer equations in 1-D, 2-D, and 3-D Cartesian coordinates, as well as in 1-D radial and 2-D cylindrical coordinates. To accommodate a diversity of geologic settings, a variety of materials can be specified within the modeling domain, including organic-rich materials, sedimentary rocks and soils, igneous and metamorphic rocks, ice bodies, borehole fluids, and other engineering materials. Porous materials are treated as a matrix of mineral and organic particles with pores spaces filled with liquid water, ice, and air. Functions describing the temperature dependence of the specific heat and thermal conductivity are built into CVPM for a wide variety of rocks and minerals, liquid water, ice, air, and other substances. For porous materials, the bulk thermal conductivity is found using a random two-phase (matrix particles, pores) relationship, while the conductivity of the pores themselves is found using a three-phase (liquid water, ice, air) mixing model. This scheme allows the bulk thermal conductivity to be determined for a wide range of porosities, water saturations ranging 0 %–100 %, and different planetary atmospheres. In addition to the lattice-vibration term (ρc_p), the volumetric heat capacity C depends on a latent-heat term proportional to the change in liquid water content with temperature ($\partial {\mathit{\varphi }}_{\mathrm{\ell }}/\partial T$). At temperatures below 0 ^∘C, the unfrozen water content is found using relationships from condensed matter physics that utilize physical quantities (i.e., particle size) rather than non-physical empirical coefficients requiring calibration. Solute and pore pressure effects are included in the unfrozen water equations. Due to the insolubility of most solutes in ice, impurities become strongly enriched in the liquid pore water as permafrost freezes. To allow for the non-ideal behavior that occurs at high solute concentrations, the water activity a_w is used to find the effect of solutes on the bulk freezing temperature of water T_f. With this approach, solute concentrations up to the eutectic point are allowed. Pore pressure effects on T_f are found in CVPM using either hydrostatic or lithostatic equations, whichever is more appropriate geologically. A radiogenic heat-production term is also included to allow simulations to extend into deep permafrost and underlying bedrock. For the current version of CVPM, liquid water velocities are assumed to be small enough that the associated advective heat flux is negligible compared to the diffusive heat flux.

Figure 18 Volume fraction of ice (ϕ_i) over a 35-year period following the creation of a 200 m wide lake. Magenta (ϕ_i=0) delineates the thaw bulb developing beneath the lake.

Numerical implementation of the governing equations is accomplished using the control-volume approach, allowing enthalpy fluxes to be exactly balanced at control-volume interfaces (e.g., at the interfaces between ice lenses, sedimentary units, bedrock, or a borehole casing). This approach was chosen because the expressions tend to be more accurate than with other methods near boundaries and where strong thermal or physical-property contrasts occur. Very large thermal-property contrasts generally occur near the water freezing point in permafrost. Despite the magnitude of the contrasts and the fact that the freezing front typically migrates over time, the numerical scheme used in CVPM remains stable as long as the stability criterium (Eq. 40) is satisfied. CVPM can be run in either explicit, implicit, or fully implicit modes.

CVPM has been designed for a wide range of scientific and engineering applications, and as an educational tool. The model is “forced” by changes in the boundary conditions at the edges of the problem domain. These conditions include user-prescribed temperatures and/or heat fluxes that are allowed to vary both spatially and temporally along the edges. The model is suitable for use at spatial scales ranging from centimeters to hundreds of kilometers and at timescales ranging from seconds to thousands of years. CVPM can be used over a broad range of depth, temperature, porosity, water saturation, and solute conditions on either the Earth or Mars. Through its modular design, CVPM can act as a stand-alone model or the physics package of a geophysical inverse scheme, or serve as a component within a larger Earth modeling system that may include vegetation, surface water, snowpack, atmospheric, or other models of varying complexity.

One of the goals of CVPM was to eliminate the empirical equations typically used to predict the unfrozen water content at temperatures below 0 ^∘C, replacing them with condensed matter relationships containing quantities that both have a physical meaning and can be determined using relatively simple laboratory techniques or from geophysical logs. The physical quantities in the resulting equations tend to occur within reasonably narrow limits for the material categories defined in CVPM. Thus, if measurements of the physical quantities are unavailable, it may be possible to estimate the behavior of permafrost in a region from a geologic description of the area.