Dual‐permeability modeling of preferential flow and snowmelt partitioning in frozen soils

The infiltrability of frozen soils modulates the partitioning of snowmelt between infiltration and runoff in cold regions. Preferential flow in macropores may enhance infiltration, but flow dynamics in frozen soil are complicated by soil heat transfer processes. We developed a dual‐permeability model that considers the interacting effects of freeze–thaw and preferential flow on infiltration and runoff generation in structured soils. This formulation was incorporated into the fully integrated groundwater–surface water model HydroGeoSphere, to represent water–ice phase change in macropores such that porewater freezing is governed by macropore–matrix heat exchange. Model performance was evaluated against laboratory experiments and synthetic test cases designed to examine the effects of preferential flow on snowmelt partitioning between infiltration, runoff, and drainage. Simulations were able to reproduce experimental observations of rapid infiltration and drainage behavior due to macropores very well, and approximated soil thaw to an acceptable degree. Simulation of measured data highlighted the importance of macropore hydraulic conductivity, as well as macropore–matrix heat and water transfer, on controlling preferential flow dynamics. Test cases replicated a range of snowmelt partitioning behavior commonly observed in frozen soils, including subsurface conditions that produce rapid infiltration and deeper drainage, the contrast between limited vs. unlimited infiltration responses to snowmelt, and the temporal evolution of runoff generation. This study demonstrates the important influence that water freezing along preferential flowpaths can have on infiltrability and runoff characteristics in frozen soils and provides a physically based description of this mechanism that links infiltration behavior to hydraulic and thermal properties of structured soils.

These interacting effects of preferential flow and soil freeze-thaw have a strong influence on hydrologic fluxes that partition snowmelt throughout a landscape, such as runoff generation, soil moisture distribution, and groundwater recharge (Espeby, 1992;Gray et al., 2001;Mohammed et al., 2019). Infiltrated meltwater that bypasses frozen soil layers also has the potential to facilitate leaching of labile nutrients and contaminants to groundwater during snowmelt (Grant et al., 2019;Holten et al., 2018), as there are limited opportunities for water and solute exchange with the frozen matrix . Similarly, as preferential flow reduces the dilution, sorption, and degradation effects of the soil matrix on groundwater quality, preferential flow in frozen soils can influence the performance of soil cover systems and infiltration-bioretention water treatment facilities in cold regions (Callaghan et al., 2017;Ding et al., 2019;Huang et al., 2018).
Physically based models that assume diffuse equilibrium flow have been widely applied to frozen soil environments (Flerchinger & Saxton, 1989;Hansson et al., 2004;Johnsson & Lundin, 1991) but, due to a lack of representation of preferential flow, are often incapable of capturing the timing and magnitude of snowmelt driven processes in macroporous soils. Dual-continuum formulations that explicitly define both macropore and matrix characteristics can describe nonequilibrium flow and transport in porous media (Šimůnek, Jarvis, van Genuchten, & Gardenas, 2003). As infiltrability and hydraulic conductivity are influenced by pore-ice content, the effects of macropore-matrix heat transfer on porewater freeze-thaw need to be considered if dual-permeability models are to be adapted for use in frozen environments (Jarvis et al., 2016;Mohammed et al., 2018).

Core Ideas
• Preferential flow in frozen soils is an important process governing snowmelt partitioning. • A dual-permeability numerical model is developed for preferential flow in frozen soils. • Simulations demonstrate the interacting effects of preferential flow and soil freeze-thaw.
Physically based representations of preferential flow in hydrological models have the potential to help understand and evaluate how these complex thermal-hydraulic processes can affect snowmelt-driven hydrologic behavior under different land use and climate conditions (Ireson et al., 2013). Accordingly, the objectives of this study were to develop and evaluate a dual-permeability description for variably saturated flow in frozen soils. The formulation was incorporated into the fully integrated surface water-groundwater model Hydro-GeoSphere (HGS; Aquanty, 2019), to include heat transfer relations describing the freezing of water within macropores Mohammed et al., 2018). Incorporating this new dual-permeability formulation into a fully integrated hydrological model allows for explicit consideration of surface runoff partitioning, as well as water flow in multiple spatial dimensions. Model performance was assessed with measured data from frozen-soil infiltration experiments. One-dimensional (1D) synthetic test cases were also performed to examine the effects of preferential flow and macropore-matrix heat transfer on runoff generation and vadose zone fluxes in frozen soils.

Conceptual model
In most unsaturated soils, macropores lack the capillarity to retain much water and will hence generally remain air filled after soil freezing, providing preferential pathways that can enhance snowmelt infiltration into and below the frost zone (Espeby, 1992;Flerchinger et al., 2005;Ishikawa et al., 2006;Mohammed et al., 2019;Stähli et al., 1996;van der Kamp et al., 2003;Watanabe & Kugisaki, 2017). Under these conditions, preferential flow dominates and most infiltration occurs through macropores, as pore-ice greatly reduces the hydraulic conductivity of the soil matrix at freezing temperatures (Demand et al., 2019;Granger et al., 1984;Holten et al., 2018;Huang et al., 2018). However, depending on the soil thermal regime, infiltrated water may be cooled by the surrounding frozen soil matrix and freeze within F I G U R E 1 (a) Conceptual model of interacting matrix and macropore water and heat transport in frozen soil by Mohammed et al. (2018); and (b) simplified conceptual model of matrix and macropore water and heat transport used for numerical model development in this study. For the simplified model, although latent heat effects (solid blue circles) are not explicitly included in the soil matrix, it is indirectly considered through the heat capacity (blue dashed circle). It should be noted that arrows work both in vertical and lateral directions within domains in all cases preferential pathways, thereby blocking flow (Demand et al., 2019;Stähli et al., 1996;Watanabe & Kugisaki, 2017). Conceptualizing preferential flow as a dual-continuum process accounts for this nonequilibrium of water pressures and temperatures between the bulk soil (matrix) and macropores during infiltration (Köhne et al., 2009;Šimůnek et al., 2003). Dual-domain models partition the soil into two interacting pore domains (macropore and matrix), which are parameterized with different flow and transport properties (Jarvis et al., 1991). Macropore-matrix interface parameters that describe mass and heat exchange between domains are also needed (Gerke & van Genuchten, 1993b). Mohammed et al. (2018) developed a dual-permeability framework for coupled water and heat transport in frozen soils ( Figure 1a). The water-retention characteristics of the macropore domain allow it to remain air filled under unsaturated conditions, enabling infiltration into frozen soil and imbibition into the matrix. The model considers effects of macroporematrix heat exchange on the freezing of infiltrating water in macropores. The conceptualization of heat transfer in the soil matrix (Figure 1a) includes (a) heat conduction due to temperature gradients, (b) advection due to flowing water, and (c) latent heat transfer due to water-ice phase change. In macropores, only advection and latent heat transfer are considered, whereas conduction is neglected due to the large thermal mass and volume of the matrix relative to macropores. Heat exchange across the macropore-matrix interface allows infiltrating water in macropores to be cooled by the surrounding matrix and freeze depending on matrix thermal conditions, potentially blocking macropores with ice. Using this framework, Larsbo et al. (2019) developed the first dualpermeability description of water and heat transport in frozen soils. The description was implemented into the MACRO model (Larsbo et al., 2005), a physically based 1D numerical model for simulating vadose zone flow and transport. However, their model was not compared against measured infiltration data, nor does it account for the possibility of higher spatial dimensions or surface flow partitioning.
It is computationally challenging to solve the Richards equation in frozen soils, due to sharp gradients at both wetting and freezing fronts (Hansson et al., 2004). In dualpermeability models, these fronts occur both within and across domains . In this study, several simplifying assumptions were used to make the model tractable, while still including relevant heat transfer processes that control the system response. Thermal conduction was assumed to be the primary heat transfer mechanism governing porewater freeze-thaw in both domains (Figure 1b). This assumes that conduction dominates in the soil matrix , that melting snow at the ground surface is isothermal (0˚C) so advective transport is negligible (Iwata et al., 2010), and that heat loss to the bulk matrix is the major mechanism governing freezing of water in macropores (Stähli et al., 1996;Watanabe & Kugisaki, 2017). Latent heat of freezing-thawing in the soil matrix is not explicitly considered, but taken into account indirectly by using an "apparent" heat capacity (limitations discussed later). The warming effect of interdomain heat exchange on the matrix due to infiltration in the macropore domain is ignored, as it is less relevant to flow dynamics compared to porewater freezing in macropores. In effect, the temperature of the bulk soil matrix and properties of the macropore-matrix interface determine whether preferential flow occurs or is reduced due to freezing of infiltrated water and pore blockage.

Numerical model
HydroGeoSphere is a three-dimensional numerical model capable of simulating coupled surface and subsurface water flow and solute or energy transport processes (Aquanty, 2019; Therrien et al., 2010). Here, we describe the existing soil matrix freeze-thaw processes in HGS (Schilling et al., 2019), along with the new dual-permeability freeze-thaw relations developed in this study. Although simulations presented here are 1D, model capabilities are extended to higher spatial dimensions (Mohammed, 2019). A full description of the model capabilities can be found in Aquanty (2019), and a list of symbols is provided in Supplemental Table S1. Subsurface flow in the dual-permeability formulation of HGS is described by the Richards equation in both the soil matrix and macropores (Gerke & van Genuchten, 1993a). In the following equations, terms describing the soil matrix domain are denoted with a subscript m, and macropore domain with subscript f. The total water content, θ t (m 3 m −3 ), in a volume of soil is given as where θ f (m 3 m −3 ) is the water content in the macropore domain, θ m (m 3 m −3 ) is the water content in the matrix domain, and w f (unitless) is the ratio of the volume of the macropore domain relative to the total soil-pore volume.
Combining the Richards equation with a fluid mass balance, assuming ice is immobile, with relative permeability as a function of total water (θ w ) and ice content (θ ice ) in each domain, yields where ρ w (= 1,000 kg m −3 ) and ρ ice (= 920 kg m −3 ) are the density of water and ice, respectively. K s (m s −1 ) is the saturated hydraulic conductivity, k r (unitless) is the rel-ative hydraulic conductivity of each domain, and ψ (m) and z (m) are the matric potential and elevation heads, respectively. Γ w (s −1 ) is the fluid transfer rate between subsurface domains. S w (unitless) is the effective liquid saturation (i.e., the volumetric liquid soil moisture normalized to porosity, ε [m 3 m −3 ] of each domain). Water retention and K s relations for both domains are described using the van Genuchten-Mualem model (van Genuchten, 1980). Relative hydraulic conductivity reduction due to pore-ice blockage is achieved by reducing the liquid water saturation by the ice content: where S w is the effective liquid water saturation, S w+ice is the total water saturation, S ice is the ice saturation, and T m (˚C) is the temperature of the soil matrix. The partitioning of water between ice and liquid phases is determined by temperature (which is a spatially and temporally variable) such that (Andersland & Ladanyi, 2004): where T freeze (0.1˚C) and T melt (−0.1˚C) are the freezing and melting temperatures, respectively, and β (= 0.42) is a dimensionless empirical fitting parameter (Rempel, 2008). This formulation of phase change in the soil matrix is not modeled using the energy balance, and matrix temperature alone determines the partitioning between ice and water (i.e., when water is added to a model cell, the temperature of the water is in instantaneous equilibrium with the cell temperature [limitations discussed below]). Fluid exchange between domains (Γ w ) is proportional to the matric potential difference between domains (using a dual-node approach): where K as (m s −1 ) and k r (unitless) are the saturated hydraulic conductivity and relatively hydraulic conductivity of the domain interface, respectively. The interface relative hydraulic conductivity is reduced by the ice content in the soil matrix. The factor has a minimum value of 10 −6 so that a small amount of unfrozen water always remains present at the interface, to minimize numerical instability. The macroporematrix interface exchange coefficient, α wd (m −2 ), is defined as (Gerke & van Genuchten, 1993b): where β d is a dimensionless geometrical shape factor (assumed 3 for a planar geometry), d (m) is the diffusion pathlength over which the water and heat exchange occurs (Gerke & van Genuchten, 1996), and γ w is a dimensionless empirical coefficient (0.4). Analytical solutions are used to calculate heat conduction and freeze-thaw in the soil matrix, thus eliminating the need to solve the fully coupled mass and energy balance and reducing computational demands when applying the model at larger scales and coupling to surface flow (Schilling et al., 2019). Temperature distribution in the soil matrix is calculated by (Schilling et al., 2019): where T m (˚C) is the temperature of the matrix, and T b (˚C) is the background temperature used as the lower boundary temperature at depth where temperature is assumed constant. Soil thermal properties are assumed to be homogeneous and constant throughout the simulation (despite changing water and ice contents), with κ m (m 2 s −1 ) as the thermal diffusivity defined as where λ m (W m −1˚C−1 ) and c m (J kg −1˚C−1 ) are the soil bulk thermal conductivity and apparent heat capacity, respectively. Equation 7 is solved subject to the following boundary conditions: where T s (˚C) is the ground surface temperature. It should be noted that this formulation of heat transport in the soil matrix does not account for latent heat effects on soil temperature, which can cause discrepancies between measured and simulated temperatures during soil freezing and thawing periods (McKenzie et al., 2007). Instead, the effects of latent heat are indirectly considered by using an "apparent" or "effective" heat capacity value to represent the effects of latent heat on reducing soil thermal diffusivity (Goodrich, 1978).

Description of porewater freeze-thaw in macropores
This section details the newly implemented freeze-thaw processes in the macropore domain, based on the conceptualization of heat transfer in Figure 1b, which does include latent heat exchange. Assuming the freezing temperature in macropores is largely unaffected by strong capillary forces due to their relatively large pore size, water freezes at 0˚C. Ice formation within the macropore domain is governed by energy exchange between domains, driven by thermal conditions in the soil matrix Mohammed et al., 2018). The rate of ice formation is determined by the difference in soil matrix temperature (T m ) and infiltrating water temperature (T in ), and a specified heat transfer coefficient defining the geometric and thermal characteristics of the macroporematrix interface: where ε f is the porosity of the macropore domain (i.e., w f ε t ), and L i ( = 335,000 J kg −1 ) is the enthalpy of fusion for water. λ m (W m −1˚C−1 ) is the thermal conductivity of the domain interface, which is assumed equal to the thermal conductivity of the matrix domain, and (α wd ) is the exchange coefficient that is defined based on the geometry of heat transfer between domains (Equation 6). If T in is assumed to be 0˚C (Stähli et al., 1996), Equation 10 reduces to Thawing of water in the macropore domain is assumed to occur coincident with matrix thaw at corresponding depths (since, once frozen, the temperature in the macropore is in equilibrium with the matrix), which implies that the freezing process described above is unidirectional. For numerical stability, any point in the macropore domain was considered to be blocked by ice and impermeable when S ice = 0.95.

Simulation of laboratory infiltration experiments
To compare the model against measured data, we simulated a subset of experiments presented in Pittman et al. (2020). Intact soil columns (50-cm length × 30-cm diam.) were frozen to −8˚C before being subjected to ponded infiltration and soil thawing conditions. The resulting infiltration and drainage were measured along with changes in soil moisture and temperature throughout the column. Details on the experimental design and results can be found in Pittman et al. (2020). Three columns had interconnected macropore networks, which produced rapid infiltration (average = 1.1 × 10 4 mm d −1 ) and drainage (average = 4.3×10 3 mm d −1 ) while the soil was still frozen. Since responses of the three columns were similar, model parameters were calibrated to the average infiltration, Macropore 0-0.5 0.9 0 10 2.0 6.0 × 10 −3 0.042 drainage, and soil thaw rates. The first 7 d of the experiments were simulated, until the soil thawed to a depth of approximately 20 cm. It should be noted that highly calibrated model parameters that would achieve an exact fit to the observed data were not the goal of these simulations. Rather, the focus was to evaluate the capability of the dual-permeability formulation to represent preferential flow and freeze-thaw processes via the model's ability to reproduce the general observations of water movement and soil thaw of the column experiments. A 1D vertical domain was used consisting of two layers representing the topsoil (0−0.25 m) and subsoil (0.25−0.5 m) ( Table 1). Measured hydraulic conductivity and soil moisture characteristic curves were used to estimate matrix hydraulic parameters (Table 1). Macroporosity (ε f = w f ε t ) of each layer was estimated from soil moisture characteristic curves as follows (Reynolds et al., 2009): pores larger than 0.1 mm were assumed to be drained at 100 mm of tension, and the volumetric water content at this tension was subtracted from the measured total porosity to estimate the macroporosity. Macropore van Genuchten-Mualem parameters could not be measured in the same fashion as matrix parameters and were estimated from literature values (Dohnal et al., 2012;Gerke & van Genuchten, 1993a), whereas K sf was estimated using Poiseuille's equation (Rosenbom et al., 2009). The most abundant macropore structures were generally < 1 mm in aperture (Pittman et al., 2020), and K sf was estimated assuming an average cylindrical macropore radius (r c ) of 0.1 mm, which resulted in K sf = 6.3 × 10 5 mm d −1 (7.3 × 10 −3 m s −1 ). K sf and K sm values were adjusted to match observations. The diffusion pathlength, d, was set to 50 mm (0.05 m) (Equation 6, Table 2), considered to be representative of a soil with a high degree of macropore-matrix exchange . Thermal diffusivity, κ m , was used as a calibration parameter to match the soil thaw rate.
A specified fluid flux was applied to the top of the domain to simulate ponding of water. As in the experiments, no runoff was allowed and all water remained ponded until it infiltrated. Unit gradient boundary conditions were applied to the base of both pore domains. Initial soil moisture conditions were estimated from pre-freezing soil moisture measurements. Recorded temperatures at the top of the column during experiments were applied as the top temperature boundary. The model was initialized with a temperature of −8˚C at the soil surface, which decreased linearly to a temperature of −6˚C, which was applied to the base of column as the lower boundary for the duration of the simulation. This was done due to the limitation of the analytical scheme for matrix temperature in Equation 7, in which the lower boundary temperature (Equation 9b) is fixed, when in the experiment it actually increased from approximately −8 to −6˚C, so using a lower boundary temperature that matched the average temperature at the base of column (−6˚C) during the experiment allowed for a better match to soil temperatures near the base of the column during the simulation.

Test Case 1: Episodic water input into frozen soil
The remaining test cases (Cases 1−3) are synthetic and intended to test the model response to different hydrological scenarios. All test cases are simulated with the same 1D vertical soil domain (Table 1). These cases are similar to test cases presented in Larsbo et al. (2019), but with a larger domain and greater amount of water input, in order to illustrate the effects of preferential flow on infiltration in addition to runoff generation and deeper drainage in frozen soil.
Test Case 1 simulates repeated water input events in frozen soil with initially air-filled macropores in order to examine the effect of ice buildup and blockage within macropores on infiltration and runoff generation. For this scenario, no matrix water flow or macropore-matrix water transfer was allowed, such that all infiltration occurs through the macropore domain and the freezing of infiltrating water in macropores is governed only by heat transfer between domains (Equation 11). A 50-cm-long soil column was simulated with a 1D vertical domain. Model domain and macropore hydraulic and macropore-matrix transfer parameters used are listed in Tables 1 and 2. The soil matrix temperature was initialized and held constant at −2˚C, by applying a temperature of −2˚C to the upper and lower temperature boundary conditions for the entire simulation. A unit gradient boundary condition was applied to base of the macropore domain. To simulate episodic water input every 30 min, a specified flux was applied at the top of the column for 6 minutes at a rate of 1,000 mm d −1 (4.17 mm total per event) and repeated for five input events. Any excess water input unable to infiltrate was partitioned to runoff.

Test Case 2: Constant water input into frozen soil
This scenario was aimed at illustrating the effect of freezing of infiltrating water on the temporal evolution of water partitioning between infiltration, runoff, and drainage. Initial and lower boundary conditions and soil hydraulic and thermal parameters in this simulation are the same as Case 1 (i.e., a constantly frozen matrix with no macropore-matrix water exchange). The simulation considered one water input event occurred for 10 h at a rate of 100 mm d −1 equaling a total input of ∼41.7 mm.

Test Case 3: Thawing of soil during constant water input event
For Test Case 3, we simulated the thawing of an initially frozen soil profile during a constant water input of 50 mm d −1 for 10 d. The model domain was initialized at −5˚C and the temperature at the top boundary was set to 5˚C for the simulation period. The domain and hydraulic parameters were the same as previous test cases, but in this case, two different heat transfer coefficients, representing low (50 mm) and high (300 mm) values of diffusion pathlengths, were used to simulate high (d = 50 mm) and low degrees (d = 300 mm) of heat transfer between macropore and matrix domains (Tables 1  and 2).
Simulated infiltration and drainage rates were controlled by the hydraulic conductivity of the macropore domain (K sf ) and macropore-matrix interface (thus, matrix K sm ). Macropore flow accounted for ∼99% of cumulative infiltration (Figure 3), which is consistent with the rapid drainage response and flow calculations from experiments (Pittman et al., 2020), as well as tracer experiments that show connected macropore networks can channel infiltrated water through frozen soil to unfrozen ground (Demand et al., 2019;Grant et al., 2019;Holten et al., 2018). Energy exchange between the soil matrix and flowing water in macropores resulted in a very small amount of the infiltrating water freezing in the macropore domain. With ponding at the soil surface, flow in the macropore domain was fast enough, relative to the freezing rate, that pore-ice formation within the macropore domain was limited. This resulted in rapid drainage, with roughly half of the water bypassing the matrix and the other half retained in the upper half (i.e., topsoil water content increasing from 0.15 to 0.38) of the matrix (Figure 2c). Increases in matrix soil moisture were mostly due to infiltrating water in the macropore domain being imbibed into the matrix via macroporematrix fluid transfer (Figure 3a), with only a small amount of vertical infiltration into the frozen matrix (Figure 3b). Higher matrix permeability and air-filled porosity of the topsoil (Table 1) enabled more macropore-matrix fluid transfer compared with the subsoil where lower matrix permeability and higher relative ice saturation limited fluid exchange. This mechanism of lateral fluid transfer during infiltration into frozen soil agrees with dye-staining experiments that found that most infiltration occurred through macropores and that dye penetration into the matrix occurred along macropore walls (Cey & Rudolph, 2009; Demand et al., 2019; Lepore F I G U R E 2 Measured and simulated cumulative (a) infiltration and (b) drainage responses in laboratory experiment (C1, C4, and C6 represent the soil column ID numbers); (c) simulated matrix total water saturations during infiltration; and (d) temperature profile showing depth of soil matrix thaw at different simulated (sim) and observed (obs) times in days (d) F I G U R E 3 (a) Simulated infiltration, drainage, and macroporematrix fluid transfer rates for laboratory experiment test case. Note the different scale of y axis for water input located on the right; (b) zoomedin view of matrix infiltration Stadler et al., 2000). Imbibition along macropore walls influences both the magnitude and depth of infiltration and depends on the matrix air-filled porosity and permeability of the macropore-matrix interface (Callaghan et al., 2017;Cey et al., 2006;Weiler, 2005).
The model was able to reasonably approximate the rate of soil thaw (Figure 2d) given the inability to explicitly represent latent heat effects in the bulk soil, as well as the inability to replicate the transience of the lower boundary temperature (Section 2.2). The thermal diffusivity (λ m /c m = 7.0×10 −8 m 2 s −1 ) in Equation 7 needed to simulate soil-thaw was lower than the range for clay-rich soils (1.0−2.5 × 10 −7 m 2 s −1 ) (Jury & Horton, 2004), which is expected when using the "apparent" heat capacity as an effective parameter value to mimic latent heat effects. Latent heat reduces the apparent heat capacity (c m ) according to Harlan (1973): is the bulk heat capacity of soil matrix. Despite the limitations of not explicitly incorporating latent heat or using the transient description given above, the model was still able to reproduce the average patterns of soil thaw. The shallow soil (10-cm depth) thawed in the first day (Figure 2d). The rate of thaw then decreased at greater depths with thawing at 20 cm occurring at 7 d, which compared well with the measured thaw depth, but with an overestimation of near-surface temperatures (Figure 2d). The ability to reproduce soil thaw at the profile scale reasonably justifies the use of simplified matrix heat transport for the purposes used here, where reproducing average rates of soil thaw and complex hydraulic behavior is more important than detailed descriptions of heat transport in the matrix, and for the flexibility it allows when upscaling to higher spatial dimensions and coupling to surface flow. However, in certain environments, these assumptions will not be acceptable and transient heat transfer relations explicitly considering latent heat effects and advection by flowing water in the soil matrix would be needed to adequately assess ground thaw mechanisms (Kane et al., 2001;Kurylyk & Watanabe, 2013). A significant, but computationally demanding, improvement to the heat transfer and freeze-thaw relations here would be to include fully coupled water and thermal transport processes in the bulk soil matrix (Larsbo et al., 2019). Overall, results show that using measured soil matrix properties (porosity, soil moisture characteristic, and K s ) and estimated macropore (macroporosity and macropore aperture) properties is effective for simulating preferential flow in frozen soils with the dual-permeability approach presented. Given that infiltration and drainage occurred so quickly in this experiment and in the absence of more detailed experimental datasets, the following synthetic test cases were designed to illustrate the effect of the freezing of water within macropores (Equation 11) on preferential flow and the partitioning of water between runoff, infiltration, and deep drainage in frozen soil.

Episodic water input in frozen soil (Test Case 1)
Test Case 1 considered infiltration of water into an airfilled macropore domain surrounded by a frozen soil matrix. This test case demonstrates the importance of infiltrationrefreezing behavior on the temporal evolution of runoff generation in frozen ground, as well as the transition of the soil profile from unlimited to limited to restricted infiltrability as described by Granger et al. (1984). During the first and second input events, all water infiltrates and drainage occurs at the base of the domain (Figure 4), showing that unlimited infiltration can occur in frozen soils with air-filled macropores. This response is similar to observations in the previous soil column experiments (Section 4.1). During the first input event, ice begins to form immediately throughout the column when water infiltrates into the macropore domain, but ice formation is insufficient to limit infiltrability below the rate of water input (Figures 5 and 6). However, with consecutive input events, a larger portion of the macropore domain becomes blocked with ice ( Figure 6), reducing infiltrability to below the rate of water input, resulting in less infiltration and more runoff in the fourth event (Figure 4). By the beginning of F I G U R E 5 Liquid and ice saturations at depths of 1-, 10-, and 30cm depths below the soil surface in the macropore domain simulated in Test Case 1 the fifth (last) input event, macropores are almost completely blocked by ice in the upper 5 cm of the column ( Figure 5 and 6), which restricts further infiltration and drainage, and all water input is partitioned to runoff (Figure 4).
These results are largely consistent with the current perception of rapid preferential flow and pore-ice formation in macropores (Holten et al., 2018;Watanabe & Kugisaki, 2017) and show that the model provides a quantitative description of the dynamics of water freezing in initially air-filled macropores, which other studies have shown are important to reduce infiltrability and generate runoff in frozen soils Stadler et al., 1997). The trend of decreasing infiltration and deep drainage over successive events as macroporosity became blocked with ice is similar to the numerical experiments of Larsbo et al. (2019) with similar boundary conditions. Likewise, the trends match the laboratory study of Grant et al. (2019) where infiltration and drainage though frozen soil columns decreased with successive irrigation events. Simulations here additionally considered the partitioning of water input into surface runoff. An analogous field situation to this scenario would be a winter period with overwinter snowmelt events . The simulated results generally agree with field studies that show freezing of infiltrated meltwater in macroporous soils over multiple midwinter events causes a progressive reduction in infiltrability and a sequential increase in runoff Thunholm et al., 1989). With consecutive input events under frozen conditions, a larger portion of the macropore Vadose Zone Journal F I G U R E 6 Simulated vertical profiles of macropore water and ice saturation for Test Case 1 showing infiltration into macropore domain with freezing at various times during the simulation period before and after water input events ranging from (a) initially air-filled and ice-free macropores at 0.001 d to (i) final state of macropores blocked with ice at 0.09 d network becomes blocked with ice ( Figures 5 and 6), and more runoff is generated in later events (Figure 4).

Constant water input in frozen soil (Test Case 2)
In Test Case 2 with constant water input, simulated infiltration is initially unlimited and begins concurrently with water input at the same rate of 100 mm d −1 (Figure 7). Drainage begins at 0.08 d and peaks at a maximum rate of 40 mm d −1 at approximately 0.1 d (Figure 7). Runoff begins at 0.12 d when ice saturation at 1-cm depth reached 0.45, reducing the infiltrability below the input rate (Figures 8 and 9). Macropores become completely blocked with ice (S ice = 0.95) at approximately 0.13 d (Figure 8), after which all water input is partitioned into runoff (Figure 7).
The dynamics of this simulation are in line with studies that argue that freezing of infiltrated water and pore-ice blockage close to the ground surface (near the water source) are needed F I G U R E 7 Input (black), infiltration (blue), runoff (purple), and drainage (green) rates simulated in Test Case 2 to reduce infiltrability before runoff (or ponding) can occur in frozen soils (Baker & Spaans, 1997;Nyberg et al., 2001;Stadler et al., 2000). The results also agree with field studies that show snowmelt in dry frozen soils with an initially high infiltrability can still generate runoff if the soil remains cold enough to cause freezing of infiltrated water and poreice blockage near the ground surface (Appels et al.

Response to water input during soil thaw (Test Case 3)
Test Case 3 examines the sensitivity of simulation results to the degree of macropore-matrix thermal energy exchange (Equation 11). Since heat transfer properties are relatively consistent in most soils, the degree of heat exchange is controlled by the geometry of the macropores and the macropore network within the soil matrix, specifically the diffusion pathlength (Equation 6), which has the dominant influence on the rate of interdomain heat transfer . At the beginning of the simulation, macropores are open and water immediately infiltrates in both the high (Case 3a, d = 50 mm) and low (Case 3b, d = 300 mm) macroporematrix heat transfer scenarios ( Figure 10). As water in the macropore domain moves downward to depths where the soil matrix is still frozen, it begins to freeze at that depth. Subsequently, as the thawing front in the matrix migrates downward, water in the macropore domain that temporarily froze subsequently melts and vertical flow continues. This repeated freeze-thaw and percolation process occurs sequentially with depth as the soil thaws and appears in Figure 11 as a shortterm increase then decrease in ice saturation at successive depths from 5 to 50 cm.
In Case 3a with a high degree of macropore-matrix heat transfer, the infiltration rate initially equals the water input rate as water infiltrates into air-filled macropores, but infiltrability is quickly reduced once macropore ice saturation increases. When ice completely blocks the macropore domain at 20-and 30-cm depth after 0.2 d (Figure 11a), infiltration and drainage cease and all input is partitioned to runoff ( Figure 10a). As the matrix subsequently thaws from the top (Figure 11a), the ice in the macropore domain begins to melt at the depths corresponding to the soil matrix thawing front. Note that there is a brief period of time between 4.8 and 5 d where the macropore domain is not fully blocked with ice at any depth (Figure 11a), which caused a brief increase in infiltration and drainage and correspondingly brief decrease in runoff (Figure 10a). Once the soil at the base of the column begins to thaw at 6.8 d, infiltration and drainage increase until all ice has melted. When ice no longer completely blocks the macropore domain at any depth (Figure 11a), input water again infiltrates and drains (Figure 10a). Water stored in the macropore domain above previously frozen layers can now flow, causing infiltration and drainage to briefly increase to a peak value of almost 80 mm d −1 , before drainage stabilizes at the water input rate (Figure 10a).
In Case 3b with low macropore-matrix heat transfer, ice saturation increases much more gradually and the macropores never get fully blocked by ice (Figure 11b). Therefore, infiltration and drainage equal the water input rate throughout the simulation (Figure 10b).
These results highlight that the degree of macroporematrix heat transfer during a water input event dictates whether macropores become blocked by the freezing of infiltrated water . In particular, it shows the range of macropore-matrix heat transfer that gives rise to contrasting unlimited and limited infiltration responses in frozen soils (Granger et al., 1984) and how these processes affect the timing and magnitude of runoff generation. Depending on the degree of macropore blockage, percolation through the frost zone may occur before soil thaws (Baker & Spaans, 1997;Mohammed et al., 2019;Stähli et al., 2004) or, alternatively, not until after complete ground thaw Pittman et al., 2020).
In all simulations (Test Cases 1−3), infiltration-runoff partitioning was most sensitive to thermal energy exchange governed by α wd in Equations 6 and 11, all else remaining equal, specifically its effect on ice formation in macropores that, in turn, controls macropore hydraulic conductivity and soil infiltrability. Simulations here expand the physical basis for conceptualizing runoff generation in frozen soils by assessing the subsurface flow and freeze-thaw processes that lead to different infiltration-runoff responses. To the best of our knowledge, the model of Larsbo et al. (2019) is the only other dual-permeability formulation for freezing of porewater in macropores, and it is thus hard to compare the suitability of the macropore-matrix heat transfer parameters used here with the literature. However, the conceptualization of the F I G U R E 9 Simulated vertical profiles of macropore water and ice saturation for Test Case 2 with infiltration in the macropore domain at various times showing (a) initially air-filled and ice-free macropores at 0.05 d, (b) the water and ice profiles when runoff begins at 0.083 d, and (c) the final state of blocked macropores when all input is partitioned to runoff at 0.125 d F I G U R E 1 0 Input (black), infiltration (blue), runoff (purple), and drainage (green) rates simulated in (a) Test Case 3a with a high rate of macropore-matrix heat transfer, and (b) Test Case 3b with low macropore-matrix heat transfer F I G U R E 1 1 Ice saturations at different depths in macropore domain simulated in (a) Test Case 3a with high macropore-matrix heat transfer, and (b) Test Case 3b with low macropore-matrix heat transfer geometry of the macropore-matrix interface over which heat transfer occurs is analogous to water and solute transfer in dual-permeability modeling studies in unfrozen soils (Köhne et al., 2009;Šimůnek et al., 2003). The first-order transfer equations are the same as those for water and solute transfer, provided that the effective hydraulic properties are replaced by the thermal properties (Equation 11). The value of α wd is related to the geometry of soil matrix aggregates, macropore areal density, and transport characteristics of the macroporematrix interface (Faúndez Urbina et al., 2019;Gerke & van Genuchten, 1993b;Jarvis et al., 1991;Nimmo, 2010). The diffusion pathlength (d) over which water, solute, and heat exchange occurs has the biggest influence on the rate of heat transfer and freezing of porewater in macropores, which governs subsequent hydraulic conductivity reduction . Conceptually, d is some average characteristic distance between macropores in a representative elementary volume of the soil profile (Faúndez Urbina et al., 2019;Gerke & van Genuchten, 1993b;Jarvis et al., 1991). As d is physically related to properties of the macropore network, it can be estimated from controlled experiments (Arora et al., 2011;Gerke & van Genuchten, 1996;Köhne & Mohanty, 2005). The results in Test Case 3 herein are very similar to Test Case 4 presented in Larsbo et al. (2019), and it should be noted that both simulation results illustrate that the effects of this interdomain transfer conceptualization on infiltration behavior can be highly counterintuitive, since it suggests a higher macropore density means that porewater in macropores freeze faster, resulting in less infiltration and more runoff (compared with a lower macropore density or longer diffusion pathlength). Given that this is one of the first models to use macroporematrix transfer relations to dictate water-ice phase change and associated reduction in hydraulic conductivity in macropores, further testing with field data on infiltration and runoff generation in frozen soils is critically needed. To this end, the terms in Equation 6 can alternatively be treated as a lumped empirical parameter to be calibrated to observed field data, which Vadose Zone Journal is a common practice (Bishop et al., 2015;Frey et al., 2012;Gerke & van Genuchten, 1996;Šimůnek et al., 2003).

CONCLUSIONS AND LIMITATIONS
We developed a dual-permeability model of preferential flow in frozen soils that considers the effects of porewater freezing in macropores. Incorporating matrix and macropore domains enabled the simulation of snowmelt infiltration into frozen soil with initially air-filled macropores, a situation commonly encountered in cold regions (Gray et al., 2001;van der Kamp et al., 2003), as well as porewater freezing in macropores and its effects on reducing infiltrability (Stähli et al., 1996;Watanabe & Kugisaki, 2017). The model includes a first-order expression to describe macropore-matrix heat transfer and the corresponding rate of water freezing in macropores . Simulations of laboratory experiments showed that the dual-permeability approach used here can reproduce the dynamics of preferential flow in frozen soils, including infiltration, drainage, and bypass flow due to macropores. The simplified heat transport description in the soil matrix was able to approximate soil thaw to an acceptable degree, which in this case meant reproducing the average rate of soil thaw. Results highlight the importance of macropore hydraulic conductivity, as well as macropore-matrix heat and water transfer on controlling preferential flow in frozen soils. Test-case simulations additionally demonstrated the important effect of porewater freezing along preferential flowpaths, and its influence on runoff generation. Simulations replicated a range of snowmelt partitioning behavior commonly observed in frozen soils, including subsurface conditions that produce rapid infiltration and drainage, the contrasting response of limited infiltration and enhanced runoff, and the evolution of runoff generation through time (Demand et al., 2019;Granger et al., 1984;Grant et al., 2019;Holten et al., 2018;Mohammed et al., 2019). Freezing of infiltrated water along preferential flowpaths was largely controlled by the geometrical and thermal characteristics of the macropore-matrix interface. This study presents the first attempt at including both frozen soil and preferential flow processes into a physically based model that can simulate coupled surface and subsurface hydrological processes, and thus there are limitations as well as opportunities for further development. The key limitation of the model presented is the simplification of heat transport and porewater freeze-thaw in the soil matrix. In addition to significantly reducing computational effort, the simplified matrix heat transfer descriptions were considered acceptable in this case where capturing the hydraulic and freezethaw behavior in macropores was the focus. However, certain aspects of frozen soil behavior cannot be simulated with the current heat transfer formulations, such as cyro-suction, transient latent heat effects in the bulk soil matrix, and advective heat fluxes in the soil matrix and macropores. However, the simplified heat transfer description enables HGS to be applied to simulate coupled surface and subsurface flows in two and three dimensions. Further research is underway to examine whether the controlling mechanisms identified at the profile scale play a similar role across larger hillslope and watershed scales. The difficulty of untangling the complex processes governing snowmelt partitioning dynamics necessitates more field and experimental data to validate these emerging models.

C O N F L I C T O F I N T E R E S T
The authors declare no conflict of interest.

A C K N O W L E D G M E N T S
The authors wish to thank Dr. Larry Bentley for helpful discussions on model development; Freda Pittman for conducting laboratory infiltration experiments; and Alberta Innovates, Alberta Environment and Parks, the Alberta Energy Regulator/Alberta Geological Survey, and the Natural Sciences and Engineering Research Council of Canada for funding support.