How Soil Freezes and Thaws at a Snow-Dominated Forest Site in the U.S.—A Synthetic Approach Using the Soil and Cold Regions Model (SCRM)

Abstract

The freeze–thaw process controls several hydrologic processes, including infiltration, runoff, and soil erosion. Simulating this process is important, particularly in cold and mountainous regions. The Soil and Cold Regions Model (SCRM) was used to simulate, study, and understand the behavior of twelve homogenous soils subject to a freeze–thaw process, based on meteorological data at a snow-dominated forest site in Laramie, WY, USA, from 2010 and 2012. The relationships of soil pore size, soil particle contact, and meteorological data were varied. Our analysis of the model compared simulations using metrics such as soil frost depth, days with ice, and maximum ice content. The model showed that the freeze–thaw process was strongest in the period with a shallow snowpack, with particle packing within the soil profile being an important factor in this process; that soil texture and water content control soil thermal properties; and that water movement towards the freezing front was especially important in fine-textured soils, where water and ice were concentrated in the upper layers. Based on these results, future research that combines a broader set of soil conditions with an extended set of field meteorology and real soil data could elucidate the influence of soil texture on the thermal properties related to soil frost.

1. Introduction

Frozen soils are an important factor that can control hydrologic processes such as meltwater infiltration, runoff timing and amount, water retention characteristics, and soil erosion [1,2,3,4]. When the ambient temperature reaches the soil’s freezing point, the freeze–thaw (FT) phenomenon plays an important role and is of interest to fields such as hydrology [5], engineering [6], agriculture [7], and forestry [8]. The FT process is especially relevant in cold regions, particularly at high elevations, where the influence of lower temperatures can result in more frequently frozen soils. These regions where seasonal soil frost is present are commonly located above 30° N–S [8]. This permeable/semipermeable ice layer can have regional effects on hydrologic response and soil thermal characteristics [9,10].

The phase change within the soil from water to ice and vice versa is controlled mainly by the movement of liquid water towards the freezing front, a process similar to water movement when the soil is drying [11,12,13,14]. Characteristics of a decrease in temperature in frozen soils can include decreased hydraulic conductivity [15] and increased water viscosity [16]. Therefore, coupling water and heat flow is necessary to predict the FT process [17].

One-dimensional (1D) numerical models have been used for decades to understand and describe soil freeze–thaw processes and soil–plant–atmosphere interactions, e.g., the CoupModel (Coupled heat and mass transfer Model for soil–plant–atmosphere systems) [18], SHAW (Simultaneous Heat and Water Model) [19], the SCRM (Soil and Cold Regions model) [20,21], and the spatially distributed ParFlow-CLM [22]. Nevertheless, the SCRM is currently the most comprehensive model, considering both meteorological data and vegetation properties.

One of the main factors controlling soil frost is the timing and amount of snowfall [23,24,25]. Freezing progresses from the ground surface downwards. Snow cover/soil (insulation ability) and cold weather intensity and timing control how deep-frozen soils occur [26,27,28]. Snowpack is an excellent thermal insulator due to its low thermal properties [29,30,31] and can insulate the soil from meteorological conditions such as cold weather [32,33]. Hence, the snowpack controls the duration and depth of the frozen layers within a soil [34], but this phenomenon is also controlled by other factors such as soil properties, precipitation patterns, and temperature.

Soil texture is an important factor controlling thermal processes in freezing soils [35,36,37]. Frozen soils control not only water infiltration into soil, but also runoff during spring [38,39]. Several studies have related freeze–thaw processes to soil texture [40,41,42,43]. The scope of previous research is restricted to a specific or local context, and does not present a broader framework required for universal analysis.

Understanding how soils behave under different weather scenarios will explain some of the hydrological processes described above. High elevation regions where soils freeze during the whole or part of a winter season [1] are dispersed worldwide, and the timing of snow melt is crucial for regional water availability. A lot of work worldwide has been dedicated to analyzing the freeze–thaw processes but focusing on specific sites and soils [44], in synthetic made soils [45], and in laboratory tests [13,46]. Despite the importance of the freeze–thaw processes, there has been little work to examine soil behavior in alpine zones using a wide range of soil types and their effect on the freeze–thaw processes. In this study, we provide a comprehensive and academic assessment of the freeze–thaw processes across a wide range of soil types in two contrasting climatic conditions. This analysis was performed by addressing the following research questions: (1) Is soil texture a major driver of the freeze–thaw process? (2) How do soil hydraulic properties function as a secondary driver of the freeze–thaw process?

We hypothesize that soil particle size, soil water content, and depth of snowpack directly influence the buildup of ice. Specifically, fine-textured soils, thinner snowpack, and higher soil water content will correlate with greater buildup of ice within the soil. This study also introduces metrics that aim to explain the behavior of different soil types in the presence of freeze–thaw processes. While these hypotheses are known, the main goal of this study was to synthesize the understanding of soil frost on different soil textures within contrasting climate input data. We intend this study to be a contribution and a guide to the development of research on the subject in areas where freezing processes have not been studied extensively, such as in South America.

2. Materials and Methods

Twelve soil textures were analyzed using the Soil in Cold Regions Model (SCRM) [20]. Section 2.1 describes the study site, its parameters of interest for the model, the sources of meteorological data, as well as the configurations and theory behind the SCRM model. Section 2.2 describes the metrics and assumptions used to evaluate the model’s performance.

2.1. Study Site, Meteorological Data, and Model Configuration

2.1.1. Study Site

A snow-dominated forest site located in the Medicine Bow National Forest (41°22′30.0″ N 106°15′30.5″ W, about 50 km west of Laramie, WY, USA) was chosen to drive the model. The site is at 3000 m elevation, south-facing, with a slope of 12.5° and is dominated by conifers, specifically Picea engelmannii Parry ex. Engelm. (Engelmann Spruce) and Abies lasiocarpa (Hook.) Nutt. (Subalpine fir). The parameters were defined in terms of vegetation as follows: average tree height = 20 m, stand stem biomass = 12 kg m⁻², maximum stand Leaf Area Index (LAI) = 5 m² m⁻², and root mass = 4.4 kg m⁻² [20].

2.1.2. Meteorological Dataset

For the model, we used meteorological data from study [20] for two hydrologic years (2010–2011 and 2011–2012). The data used in the model are described in [20,47]. The meteorological data were derived from a GLEES (Glacier Lakes Ecosystems Experiments Sites) tower from 5 to 15 min intervals. Rainfall (and snow) was collected from a rain gauge also at GLEES. The data utilized in the model were relative humidity (%), air temperature (°C), wind speed (m s⁻¹), atmospheric pressure (Pa), precipitation (m), cloud cover fraction (dimensionless), atmospheric turbidity (dimensionless), leaf area index (m² m⁻²), and stem area index (m² m⁻²). Additionally, the wavelength (m), bandwidth (m), extraterrestrial radiation (J m⁻² sm⁻¹), water vapor absorption (1 m⁻¹), ozone absorption (1 m⁻¹), and miscellaneous gas optical depths were required to calculate incoming solar radiation at the Earth’s surface.

The monthly average precipitation and temperature distribution (Figure 1) in the 2010–2011 hydrologic year were characterized by temperatures ranging from −33.5 °C to 21.7 °C (0.6 °C average) and annual total precipitation of 2128 mm. The 2011–2012 hydrologic year was warmer and with less precipitation than 2010–2011. The temperature ranged from −23.24 °C to 23.39 °C (2.25 °C average) and annual total precipitation was 989 mm.

2.1.3. Model Configuration

The Soil and Cold Regions Model (SCRM) is a one-dimensional numerical model that couples water flow and heat transport in soil and snow using all three phases of water: liquid, vapor, and ice [20,21]. The model is driven by precipitation and the surface energy balance at the top boundary and free drainage water flow, and zero temperature gradient for heat transport at the bottom. SCRM is physics-based and requires few empirical parameters. The model can deal with canopy–ground surface energy balance and root water uptake. The SCRM model is based on five main physical processes:

• Snow water flow, which depends on gravity (water) and temperature gradient (vapor).
• Soil water flow, which depends on the pressure head gradient, temperature gradient, and gravity (ignored for water vapor).
• Snow–soil heat transport, which is based on the convection and conduction. The model represents bulk snow–soil heat capacity, latent heat of vaporization, latent heat of fusion, snow–soil thermal conductivity, liquid water heat capacity, water vapor heat capacity, liquid water gravity, and water vapor flux.
• Soil–plant–atmosphere transfer, which is based on the canopy and ground surface energy balance. Here, evaporation is in terms of liquid and vapor water flow (at the upper half of the soil profile), and the water vapor heat transfer at the canopy.
• Root water uptake, which is a sink term within the soil water flow equation. This sink term is based on macroscopic root analysis.

The model does not include preferential flow, 2D or 3D processes, solute concentration changes, or changes in organic matter content. The main governing equations are described below and the mathematical derivation of the equations describing the entire process can be found in [20,21].

Soil water flow and snow–soil heat transport are described by Equations (1) and (2), respectively:

$$\frac{\partial {\theta }_w}{\partial t}+\frac{1}{{\rho }_w}\frac{\partial \left(H_r{\rho }_{vs}{\theta }_a\right)}{\partial t}+\frac{{\rho }_i}{{\rho }_w}\frac{\partial {\theta }_i}{\partial t}=\frac{\partial }{\partial z}\left(K_{wh}\frac{\partial h}{\partial z}+K_{wh}+K_{wT}\frac{\partial T}{\partial z}+K_{vh}\frac{\partial h}{\partial z}+K_{vT}\frac{\partial T}{\partial z}\right)-S_w$$

(1)

$$\frac{\partial \left(CT\right)}{\partial y}+{\gamma }_v\frac{\partial \left(H_r{\rho }_{vs}{\theta }_a\right)}{\partial t}-{\gamma }_f{\rho }_i\frac{\partial {\theta }_i}{\partial t}=\frac{\partial }{\partial z}\left(\kappa \frac{\partial T}{\partial z}-C_w{q}_{w}T-c_v{\rho }_w{q}_{v}T-{\gamma }_v{\rho }_w{q}_v\right)-C_w{S}_{w}T$$

(2)

where t is time (s), ${\theta }_w$ is liquid water content (m³ m⁻³), $H_r$ is soil air relative humidity (dimensionless),${\rho }_{vs}$is saturated water vapor density (kg m⁻³),${\theta }_a$ is air content (m³ m⁻³),${\rho }_i$ is ice density (kg m⁻³),${\rho }_w$is liquid water density (kg m⁻³),${\theta }_i$is ice content (m³ m⁻³), z is vertical coordinate (m),$K_{wh}$is isothermal liquid water hydraulic conductivity (m s⁻¹), h is soil water pressure head (m),$K_{wT}$is thermal liquid water hydraulic conductivity (m² s⁻¹ K⁻¹),$K_{vh}$is isothermal water vapor hydraulic conductivity (m s⁻¹), $K_{vT}$is thermal water vapor hydraulic conductivity (m² s⁻¹ K⁻¹), $S_w$is root water uptake sink term (s⁻¹), C is the soil heat capacity (J m⁻³ K⁻¹), T is temperature (°C),${\gamma }_v$is latent heat of vaporization (J kg⁻¹),${\gamma }_f$is latent heat of fusion (J kg⁻¹), κ is soil thermal conductivity (J m⁻¹ s⁻¹ K⁻¹),$C_w$is liquid water heat capacity (J m⁻³ K⁻¹), $q_w$is liquid water flux (m s⁻¹), $c_v$is water vapor heat capacity (J kg⁻¹ K⁻¹), and $q_v$is water vapor flux (expressed as an equivalent liquid water flux in m s⁻¹). Note that $H_r$ = 1 for snow and 0 < $H_r$ < 1 for soil.

The energy balance equation can be expressed as:

$$\frac{\partial \left(CT\right)}{\partial t}+{\gamma }_v{\rho }_{vs}{\theta }_a\frac{d{H}_r}{dt}\frac{\partial T}{\partial t}+{\gamma }_v{H}_r\frac{\partial T}{\partial t}+{\gamma }_v{H}_r{\rho }_{vs}\frac{\partial {\theta }_a}{\partial t}+{\gamma }_f{\rho }_i\frac{d{\theta }_w}{dT}\frac{\partial T}{\partial t}=\frac{\partial }{\partial z}\left(\kappa \frac{\partial T}{\partial z}-C_w{q}_{w}T-c_v{\rho }_w{q}_{v}T-{\gamma }_v{\rho }_w{q}_v\right)-C_w{S}_{w}T$$

(3)

In Equation (3), it should be noted that $\frac{d{\theta }_a}{dT}$ is set to zero and T is solved according to [48].

The model profile configuration holds 81 nodes ranging from 0.01 to 10.0 m deep. The van Genuchten–Mualem (vG) parameterization model was used to calculate the vG parameters for 12 soils (Figure 2 and

Table 1. Soil texture distribution and vG parameters. Each soil number matches those in Figure 2.

Soil #	Clay %	Sand %	θr (m3 m−3)	θs (m3 m−3)	Alpha (1 m−1)	n	Ksat (m s−1)
1	5	90	0.051535	0.376943	3.3207	2.5032	3.06 × 10−5
2	5	85	0.044628	0.382347	3.8190	2.0024	1.77 × 10−5
3	10	65	0.041513	0.386307	3.0492	1.4038	5.77 × 10−6
4	30	60	0.071651	0.385478	2.7089	1.2740	2.28 × 10−6
5	45	50	0.084738	0.408392	2.9633	1.1924	2.97 × 10−6
6	20	40	0.062679	0.406259	0.9704	1.4966	1.89 × 10−6
7	35	30	0.084063	0.443478	1.2939	1.3892	1.52 × 10−6
8	60	20	0.097124	0.485228	2.1038	1.2058	3.01 × 10−6
9	20	20	0.069263	0.4327265	0.4932	1.6382	2.63 × 10−6
10	35	10	0.09205	0.479122	0.9078	1.4813	2.25 × 10−6
11	50	5	0.103863	0.509586	1.4156	1.3335	2.93 × 10−6
12	10	5	0.060647	0.486972	0.6897	1.6471	4.16 × 10−6

). The ideal soils represent the basic 12 soil textures: sandy, loamy sand, sandy clay loam, sandy loam, sandy clay, loam, clay loam, clay, silty loam, silty clay loam, silty clay, and silt. The vG parameters were extracted using the Pedotransfer Function ROSETTA [49] and soil textures (sand, clay, and silt %) as input. Each texture was chosen randomly but setting up the criteria to be closer to the center of the soil type in the USDA soil triangle (red stars in Figure 2). It should be noted that at the model soil profile, the first 0.6 m of soil comprises layers 1 to 60 (0.01 m each), while the final (bottom) describes the remaining 9.4 m (node 61 to 81).

2.2. Data Processing and Metrics Definition

The initial soil conditions (water content and temperature per node) were measured by [20] and then used as initial conditions in the model. The model includes two logic restrictions when it calculates the soil water pressure head, which are θ ≥ ${\theta }_r$and ${\displaystyle \sum }\theta +{\theta }_i<n\left({\theta }_s\right)$; both inequalities are physically based. Using the soil matric potential (MP) in the following way ensured that both statements are false. First, the initial water content was translated to pressure head values according to the van Genuchten equation and solving for suction pressure or matric potential (h).

Several metrics were developed to differentiate the behavior of the 12 soil textures:

(i)

Freeze–thaw cycles (FTCs) within a soil is when the ice water content (θ_i, in m³ m⁻³) goes above zero θ_i and then comes back again to zero. It was calculated as the number of annual cycles for each of the 12 soils.

(ii)

Maximum ice content measures the maximum ice content within the node. The first (0.01 m) and second (0.02 m) nodes were used to compute this metric because the maximum ice content can be found at these two depths.

(iii)

Total elapsed time (days/hours) where each soil had frost (i.e., days with ice). This value was computed within the first 5 nodes (0.01 to 0.05 m depth). It is important to note that this metric accounts for both continuous and discontinuous ice content.

(iv)

Temporal ice mass was calculated from the area under the curve of the ice mass time series, using the trapezoidal approximation method.

(v)

Maximum frost depth in each soil texture was defined as the last non-zero value within the soil profile.

3. Results

3.1. Snowfall Amount and Timing

In terms of snowfall, the first snowfall event differed between hydrologic years. In 2011–2012, it happened 20–30 days before that in 2010–2011, most likely because temperatures were lower in the first 5 days of the 2011–2012 hydrologic year compared to the previous hydrologic year. The snowpack behaviors in both hydrologic years were different. The 2010–2011 hydrologic year had a maximum snow depth of 4.64 m, whereas the 2011–2012 hydrologic year had a maximum snow depth of 1.99 m. The model did not show important differences in snowpack and Snow Water Equivalent (SWE) areas (days*m) between soil textures (

Table 2. Snow height and snow water equivalent (SWE) areas for both hydrologic years and all soils.

Hydrologic Year	2010–2011		2011–2012
Soil Texture	Snow Height (days*m)	Snow Water Equivalent (days*m)	Snow Height (days*m)	Snow Water Equivalent (days*m)
Sand	598.751	204.735	171.925	47.175
Loamy Sand	597.284	204.591	170.475	46.391
Sandy Clay Loam	588.703	201.180	165.708	44.876
Sandy Loam	585.536	199.899	163.761	44.876
Sandy Clay	584.888	199.552	162.717	44.410
Loam	587.246	200.521	164.303	44.615
Clay Loam	587.706	200.769	165.233	45.278
Clay	587.829	200.830	164.972	44.663
Silty Loam	589.067	201.323	165.693	45.426
Silty Clay Loam	590.139	201.706	166.422	45.700
Silty Clay	590.568	201.851	166.733	45.814
Silt	592.179	202.648	166.996	45.641
Difference	14.903	5.693	9.676	2.882

). The difference in snow depth area was 14 days*m and 9 days*m and a SWE difference of 5.6 and 2.8 days*m for 2010–2011 and 2011–2012 hydrologic years, respectively.

Minimal differences were found among each soil texture in terms of snow depth and SWE, as calculated by the model. In the 2010–2011 hydrologic year, the maximum snow depth difference was 4.4 cm and 1.7 cm for Snow Water Equivalent (SWE), while the 2011–2012 hydrologic year saw 3.7 cm for snow depth and 0.97 cm for SWE. Nevertheless, substantial differences were seen in the amount of snow, which is primarily due to higher precipitation levels and lower temperatures of the 2010–2011 hydrologic year.

3.2. Freeze–Thaw Metrics

3.2.1. Freeze and Thaw Cycles (FTCs)

The FTCs in the first 5 cm within the soil profile show significant variability across soil types (

Table 3. Freeze and thaw cycles within first 5 cm depth (first five soil nodes) and maximum ice content in both hydrologic years (Hy) analyzed.

	Soil Type	2010–2011 Hydrologic Year (Hy)					2011–2012 Hydrologic Year (Hy)
		Soil Depth (cm)					Soil Depth (cm)
		1	2	3	4	5	1	2	3	4	5
Freeze–Thaw Cycles	Sandy	1	1	1	1	1	2	1	1	1	3
	Loamy Sand	1	1	1	2	2	3	2	2	2	3
	Sandy Clay Loam	2	3	3	1	1	5	3	5	5	4
	Sandy Loam	2	2	2	1	1	6	6	5	4	1
	Sandy Clay	2	2	2	1	1	5	5	4	4	3
	Loam	3	2	2	1	1	6	5	5	6	4
	Clay Loam	3	2	2	1	1	6	4	4	5	3
	Clay	4	2	2	2	1	5	5	6	4	5
	Silty Loam	6	6	6	5	6	4	3	2	3	2
	Silty Clay Loam	2	3	2	1	1	8	7	7	4	4
	Silty Clay	3	3	2	1	1	5	5	5	3	7
	Silt	7	7	7	5	6	3	3	3	2	2
Maximum Ice Content	Sandy	0.325	0.325	0.241	0.18	0.145	0.214	0.177	0.152	0.154	0.108
	Loamy Sand	0.338	0.338	0.252	0.153	0.136	0.288	0.337	0.221	0.197	0.164
	Sandy Clay Loam	0.299	0.225	0.21	0.201	0.191	0.311	0.307	0.306	0.303	0.303
	Sandy Loam	0.251	0.223	0.215	0.211	0.162	0.253	0.148	0.138	0.125	0.115
	Sandy Clay	0.23	0.219	0.216	0.206	0.194	0.216	0.196	0.191	0.188	0.188
	Loam	0.289	0.21	0.192	0.183	0.175	0.304	0.302	0.301	0.296	0.292
	Clay Loam	0.318	0.248	0.241	0.2	0.188	0.249	0.179	0.165	0.152	0.222
	Clay	0.282	0.268	0.257	0.244	0.227	0.236	0.27	0.266	0.265	0.26
	Silty Loam	0.304	0.208	0.206	0.173	0.17	0.216	0.179	0.163	0.144	0.119
	Silty Clay Loam	0.356	0.27	0.25	0.223	0.207	0.261	0.228	0.203	0.184	0.165
	Silty Clay	0.331	0.294	0.277	0.252	0.212	0.272	0.241	0.204	0.186	0.168
	Silt	0.376	0.268	0.255	0.243	0.236	0.398	0.398	0.393	0.4	0.417

). Coarse soils such as sand and loamy sand had between one and three cycles in both periods. These values were the lowest within all 12 soil types, whereas the fine-textured soils such as clay and silt had a higher number of cycles.

3.2.2. Maximum Ice Content

The maximum ice content is typically between 1 and 2 cm depth (Table 3). In general, the 2010–2011 hydrologic year holds the highest values of max ice content within all soil textures. In regard to soil texture ice content, silt texture held higher ice content at the upper soil in both hydrologic years.

A special case is the silt soil in the 2011–2012 hydrologic year, where the ice content goes up within the soil profile. This situation occurred in the spring season when the snow begins to melt, between January and late May, and a water refreezing event occurred, yielding infiltration. At 5 cm depth, a peak of water content occurs and remains constant until 40 cm depth; therefore, in response to cold temperatures, water in the upper part of the soil profile freezes.

3.2.3. Days with Ice

Similar to the FTCs, colder temperatures translate to more days with ice (

Table 4. Days with ice and temporal ice mass metrics for both hydrologic years analyzed within the first 5 cm depth (first five soil nodes).

	Soil Type	2010–2011 Hydrologic Year					2011–2012 Hydrologic Year
		Soil Depth (cm)					Soil Depth (cm)
		1	2	3	4	5	1	2	3	4	5
Days with Ice	Sandy	240.58	244.17	257.67	255.75	254.58	151.00	200.17	198.92	197.33	189.08
	Loamy Sand	224.50	257.83	256.50	221.17	193.17	158.58	198.58	196.92	195.25	181.92
	Sandy Clay Loam	169.33	134.92	120.50	114.08	113.75	168.67	167.33	160.58	148.00	141.17
	Sandy Loam	171.08	165.08	158.67	152.83	150.33	173.25	167.92	160.17	129.92	118.83
	Sandy Clay	169.83	164.92	159.25	152.75	150.25	150.92	141.08	127.58	107.67	94.08
	Loam	164.08	131.00	128.33	127.58	127.58	165.75	163.75	151.17	139.50	131.25
	Clay Loam	177.75	154.83	146.58	132.58	133.50	154.67	148.83	137.75	124.42	128.83
	Clay	183.33	168.58	164.17	153.83	149.75	157.83	161.17	145.00	128.92	110.17
	Silty Loam	204.17	182.42	166.50	159.83	151.83	184.58	180.83	168.58	157.92	149.17
	Silty Clay Loam	202.58	168.08	142.17	131.25	132.42	185.08	180.83	166.50	152.50	146.17
	Silty Clay	219.75	195.17	167.92	153.92	150.58	185.33	180.83	163.17	151.08	136.67
	Silt	228.67	221.33	198.92	168.92	156.00	176.58	172.75	173.92	159.42	153.25
Temporal Ice Mass	Sandy	62.18	66.73	53.32	40.72	31.25	28.63	33.19	27.88	23.49	16.58
	Loamy Sand	62.93	73	49.63	26.47	21.85	35.50	43.73	36.59	31.41	24.60
	Sandy Clay Loam	26	18.3	16.03	14.62	13.76	21.20	18.67	15.65	13.44	12.56
	Sandy Loam	29.02	25.19	22.58	22.82	17.73	20.54	12.28	10.90	6.65	5.92
	Sandy Clay	30.43	28.01	24.93	23.57	21.03	11.76	9.86	7.55	5.72	4.61
	Loam	23.31	17.79	15.34	14.62	14.04	19.08	16.20	13.19	11.54	9.76
	Clay Loam	30.5	25.72	22.27	17.18	15.63	18.73	14.64	11.60	8.82	10.15
	Clay	37.94	34.12	31.22	27.56	23.96	19.23	15.30	12.91	8.63	6.49
	Silty Loam	25.43	23.05	20.62	18.43	16.44	23.35	20.24	16.51	14.48	12.18
	Silty Clay Loam	35.22	27.39	22.13	18.49	17.1	26.82	21.26	16.79	14.22	11.44
	Silty Clay	40.32	34.78	30.63	26.07	21.47	27.87	22.57	16.42	13.04	10.22
	Silt	37.11	31.92	28.11	24.87	22.66	32.87	28.72	26.40	25.11	22.97

). The very coarse and very fine-textured soils resulted in the most days with ice (Table 4; sand, loamy sand, and silt family soils). Interestingly, in the coarse soils, the highest volumetric ice content values were between 2 and 3 cm depth. This contrasts with the other soils, where freezing occurs within 1 cm depth. Loamy sand soil had the most days with ice during both periods due to its water content, as compared with sandy soil.

3.2.4. Temporal Ice Mass

Unlike the other metrics, temporal ice mass becomes evident within the first 5 cm depth (Table 4). The model computed the total, average, and maximum temporal ice mass per soil profile, and the mass ratio between the total and the maximum ice content. The temporal ice mass is concentrated between 1 and 2 cm depth, and the depth of the largest temporal ice areas is around 8 cm. The coarse soils have the greatest temporal ice mass in the upper layers, which is caused by water suction from lower layers and the capacity to hold more water than the other soils (low initial soil water content). On the contrary, fine-textured soils have a diminished ability to suck up water due to less pore availability and a higher initial soil water content.

3.2.5. Frost Depth

Coarse textured soils, such as sand and loamy sand, had the deepest soil frost among the soils. The fine-textured soils (sandy loam, sandy clay, and clay) had a shallower frost depth (

Table 5. Frost depth for both hydrologic years analyzed within first 5 cm.

Soil Type	2010–2011 Hy	2011–2012 Hy
	Frost Depth (cm)
Sandy	129	112
Loamy Sand	76	129
Sandy Clay Loam	35	35
Sandy Loam	29	15
Sandy Clay	14	19
Loam	35	27
Clay Loam	30	29
Clay	14	36
Silty Loam	41	26
Silty Clay Loam	31	22
Silty Clay	28	20
Silt	37	24

Hy: hydrologic year.

). In general, the colder period (2010–2011) has deeper frost depth than the warmer period (2011–2012). However, some exceptions were found. Sandy clay and clay went down to only 14 cm depth, and they had some of the lowest FTCs in the 2010–2011 period, yet the same soils in the 2011–2012 period had a higher frost depth. One explanation is the shallow snowpack and decreased insulation capacity. This result is noticeable in the FTCs, where more cycles were found in the 2011–2012 hydrologic year.

In general, maximum frost depth in 2010–2011 was at node 70 (1.12 m depth), whereas in 2011–2012, it occurred at node 68 (0.90 m). Ice was concentrated in the topsoil profile within the first 5 to 6 cm depth in both simulations.

4. Discussion

Variations in soil texture (i.e., pore size) and its hydraulic properties (i.e., soil water content) are the principal factors governing soil behavior within the freeze–thaw processes [43], and these factors also have a major influence on infiltration processes [13]. Thus, coarse soils versus fine-textured soil behave differently under the same meteorological conditions. Fine-textured soils, due to capillary forces, hold more water at subzero temperatures [50], whereas coarse soils have weaker capillary forces that allow more water to freeze within the profile. The present study provides a comprehensive analytical evaluation of the principal factors governing the freeze–thaw process based on a broad range of soil textures.

4.1. Model General Results

In general, the model has the ability to reproduce the main processes in the soil regarding the freeze–thaw process such as ice/water coexistence [51,52], water flow in soil with ice [53], and water movement towards the freezing front [17]. The model reproduces the coexistence between ice and water, which means that even subzero temperatures are not enough to freeze all water present in the soil.

For instance, in the 2010–2011 hydrologic year, the timing of when water freezes and soil ice content increases can be observed in the sand and clay soil types and is especially clear between January and July, where soil water content never decreases to zero, i.e., the residual water content remains unfrozen. This result can be explained because the model uses the Clausius–Clapeyron equation (liquid vapor equilibrium), from which it is possible to calculate the water potential using the temperature gradient [54]. This equation can also be used to calculate the liquid water content and ice content [55].

4.2. Freeze–Thaw Metrics

In general, when the soil is not 100% frozen, there will be some water flow upwards or downwards. The soil water content and pore size drives how much water can be frozen. Hence, in fine-textured soils, a higher energy input is necessary to freeze the soil because of the greater matric potential required; therefore, less liquid water will be frozen. These findings match those in study [53], which found that the middle section of the pore (represented as a capillary bundle) freezes, allowing water flow around the periphery until very low temperatures completely freeze the water within the pore. These results do not agree because of a shared equation in the model, but because of the physical basis of the model that allows water flowing through available pore spaces within the soil when it is not completely frozen.

All metrics results were normalized between 0 and 1. Each soil texture within the USDA triangle has a circle which represents the scaled result (Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7). A bigger circle (value of 1) indicates the highest value in that metric (i.e., a big circle in the FT cycles metric indicates a soil with several cycles within the year).

4.2.1. Freeze–Thaw Cycles (FTC)

Coarse soils such as sand are not as susceptible to freezing [12,42] as fine-textured soils such as loams, clays, and silt soils, which are more susceptible to freezing [56]. In fine-textured soil, if the freezing process is slow, the unfrozen moisture may move towards the freezing front, causing ice to accumulate in the upper layers [57]. When an FT process interacts with the soil, depending on the different soil textures, the soil structure might change. Fine-grained soil such as silt and clay are more vulnerable to FT cycles in their aggregate stability, while coarse-grained soils are more stable [3,58,59]. These behaviors are demonstrated in the model simulations (Table 2) where coarse soils (sand and loamy sand) have a lower number of FT cycles versus fine-textured soils such as clay and silt, which have higher fluctuations of FT cycles. This behavior is related to two processes: (i) the snowpack thickness (related to the air temperature), and (ii) the soil thermal conductivity.

The snowpack thickness and air temperature are important parts of the system, as they regulate system insulation capacity. A shallow snowpack in 2011–2012 coupled with high temperatures had less insulation capacity than a thicker snowpack and low temperature environment. It is clear that the 2011–2012 hydrologic year had higher temperatures in the winter and spring seasons than the 2010–2011 hydrologic year (Figure 1), which could explain the lower incidence of freeze–thaw cycles (e.g., warm environment, less ice; see normalized FTCs metric in Figure 3). Differences in the FT cycles between water years are apparent in the soil ice content for silt soil (1 cm depth). The 2010–2011 hydrologic year was more stable during winter, and it had more cycles when snowmelt and a shallow snow cover started to allow freeze–thaw processes in the soil during the spring season. In contrast, in the 2011–2012 hydrologic year, the ice content was quite unstable until the snowpack melted due to higher temperatures. For instance, fine-textured soils such as silt can have several freeze–thaw cycles, allowing water to infiltrate the soil below a thick snow cover. Under a shallower snow cover (2011–2012 hydrologic year), silt soil will have more FT cycles because of decreased insulation capability of the snow cover, which will produce more ice but also more thaw processes.

Higher water content can also enhance heat flow in soils, except for sand [60,61,62]. One possible explanation for this exception could be that the lower water content, particle contact (due to grain size and air content), and high permeability of sand do not allow many FT cycles, allowing for a roughly constant ice content. Silt soil has the lowest thermal conductivity among the four principal soils, but higher water content, better particle contact, low permeability, and higher water retention (due to pore size), which enhance heat flow.

Figure 3. Freeze–thaw cycles for 1 cm depth, 2010–2011 and 2011–2012 hydrologic years. Results were normalized ranging from 0 to 1; each circle represents the scaled result. A bigger circle indicates the highest value in the metric.

Figure 3. Freeze–thaw cycles for 1 cm depth, 2010–2011 and 2011–2012 hydrologic years. Results were normalized ranging from 0 to 1; each circle represents the scaled result. A bigger circle indicates the highest value in the metric.

4.2.2. Maximum Ice Content

Coarse soils such as sand and loamy sand have less thermal conductivity because of the greater air content, lower soil water content, and less particle contact, causing a decrease in both convection and conduction. The particle–particle contact is an efficient heat transfer process [63] which make these soils thermal insulators, and as a result, less ice builds up within coarser soils.

Our results from the maximum ice content (example for all soils at 1 cm depth in Figure 4) metrics match the finer textured soil behavior, experiencing greater convection from the higher water content and greater conduction due to particle–particle contact. Fine-textured soils yield less ice due to efficient heat transfer processes, which causes more FT cycles to occur and increases ice production [64].

Figure 4. Maximum ice content, 2010–2011 and 2011–2012 hydrologic years. Results were normalized ranging from 0 to 1; each circle represents the scaled result. A bigger circle indicates the highest value in the metric.

Figure 4. Maximum ice content, 2010–2011 and 2011–2012 hydrologic years. Results were normalized ranging from 0 to 1; each circle represents the scaled result. A bigger circle indicates the highest value in the metric.

4.2.3. Days with Ice

Ambient air temperatures are a crucial factor in ice formation. Colder temperatures translate to more days with ice. In the first 2 cm within the soil profile, the snowpack is insulating the soil against warm air temperature; thus, a shallow snowpack will result in less ice formation and fewer days with ice, although the first 5 cm within the soil profile reflect how low air temperatures drive the soil ice content for sand, clay, and silt soils (Table 4). The 2010–2011 period had more snow and lower temperatures than the 2011–2012 period, which caused significantly more days with ice.

As a general comment, it seems that ice concentration in the upper layers may be due to water movement upwards due to pressure and temperature gradients. It may also be controlled by the differences in water content (silt and silty clay) and pore size (sand and loamy sand), i.e., bigger pores do not need higher pressure values to hold water.

For example, in study [20], the soil layer was a mix that contained a majority of sandy soil at the bottom and small layers of loamy soils at the top profile, and this soil had the fewest days with ice (e.g., 139 days at 1 cm depth in 2010–2017 period). These results are probably due to the saturated and residual water contents. Recalling the van Genuchten parameters, in [20], soils had the lowest saturated water content among the soils and zero residual water content. With less water available in the soil profile, there will likely be less freezing.

The days with ice metric (example for all soils at 1 cm depth in Figure 5) showed little difference between all 12 soils. This metric is related to two processes that are difficult to distinguish from one another within the model context: (i) fine-textured soils’ water content will increase the days with ice due to more water availability, and (ii) the diminished thermal conductivity of coarse soils can allow more ice formation and a slower or reduced thawing process.

Figure 5. Days with ice, 2010–2011 and 2011–2012 hydrologic years. Results were normalized ranging from 0 to 1; each circle represents the scaled result. A bigger circle indicates the highest value in the metric.

Figure 5. Days with ice, 2010–2011 and 2011–2012 hydrologic years. Results were normalized ranging from 0 to 1; each circle represents the scaled result. A bigger circle indicates the highest value in the metric.

4.2.4. Temporal Ice Mass

The temporal ice mass is concentrated between 1 and 2 cm depth, and the depth of the largest temporal ice areas is around 8 cm. The coarse soils have the greatest temporal ice mass in the upper layers, which is caused by water suction from lower layers and the capacity to hold more water than the other soils (low initial soil water content). On the contrary, fine-textured soils have a diminished ability to suck water up due to less pore availability and higher initial soil water content (example at 1 cm depth in Figure 6).

Figure 6. Temporal ice mass, 2010–2011 and 2011–2012 hydrologic years. Results were normalized ranging from 0 to 1; each circle represents the scaled result. A bigger circle indicates the highest value in the metric.

Figure 6. Temporal ice mass, 2010–2011 and 2011–2012 hydrologic years. Results were normalized ranging from 0 to 1; each circle represents the scaled result. A bigger circle indicates the highest value in the metric.

4.2.5. Frost Depth

In general, the colder period (2010–2011) had deeper frost depth (Table 2 and Figure 7) than the warmer period (2011–2012). However, some exceptions were found. Sandy clay and clay froze to 14 cm depth, and they each had some of the lowest FT cycles in the 2010–2011 period, yet the same soils in the 2011–2012 period had a higher frost depth. One explanation is the shallow snowpack and thus a decreased insulation capacity. This result is noticeable in the FT cycles as well, where more cycles were found in the 2011–2012 hydrologic year.

Figure 7. Frost depth, 2010–2011 and 2011–2012 hydrologic years. Results were normalized ranging from 0 to 1; each circle represents the scaled result. A bigger circle indicates the highest value in the metric.

Figure 7. Frost depth, 2010–2011 and 2011–2012 hydrologic years. Results were normalized ranging from 0 to 1; each circle represents the scaled result. A bigger circle indicates the highest value in the metric.

The thermal conductivity of coarse and fine-textured soils is different, and consequently, how the heat is transferred also varies. Hence, in a warmer environment, such as in 2011–2012, with a shallow snowpack, the temperature oscillation had a greater effect on these soils, and earlier snow accumulation will result in decreasing soil frost depth [65]. This delay can also lead to soil frost, which makes the soil less permeable for the spring snowmelt, as observed in 2011–2012, but thick snow cover appears to diminish this effect (e.g., [66,67]).

Frost depth can be explained then by (i) the soil permeability, (ii) the soil susceptibility to frost, (iii) the soil pore size, and (iv) the soil water content. Therefore, in fine-textured soils (clay and silt), it is more difficult for water to migrate towards the freezing front due to low permeability, which results in less ice formation [68]. Clay soils are better insulators than silt and sandy soils and more susceptible to forming ice layers due to their higher water content, which fills in the voids within the soil matrix [69]. Contrary to the behavior in clays, in coarse, sandy soils, water migrates towards the freezing front through the available void spaces [70]; thus, sand will likely have a greater frost depth than other soil textures. Moreover, because sand has more void spaces and lowest water content, it will have more air-filled porosity, and thus the sand will have lower values of thermal conductivity, reflecting the low thermal conductivity of air. The model did not account for these low values in the calculation. Instead, by using the ice and soil conductivities, which are higher than the thermal conductivity of air, the simulation erroneously showed a high thermal conductivity for sand.

In general terms, coarse soils such as sand and loamy sand have less thermal conductivity because of the greater air content, lower soil water content, and less particle contact, which causes a decrease in both convection and conduction. The particle–particle contact is an efficient heat transfer process [63], which make these soils thermal insulators, and as a result, less ice builds up within coarser soils.

4.3. Final Remarks

Thermal conductivity (λ) is enhanced by particle contact and the increasing water content in fine-textured soils, whereas in sandy and loamy sand coarse soils, low water content and high air content will reduce the thermal conductivity. Additionally, the thermal conductivity depends on the water phase evolution [71].

Two assumptions made by [72,73] were contrasted regarding soil thermal conductivity and soil texture. Ref. [72] mentioned that when the thermal conductivity and grain size decrease, this might be due to the presence of more particles within a soil; hence, more thermal resistance would exist among the particles, which seems to be contrary to our results. However, Ref. [74] mentioned that soils with low thermal conductivities might experience larger gradients in surface temperature, which agrees with our fine-textured results.

In contrast, Ref. [73] mentioned when the contact between soil particles increases, heat transfer is more efficient. Therefore, in fine-sized particle soils (silt and clay), the particle contact is more efficient than coarse soils. An example of this is the variability observed in finer textured soil within the freeze–thaw cycles. Ultimately, a primary driver is that the water content plays an important role in facilitating thermal exchange between soil particles, as our results also show.

The thickness of snowpack plays an important role as a soil insulator affecting soil temperature and energy fluxes into the subsurface [75,76,77], the temperature becoming another factor governing freeze–thaw process. A shallow snowpack allows more ice formation and less infiltration within the soil profile [78] and may increase the freeze–thaw cycles within the soil [79,80,81,82]. By contrast, a thicker snowpack insulates the soil, thereby maintaining a more moderate temperature and less building up of ice. In summary, thin snowpack during winter results in low soil temperatures, more FT cycles, and deeper and extensive soil freezing [83].

In vegetated areas with shallower snowpack, the uptake of water by roots has been shown to influence the soil water content. This is especially true in early spring, when the growing season begins after snowmelt. Plants can access root water and photosynthesize sooner due to a thinner snowpack and more incoming light. A shift in the snowfall (late snowpack) may affect the soil frost duration and depth [84,85], and soil freezing could lead to fine root mortality, affecting a tree’s development [77]. For instance, the simulated shallow snow cover in the 2011–2012 hydrologic year shows decreasing water content in silt and clay soil types. This behavior can be explained by (i) refreezing water from snow melt, (ii) increasing incoming radiation due to a thin snowpack, and (iii) root water uptake. The activity of plants beneath the snowpack can also disturb the frozen soil and allow snow melt [30]; however, this activity depends on the thickness of snowpack because plants need light to grow. The amount of available unfrozen water dictates plant uptake in late winter to early spring [86,87]. Because the 2011–2012 hydrologic year had more unfrozen water, it is possible that the decrease in soil water content is due to root water uptake stimulated by an increase in incoming light.

5. Conclusions

In two contrasting climatic years, we used a 1D numerical model (SCRM) coupling heat transport and water flow to synthesize the understanding of soil frost on 12 different soil textures at a snow-dominated site. Model simulation on these soils agreed with the knowledge about soil freezing–thawing processes. As we hypothesized, soil texture and soil hydraulic properties, such as water content, are the primary drivers that control the freeze–thaw process in soil. However, weather is another important factor affecting this process. Our results show that coarse soils are less conductive and subject to less convective processes (i.e., a good insulator) which combine to prevent ice formation within the soil profile. In contrast, the high particle–particle contact in fine-textured soil results in greater conductive heat transfer, and the high water content enhances thermal conductivity and convection processes. The metrics accurately described all these soil textures variations in the freeze–thaw process.

The freeze–thaw process is strongest in the period with the late snowpack according to the metrics described. Therefore, in the presence of a late snowfall and subzero air temperatures, the soil will freeze easily. Our results have shown that late snowfall will lead to more days with ice and a greater concentration and depth (especially in coarse textured soils). These results are due to the absence of the insulating capacity of snowpack.

Despite the inherent limitations of a synthetic simulation without real field data to compare, our study provides a methodological framework and metric references for future studies addressing and inferring freeze–thaw processes in soil in a snow-dominated forest site. This method is particularly important for understanding the implications of climate change in high mountain regions, where increasing temperatures will greatly affect rainfall patterns and the snowpack distribution [87,88,89].