The diffusion equation is used to simulate the transport of energy in the soil.
where T is the soil temperature, C is the volumetric heat capacity [J/(m3K)], θ is the volumetric soil moisture content [m3/m3], K is the bulk thermal conductivity [W/(mK)], Lf is the latent heat of fusion (J/Kg), rho_ice is the ice density, θ_ice is the soil ice content [m3/m3], t is the time, and z is the depth. The last term on the right side of Eq. (1) represents the latent heat of released/consumed during the phase change. The soil thermal conductivity follows the parameterization of Peters-Lidard (Ref) given by:
where Ke , Ksat, and Kdry are the Kersten number, saturated, and dry soil thermal conductivities, respectively, and are defined in Peters-Lidard et al., 1998. The effective volumetric heat capacity, C, is calculated based on the respective fraction of each component (water, ice, air, and rock):
here, θ, θ_w, and θ_s are the total water content, the liquid water content and the maximum soil moisture content (porosity), respectively. The parameters Cw, Cice, Cair, and Csoil are the volumetric heat capacities of water, ice, air, and soil, respectively. The diffusion equation (Eq. 1) is discretized using the a fully implicit finite difference scheme and thus unconditionally stable.
Freezing-point depression equation: Frozen soils not only impact the soil thermal and hydrological properties, but its presence strongly affects the partition of precipitation into surface runoff and infiltration leading to a strongly coupled surface and subsurface system. Due to the formation of impermeable soil over the winter period, infiltration of water significantly reduces during the spring freshet generating high peaks in the surface discharge.
In unfrozen conditions, the liquid water content is equal to the total water content in the soil. However, when soil freezes the liquid water content also depends on the soil temperature. The soil matric potential, Ψ_s [m], under the assumption that soil water potential remains in equilibrium with vapor pressure over pure ice when ice is present in the soil, and neglecting soil water osmotic potential, is given by:
where T and To are the soil temperature and freezing temperatures, respectively (Ref), g is the acceleration of gravity [m/s^2] and Lf is the latent heat of fusion [J/Kg]. The soil matric potential, psi [m], as a function of soil moisture content is given by (ref: Clapp-Hornberger):
where Ψ_s, θ_s, and θ_l are saturated soil matric potential [m], saturated soil moisture [-] (porosity), and soil liquid moisture [-], respectively. The exponent b is the Clapp-Hornberger parameter. According to the ‘freezing equals drying’ approximation (Ref), we equate Eq. (4) and Eq. (5) to obtain an expression for liquid water content:
Equation (6) is called the “freezing-point depression equation” (Refs) and gives the maximum amount of liquid water (unfrozen soil moisture content) that can exist below the subfreezing temperature. The frozen soil moisture content (i, ice fraction) is given by: θ_ice = θ_t- θ_l , where t is the total water content given by soil water retention curve:
