Three-phase numerical model for subsurface hydrology in permafrost-affected regions (PFLOTRAN-ICE v1.0)

Degradation of near-surface permafrost due to changes in the climate is expected to impact the hydrological, ecological and biogeochemical responses of the Arctic tundra. From a hydrological perspective, it is important to understand the movement of the various phases of water (gas, liquid and ice) during the freezing and thawing of near-surface soils. We present a new non-isothermal, single-component (water), three-phase formulation that treats air as an inactive component. This single component model works well and produces similar results to a more complete and computationally demanding two-component (air, water) formulation, and is able to reproduce results of previously published laboratory experiments. A proof-of-concept implementation in the massively parallel subsurface flow and reactive transport code PFLOTRAN is summarized, and parallel performance of that implementation is demonstrated. When water vapor diffusion is considered, a large effect on soil moisture dynamics is seen, which is due to dependence of thermal conductivity on ice content. A large three-dimensional simulation (with around 6 million degrees of freedom) of seasonal freezing and thawing is also presented.


Introduction
The Arctic and sub-Arctic regions of the Earth are warming at a rate significantly faster than the rest of the planet (Turner et al., 2007;Hansen et al., 1999) and are experiencing environmental change at a rapid pace.Permafrost occupies nearly one-quarter of the landmass of the Northern Hemisphere and contains approximately 1600 Gt of organic carbon (Tarnocai et al., 2009).This carbon is potentially available to be released to the atmosphere, thus driving further climate change.However, the timing, rate and chemical form of future carbon releases to the atmosphere are highly uncertain.Much of the uncertainty about mobilization of thawed carbon derives from uncertainty in future soil moisture conditions after anticipated reorganization of permafrost-affected landscapes through permafrost degradation, thaw-induced subsidence and hydrologic processes.In order to predict the amount of carbon released to the atmosphere as well as other adverse effects of permafrost degradation, it is therefore important to have the capability to simulate hydrologic response of permafrost-affected regions to an increase in the mean annual temperatures (Kane et al., 1991;Lunardini, 1996;Osterkamp and Romanovsky, 1999;Schuur et al., 2008).Additionally, it is also important to have simulation capability for assessing vulnerability of structures in cold regions where thawing of permafrost can lead to soil consolidation causing considerable damage.
Published by Copernicus Publications on behalf of the European Geosciences Union.

S. Karra et al.: Three-phase numerical model for subsurface hydrology
These models solve for temperature and ice content using extensions to Richards equation and an equilibrium closure relationship between unfrozen water content and temperature.Closure relationships were empirical in many of these models although closure relations that combine thermodynamic constraints with the unfrozen soil moisture characteristic curves are available and were used in some of the models.All of the aforementioned models focused on small spatial scales (typically on the order of centimeters) and short time frames (hours to days) and many focused on validation against laboratory data from one-dimensional freezing soil columns.These previous studies thus provide much of the prerequisite understanding for predictive capability for cold-region hydrology, and are thus an important complementary approach to permafrost parameterizations in Earth system models (e.g., Lawrence et al., 2011).
To understand the evolution of cold-region hydrologic systems, simulations must address spatial scales of tens of meters to kilometers and multi-decadal time frames.Capability to model at those scales is currently limited.The SUTRA-ICE code (McKenzie et al., 2007) was developed specifically for this purpose and has been used in a number of applications involving groundwater systems that are fully saturated with a combination of ice and liquid.Frampton et al. (2011) used the MarsFlo (Painter, 2011) code to model multidecadal response of a hypothetical partially saturated hydrologic system to a warming trend at the hillslope scale.
Although both SUTRA-ICE and MarsFlo have been used to model application-relevant scales, neither has sufficient flexibility to form a general modeling capability.That a massively parallel implementation is needed follows from the nature of the subsurface thermal hydrologic modeling problem: (a) conditions in the active layer are highly dynamic and respond to seasonal temperature and infiltration forcing, which necessitates a relatively small time step (hours to days); (b) simulations need to span time periods of multiple decades; (c) subsurface thermal hydrology will eventually form one component in a larger multi-process modeling capability (Painter et al., 2012); (d) flows in the active layer may be sensitive to microtopography, thus requiring relatively fine spatial resolution (centimeters to tens of centimeters).
As far as we are aware, there are currently no simulations tools capable of representing the entire range of processes required for modeling hydrology of permafrost-affected regions at the appropriate scale.This is true even if the mechanical and surface processes are neglected and the focus is exclusively on subsurface thermal hydrology.The SUTRA-ICE code is limited to the situation where the pore space is filled with a combination of liquid and ice (i.e., no gas phase) and is thus not appropriate for modeling the dynamics of the active layer, the uppermost layer of soil that freezes and thaws and often partially drains on an annual basis.Water flows in a deepening and partially draining active layer have been identified (Painter et al., 2012) as a key response of permafrost affected regions to warming temperatures, which is why full three-phase capability is a key requirement.Mars-Flo meets that requirement; it is two-component (air, water), uses general three-dimensional and fully unstructured grids, and is capable of representing all possible combinations of the ice, liquid and gas phases in Earth-and Marsrelevant conditions.The generality of MarsFlo was required for the Mars applications that it was originally designed for, which exhibited ice-liquid-gas, ice-liquid, ice-gas, liquidgas, ice-only, liquid-only and gas-only conditions in a single simulation (Grimm and Painter, 2009).However, the general two-component capability is computationally demanding and likely not required for Earth permafrost applications.In addition, neither MarsFlo nor SUTRA-ICE is capable of using massively parallel computing hardware.
Although general requirements for modeling subsurface hydrology in permafrost-affected regions are clear, the computationally demanding nature of the three-phase thermal hydrology simulations places a premium on fine-tuning the process representations as well as the software implementation.This paper addresses the details of process representations and summarizes parallel implementation of a threephase model (PFLOTRAN-ICE version 1.0) for use in projecting hydrologic response of degrading permafrost.Specific questions addressed here include the choice between a Richards-like formulation with passive gas phase and a full two-component formulation, building confidence of the Richards-like formulation by comparing with column experiments, appropriateness of neglecting vapor-phase diffusion, parallel implementation and model initialization strategies.
The balance equations and the solution methodology are described in Sect. 2. Results from a one-dimensional horizontal problem are compared to experimental data in Sect.3. Comparison between the current model and two-component air-water multi-phase model based on Painter (2011) is performed in Sect. 4. The effect of water vapor diffusion on freezing is addressed in Sect. 5. Three-dimensional simulations using the numerical model presented in this paper are shown in Sect.6. Final remarks are provided in Sect.7.

Balance equations
In this formulation, we do not track the movement of air, and hence we do not consider the mass balance for air.With that approximation, which is equivalent to the approximations that lead to Richards equation, balance equations for water and energy are required.The balance of mass and energy for the water component that can be in three phases (liquid, gas, ice) are given by (Painter, 2011) where the subscripts l, i, g denote the liquid, ice and gas phases respectivley; φ is the porosity; s α (α = i, l, g), is the saturation index of the αth phase; η α (α = i, l, g) is the molar density of the αth phase; X α w (α = i, l, g) is the mole fraction of H 2 O in the αth phase; τ g is the tortuosity of the gas phase; D g is the diffusion coefficient in the gas phase; T is the temperature (assuming all the phases and the soil are in thermal equilibrium); c r is the specific heat of the soil; ρ r is the density of the soil; U α (α = i, l, g) is the molar internal energy of the αth phase; H α (α = l, g) is the molar enthalpy of the αth phase; Q e is the heat source; ∇() is the gradient operator and ∇ • () is the divergence operator.
The Darcy velocity for the gas and liquid phases given as follows: where k is the absolute permeability; k rα (α = l, g) is the relative permeability of the αth phase; ρ g , ρ l are the mass densities of the gas and liquid phases; Q w is the mass source of H 2 O; µ α (α = l, g) is the viscosity of the αth phase; p α (α = l, g) is the partial pressure of the α-th phase; g is acceleration due to gravity and z is the vertical distance from a reference datum.In Eq. (1) we account for the compressibility of ice and liquid.Additionally, compressibility of the pore space is accounted for through a relation between porosity and liquid pressure as follows where φ 0 is the porosity of the undeformed soil, p ref is a reference liquid pressure set to 1 atm, c is a compressibility coefficient that is set to 10 −7 Pa −1 .Constraint on the saturations of the various phases of water is given by Furthermore, neglecting the amount of air in liquid and ice phases, we have where X β a (β = l, i) is the mole fraction of air in βth phase, and so Eqs.( 1) and (2), based on the assumption that p g is hydrostatic (i.e., p g = (p g ) 0 − ρ g g z; (p g ) 0 is 1 atm), reduce to In the above formulation, temperature and liquid pressure are chosen to be primary variables.With this approach, one does not have to change the primary variables with a change of phase; such a method, also known as variable switching, is typically used in multi-component, multi-phase systems (e.g., Painter, 2011).

Constitutive relations
In addition to the previously described balance equations, constitutive relations are required to model non-isothermal, multi-phase flow of water.Relations for mole fraction of water vapor, saturations of the phases, thermal conductivity, relative permeability and water vapor diffusion coefficient are specified in this section.
The mole fraction of water in vapor phase is given by the relation, where p v is the vapor pressure and p g is the gas pressure (since we are interested in near-surface regions, for our calculations we shall assume that p g = 1 atm throughout the domain).Assuming thermal equilibrium among the ice, liquid and vapor phases, vapor pressure is calculated using Kelvin's relation (Edlefsen and Anderson, 1943) which includes vapor pressure lowering due to capillary effect as follows where P sat is the saturated vapor pressure, P cgl is the liquidgas capillary pressure, given by P cgl = p g − p l and R is the ideal gas constant.Empirical relations for saturated vapor pressure are used for both above and below freezing conditions (i.e., T = 273.15K).The gas molar density η g is calculated using ideal gas law.
To calculate the partitioning of ice, liquid and vapor phases, at a known temperature and liquid pressure, the following two relations are solved simultaneously for s l and s i (Painter and Karra, 2014):

S. Karra et al.: Three-phase numerical model for subsurface hydrology
Here, S * is the retention curve for unfrozen liquid-gas phases.In these equations, h 0 iw is the heat of fusion of ice at 273.15 K, ρ i is the mass density of ice, ϑ = T − T 0 T 0 and T 0 = 273.15K. Equation ( 9a) is derived assuming that ice can be treated as a solid for the purposes of relating capillary pressure and phase saturations, so the remaining pore space is divided into vapor and liquid phases using the retention curve for unfrozen liquid-vapor.The second relation in Eq. ( 9b) is derived as follows: the first term in the square brackets is the capillary pressure between ice-liquid phases, when gas phase is absent (see Painter and Karra, 2014), and the second term is the addition to the ice-liquid capillary pressure due to the presence of the gas phase.Equations ( 9a) and ( 9b) are derived assuming no freezing-point depression.Furthermore, it has been shown that Painter and Karra (2014) the results based on generalizations of Eqs. ( 9a) and (9b) match well with the experimental results for liquid water content as a function of temperature for different total water content values as measured in Watanabe and Wake ( 2009) and Wen et al. (2012).Although the constitutive equations for calculating the saturations of ice, water and vapor are implicit in nature, closed-form expressions for the derivatives of these saturations with respect to temperature and liquid pressure can be derived, as shown in Appendix A. These derivatives are used for Jacobian evaluation when the partial differential Eq. ( 6) are solved using temperature (T ) and liquid pressure (p l ) as the primary unknown variables.
For S * , we use van Genuchten's model (van Genuchten, 1980), as follows: with the Mualem model (Mualem, 1976) for the relative permeability of liquid water, where λ, α are parameters, with γ = 1 1 − λ .Equation ( 10) assumes that the residual saturation is zero, although a non-zero residual saturation can be easily account for.Additionally, note that from Eq. ( 10), S * is non-zero for finite values of P c .This ensures that complete dry-out does not occur, and that liquid (even if the liquid saturation is very small) is present at all times.
The thermal conductivity for the frozen soil is chosen to be (Painter, 2011) where κ wet,f , κ wet,u are the liquid-and ice-saturated thermal conductivities, κ dry is the dry thermal conductivity, Ke f , Ke u are the Kersten numbers in frozen and unfrozen conditions and are assumed to be related to the ice and liquid saturations by power law relations as follows with α f , α u being the power law coefficients.The gas diffusion coefficient D g is assumed to depend on temperature and pressure as follows: where D 0 g is the reference diffusion coefficient at some reference temperature, T ref , and pressure P ref .

Description of PFLOTRAN
PFLOTRAN (Lichtner et al., 2013) is a massively parallel multi-phase, multi-component, surface-subsurface flow, geomechanics and reactive transport code.The continuum mass, energy (or flow equations) for multiple components including water, supercritical CO 2 , are sequentially coupled to the reactive chemistry equations for a network of geochemical components.The continuum partial differential equations are spatially discretized (for both structured and unstructured grids) using a finite volume technique, and backward Euler scheme is used for time discretization.The discretized equations at each implicit time step reduce to a set of nonlinear algebraic equations which are iteratively solved using an inexact Newton-Krylov method in PFLOTRAN.PFLO-TRAN is written in modular, object-oriented Fortran9X.Parallelization is done using domain decomposition by implementing the PETSc toolkit which takes care of communication between processor core along with providing solvers for the linear and non-linear equations.PFLOTRAN performs parallel I/O via both collective HDF5 calls and direct MPI-IO calls inside the PETSc routines.Information regarding PFLOTRAN installation and user documentation can be obtained from http://www.pflotran.organd PFLOTRAN-ICE v1.0 which has the ice module can be dowloaded freely from https://bitbucket.org/satkarra/pflotran-ice-release-v1.0.

Solution methodology
The system (Eqs.6a and 6b) can be written in the form (assuming no source/sink) where A, F are the accumulation and flux terms.Equation (15) is discretized using finite volume method with backward Euler temporal discretization, to obtain the following form: where the superscript i denotes the time step, the subscript n denotes the cell n, and n being the neighboring cell to cell n, A nn denotes the area of the interface between the cells n and n , V n is the volume of the cell n, A n is the accumulation term in the nth cell and F nn is the flux term across the interface between the cells n and n .Finite difference is used to calculate the gradients in the flux term, and the material properties in F nn is based on the intercell average.The gradient terms in F are discretized using a two-point flux approximation between the neighboring grid cells.This requires the flux to be orthogonal to the face common to the grid cells.The two discretized set of equations are solved in a fully coupled fashion using inexact Newton-Krylov method.The calculations shown in Appendix A are used for the calculation of Jacobian (needed in the Newton-Krylov method), namely, for derivatives of A n and F nn with respect to p l and T .The intercell averages for the flux terms are calculated as follows: k is harmonic averaged, (τ φ) is harmonic averaged, η l , ρ l , k r are upwinded, κ is harmonic averaged.For the gas diffusion term, the coefficient for the gradient (φ s g τ g η g D g ) is chosen based on values from the cell with smaller X w g (see Painter, 2011, for the reasoning behind this averaging), η g is calculated using the temperature in the cell assuming ideal gas law and by assuming that the gas pressure is 1 atm.
The above solution methodology is used to implement the governing equations in Sect.2.1 in a massively parallel fashion in PFLOTRAN.At each grid cell, in addition to solving for liquid pressure and temperature, an inert tracer concentration is also solved for.This involves solving an additional advection-diffusion equation, although this concentration is not used in this paper.Thus for each grid cell, three degrees of freedom (or unknowns) are solved for.Figure 1 shows the parallel scalability in terms of strong scaling (that is, for a fixed problem size, the solution time is measured as a function of number of processor cores) of the implementation in PFLOTRAN without any input or output.Two cases with 3 million and 12 million degrees of freedom (dofs) are considered.For the 3 million case, the code scales up to 1024 processor cores with about 2929 dofs per processor core, while for the 12 million case, it scales up to 2048 processor cores with approximately 5850 dofs per processor core.The degradation in the performance beyond 1024 for 3 million dofs and 2048 processor cores for 12 million dofs is attributed to decrease in work load per processor core, thus increasing the parallel communication time over computation time.Additionally, it is also known that increase in the processor cores increases the number of linear solver iterations needed for convergence (Hammond et al., 2014); this adds to the degradation of performance at larger processor core count.To evaluate the performance of the implementation on small clusters that are more easily available to researchers, we looked at the strong scaling on a computer with 32 floating point units.Four cases of varied degrees of freedom -12 000, 24 000, 48 000 and 96 000 -were considered.Each of the four cases seem to scale reasonably up to 8, 16, 16 and 32 processor cores, respectively, with up to approximately 1500, 1500, 3000 and 3000 dofs per processor core.The degradation beyond these processor cores is due to larger communication compared to computation.

Comparison with experimental data
In this section, we shall compare the numerical results against the experimental data from Jame and Norum (1980) for a partially saturated porous medium.The Jame and Norum experimental setup was as follows: a 30 cm long horizontal tube with partially saturated no.40 silica flour sealed at the ends was used.The sample was initially unfrozen and the temperature at one end was lowered while maintaining the temperature at the other end at the initial temperature.Total water content (ice plus water) was measured at different times using gamma ray attenuation.Results from three tests were reported in Jame and Norum (1980).In the first test, the sample had a water content of 15.6 % (by dry weight), with an initial temperature of 20 • C, and the temperature at the cold end set to −10 • C. For the second test, a water content of 15 %, an initial temperature of 5 • C, and a cold end temperature of −5 • C, was used.Finally, in the third test, a water content of 9.5 %, an initial temperature of 5 • C, and a cold end temperature of −5 • C, was used.Figure 2 shows experimental and numerical results for the water content (by dry weight) as a function of position for 6, 24 and 72 h.The temperature profiles are also compared at the three instances in time.The same set of parameters, listed in Fig. 2, were used for all the three tests.A good comparison can be seen between the numerical results and experiments for both water content as well as the temperature.The reason behind the differences seen in the water content at the cold end of the tube is unknown but similar differences have been seen previously by others as well (Jame and Norum, 1980; White and Oostrom, 2006;Painter, 2011).

Our approach vs. two-component approach
In this section, two configurations are considered to compare the results from the current approach with the twocomponent (air-water) approach based on Painter (2011).

1-D horizontal domain
First, we shall consider the one-dimensional horizontal experiment by Jame and Norum (1980) discussed in Sect.3. The comparison between PFLOTRAN and the results from a two-component approach are shown in Fig. 3.The initial and boundary conditions used are same as the ones used in the experiments, discussed in Sect.3. Overall a good match can be seen with minor differences in the solution at the boundaries and at the freezing front.This demonstrates that the single-component Richards model is adequate for this application.

2-D domain
In the one-dimensional simulations summarized in Sect.4.1, the single-component model gave very similar results to the more complete two-component model (Painter, 2011) that accounts for advective transport of water vapor.However, a comparison between the two models in a one-dimensional configuration is not very demanding because excursions in gas-phase pressure, which are neglected in the Richardsbased model but may occur during freeze up in a twocomponent model, are not able to induce significant advective transport of water vapor in one-dimensional configurations.Numerical experiments in a two-dimensional configuration provide a more sensitive test of the adequacy of the single-component model.
The domain for the two-dimensional tests is rectangular with depth of 20 m and horizontal extent of 50 m.The regular grid spacing is 1 m in the horizontal and 0.2 m in the vertical.The initial conditions, thermal properties, and mean annual surface temperature are selected to cause the active layer depth to be at approximately 1 m with fully saturated frozen soil below that depth.No flow conditions are applied on the left and bottom boundary.The top is specified as an infiltration boundary (specified infiltration rate and temperature in the single component model; specified infiltration rate, temperature and gas pressure in the two-component model).A cyclic temperature condition representing seasonal variations is applied at the top.Infiltration is applied when the temperature is above freezing; no infiltration is applied when the temperature is below freezing.The temperature on the boundary on the right face is held at 2 • C between depths of 1 and 2 m, mimicking a talik.The boundary condition for flow in that region of the right boundary corresponds to a seepage face.
The boundary and initial conditions in this twodimensional simulation are designed to cause a shallow perched aquifer to form in the active layer during summer.Water then flows toward the right seepage face.The simulations are designed to test whether gas pressure excursions induced by soil freezing in fall, which are not represented in the single-phase passive gas model, will enhance lateral water and vapor flow.Comparisons between the single-component, passive-gas model and the more complete two-component   2011).The curves are liquid saturation vs. horizontal distance (the talik is on the right) at depths of 10, 30, 50 and 70 cm.The twocomponent model does show more lateral movement, but the differences are quite small (note the narrow range on the y axis).For these and similar comparisons, it can be concluded that the single-component, passive-gas approximation is adequate for the purposes of modeling water dynamics in Earth permafrost.This is in contrast to applications involving the hydrologic system of Mars, which were found to be sensitive to advective transport of water in the vapor phase (Grimm and Painter, 2009).

Effect of vapor diffusion
To study the effect of vapor diffusion on the formation and evolution of permafrost, a one-dimensional vertical column of height 30 m was considered.The domain was initialized with a water table at a height of 15 m and a temperature of 1 • C. A geothermal heat flux of 100 mW m −2 was applied along with a no flow boundary condition at the bottom of the domain.A temperature of −5 • C was applied at the top with no infiltration.The simulation was run to 3000 year.The temperature and ice saturation profiles for cases with and without vapor diffusion are shown in Figs. 5 and 6.For the case without vapor diffusion, as the temperature in the vadoze zone between z = 15 and z = 20 dropped below freezing, the vapor converted into ice, and a thin ice layer starts to form.The position and thickness of the ice layer does not change significantly as a very small increase in the ice content is seen.On the other hand, for the case with diffusion, the thickness of the ice layer increases with time.Also, the fraction of ice in this layer can be seen to increase significantly.This is due to two mechanisms: the first being that the vapor layer below the ice layer diffuses to the bottom of the ice layer which is cooler as seen in Fig. 6b, and second that a feedback from soil thermal conductivity causes further decrease in temperature, which in turn increases ice layer thickness as well as ice content.This feedback from soil thermal conductivity is primarily due to its dependence on ice saturation.Furthermore, for the case with diffusion, as seen in Fig. 6b the diffusion of vapor to a cooler region of the domain causes the height of the water table to decrease.Note that in both cases the temperature at the bottom of the domain increases due to the geothermal flux.winter the soil is completely filled with ice and as the temperature on the top region warms in spring, the ice in the top melts.In summer, the ice melts to a depth of around 0.8-1 m.As the top temperature cools down in the fall season, the ice layer starts to freeze from the top to an essentially completely frozen state in winter.Reasonably high amounts of gas are seen in the top layer during spring, summer and early fall seasons with the peak being in summer.Previous models such as SUTRA-ICE cannot capture this effect since gas is not tracked in their formulation.Figure 8b also clearly shows the formation of ice in the topmost part of the domain in early fall.Additionally, a point to be noted from this simulation is that even though ice was initially present in the entire domain as an initial condition, a cyclic profile was reached fairly quickly (in about 5 year) and then the active layer generally followed the surface topography.

Discussion and conclusions
Numerical models are increasingly being used to help understand how subsurface hydrology in permafrost-affected regions will respond to increasing air temperatures and changes in precipitation.Such models generally fall into two classes.One class focuses on groundwater systems at large scales with approximate treatment of active layer and intrapermafrost physics (e.g., McKenzie et al., 2007;Bense et al., 2009;Bosson et al., 2013;Vidstrand et al., 2013;Grenier et al., 2013;McKenzie and Voss, 2013).The second class includes more realistic descriptions of water dynamics in the active layer, including the effects of non-zero gas content (e.g., Painter, 2011;White, 1995).However, those models have generally been limited to relatively small scales (e.g., the column scale) and one spatial dimension because of computational demands of the three-phase models.The implementation described here takes advantage of highly scalable parallel subsurface multi-physics capability in PFLOTRAN (Lichtner et al., 2013); thus, enabling an important class of applications involving degradation of ice-wedge polygon bogs that require both three-phase physics and relatively large domain sizes (Painter et al., 2012).The implementation described here represents a singlecomponent (water substance) partitioned over three-phases (ice, liquid, vapor) coupled with an energy balance equation.The single-component multi-phase formulation gives nearly identical results to the more complete two-component formulation (Painter, 2011) for applications of interest.Thus, the less demanding single-component model is preferred for applications involving hydrology of Earth permafrost.However, Mars applications (e.g., Grimm and Painter, 2009) will generally require the two-component model.
Successful comparisons with laboratory freezing-column experiments build confidence in both the numerical implementation and the constitutive model (Painter and Karra, 2014) for partitioning among ice, liquid and gas phases.In the constitutive model used here, the partitioning among the three-phases follows from information about the soil water characteristic curve in unfrozen conditions.This is preferable to purely empirical freezing curves, as those empirical freezing curves would need to be developed anew for each application in contrast to the soil water characteristic curve, which may be estimated from information about soil texture.
Although the gas phase is passive in the implementation described here, as it is in Richards equation, diffusion of water vapor is included.In our one-dimensional simulations of Sect.5, vapor diffusion had a surprisingly large effect on the subsurface soil moisture dynamics in unsaturated conditions.The sensitivity to the vapor diffusion process results partially from a dependence of the thermal conductivity model on ice content.As vapor diffuses to cold regions and cold traps as ice, the thermal conductivity increases, which decreases the soil temperature during winter and further increases vapor cold trapping.However, the vapor diffusion model used here is approximate.Further evaluation of the importance of vapor diffusion for Arctic soils using better vapor diffusion models (e.g., Webb and Ho, 1998) is thus needed.
Our code (PFLOTRAN-ICE v1.0) can be used for length scales ranging from the order of a representative element volume (Bear, 2013) to kilometers, where resolving the physics in more detail is important.Existing permafrost thaw representation (Lawrence et al., 2011) in codes such as CLM (Oleson et al., 2010), for instance, can be used for even larger domains of the order hundreds of kilometers but such codes assume a one-dimensional vertical flow and heat transport, and do not account for lateral flow.For example, if one is interested in looking at the hydrology in a polygonal region or a small set of such regions, CLM cannot model the flow between troughs and centers while our code can be used to resolve the non-isothermal flow in these regions.Furthermore, our code can be used for comparing with larger-scale codes such as CLM to build confidence in these codes which have simplified physics.
The work described here focuses only on the subsurface hydrology without consideration of surface flows.As Painter et al. ( 2012) discuss, a comprehensive modeling capability for hydrology in permafrost-affected regions will also require representation of surface flow, surface energy balance, and evolution of topography caused by thawing of permafrost and melting of ground ice.Those important couplings will be addressed in the future.For known mass and energy fluxes (m 0 , e 0 ), Eq. ( C4b) can be used to solve for temperature (T ) as a function of z.Using this temperature profile, liquid pressure (p l ) can be then evaluated using Eq.(C4a).Once p l , T are known as functions of z, liquid, ice and water vapor saturations can be evaluated using Eq. ( 9).

Figure 1 .
Figure 1.Strong scaling performance of PFLOTRAN using Jaguar Cray XK6 supercomputer at Oak Ridge National Laboratory for the non-isothermal, multi-phase (ice, vapor and liquid) sub-surface water flow problem (no I/O).Domain sizes with 3 million and 12 million dofs are considered.For the 3 million dofs case, the code scales up to 1024 processor cores with about 2929 dofs per processor core, and for the 12 million dofs case, it scales up to processor cores with approximately 5850 dofs per processor core.

Figure 7 .
Figure 7. Three-dimensional domain based on surface topography measured using lidar from Barrow, AK (C.Tweedie, personal communication, 2012).The size of the domain in the horizontal plane is 25 m × 25 m and the height variation is between 4.2-4.6 m.

Figure 9 .
Figure 9. Ice thawing and freezing with seasonal surface temperature variation (continued).Ice and gas saturations for summer and early fall shown here.For the values of parameters used see Fig. 8.

Figure 10 .
Figure 10.Ice thawing and freezing with seasonal surface temperature variation (continued).Ice and gas saturations for peak fall shown here.For the values of parameters used see Fig. 8.

S.
Karra et al.: Three-phase numerical model for subsurface hydrology Integrating Eqs.(C1a) and (C1b), we get v l η l − φs g τ g η g D e 0 are constant mass and energy fluxes.The mole fraction of water vapor X g w can be calculated using (without including the lowering factor due to capillary effects)