Coupled cellular automata for frozen soil processes

Heat and water movement in variably saturated freezing soils is a strongly coupled phenomenon. The coupling is a result of the effects of sub-zero temperature on soil water potential, heat carried by water moving under pressure gradients, and dependency of soil thermal and hydraulic properties on soil water content. This study presents a one-dimensional cellular automata (direct solving) model to simulate coupled heat and water transport with phase change in variably saturated soils. The model is based on first-order mass and energy conservation principles. The water and energy fluxes are calculated using first-order empirical forms of Buckingham–Darcy’s law and Fourier’s heat law respectively. The liquid–ice phase change is handled by integrating along an experimentally determined soil freezing curve (unfrozen water content and temperature relationship) obviating the use of the apparent heat capacity term. This approach highlights a further subtle form of coupling in which heat carried by water perturbs the water content–temperature equilibrium and exchange energy flux is used to maintain the equilibrium rather than affect the temperature change. The model is successfully tested against analytical and experimental solutions. Setting up a highly non-linear coupled soil physics problem with a physically based approach provides intuitive insights into an otherwise complex phenomenon.


Introduction
Variably saturated soils in northern latitudes undergo repeated freeze-thaw cycles.Freezing reduces soil water potential considerably because soil retains unfrozen water (Dash et al., 1995).The resulting steep hydraulic gradients move considerable amounts of water upward from deeper warmer soil layers which accumulates behind the freezing front.The resulting redistribution of water alters soil thermal and hydraulic properties, and transports heat from one soil zone to another.As water freezes into ice, the latent heat maintains soil temperatures close to 0 • C for long periods of time.The water and energy redistribution has significant implications for regional hydrology, infrastructure and agriculture.Understanding the physics behind this complex coupling remains an active area of research.Field studies have been widely used to better understand the mecha-nism of these thermohydraulic cycles (e.g.Hayashi et al., 2007).Innovative column studies under controlled laboratory settings have allowed for further insights by isolating the effects of factors that drive soil freezing and thawing, a separation impossible to achieve in the field (e.g.Nagare et al., 2012).Mathematical models, describing the mechanism of water and heat movement in variably saturated freezing soils, have been developed to complement these observational studies.Analytical solutions of freezing and thawing front movements have been developed and applied (e.g.Stefan, 1889;Hayashi et al., 2007) and numerical models have replicated the freezing-induced water redistribution with reasonable success (e.g.Hansson et al., 2004).The optimization of existing numerical modelling approaches also remains an active area of research.For example, improvements to numerical solving techniques to address sharp changes in soil properties, especially behind freezing and thawing fronts, Published by Copernicus Publications on behalf of the European Geosciences Union.and during special conditions such as infiltration into frozen soils have been reported recently (e.g.Dall'Amico et al., 2011).
Although the coupling of heat and water movement in variably saturated freezing soils is complex, fundamental laws of heat and water movement coupled with principles of energy and mass conservation are able to explain the physics to a large extent.There is a paradigm shift in modelling of water movement in variably saturated soils using physically based approaches.For example, HydroGeoSphere and Parflow (Brunner and Simmons, 2012;Kollet and Maxwell, 2006) are examples of codes that explicitly use Richard's equation to model subsurface flow.Thus, the use of derived terms such as specific yield is not required.Mendicino et al. (2006) reported a three-dimensional (3-D) CA (cellular automata; direct solving) model to simulate moisture transfer in the unsaturated zone.Cervarolo et al. (2010) extended the application of this CA model by coupling it with a surface-vegetation-atmosphere-transfer scheme to simulate water and energy flow dynamics.Direct solving allows for unstructured grids while describing the coupled processes based on first-order equations.Use of discrete first-order formulations allow one to relax the smoothness requirements for the numerical solutions being sought.This has advantages, particularly in large-scale models, wherein use of relatively coarse spatial discretization may be feasible.Therefore, it is important to expand the application of direct solving to further complicated unsaturated soil processes.
This study presents a coupled CA model to simulate heat and water transfer in variably saturated freezing soils.The system is modelled in terms of the empirically observed heat and mass balance equations (Fourier's heat law and Buckingham-Darcy equation) and using energy and mass conservation principles.The liquid-ice phase change is handled with a total energy balance including sensible and latent heat components.In a two-step approach similar to that of Engelmark and Svensson (1993), the phase change is brought about by the residual energy after sensible heat removal has dropped the temperature of the system below freezing point.Knowing the amount of water that can freeze, the change in soil temperature is then modelled by integrating along the soil freezing curve.To our knowledge, coupled cellular automata have not yet been used to explore simultaneous heat and water transport in frozen variably saturated porous media.The model was validated against the analytical solutions of (1) the heat conduction problem (Churchill, 1972), (2) steady state convective and conductive heat transport in unfrozen soils (Stallman, 1965), (3) unilateral freezing of a semi-infinite region (Lunardini, 1985), and (4) the experimental results of freezing-induced water redistribution in soils (Mizoguchi, 1990).

Cellular automata
Cellular automata were first described by von Neumann in 1948 (see von Neumann andBurks, 1966).The CA describe the global evolution of a system in space and time based on a predefined set of local rules (transition rules).Cellular automata are able to capture the essential features of complex self-organizing cooperative behaviour observed in real systems (Ilachinski, 2001).The basic premise involved in CA modelling of natural systems is the assumption that any heterogeneity in the material properties of a physical system is scale dependent and there exists a length scale for any system at which material properties become homogeneous (Hutt and Neff, 2001).This length scale characterizes the construction of the spatial grid cells (elementary cells) or units of the system.There is no restriction on the shape or size of the cell with the only requirement being internal homogeneity in material properties in each cell.One can then recreate the spatial description of the entire system by simple repetitions of the elementary cells.The local transition rules are results of empirical observations and are not dependant on the scale of homogeneity in space and time.The basic assumption in traditional differential equation solutions is of continuity in space and time.The discretization in models based on traditional numerical methods needs to be over grid spacing much smaller than the smallest length scale of the heterogeneous properties making solutions computationally very expensive.The CA approach is not limited by this requirement and is better suited to simulate spatially large systems at any resolution, if the homogeneity criteria at elementary cell level are satisfied (Ilachinski, 2001;Parsons and Fonstad, 2007).In fact, in many highly non-linear physical systems such as those describing critical phase transitions in thermodynamics and the statistical mechanical theory of ferromagnetism, discrete schemes such as cellular automata are the only simulation procedures (Hoekstra et al., 2010).
On the contrary, explicit schemes like CA are not unconditionally convergent and hence given a fixed space discretization, the time discretization cannot be arbitrarily chosen.Another limitation of the CA approach was thought to be the need for synchronous updating of all cells for accurate simulations.However, CA models can be made asynchronous and can be more robust and error resistant than a synchronous equivalent (Hoekstra et al., 2010).
The following section (Sect.2.1) describes a 1-D CA in simplified, but precise mathematical terms.It is then explained with an example of heat flow (without phase change) in a hypothetical soil column subjected to a time varying temperature boundary condition.

Mathematical description
Let S t i be a discrete state variable which describes the state of the ith cell at time step t.If one assumes that an order of N elementary repetitions of the unit cell describe the system spatially, then the complete macroscopic state of the system is described by the ordered Cartesian product S t 1 ⊗ S , where G t is the global state variable of the system defining the physical state of the system at time t.Given this algebra of the system, G t+1 is given by where The quantity r is generally called the radius of interaction and defines the spatial extent on which interactions occur on the local scale.In the case of the 1-D CA, the only choice of neighbourhood which is physically viable is the standard von Neumann neighbourhood (Fig. 1).

Physical description based on a heat flow problem in a hypothetical soil column
Let us consider the CA simulation of heat flow in a soil column of length L c and a constant cross-sectional area.The temperature change in the column is driven by a time varying temperature boundary condition applied at the top.It is assumed that no physical variation in the soil properties exist in the column at length intervals smaller than x.Each cell in the 1-D CA model can therefore be assumed to be of length x.Therefore, the column can be discretized using L c / x elementary cells.To simulate the spatio-temporal evolution of soil temperature in the column, an initial temperature for each elementary cell has to be set.To study the behaviour of the soil column under external driving (time varying temperature), a fictitious cell is introduced at the top and/or the bottom of the soil column and subjected to time varying temperatures.The transition rules need to be defined now.Once the transition rules of heat exchange between neighbours are defined, the fictitious boundary cells interact with the top and/or bottom cells of the soil column as any other internal cell based on the prescribed rules and the predefined temperature time series.Although the same set of rules govern interaction among all cells of the column, heat exchange cannot affect the temperature of the fictitious cells as that would corrupt the boundary conditions.This is handled by assigning infinite specific heats to the fictitious cells.This allows evolution of the internal cells and the boundary cells according to the same mathematical rules/empirical equations.The preceding mathematical description of the CA algebra is based on the assumption that the state variable defining each cell is discrete in space and time.But soil temperatures are considered to be continuous in space and time.The continuous description of the soil temperature can be adapted to the CA scheme by considering small time intervals over which the temperature variations are not of interest and hence for all practical purposes can be assumed constant.Conditions for convergence of the numerical temperature profile set an upper limit on the size of this time interval for a given value of x.Therefore, once the length scale of homogeneity x in the system and the local update rules have been ascertained, the CA is ready for simulation under the given initial and boundary conditions.Equations (2) and (3) (Sect.3), applied sequentially, would be the local update rules for this simple case of heat flow in a soil column (without phase change) driven by time varying temperatures at the top.
The meaning of the terms used in the mathematical description of CA can now be explained with respect to the heat flow simulation for the hypothetical soil column: S t i is the temperature of the ith cell at time t, r = 1, φ is a sequential application of Eqs. ( 2) and (3) describing heat loss/gain by a cell due to temperature gradients with its two nearest neighbours and temperature change due to the heat loss/gain respectively, and χ is the identity mapping.

Coupled heat and water transport in variably saturated soils
The algorithm developed for this study simultaneously solves the heat and water mass conservation in the same time step.
The implementation is based on the assumption of nearest neighbour interactions, i.e. r = 1.The one-dimensional conductive heat transport in variably saturated soils can be given by the heat balance equation where subscripts i and ζ refer to the cell and its active neighbours, q h is the net heat flux (J s −1 m −2 ) for the ith cell, T is cell temperature ( • C), λ i,ζ is average effective thermal conductivity of the region between the ith and the ζ th cells (J s −1 m −1 • C −1 ), and l i,ζ is the distance between the centres of the ith and the ζ th cells (m).Effective thermal conductivity can be calculated using one of the popular mixing models (e.g.Johansen, 1975;Campbell, 1985).The empirical relationship between heat flux from Eq. ( 2) and the resulting change in cell temperature ( where l i is the length of the cell (m) and C i (J m −3 • C −1 ) is the effective volumetric heat capacity of the cell such that where θ is volumetric fraction (m 3 m −3 ) and subscripts w, ice, s, and a represent water, ice, soil solids and air fractions.
The mass conservation equation in 1-D can be written as where ρ is density (kg m −3 ), is the total volumetric water content (m 3 m −3 ), q w is the Buckingham-Darcy flux (m s −1 ), and S s is sink/source term.In unfrozen soils, θ ice = 0 and = θ w .Buckingham-Darcy's equation is used to describe the flow of water under hydraulic head gradients wherein it is recognized that the soil matric potential (ψ) and hydraulic conductivity (k) are functions of liquid water content (θ w ).The dependency of ψ and k on θ w can be expressed as a constitutive relationship.The constitutive relationships proposed by Mualem-van Genuchten (van Genuchten, 1980) defining ψ(θ w ) and k(θ w ) are used in this study: where θ res (m 3 m −3 ) is the residual liquid water content, η (m 3 m −3 ) is total porosity, K s (m s −1 ) is the saturated hydraulic conductivity, and α (m −1 ), n and m are equation constants such that m = 1 − 1/n.For an elementary cell in a 1-D CA model, the Buckingham-Darcy flux in its simplest form can be written as where all subscripts have the same meaning as introduced so far, z is the cell elevation and k represents the average hydraulic conductivity of the region between the ith and the ζ th cells.In this study, phase change and associated temperature change is brought about by integrating along a soil freezing curve (SFC).SFCs can be defined because the liquid water content in frozen soils must have a fixed value for each temperature at which the liquid and ice phases are in equilibrium, regardless of the amount of ice present (Low et al., 1968).Soil freezing curves for different types of soils developed from field and laboratory observations between liquid water content and soil temperature have been widely reported (e.g.Anderson and Morgenstern, 1973;Stähli and Stadler, 1997).Van Genuchten's model can be used to define a SFC (Eq.7), wherein ψ(θ w ) is replaced with T (θ w ), and n, m and α ( • C −1 ) are equation constants.

The coupled CA model
Figure 2 shows a flow chart describing the algorithm driving the coupled CA code.The code was written in MATLAB ® .Let the superscript t denote the present time step and subscript i be the spatial index across the grid where each node represents centres of the cell.The thermal conduction and hydraulic conduction modules represent two different algorithms that calculate the net heat (q h,i ) and water (q w,i ) fluxes respectively across the ith cell.In essence, the thermal conduction and hydraulic conduction codes run simultaneously and are not affected by each other in the same time step.However, the processes are not independent and are coupled through updating of model parameters and state variables at the end of each time step.Hydraulic conduction is achieved by applying Eq. ( 10) to each elementary cell using the hydraulic gradients between it and its immediate neighbours (r = 1).Similarly, Eq. ( 2) is used to calculate the heat flux between each elementary cell and its immediate neighbours using the corresponding thermal gradients.The change in mass due to the flux q w,i is used to obtain change in pressure head ( ψ i = ψ t+ t i − ψ t i ) from the ψ(θ w ) relationship.The updated value of total water content is then used to update the volumetric heat capacity C i (Eq.4).The updated value of C i is used as an input to the energy balance module along with computed heat flux q h,i .This represents the first stage of coupling between hydraulic and thermal processes.The energy balance module computes the total change in ice and water content due to phase change, and the total temperature change ( T i ) due to a combination of thermal conduction and phase change.The energy balance module is explained using an example of a system wherein the soil temperature is dropping and phase change may take place if cell temperature drops below the freezing point of pure water (T fw = 0 • C).Inside the energy balance module, the change in temperature ( T i ) is calculated using Eq.(3) and values of C i and Q h,i assuming that only thermal conduction takes place.If the computed T i for a given cell is such that T t+ t i ≥ T fw , then water cannot freeze; cell temperatures are updated without phase change and the code moves into the next time step.In the approach of this study, phase change and associated temperature change can occur if and only if the present cell temperature (T i ) and water content (θ w,i ) represent a point on the SFC.This point along the SFC (Fig. 3) is defined here as the critical state point (T crit , θ wcrit ).If T i gives T t+ t i < T fw for any cell, then freezing point depression along the SFC accounts for change in temperature due to freeze-thaw.The freezing point depression or T crit is defined for the cell by comparing the cell θ w,i with the SFC.However, the coupled nature of heat and water transport in soils perturbs the critical state from time to time, e.g. when freezing induces water movement towards the freezing front or infiltration into frozen soil leads to accumulation or removal of extra water from any cell.In such a case, Q h,i needs to be used to bring the cell to the critical state.This may require thermal conduction without phase change (T crit > T i ) or freezing of water without temperature change (T crit < T i ).This process gives us an additional change in temperature or water con- tent which is purely due to the fact that the additional water accumulation disturbs the critical state.This is another and a subtle form of coupling between heat and water flow.Because of the above consideration to perturbation of critical state caused by additional water added/removed from a cell, infiltration into frozen soils during the over-winter or spring melt events need no further modifications to the process of water and heat balance.
If Q h,i is such that a cell can reach critical state and still additional heat needs to be removed, then this additional heat (Q res,i ) removal leads to the freezing of water.The freezing of water leads to change in the temperature of the cell such that where L f is the latent heat of fusion (334 000 J kg −1 ) and T new,i is the new temperature of the cell (Fig. 3).If the change in water content due to freezing is such that θ w,i = θ res , then no further freezing of water can take place and Q res,i is used to decrease the temperature of the cell using Eq. ( 3) and the updated value of C i (i.e. after accounting for change in C i due to phase change).The soil thawing case is exactly similar as described above; the only dissimilarity is that a different SFC may be used if hysteric effects are observed in SFC paths as observed in studies by Quinton and Hayashi (2008), and Smerdon and Mendoza (2010).If the cell temperature is above freezing, then the matric potential is calculated using Eq. ( 7).For cell temperatures below freezing point, the water pressure (matric potential) can be determined by the general-R.M. Nagare et al.: Coupled cellular automata for frozen soil processes ized Clausis-Clapeyron equation by assuming zero ice gauge pressure: where g is acceleration due to gravity (9.81 m s −2 ).At the end of the energy balance calculations, temperatures of all the cells are updated using the T i computed in the energy balance module.Water content for each cell is updated by considering the change due to freeze/thaw inside the energy balance module and q w,i .Hydraulic conductivity of each cell is updated (Eq.8) using the final updated values of water content.Pressure and total heads in each cell are updated considering water movement (Eq.7) and freezing/thawing (Eq.12).The volumetric heat capacity of each cell is updated one more time (Eq.4) to incorporate the changes due to freeze/thaw inside the energy balance module.Thermal conductivity of each cell is updated using a mixing model (e.g.Johansen, 1975).This completes all the necessary updates and the model is ready for computations of the next time step.
The CA scheme described here is not unconditionally convergent.Hence, the size of the time step cannot be arbitrarily chosen.In our implementation of the CA model, adaptive time stepping has been achieved following the convergence analysis reported in Appendix A.

Heat transfer by pure conduction
The ability of the CA model to simulate pure conduction under hydrostatic conditions was tested by comparison to the analytical solution of one-dimensional heat conduction in a finite domain given by Churchill (1972).A soil column with total length (L c ) of 4 m was assumed to have different initial temperatures in its upper (T u = 10 • C) and lower (T l = 20 • C) halves (Fig. 4).The system is hydrostatic at all times and there is no flow.At the interface, heat conduction due to the temperature gradient will occur until the entire domain reaches an average steady state temperature of 15 • C. The analytical solution given by Churchill (1972) can be expressed as The parameters used in analytical examples for Churchill (1972), and the CA code are given in Table 1.There is excellent agreement between the analytical solution and the CA simulation (Fig. 4).

Heat transfer by conduction and convection
Stallman's analytical solution (1965) to the subsurface temperature profile in a semi-infinite porous medium in response to a sinusoidal surface temperature provides a test of the CA model's ability to simulate one-dimensional heat convection and conduction in response to a time varying Dirichlet boundary.
Given the temperature variation at the ground surface described by the temperature variation with depth is given by , where A is the amplitude of temperature variation ( • C), T surf is the average surface temperature over a period of τ (s), T ∞ is the initial temperature of the soil column and temperature at infinite depth, and q f is the specific flux through the column top.
The parameters used in analytical examples for Stallman (1965), and the CA code are given in Table 2.The coupled CA code is able to simulate the temperature evolution due to conductive and convective heat transfer as seen from the excellent agreement with the analytical solution (Fig. 5).Lunardini (1985) presented an exact analytical solution for propagation of subfreezing temperatures in a semi-infinite, Table 1.Simulation parameters for heat conduction problems.Analytical solution for this example is given by Eq. ( 13) as per Churchill (1972).

Symbol Parameter
Value initially unfrozen soil column with time t.The soil column is divided into three zones (Fig. 6a) where zone 1 is fully frozen with no unfrozen water; zone 2 is "mushy" with both ice and water; and zone 3 is fully thawed.The Lunardini (1985) solution as described by McKenzie et al. ( 2007) is given by following set of equations:  (Stallman, 1965) and coupled CA model steady state solutions for conductive and convective heat transfer.The soil column in this example is infinitely long, initially at 20 • C, and the upper surface is subjected to a sinusoidal temperature with amplitude of 5 • C and period of 24 h.
where T 1 , T 2 and T 3 are the temperatures at distance x from the temperature boundary for zones 1, 2, and 3 respectively; T 0 , T m , T f , and T s are the temperatures of the initial conditions, the solidus, the liquidus, and the boundary respectively; D 1 and D 3 are the thermal diffusivities for zones 1 and 3, defined as λ 1 /C 1 and λ 3 /C 3 where C 1 and C 3 , and λ 1 and λ 4 are the volumetric bulk-heat capacities (J m −3 • C −1 ) and bulk thermal conductivities (J s −1 m −1 • C −1 ) respectively of the two zones.The thermal diffusivity of zone 2 is assumed to be constant across the transition region, and the thermal where γ d is the dry unit density of soil solids, and ξ = ξ 1 − ξ 3 where ξ 1 and ξ 3 are the ratio of unfrozen water to soil solids in zones 1 and 3 respectively.For a time t in the region from 0 ≤ x ≤ X 1 (t) the temperature is T 1 and X 1 (t) is given by and from X 1 (t) ≤ x ≤ X(t) the temperature is T 2 where X(t) is given by and for x ≥ X(t) the temperature is T 3 .The unknowns, ϑ and γ , are obtained from the solution of the following two simultaneous equations: .
The verification example based on Lunardini's (1985) analytical solution used in this study is the same as that used  Lunardini's (1985) three-zone problem.Equations ( 18), ( 19), and (20) are used to predict temperatures in the completely frozen zone (no phase change and sensible heat only), mushy zone (phase change and latent heat + sensible heat), and unfrozen zone (sensible heat only) respectively.(b) Linear freezing function used to predict unfrozen water contents for two cases used in this study (T m = −1 • C and by McKenzie et al. (2007).Lunardini (1985) assumed the bulk-volumetric heat capacities of the three zones, and thermal conductivities in each zone, to be constant.It was also assumed for the sake of the analytical solution that the unfrozen water varies linearly with temperature.As stated by Lunardini (1985), if unfrozen water varies linearly with temperature then an exact solution may be found for a three-zone problem.Although this will be a poor representation of a real soil system, it will constitute a valuable check for approxi-Table 3. Simulation parameters for predicting the subsurface temperature profile with phase change in a three-zone semi-infinite porous medium.The analytical solution to this one-dimensional problem with sensible and latent heat zones is given by Eqs. ( 18)-( 25) as per Lundardini (1985).

Symbol Parameter
Value mate solution methods.The linear freezing function used in this study is shown in Fig. 6b and the parameters used in Lunardini's analytical solution are given in Table 3.The excellent agreement between the analytical solution and coupled CA model simulations (Fig. 7a, b) for two different cases of T m shows that the model is able to perfectly simulate the process of heat conduction with phase change.
6 Comparison with experimental data Hansson et al. (2004) describe laboratory experiments of Mizoguchi (1990) in which freezing-induced water redistribution in 20 cm long Kanagawa sandy loam columns was observed.The coupled CA code was used to model the experiment as a validation test for simulation of frost-induced water redistribution in unsaturated soils.Four identical cylinders, 8 cm in diameter and 20 cm long, were packed to a bulk density of 1300 kg m −3 resulting into a total porosity of 0.535 m 3 m −3 .The columns were thermally insulated from all sides except the tops and brought to uniform temperature (6.7 • C) and volumetric water content (0.33 m 3 m −3 ).The tops of three cylinders were exposed to a circulating fluid at −6 • C. One cylinder at a time was removed from the freezing apparatus and sliced into 1 cm thick slices after 12, 24, and 50 h.Each slice was oven dried to obtain the total water content (liquid water + ice).The fourth cylinder was used to precisely determine the initial condition.The freezing-induced water redistribution observed in these experiments was simulated using the coupled CA code.Parameters used were sat-  (Lunardini, 1985) and coupled CA model solutions for heat transfer with phase change.Lunardini's (1985) solution is shown and compared with the CA simulation for two cases:  3, Fig. 6).
urated hydraulic conductivity of 3.2 × 10 −6 m s −1 and van Genuchten parameters α = 1.11 m −1 and n = 1.48.The hydraulic conductivity of the cells with ice was reduced using an impedance factor of 2. The thermal conductivity formulation of Campbell (1985) as modified and applied by Hansson et al. (2004) was used.In their simulations of the Mizoguchi (1990) (Mizoguchi, 1990, as cited by Hansson, 2004) and coupled CA model results: (a) 12, (b) 24, and (c) 50 h.
the model using a heat flux boundary at the top and bottom of the columns.The heat flux at the surface and bottom was controlled by heat conductance terms multiplied by the difference between the surface and ambient, and bottom and ambient temperatures respectively.boundary conditions were used in the CA simulations.The value of heat conductance at the surface was allowed to decrease non-linearly as a function of the surface temperature squared using the values reported by Hansson et al. (2004).The heat conductance coefficient of 1.5 J s −1 m −2 • C −1 was used to simulate heat loss through the bottom.Hansson and Lundin (2006) observed that the four soil cores used in the experiment performed by Mizoguchi (1990) were quite similar in terms of saturated hydraulic conductivity, but probably less so in terms of the water-holding properties where more significant differences were to be expected.Such differences in waterholding capacity would result in significant differences in unsaturated hydraulic conductivities of the columns at dif-ferent times during the freezing experiments.The simulated values of total water content agree very well with the experimental values (Fig. 8).The region with a sharp drop in the water content indicates the position of the freezing front.
There is clear freezing-induced water redistribution, which is one of the principal phenomena for freezing porous media and is well represented in the coupled CA simulations.Mizoguchi's experiments have been used by a number of researchers for validation of numerical codes (e.g.Hansson et al., 2004;Painter, 2011;Daanen et al., 2007).The CA simulation shows a comparable or improved simulation for total water content as well as for the sharp transition at the freezing front.

Conclusions
The study provides an example of application of direct solving to simulate highly non-linear processes in variably saturated soils.The modelling used a one-dimensional cellular automata (CA) structure wherein two cellular automata models simulate water and heat flow separately and are coupled through an energy balance module.First-order empirical laws in conjunction with energy and mass conservation principles are shown to be successful in describing the tightly coupled nature of the heat and water transfer.In addition, use of an observed soil freezing curve (SFC) is shown to obliviate the use of non-physical terms such as apparent heat capacity and provide insights into a further subtle mode of coupling.This approach of coupling and use of SFC is easy to understand and follow from a physical point of view and straight forward to implement in a code.The results were successfully verified against analytical solutions of heat flow due to pure conduction, conduction with convection, and conduction with phase change.In addition, freezing-induced water redistribution was successfully verified against experimental data.

Appendix A: Convergence analysis
The CA scheme described in this paper is not unconditionally convergent.Hence, the size of the time step cannot be arbitrarily chosen.In this section we present a detailed evaluation of the convergence criteria of our code to address the choice of the time step.The heat and flow convergence criteria are derived one after another.We start with the heat balance portion.The local energy balance is the basic principle used in our approach.This is imposed by ensuring flux continuity of heat.The local heat balance is described by Eqs. ( 1) and ( 2) and the freezethaw effect.For a 1-D CA application, assuming r = 1, this can be written as where l i is the uniform cell size and λ i,i+1 and λ i,i−1 are the average effective thermal conductivity of the region between the ith, and the i + 1th and i − 1th cells respectively.The second term on the left-hand side of equation is the contribution of freeze-thaw to the thermal energy conservation.
T t i = T t i + e t i is some approximation of the exact solution for temperature T t i at time t and cell index i given an approximation error e t i .Similarly, θ t wi = θ t wi + e t i is an approximation, subject to the discretization error e t i , of the exact solution for the volumetric fraction of water θ t wi .It is useful to rewrite the equation in the following simpler form in order to decouple the thermal and hydraulic processes in terms of known parameters: quantity in parentheses in the second term of the equation can be approximated as ρ w L f dθ w /dT | T =T t i , where dθ w /dT | T =T t i is the slope of the soil freezing curve at T = T t i , a known quantity.Finally, we introduce the term as apparent heat capacity, and rewrite Eq. (A1) as Rearranging the terms, we obtain where λi = λ i,i−1 + λ i,i+1 .Replacing all the error terms by the maximum absolute error term, defined as All coefficients of error terms on the right-hand side of Eq. (A4) are either positive or zero.Given this, the upper bound on the error at time t + t, defined as where f (T , λ) is the term in squared brackets in Eq. ( A4) and F = max {f (T , λ)}.Therefore as long as Eq. ( A5) is satisfied, the error always has an upper bound controlled solely by the discretization error.This is the condition for stability.But, because C i is a function of time, an adaptive time stepping scheme would be well suited to solve the problem.The adaptive time stepper would need to satisfy Eq. (A5) at each time step.As long as the thermal energy balance component of our CA algorithm obeys the time stepping-spatial discretization relationship in Eq. (A5) it remains stable.For such time step control, using the Lax-Richtmeyer equivalence theorem, one only needs to show that the thermal module represents a consistent numerical approximation to the full diffusion equation (including C i to account for the freeze-thaw effect) in order to prove convergence of our method.To do this we note the following recurrence relations: It is worth noting that here we have assumed a constant value of F through all time steps.In the following we argue that this does not affect the generality of the convergence analysis that follows next.Clearly, if the only source of error in our approximate solution is the discretization of a continuous process, then our initial values must be error-free; i.e.E 0 = 0. Therefore, For lim t, l i → 0, we have the cluster of terms within the square brackets converge to the expression As T is an exact solution of the above diffusion-equation form, we must have the terms within square bracket converge to 0 as lim t, l i → 0. This argument for the boundedness of F as t, l i → 0 holds at each time step and, hence, would have led to the same conclusion if we would have used a time variable maximum value of f (T , λ) in Eq. (A8).Therefore, in general, at any time t This formally shows that our numerical algorithm, with time stepping satisfying Eq. (A5), is consistent and hence follows the convergence of the thermal module.
We can construct a similar convergence analysis for the hydraulic module.But we will approach this problem from the continuum version of the modified Richard's equation for variably saturated flow for the sake of brevity.The modified Richard's equation for variably saturated flow can be written as (in the absence of a source term) The left-hand side of Eq. (A12) follows from the continuum version of the first term on the left-hand side of Eq. ( 5) where Eq. ( 6) has been used to eliminate .The term on the righthand side is the Darcy flux, introduced as the continuum version of Eq. ( 10) where the total head H = ψ + z.The effect of freeze-thaw on the total head can be accounted for as a Clausis-Clapeyron process as given in Eq. ( 12).To make this clear, we rewrite Eq. (A12) as follows: We can make use of the following relations to eliminate temperature from Eq. (A13): where C w = ∂θ w /∂H | H =H (t) is the local slope of the soil retention curve which can be derived from Eq. ( 7).Equation (A16) now has the same form as the expression in Eq. (A10).It is immediately clear that, if one would have followed the full formal arguments as outlined for the thermal module, the condition for stability of the variably saturated flow dynamics part of our algorithm is of the form where C w = C w (1−ρ ice /ρ w ) and ki = k i,i−1 +k i,i+1 .We refer the reader to Mendicino et al. (2006) for a thorough convergence analysis of the water flow problem.The formal argument is exactly equivalent to that presented by us for the heat flow problem.Combining Eqs.(A5) and (A17), the following condition gives stability and convergence conditions for the overall CA problem: )

Figure 1 .
Figure 1.One-dimensional cellular automata grids based on the von Neumann neighbourhood concept.How many neighbours (grey cells) interact with an active cell (black) is controlled by the indicial radius (r).

Figure 2 .
Figure2.Flow chart describing the algorithm driving the coupled CA code.Subscripts TC, HC and FT refer to changes in physical quantities due to thermal conduction, hydraulic conduction and freeze-thaw processes respectively.Hydraulic conduction and thermal conduction are two different CA codes coupled through updating of volumetric heat capacity and the freeze-thaw module to simulate the simultaneous heat and water movement in soils.Corresponding equations or sections containing module description are shown in red text in squared brackets.

Figure 3 .
Figure 3. Graphical description of the phase change approach used in this study.The curve is a soil freezing curve for a hypothetical soil.The change in water content (dθ w ) due to Q res,i is used to determine T new by integrating along the SFC (Eq.11).

Figure 4 .
Figure 4. Comparison between the analytical solution given by Churchill (1972) and coupled cellular automata model simulation for a perfectly thermally insulated 4 m long soil column.Lines represent the analytical solution and symbols represent the CA solution for time points as shown in the legend.The initial temperature distribution is shown on the right.

Figure 5 .
Figure 5.Comparison between the analytical(Stallman, 1965) and coupled CA model steady state solutions for conductive and convective heat transfer.The soil column in this example is infinitely long, initially at 20 • C, and the upper surface is subjected to a sinusoidal temperature with amplitude of 5 • C and period of 24 h.

Figure 6 7 Figure 6 .
Figure 6 (a) Diagram showing the setting of Lunardini (1985) three zone problem.Equations 2 18, 19, and 20 are used to predict temperatures in completely frozen zone (no phase change 3 and sensible heat only), mushy zone (phase change and latent heat + sensible heat), and 4 unfrozen zone (sensible heat only) respectively.(b) Linear freezing function used to predict 5 unfrozen water contents for two cases used in this study (Tm = -1°C and Tm = -4°C).6 7

Figure 7 .
Figure 7.Comparison between analytical solution of heat flow with phase change(Lunardini, 1985) and coupled CA model solutions for heat transfer with phase change.Lunardini's (1985) solution is shown and compared with the CA simulation for two cases: (a) T m = −1 • C and (b) T m = −4 • C (Table3, Fig.6).
Figure 7.Comparison between analytical solution of heat flow with phase change(Lunardini, 1985) and coupled CA model solutions for heat transfer with phase change.Lunardini's (1985) solution is shown and compared with the CA simulation for two cases: (a) T m = −1 • C and (b) T m = −4 • C (Table3, Fig.6).
t 2 ⊗ . . .⊗ S t i ⊗ . . .⊗ S t N at time t.Let a local transition rule φ be defined on a neighbourhood of spatial indi-

Table 2 .
Stallman (1965)meters for predicting the subsurface temperature profile in a semi-infinite porous medium in response to a sinusoidal surface temperature.The analytical solution to this one-dimensional heat convection and conduction problem in response to a time varying Dirichlet boundary is given by Eqs.(14)-(17) as perStallman (1965).