Comparison of algorithms and parameterisations for infiltration into organic-covered permafrost soils

Infiltration into frozen and unfrozen soils is critical in hydrology, controlling active layer soil water dynamics and influencing runoff. Few Land Surface Models (LSMs) and Hydrological Models (HMs) have been developed, adapted or tested for frozen conditions and permafrost soils. Considering the vast geographical area influenced by freeze/thaw processes and permafrost, and the rapid environmental change observed worldwide in these regions, a need exists to improve models to better represent their hydrology. In this study, various infiltration algorithms and parameterisation methods, which are commonly employed in current LSMs and HMs were tested against detailed measurements at three sites in Canada’s discontinuous permafrost region with organic soil depths ranging from 0.02 to 3 m. Field data from two consecutive years were used to calibrate and evaluate the infiltration algorithms and parameterisations. Important conclusions include: (1) the single most important factor that controls the infiltration at permafrost sites is ground thaw depth, (2) differences among the simulated infiltration by different algorithms and parameterisations were only found when the ground was frozen or during the initial fast thawing stages, but not after ground thaw reaches a critical depth of 15 to 30 cm, (3) despite similarities in simulated total infiltration after ground thaw reaches the critical depth, the choice of algorithm influenced the distribution of water among the soil layers, and (4) the ice impedance factor for hydraulic conductivity, which is commonly used in LSMs and HMs, may not be necessary once the water potential driven frozen Correspondence to: Y. Zhang (yinsuozhang@carleton.ca) soil parameterisation is employed. Results from this work provide guidelines that can be directly implemented in LSMs and HMs to improve their application in organic covered permafrost soils.


Introduction
Infiltration of snowmelt or rain into frozen ground or the unfrozen active layer is a critical hydrological process in permafrost regions (Woo, 1986) and its simulation is a key component in almost all process-based Land Surface Models (LSMs) (e.g. Bonan, 1991;Verseghy, 1991;Gusev and Nasonova, 2003;Dai et al., 2003) and Hydrological Models (HMs) (e.g. Beven and Kirkby, 1979;Liang et al., 1994;Yang and Niu, 2003;Pomeroy et al., 2007;Peckham, 2008). However, mathematically quantifying infiltration has always been a challenge (Smith et al., 2002), due mainly to the heterogeneity of most natural soils and highly dynamic changes of soil water status and hydraulic properties during infiltration. Those difficulties become extreme in permafrost environments due to ground thawing/freezing processes and a surface organic layer that frequently mantles permafrost terrain (Kane and Stein, 1983;Kane and Chacho, 1990;Slater et al., 1998). Soil hydraulic properties change rapidly during infiltration or abruptly within infiltration depth between: (1) frozen and unfrozen states (Burt and Williams, 1976;Kane and Stein, 730 Y. Zhang et al.: Cold region infiltration algorithms Kane, 1979;Quinton et al., 2005), which may cause convergence problems for numerical infiltration schemes (e.g.  or violate the assumptions for many analytical infiltration schemes (e.g. Green and Ampt, 1911). Other complicating factors for infiltration in permafrost terrain include macropore-induced preferential flow (Mackay, 1983), ineffective pore spaces in organic soils , and hysteresis effects during thawing and freezing (Horiguchi and Miller, 1980). Mathematical representation of infiltration into permafrost soils is poorly developed compared to those in non-permafrost soils (Kane and Chacho, 1990;Luo et al., 2003).
Most early LSMs and HMs do not have an explicit frozen soil scheme (Luo et al., 2003;Zhang et al., 2008). The influence of frozen soil on infiltration and runoff is typically treated with a few simple assumptions such as: (1) liquid soil moisture remains at zero or a small constant value once the soil temperature passes below 0 • C (e.g. Verseghy., 1991;Dai et al., 2003), and (2) hydraulic conductivity becomes zero once frozen (e.g. Bonan, 1991;Verseghy, 1991;Dai et al., 2003). Recent improvements of frozen soil processes in LSMs and HMs include: (1) soil ice content is explicitly represented as a diagnostic variable (e.g. Cherkauer et al., 2003;Niu and Yang, 2006;Nicolsky et al., 2007), (2) thawing/freezing depth is recognized as a controlling factor for infiltration/runoff and is dynamically simulated with improved algorithms, parameterisations and model configurations (e.g. Slater et al., 1998;Kuchment et al., 2000;Yi et al., 2006;Zhang et al., 2008), (3) variable unfrozen water content is parameterized using relationships with subfreezing soil temperature (e.g. Li and Koike, 2003;Zhang et al., 2008), and (4) frozen soil infiltration is allowed based upon soil ice content or subfreezing soil temperature (Niu and Yang, 2006;Nicolsky et al., 2007;Pomeroy et al., 2007). In some cold region hydrological models (e.g. Flerchinger and Saxton, 1989;Zhang et al., 2000;Pomeroy et al., 2007), infiltration schemes for frozen soil have been explicitly designed. However, their algorithms vary widely from first order empirical estimation (e.g. Gray et al., 1985) to complex numerical solutions of the simultaneously coupled thermal and moisture transfer equations with phase changes (e.g. Tao and Gray, 1994;. Some models even provide multiple options for infiltration simulations during different infiltration stages or different site conditions (e.g. Pomeroy et al., 2007;Peckham, 2008). Testing and comparison of infiltration schemes were only found for mineral soil conditions (e.g. Slater et al., 1998;Cherkauer and Lettenmaier, 1999;Cherkauer et al., 2003;Niu et al., 2005), and many of them only dealt with homogeneous soils (e.g. Flerchinger et al., 1988;Boike et al., 1998;Mishra et al., 2003;Chahinian et al., 2005). The validation of infiltration simulations in organic covered permafrost soils is extremely scarce due to the limited quantity and quality of field data in such regions.
In this study, we present a comprehensive review of infiltration algorithms and parameterisations and evaluate their applicability for organic-covered permafrost soils. Selected algorithms and parameterisations are evaluated using field data obtained from three organic-covered sites in Canada's discontinuous permafrost region. The overall objective is to provide guidelines for the implementation of appropriate infiltration algorithms/parameterisations in LSMs and HMs to improve their performance in permafrost regions.

Infiltration algorithms
Infiltration of surface water is controlled by many factors, including soil depth and its texture profile, soil hydraulic properties and water status, water supply intensity and patterns, infiltration time, and thawing/freezing depth. Efforts have been made to numerically solve the water transfer equation or its coupled form with the heat transfer equation for non-uniform unfrozen soil infiltration (e.g. Celia et al., 1990;Ross, 1990;Šimůnek et al., 2005), or uniform frozen soil infiltration (e.g. Harlan, 1973;Guymon and Luthin, 1974;Tao and Gray, 1994;Hansson et al., 2004), but to the best of our knowledge, no successful application exists for infiltration problems in non-uniform soil with thawing/freezing process involved. Moreover, since most of the infiltration events in cold environments involve large volumes of water flux in a short period, extremely fine resolutions in temporal (seconds or less) and in spatial (centimeters or less) domains are required to achieve stable numerical solutions (Jame and Norum, 1980;Tao and Gray, 1994;Smith et al., 2002), imposing considerable computational expense. Consequently, operational LSMs and HMs rarely utilize numerical schemes for infiltration. Typically, infiltration is separately calculated using conceptual, empirical or analytical methods and added as a source term to the numerical scheme, which is used to calculate the heat transfer and water redistribution within the vadose zone (e.g. SHAW, Flerchinger and Saxton, 1989;CLASS, Verseghy, 1991;CLM3.5, Oleson et al., 2008). Considering this, numerical infiltration schemes are excluded from this study. Table 1 lists equations of infiltration algorithms and parameterisations referenced in this study. Table 2 summarises most infiltration schemes in current LSMs and HMs that are applicable for soils involving thawing/freezing. The conceptual models are typically developed for extreme soil conditions. For example, examining snowmelt infiltration in frozen prairie soils, Granger et al. (1984) and Gray et al. (1985) grouped infiltration patterns into three broad categories: (1) restricted; for soils with impermeable surface layers such as ice lenses, (2) unlimited; for soils with a high percentage of air-filled macropores, and (3) limited; for soils   Bonan, 1991;Verseghy, 1991;Dai et al., 2003); CHRM (Gray et al., 1985;Pomeroy et al., 2007) Unlimited infiltration All surface water infiltrate (b) CHRM (Gray et al., 1985;Pomeroy et al., 2007) Instantaneous infiltration All surface water instantaneously percolate to the supra-permafrost water table (b), (d), and (e) ARHYTHM (Zhang et al., 2000) Gray's empirical infiltration Snow melt infiltration, Eq.
in between the first two categories. For the limited condition, Eq.
(1) was proposed to estimate infiltration amount for the entire snowmelt season.  developed a semi-empirical parametric infiltration scheme (Eq. 2) for uniform frozen soil, which relates the time dependant infiltration to surface and initial saturation ratios, saturated hydraulic conductivity and soil temperature. Another category of semi-empirical schemes is termed "distributed", which re-lates probability distributions of infiltration capacity to a certain topographic index or spatial distribution of soil moisture (e.g. Beven and Kirkby, 1979;Liang et al., 1994;Niu et al., 2005). Distributed infiltration algorithms are mainly developed for regional or global applications (e.g. TOPMODEL, Beven and Kirkby, 1979;VIC, Cherkauer et al., 2003;VISA, Yang and Niu, 2003), thus excluded from this point-based model comparison study.
Analytical algorithms are exact solutions of the water transfer equation (e.g. Richards' equation) under specific soil conditions and water supply patterns. Despite their limiting assumptions (Table 2), analytical algorithms are the most frequently employed algorithms in LSMs and HMs, due to their solid physical base and ability to obtain parameters through field measurements, texture associations (Clapp and Hornberger, 1978) or pedo-transfer functions (Wösten, 1999;Wagner et al., 2001). The most frequently cited Green and Ampt (1911) algorithm (Eq. 4) has the following assumptions: (1) uniform soil extending to half infinite plane, (2) uniform antecedent water content, (3) constant head ponding at the surface, and (4) a piston-like sharp wetting front. Among the numerous efforts (e.g. Bouwer, 1969;Smith et al., 2002;Chu and Marino, 2005;Talbot and Ogden, 2008) in relaxing these assumptions, the twostage Mein and Larson (1973) infiltration scheme simulates both the pre-ponding and ponded infiltration of steady rain into uniform soil (Eq. 5), while the two-stage Smith and Parlange (1978) scheme allows for variable rainfall rates (Eq. 6). Flerchinger et al. (1988) modified Green and Ampt algorithm to simulate infiltration into layered non-uniform soil and employed it in the Simultaneous Heat and Water (SHAW) model (Flerchinger and Saxton, 1989).

Essential parameters for simulating infiltration
Practically all of the parameters in Eqs. (1-6) can be determined from basic soil hydraulic properties (i.e. θ s , K s , ψ 0 ) and characteristics (i.e. relationships among water potential, water content and hydraulic conductivity). In the context of LSMs and HMs, basic hydraulic properties and characteristics are typically associated with texture classes (e.g. Clapp and Hornberger, 1978;Letts et al., 2000) or other soil attributes such as bulk density or grain-size fraction (e.g. Wösten, 1999;Wagner et al., 2001). For frozen soils, parameterisations of unfrozen water content and ice impedance to hydraulic conductivity are also crucial to infiltration simulation (Kane and Stein, 1983;Kane and Chacho, 1990;Slater et al., 1998;Quinton et al., 2008).

Soil hydraulic properties
Field and laboratory measurements of θ 0 , θ s , ψ 0 and K s for permafrost soil are possible (e.g. Burt and Williams, 1976;Dingman, 2002;McCauley et al., 2002;Hayashi and Quinton, 2004;Carey et al., 2007;Quinton et al., 2008), yet not feasible for LSM and HM applications. The most frequently employed texture associations for hydraulic properties in LSMs and HMs were those compiled by Clapp and Hornberger (1978) for 11 classes of mineral soils and by Letts et al. (2000) for 3 classes of organic soils. Both compilations contain only unfrozen soil samples from nonpermafrost regions. For comparison, Table 3 presents a summary of hydraulic properties obtained from permafrost soils and those values from Clapp and Hornberger (1978) and Letts et al. (2000). To our knowledge, no values of K s for frozen organic soil have been reported. LSMs and HMs typically apply the same value to θ 0 and θ s (e.g. Verseghy, 1991;Niu and Yang, 2006). While this could be true for many mineral soils, smaller θ s vs. θ 0 values were frequently observed for organic soils in permafrost sites Quinton et al., 2008), due mainly to the dead-end or self-closed pores (Hoag and Price, 1997). During soil freezing, effective pore space is lowered due to the presence of ice, which blocks pores and therefore reduces both water storage capacity and conductivity. Frozen soil K s exhibits extremely large temperature dependence in the small temperature range just below freezing (Table 3). For example, for a small temperature increase from −0.26 • C to 0 • C, Burt and Williams (1976) observed K s increase almost 8 orders of magnitude for a fine sand. K s of organic soil shows a strong dependence on its state of decomposition, or more apparently, on soil depth (Z). Quinton et al. (2008) developed a simple relationship between K s and Z (Eq. 7) based on field measurements at three organic covered permafrost sites in Canada. ψ 0 is normally estimated from the soil characteristic curve (Clapp and Hornberger, 1978;van Genuchten, 1980). Most reported ψ 0 for organic soil is close to −0.01 m (Letts et al., 2000), which is much higher (closer to zero) than those for mineral soils (Table 3).
Soil type (frozen temperature)    (17)(18)(19) are widely utilized in LSMs and HMs (e.g. Flerchinger and Saxton, 1989;Hannson et al. 2004;Niu and Yang, 2006;Soulis and Seglenieks, 2008). Zhang et al. (2008) evaluated three types of unfrozen water formulations against field measurements at four permafrost sites, of which three had organic cover, and showed that all three methods could represent the field measurements reasonably well if appropriate parameters were chosen. Among the three, a water potential-freezing point depression equation (Eq. 20) (Cary and Mayland, 1972), was frequently chosen by models with coupled thermal and hydrological simulations (Flerchinger and Saxton, 1989;Cherkauer and Lettenmaier, 1999;Koren et al., 1999;Niu and Yang, 2006). The freezing-point (T f ) normally has a value of 0 • C, but could be slightly below zero for many clayey soils and some organic soils (Koopmans and Miller, 1966;Osterkamp, 1987;Quinton et al., 2005).

Site descriptions
All three field sites are located in Canada's discontinuous permafrost regions above 60 • N latitude (Table 4). The Scotty Creek peat plateau site (referred to as SC P hereafter) is located in a wetland-dominated region near Fort Simpson, Northwest Territories. The two other sites: a boreal forest site (WC F hereafter) and an alpine tundra site (WC A here-after), are located within the Wolf Creek Research Basin, Yukon Territory. SC P is a peat plateau that rises 0.9 m above a surrounding wetland, and is underlain by permafrost with an active layer ∼0.7 m deep. Vegetation is predominantly open canopy black spruce (Picea mariana) mixed with some northern shrubs and lichen and moss on the forest floor. WC F site has closed-canopy white spruce (picea glauca) mixed with other spruce, pine and poplar species. The understory consists of a wide range of shrub species with feather moss and grasses at the surface. WC A is situated on a windswept ridge with sparse vegetation of mosses, lichens, grasses and occasional patches of short shrubs. All three sites had various organic cover depths over mineral horizons (Table 4). The organic soil at SC P has two distinct layers with the upper 0.1-0.15 m consisting of living plants and lightly decomposed organic materials overlying peat in a more advanced state of decomposition. Below the peat is clay to silt-clay soil with very low permeability (Hayashi et al., 2007 (Carey and Woo, 2001;Wright et al., 2008).

Field measurements and water balance components
Field measurement periods extending from late March or early April to end of August were chosen to conduct model tests at all three sites. Infiltration scenarios include snowmelt infiltration into frozen or thawing ground, and rainfall infiltration into the thawed active layer. Two seasons from each site, i.e. 2004 and 2005 for SC P, and 1998 and 1999 for WC F and WC A, were selected based on the availability of field data. Data from 1998 and 2004 were used for calibration of unknown parameters and initial conditions, while data from 1999 and 2005 were used for model validation. Details of field measurements and the methods to quantify the water balance components can be found in Hayashi et al. (2007) and Wright et al. (2008) for Scotty Creek site, and in Pomeroy and Granger (1999) and Janowicz (2000) for Wolf Creek sites.

Snowmelt (M sn )
Daily M sn was calculated from the difference in successive daily values of snow water equivalent (SWE). This approach assumes that sublimation and evaporation from the melting snow are negligible. SWE at SC P was directly measured at 5 m intervals along a 41 m transect, while SWE at WC F and WC A were determined by daily snow depth measurements via ultrasonic depth sensors (Campbell UDG01) and snow density sampled at variable times at 25 m intervals along a 625 m transect. Average SWE values along the transects were used in this study.

Soil temperature (T ) and thaw depth (Z T )
Soil temperatures at various depths were continuously recorded at all three sites.  Zhang et al. (2008). Ground surface temperature (T s ) at SC P was directly measured by thermistor under snow cover and by infrared sensor once snowfree. T s at WC F and WC A is estimated from air temperature during snowfree period and from temperatures measured at 0.1 m above and 0.05 m below the ground surface during snowcover period. Wright et al. (2008) calibrated the coefficient C 6 in Priestley-Taylor ET Eq. (21) for three different land-cover types at SC P site with lysimeter measured ET data. The average C 6 value for the peat plateau ground surface ranged from 0.68-0.91. Here, a value of 0.82 gave the best soil moisture simulations during the calibration period (2004), thus adopted for model testing period (2005). At Wolf Creek (WC F and WC P), an ET estimation method (Eq. 22), developed by Granger and Gray (1989), was implemented by , and thus adopted in this study.

Soil water content
Daily liquid soil water contents (θ l ) were measured at all three sites using site-calibrated TDR or water content reflectometer (CS-615) probes throughout the study period. The At WC F and WC A, changes of total soil water content (ice + liquid water, θ T ) were monitored during snowmelt seasons using twin-probe gamma attenuation techniques as described by Gray and Granger (1986). The maximum monitoring depths were 1.2 m at WC F and 0.8 m at WC A. Absolute values of θ T were estimated from the changes with respect to a reference measurement the previous fall or subsequent summer when the active layer is thawed. θ T was not regularly measured at SC P site, however, values for θ T prior to snow melt were estimated from two θ T profile measurements, one by TDR in the fall of 2002 immediately before the freeze-up and another by two 0.7 m deep frozen peat cores sampled near the study pit on 6 April 2003 (Hayashi et al., 2007), using the procedure described in Wright et al. (2008).

Infiltration and runoff
The cumulative infiltration (CINF i ) of snowmelt (M sn ) and/or rainfall (R) is estimated from other measured water balance components as listed in Eqs. (23-24). Equation (23) is used for SC P and Eq. (24) is used for WC F and WC A. Equation (23) assumes that the freezing of infiltrated liquid water during the melt period is negligible, thus the melted soil water (M sw ) could be estimated from changing Z T between two time steps. Runoff (CROF) is estimated by the difference between water input (R + M sn ) and infiltration (Eq. 25).

Fig. 1.
Measured soil water content-pressure head relationships and best fitting curves with three parameterisation methods for several organic soils at a peat plateau of Scotty Creek and a north-facing slope of Grange Creek.

Soil properties and hydraulic parameters
One or more soil pits were excavated at each site to determine soil texture profiles . Additional soil cores were taken to determine basic physical properties such as fractions of sand, silt, clay and organic, bulk density (ρ b ) and total porosity (θ 0 ) Quinton et al., 2008). Field measurements of saturated hydraulic conductivity, soil water retention curves and unsaturated hydraulic conductivity curves were conducted for several organic soils at Scotty Creek and Wolf Creek watersheds (Carey and Woo, 2001;Hayashi and Quinton, 2004;Quinton et al., 2005;Carey et al., 2007;Quinton et al., 2008). Parameters for BC-Para, VG-Para and CH-Para were derived by curve-fitting to the measured data ( Figs. 1 and 2). Table 4 lists soil texture profiles and corresponding parameters at the sites. Table 5 lists the infiltration algorithms and parameterisations evaluated in this study. The algorithms/parameterisations are chosen as: (1) used currently in LSMs and HMs; (2) designed (or modified) to be used in organic covered permafrost soils, and (3) required parameters and inputs can be achieved/calibrated from the observations and data described above.  (12) and (14) with parameters from Fig. 1b.

Infiltration algorithms, parameterisations and evaluation methods
The Simultaneous Heat and Water (SHAW) Model (Flerchinger and Saxton, 1989;Flerchinger, 2000) was selected as a common platform to host most of the algorithms and parameterisations listed in Table 5 to allow reasonable comparisons. Infiltration by semi-empirical algorithms (GRAY-IN and ZHAO-IN) is calculated independently. SHAW was developed to simulate heat, water and solute transfers for soils experiencing freezing and thawing (Flerchinger and Saxton, 1989). The original SHAW model consists of full balances of energy, water and solutes within a one-dimensional profile including layers of plant canopy, snow, residue and soil (Flerchinger, 2000). Its modularized coding structure makes it effective to add, to disable or to modify individual processes or parameterisations. In this study, to focus on infiltration and reduce uncertainties, most processes that do not directly influence infiltration, such as canopy process, snow process, surface energy balance and solute transfer were disabled, and only the soil thermal and moisture transfers including thawing/freezing and infiltration/runoff are simulated. This reduces the vertical profile to only the organic and mineral soil layers. The coupled soil temperature and moisture (ice, liquid and vapor) transfer equations (Eqs. 26 and 27) were iteratively solved with a finite difference scheme. The numerical scheme quantifies the depth of soil thawing/freezing and the soil moisture redistribution among soil layers by non-infiltration processes. The thawing/freezing process is simulated by an apparent heat capacity parameterization, which has been proved effective in permafrost regions by Zhang et al. (2008). The upper boundary conditions were supplied by ground surface temperature (T s ) and evapotranspiration (ET) described above. Snowmelt and rainfall water were supplied to ground surface, but their infiltration into the soil profile is simulated by a separate module, i.e. the modified Green-Ampt scheme (GA-SHAW) for layered soils (Flerchinger et al., 1988). Zero water and heat fluxes are assumed at the lower boundary (5 m soil depth), which is adequate for short-term simulations in this study (Zhang et al., 2008). The original SHAW code uses CH-Para for soil hydraulic parameterisation and LN-Ice as ice impedance factor (Table 5). ML-CLASS is an infiltration module taken from version 3.4 of Canadian Land Surface Scheme (Verseghy, 2009), which uses Mein and Larson (1973) as infiltration algorithm, CH-Para as soil hydraulic parameterisation, SQ-Ice as ice impedance factor (Soulis and Seglenieks, 2008). IT-TOPO is coded as another optional infiltration algorithm following the principles described in Zhang et al. (2000). Parameterisations listed in Table 5, but not in the original SHAW model, were also coded into SHAW as alternative modules for comparison. The modified SHAW model is fed with daily surface forcing variables including T s , R, Msn and ET. Much smaller and dynamic time step (half-hourly or smaller) is used internally to ensure convergence of the numerical scheme. A 16-layer soil vertical resolution is used at all three sites. The layer depths are 0.05 m for top two layers, 0.1 m for 0.1-0.8 m depth and progressively increasing for deeper layers until the simulated soil bottom at 5 m. Site specific parameters (e.g. coefficients to quantify T s and ET from known meteorological variables) and unknown initial conditions (e.g. initial soil temperatures and moisture profiles below the observation depth) were achieved by fitting the simulated diagnostic variables to their observed values during calibration years (1998 and 2004). The principal diagnostic variables are thawing depth, cumulative infiltration/runoff and soil liquid water content corresponding to the measurement depths. Since multiple unknowns and multiple diagnostic variables are involved, an iterative procedure is performed until all the diagnostic variables achieve their best fitting results. Similar procedure has been performed in evaluating the thawing/freezing simulations in Zhang et al. (2008). The calibrated parameters and conditions at the end of calibration periods are then used to quantify the required inputs and initiation conditions during the model validations years (1999 and 2005). A common set of inputs and initial conditions are used for all the model validation runs with different infiltration algorithms to ensure valid comparison of the algorithms.

Soil hydraulic parameterisations
The most important soil hydraulic parameterisations for infiltration/redistribution are the water retention curve (water potential vs. water content) and hydraulic conductivity curve (hydraulic conductivity vs. water potential or content). All three commonly used methods in Table 5 are able to fit observed soil water retention curves in moderate soil moisture ranges for several organic soils (Fig. 1). Upon approaching saturation, both CH-Para and BC-Para calculate values above θ s , which have to be capped by θ s . When liquid water declines under frozen conditions with very low water potential, CH-Para gives liquid water content values below θ r . To counter this, many LSMs and HMs assume a minimum value for liquid water content (e.g. CLASS, Verseghy, 1991;SHAW, Flerchinger and Saxton, 1989). The discontinuity of CH-Para and BC-Para for saturated or extremely dry (frozen) conditions may result in numerical convergence problems for the moisture transfer equations (e.g. Eq. 27). Alternative treatments such as the water balance method, or explicit solutions must be used when soil moisture approaches θ s or θ r (Flerchinger, 2000). Theoretically, VG-Para is more suitable for numerical water transfer models due to its smoothness over the entire soil moisture range. However, its current ap-plication in operational LSMs and HMs is limited due to poor parameter availability of many soil types. Figure 2 shows unsaturated hydraulic conductivity K values observed by Carey et al. (2007) and simulated by the three parameterisation methods using an estimated K s value of 2.5×10 −4 m s −1 and parameters in Fig. 1b. All three methods gave similar K values in normal pressure head ranges, except under saturated conditions when pressure head reaches zero. In this case, K values calculated by CH-Para and BC-Para have to be capped by K s . Although only in a small pressure head range, observed K values generally match calculated values.

Parameterisation of unfrozen water content
In this study, unfrozen water content is calculated by a water potential-freezing point depression equation (Eq. 20), combined with the reversed form of one of the three water retention equations (Eq. 9,or 11,or 13). Values for ψ 0 are taken from Letts et al. (2000) for organic soils and from Clapp and Hornberger (1978) for mineral soils, based on corresponding texture class. Similar to the performance in water retention simulations (Fig. 1), the discrepancies in simulated unfrozen water curves by the three water retention equations were mostly found when liquid water content reaches maximum or minimum values due to soil temperature changes (Fig. 3). Those errors could be easily corrected by bounding the calculated unfrozen water content with observed θ r and θ s . Although from two different data sets, parameters obtained in Figs. 1 and 3 are similar for the same soil. For example, the parameters in Fig. 1c worked equally well in Fig. 3a. Since unfrozen water and soil temperature are easier to measure in permafrost soil than the soil water potential, these datasets can be an effective alternative to derive soil hydraulic parameters from traditional pressure head measurements, as demonstrated by Spaans and Baker (1996) and Flerchinger et al. (2006). Most of the soil hydraulic parameters in Table 4 are derived by this method.

Reduction of hydraulic conductivity (K) due to soil freezing
During soil freezing, two effects act to reduce K. First, the reduction of liquid water content will lower the water potential ψ (Eq. 9, or 11, or 13), reducing K in a similar manner as soil drying (Eq. 10, or 12, or 14). Second, an impedance factor due to the presence of ice is applied to K (e.g. Eq. 17,or 18,or 19). Figure 4 illustrates the changes of ice impedance factors (Fig. 4a) and hydraulic conductivity with increasing soil ice fraction (Fig. 4b). Soil parameters used are the same as in Fig. 3d and K s is set for a typical loam as outlined in Clapp and Hornberger (1978). The drying effect alone reduces K to similar orders of magnitude as those reported for frozen soils in Table 3 (Fig. 4b). Although further reduction by an impedance factor such as Eq. (17) with a C 5 value as high as 10 is noticeable, restrictions imposed by impedance equations currently employed in LSMs and HMs are relatively small compared to the effect of decreased water potential (Eq. 10). K-CLASS in Fig. 4b is the K parameterisation of frozen soil in ML-CLASS; it applies f imp,2 (Eq. 17) to Eq. (12) and uses effective pore space (θ s − θ i ) instead of θ s . Using effective pore space makes θ l /(θ s − θ i ) in Eq. (12) always equal to 1.0 when total soil water content is θ s , thus the only reduction to K s is f imp,2 . This treatment underestimates the magnitude of K reduction shown in Table 3.

Sensitivity tests
To reduce the redundancy in evaluation procedure caused by the numerous possible combinations of algorithms and parameterisations, sensitivity tests of algorithms, parameterisa-tions and some key input variables were performed for infiltration simulation. Daily infiltration outputs were cumulated over three ground thawing stages: i.e. frozen, thawing and thawed (Table 6). Frozen and thawing stages are separated by the date at which ground thawing starts, while thawing and thawed stages are separated by the date at which ground thaw reached a prescribed depth (0.4 m at SC P and WC F, and 0.15 for WC A). A baseline run with specified configurations was first conducted with data from each of the three sites (Table 6). Test runs with only one changed attribute (algorithm/parameterisation/input) were followed and results summarized in Table 6. For comparison and analysis, some observed water inputs and infiltration, and snowmelt infiltration calculated by the semi-empirical methods were also included in Table 6.  The sensitivity of the simulated infiltration to the changes of algorithms, parameterisations or inputs only occurred during soil frozen and thawing stages at all sites (Table 6). Once the ground thawed to a certain depth, all water infiltrated into the soil regardless the configured algorithms, parameterisations or inputs. The three soil hydraulic parameterisations (CH-Para, BC-Para and VG-Para) had little influence on the simulations during all three stages. In this study,  (Table 4) were deliberately chosen to represent the hydraulic curves (Figs. 1-3). For CH-Para and BC-Para, the soil water content values were capped by θ s and θ r during saturated or extremely dry or frozen conditions. The various ice impedance factors gave marginal differences for the simulated frozen soil infiltration in most of the tested cases, even disabling this factor (f imp ≡ 1) did not show obvious increase in frozen soil infiltration at all three sites. This further confirmed the result that once Eq. (20) was employed, ice impedance factors such as Eqs. (16-18) may not be necessary as an extra restriction on K, as also indicated by Fig. 4. Infiltration simulation during ground frozen and thawing stages is very sensitive to ground surface temperature, mainly through its influence on thawing development. Therefore it is cru-cial for LSMs and HMs to obtain or to simulate accurate T s in order to better simulate infiltration/runoff during soil frozen or thawing stages. Changes of ET did not influence the infiltration amount of any stage at any site, but did have a large influence on soil water content (results not shown). The insensitivity of ET on infiltration is due to the fact that the surface organic soil has a large conductivity and water holding capacity once thawed, which normally exceeds the water supplying rate from the surface, and thus does not respond to the soil water content changes caused by ET.

Comparison of infiltration algorithms
Based on the sensitivity tests, one fixed set of parameterisations (i.e. CH-Para for soil water potential and hydraulic conductivity, Eq. (20) and Eq. (11) for unfrozen water content and Eq. (18) for ice impedance), was adopted during the comparison of the three analytical algorithms (GA-SHAW, ML-CLASS and IT-TOPO). Except for its infiltration algorithm, another difference between ML-CLASS and the other two algorithms was ML-CLASS uses effective pore space (θ s − θ i ) instead of θ s for frozen soil while others retain θ s in all conditions. The simulation statistics in the validation years compare well to those achieved during the calibration years, indicating the robustness of the parameters (Table 7).
Detailed comparison of the three analytical algorithms at the three test sites during validation years are presented in Figs. 5, 6 and 7. Frozen soil infiltration was observed at both SC P and WC A, but not at WC F, where very little snowmelt water was supplied to the frozen ground (Table 6, Fig. 6a). IT-TOPO does not allow frozen soil infiltration, while both GA-SHAW and ML-CLASS algorithms simulated infiltration into frozen ground. The frozen soil infiltration simulated by ML-CLASS was approximately double the amount simulated by GA-SHAW at both SC P and WC A sites, due mainly to the effective pore space applied in ML-CLASS, which increased the hydraulic conductivity in Eq. (12) (Fig. 4). Although ML-CLASS did provide a closer value to the observation-based estimation of the total frozen soil infiltration at SC P (Table 6), GA-SHAW gave better cumulative patterns at all three sites (Figs. 5c, 6c and 7b). Liquid water content during frozen stages was mainly controlled by soil temperature, and no differences were found among the runs with different algorithms despite the large differences in simulated infiltration.
No notable difference was found for simulated ground thawing among the three infiltration algorithms. The observed ground thawing patterns were well simulated at SC P and WC A, but underestimated at WC F during the early stages. At WC F, the observed thawing started only 4 days after snowmelt began, and thawed to 0.15 m depth in 5 days, while 0.1 m of snow still remained on the ground surface. As the ground surface temperature directly beneath the snowcover can not be above zero, the model is incapable of predicting this thawing. A potential cause of this thaw is the advection of heat and water from snow-free areas nearby which cannot be accounted for in a 1-D model. The differences between the simulated infiltration by GA-SHAW and ML CLASS are small during thawing stages at all three sites and are well comparable to observed values (Table 6 and Figs. 5c, 6c and 7b). No observation was available for WC A during thawing, but most of the input water infiltrated into ground as simulated by GA-SHAW and ML-CLASS. IT-TOPO gave the smallest infiltration volumes during thawing stages as it does not allow water to infiltrate into frozen ground. The liquid soil water content during the thawing stage is controlled both by infiltration water and the liquid water released as ice melts. The simulations by GA-SHAW and ML-CLASS were similar at all three sites during thawing stages, with better results at SC P and WC A than at WC F. The weak soil water simulation at WC F was caused by the poor estimation of thawing development. Simulated liquid soil water content by IT-TOPO was lower than the other two algorithms at WC F and WC A, but similar at SC P as simulated infiltration at this site was close to observed values.
The rapid and complete infiltration of surface water during thawed stages simulated by all three algorithms was achieved by the large hydraulic conductivity and water holding capacity of the surface organic soil. While the same infiltration amounts were simulated during thawed stages (Table 6), In general, during the thawed stage, GA-SHAW performed best compared with the observed liquid water at all sites, with simulations at SC P and WC A better than at WC F ( Table 7). The poor simulation at WC F may be attributed to the delayed soil thawing that altered soil water status in the early stages of thaw (Fig. 6f, g and h). In Table 6, GRAY-IN gave the total amount of snowmelt water infiltration, which included all the infiltration during frozen stage and some of the infiltration during thawing. The value calculated by GRAY-IN for WC F is close to the observation, but at SC P infiltration is poorly underestimated. Although no comparable observation was available at WC A, the snowmelt infiltration calculated by GRAY-IN was much smaller than the total infiltration during frozen and thawing stages calculated by GA-SHAW and ML-CLASS. ZHAO-IN overestimated the infiltration during frozen stage at SC P, but gave comparable results at WC F and WC A with both observation and the results from GA-SHAW and ML CLASS. These results suggest that the parameters of empirical and parametric algorithms are highly site dependent and must be calibrated accordingly.

Conclusions
Based on various field measurements at three discontinuous permafrost sites, this study evaluated five infiltration algorithms, three soil hydraulic property parameterisations, various ice impedance schemes on frozen soil hydraulic conductivity, and their influences and sensitivity on water infiltration into organic covered permafrost soils. The following conclusions are presented to provide guidelines to improve the infiltration schemes and parameterisations of current LSMs and HMs, for their applications in permafrost environments. Some limitations of this study are also addressed.
1. The single most important factor controlling infiltration into permafrost soils is ground thaw status. The infiltration during the soil frozen stage is largely controlled by soil ice content, while thaw depth controls the infiltration during the thawing stage. Once the ground thawed to certain depth (i.e. 15-30 cm in this study), infiltration became "unlimited" as described in Gray et al. (1985).
2. The performance of the semi-empirical infiltration algorithms (GRAY-IN and ZHAO-IN) varied among the three sites, indicating that they require site-specific parameter calibration, limiting their applications in HMs and LSMs.
3. The conceptual instantaneous infiltration algorithm (IT-TOPO) restricts infiltration during the frozen stage and underestimates infiltration during thawing stage. Even though the simulated infiltration during the thawed stage is the same as other algorithms, its water redistribution scheme typically underestimates soil water content in the upper thawed layers. Consequently, IT-TOPO is not recommended for applications in organiccovered permafrost soils.
4. The two analytical algorithms (GA-SHAW and ML-CLASS) gave similar infiltration simulations during most stages and site conditions. Some differences exhibited during the ground frozen and thawing stages were caused by their different parameterisation of frozen soil hydraulic conductivity (K), rather than the infiltration algorithms used: i.e. modified Green-Ampt and Mein-Larson for layered soil. This study recommends both algorithms for infiltration simulation at organic cover permafrost sites, but the parameterisation of frozen ground K in current CLASS frozen soil module (Verseghy, 2009), particularly the introduction of an effective pore space (Soulis and Seglenieks, 2008), may overestimate the frozen soil K.
5. With properly chosen parameters, the three soil hydraulic property parameterisations, i.e. Clapp and Hornberger (CH-Para), Brooks and Corey (BC-Para), and van Genuchten (VG-Para), achieve similar soil water retention and hydraulic conductivity curves for normal soil moisture ranges. However, the calculated soil liquid water content and hydraulic conductivity by CH-Para and BC-Para should be bounded by maximum and minimum values during saturation or frozen conditions. This treatment could cause convergence problems for infiltration schemes coupled with numerical moisture redistribution schemes. While VG-Para gives smooth parameterisation curves for all soil moisture range, its application is restricted by the general availability of parameters. Paired measurement data of unfrozen water content (θ l ) and subfreezing soil temperature (T s ) could be used to derive soil hydraulic parameters by fitting θ l − T s relationships derived from Eq. (20) and the soil water retention equations.
6. Only by applying Eq. (20) to unsaturated hydraulic conductivity (K) equations (Eq. 10, 12 or 14) could realistic simulation of the reduction of K due to soil frozen be achieved in the range of observations. Further reduction by various ice impedance factors as employed in many land surface and hydrological models may not be necessary in organic-covered permafrost soils.
7. Sensitivity tests indicate that simulated infiltration is sensitive to algorithms, parameterisations and input changes only during the frozen and fast thawing stages. The infiltration at these stages is sensitive to ground surface temperature but not to evapotranspiration.
8. All three sites in this study are located in relative flat areas, reducing the potential for lateral flow. However, the abnormal ground thaw under snowcover at WC F site (Fig. 6b) was most likely caused by the advection of heat and/or water from snow-free patches. Slopes are common feature in permafrost terrain and cannot be omitted in operational LSMs or HMs for infiltration and runoff simulations.
9. The preferential flow suggested by many permafrost infiltration studies (e.g. Mackay, 1983) can be represented by high hydraulic conductivity values parameterised for surface organic layers. No additional algorithm is necessary to account for preferential flow at these sites.