Interaction of soil water and groundwater during the freezing–thawing cycle: field observations and numerical modeling

Freezing-induced groundwater-level decline is widely observed in regions with a shallow water table, but many existing studies on freezing-induced groundwater migration do not account for freezing-induced water-level fluctuations. Here, by combining detailed field observations of liquid soil water content and groundwater-level fluctuations at a site in the Ordos Plateau, China, and numerical modeling, we showed that the interaction of soil water and groundwater dynamics was controlled by wintertime atmospheric conditions and topographically driven lateral groundwater inflow. With an initial water table depth of 120 cm and a lateral groundwater inflow rate of 1.03 mm d⁻¹, the observed freezing and thawing-induced fluctuations of soil water content and groundwater level are well reproduced. By calculating the budget of groundwater, the mean upward flux of freezing-induced groundwater loss is 1.46 mm d⁻¹ for 93 d, while the mean flux of thawing-induced groundwater recharge is as high as 3.94 mm d⁻¹ for 32 d. These results could be useful for local water resources management when encountering seasonally frozen soils and for future studies on two- or three-dimensional transient groundwater flow in semi-arid and seasonally frozen regions. By comparing models under a series of conditions, we found the magnitude of freezing-induced groundwater loss decreases with initial water table depth and increases with the rate of groundwater inflow. We also found a fixed-head lower boundary condition would overestimate freezing-induced groundwater migration when the water table depth is shallow. Therefore, an accurate characterization of freezing-induced water table decline is critical to quantifying the contribution of groundwater to hydrological and ecological processes in cold regions.

1 Introduction

Frozen soils, which can be divided into permafrost and seasonally frozen soils, have great impacts on many hydrological, hydrogeological, and ecological processes (Nelson, 2003; Kurylyk et al., 2014; Schuur et al., 2015; Walvoord and Kurylyk, 2016; Evans et al., 2018). In seasonally frozen soils, the behaviors of subsurface water flow and storage are distinct from those in colder regions with permafrost as well as in warmer regions without frozen soil (Ireson et al., 2013). Understanding the effects of soil freezing and thawing on subsurface water flow and storage is a necessary path for better water resources management in seasonally frozen regions (Yu et al., 2020).

By blocking pores and reducing hydraulic conductivity, soil freezing could lead to decreased infiltration (van der Kamp et al., 2003; Iwata et al., 2008; Demand et al., 2019) and limited surface evaporation (Kaneko et al., 2006; Wu et al., 2016). As a result of cryosuction generated by soil freezing (Williams and Smith, 1989; Hohmann, 1997; Yu et al., 2018; Zheng et al., 2021), water migrates from the unfrozen zone to the freezing front. By assuming that freezing-induced water migration is restricted within the shallow part of the unsaturated zone, many numerical studies deployed a free-drainage lower boundary condition and examined the effects of freezing-induced water redistribution on infiltration and runoff (Cherkauer and Lettenmaier, 1999; Okkonen et al., 2017; Zhang et al., 2017). In fact, when the water table is shallow, groundwater in the saturated zone can also be migrated to the freezing front (Harlan, 1973; Shoop and Bigl, 1997; Stähli et al., 1999; Hansson and Lundin, 2006; Ireson et al., 2013; Chen et al., 2019). Although there have been many field observations of freezing-induced groundwater-level decline since the 1950s (Drescher, 1955; Schneider, 1961; Vinnikov et al., 1996; Daniel and Staricka, 2000; van der Kamp et al., 2003; Ireson et al., 2013; Zhang et al., 2019), numerical research that satisfactorily considers both freezing-induced water gain in the frozen zone and freezing-induced groundwater-level decline in the saturated zone is limited.

Freezing-induced water migration from the saturated zone to the freezing front could strengthen the possibility of frost heave (Chamberlain, 1981; Bronfenbrener and Bronfenbrener, 2010; Rui et al., 2019). During the thawing stage, accumulation of thawed water above the frozen zone would accelerate soil evaporation (Fetzer et al., 2017; Vanderborgh et al., 2017; Li et al., 2020), which is critical to soil salinization (Liu et al., 2009; Lopez et al., 2010; Bechtold et al., 2011). Therefore, the amount of freezing-induced water migration from the saturated zone to the frozen zone is key to understanding the effects of soil freezing on engineering problems and ecological processes (Yu et al., 2020). In some previous experimental or numerical studies, freezing-induced water migration from groundwater to the frozen zone was obtained by assuming a fixed water table (Shoop and Bigl, 1997; Hansson and Lundin, 2006; Alkhaier et al., 2012; Chen et al., 2019). It was found that a fixed-head lower boundary condition, even if the fixed head is close enough to the mean water table, would overestimate groundwater loss to support evapotranspiration because such a boundary condition implies that loss of groundwater can be replenished instantaneously (Zhu et al., 2010). However, the degree of overestimation of freezing-induced water migration using a fixed-head lower boundary condition remains unknown.

As initially proposed by Hubert (1940) and Tóth (1962), lateral groundwater flow is ubiquitous in regions with undulating topography and water table. Although wintertime water table fluctuations caused by freezing and thawing have been widely recognized (Willis et al., 1964; van der Kamp et al., 2003; Zhang et al., 2019), the simultaneous contribution of lateral groundwater flow and soil freezing–thawing to wintertime water table fluctuations has received little attention. Ireson et al. (2013) found the wintertime water table recession in their two study sites in the Canadian Prairies was partly due to lateral groundwater outflow, and Jiang et al. (2017) found the pattern of wintertime water table fall–rise in their study site in the Ordos Plateau, China, was also influenced by lateral groundwater inflow. Therefore, the wintertime water table dynamics is a combined product by atmospheric conditions and lateral groundwater inflow–outflow, instead of a prior known fixed boundary condition.

To obtain a clear understanding of wintertime water migration and water budget in regions with shallow groundwater, it is appealing to numerically combine hydrological processes in both saturated and unsaturated zones, including topographically driven lateral flow in the saturated zone, freezing-induced water migration from the saturated to the unsaturated zone, and the thawing-induced water movement. In the current study, the field site reported in Jiang et al. (2017), with a shallow water table and lateral groundwater inflow, is examined as an example. We also use a series of scenarios with different water table depths and lateral groundwater inflow rates to numerically investigate their controls on water table fluctuations and freezing-induced water gain in the frozen zone.

2 Methods

2.1 Study area and conceptual model

The Ordos Plateau in northwestern China has a semi-arid and seasonally frozen climate and is well known for the occurrence of topography-driven regional groundwater flow (Hou et al., 2010; Jiang et al., 2018). The main aquifer is the thick Cretaceous sandstone, which is overlain by a thin layer of sand. Due to the spatially undulating topography and low permeability of the sandstone aquifer, water table undulation is a subdued replica of the topography. On the catchment scale, groundwater generally flows laterally from topographic highs to topographic lows, with a downward component of groundwater flow in the recharge area and an upward component in the discharge area. In the discharge area, as a result of the contribution of lateral groundwater flow from topographic highs towards the lake, the water table is shallow near the lake. The Wudu Lake catchment (Wang et al., 2015; Zhao et al., 2020) is one of the small catchments located to the southeast of the first-order water divide of the Ordos Plateau. Based on measurements of 53 wells in the catchment (the location of wells can be found in Jiang et al., 2018), water table depth increases from less than 0.5 m in topographic lows to 26.1 m near the divide. In topographic lows with water table ranging from around 1 m to around 2 m, the typical vegetation is Achnatherum splendens.

Figure 1 The conceptual model showing main hydrological and hydrogeological processes in semi-arid regions during the freezing–thawing cycle.

The Otak national meteorological station, which is around 35 km away from the center of the Wudu Lake catchment, is the nearest meteorological station with long-term meteorological observations. From 1955 to 2016, the annual mean precipitation is 265.0 mm, including an annual mean snowfall of 16.0 mm, and the annual mean pan evaporation is 1370.0 mm. Note that the amount of snowfall is determined by weighting the water equivalent of snowfall, and the accuracy of precipitation and pan evaporation measurements is ± 0.1 mm. July has the largest monthly mean temperature equaling 22.5 ^∘C, while January has the lowest monthly mean temperature equaling −10.4 ^∘C. The period from late November to late March corresponds to the freezing–thawing cycle. As shown in Jiang et al. (2018), freezing-induced groundwater-level decline occurs in the monitoring well DK2, located in topographic lows, but does not occur in the monitoring well DK1, located in topographic highs where the water table is deep.

In the Wudu Lake catchment, the occurrence of numerous flowing wells in topographic lows indicates an upward hydraulic gradient (Wang et al., 2015), which also induces groundwater inflow into the shallow part of the aquifer. At the monitoring site, lateral groundwater inflow had been found to be an indispensable component of groundwater budget (Jiang et al., 2017). Figure 1 shows the control of main hydrological and hydrogeological processes on groundwater level during the freezing–thawing cycle. The processes leading to water-level decline include water loss induced by freezing and evaporation, and the processes leading to water-level rise are lateral groundwater inflow throughout the whole period and infiltration of thawed water during the thawing stage.

2.2 Field measurement

To examine how atmospheric conditions control groundwater dynamics in topographic lows, we set up a monitoring profile of soil water contents and soil temperatures adjacent to the groundwater-level monitoring well DK2. Soil temperatures and liquid water contents in the Quaternary sands are recorded at eight depths (at around 10, 20, 30, 50, 70, 90, 110, and 150 cm below surface) using 5TM sensors produced by Decagon Devices. The 5TM sensors, which measure the soil dielectric permittivity to represent liquid soil water content (Zheng et al., 2015), have an accuracy of around ± 2 % volumetric water content. To ensure the representativeness of measurement, we performed site-specific calibration by establishing the correlation between liquid water content measured by the 5TM sensors and that measured by the gravimetric method. It has been reported that 5TM sensors are accurate enough to measure the liquid water content, even if soil is partially frozen (Yang et al., 2013; Zheng et al., 2017; Xue et al., 2020).

Table 1 The measured soil texture and bulk density and calibrated hydraulic parameters.

To obtain the bulk density, 12 soil samples were collected by the cutting-ring method, with a maximum depth of 120 cm and a resolution of 10 cm (Table 1). We also collected soil samples from the 12 depths for measuring soil particle size using the Mastersizer 2000 instrument (Malvern Instruments, England). The saturated hydraulic conductivity of soil samples from selected layers was measured by the HYPROP (https://www.metergroup.com/environment/products/hyprop-2/, last access: 20 July 2021). The soil samples with low clay content above 70 and below 100 cm were measured to be around 43.2 m d⁻¹, while soil with higher clay content at depths ranging between 70 and 100 cm is measured to be 1.9 m d⁻¹.

Figure 2 The contour map showing the evolution of soil temperature at the study site during the freezing–thawing cycle. Also shown is the monitored fluctuating groundwater level.

In this study, the data collected during 28 November 2015 and 1 April 2016 are used as an example. As shown in Fig. 2, water table depth increases from 115 cm on 28 November 2015 to 143 cm on 29 January 2016, which corresponds to the stage with an increasing frost depth.

2.3 The SHAW model

By assuming that the mechanisms of water transport in partially frozen soils are similar to those in unsaturated soils, Harlan (1973) pioneered modeling studies of heat–fluid transport with freezing–thawing by considering the co-existence of a frozen zone, an unsaturated unfrozen zone, and a saturated unfrozen zone. The Simultaneous Heat and Water (SHAW) model (Flerchinger and Saxton, 1989) is one of the most robust Harlan-type models for freezing-induced water migration in the one-dimensional (1D) domain (Hayhoe, 1994; DeGaetano et al., 2001; Kahimba et al., 2009; Li et al., 2010; Zhang et al., 2017). A 1D soil column is justified to account for water redistribution induced by freezing and thawing because water flow in the unsaturated zone is predominantly vertical (Stephens, 1996; Romano et al., 1998; van Dam and Feddes, 2000). The contributions of liquid water, ice, and vapor contents have been considered in the water balance equation, and the contributions of heat conduction, phase change, liquid flow, and vapor gas diffusion have been considered in the energy balance equation. By adding a source term representing the flux of groundwater inflow–outflow in the water balance equation, it has the ability to simultaneously account for the contributions of lateral groundwater flow and freezing–thawing in the soil column. Therefore, the SHAW model is suitable to couple the groundwater dynamics with soil water dynamics during the freezing–thawing cycle, as observed in our field site.

Considering convective heat transfer by liquid, latent heat transfer by ice, and latent heat transfer by vapor, the energy balance equation in the soil matrix is expressed as

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& C_{\mathrm{s}}{\displaystyle \frac{\partial T}{\partial t}}-{\mathit{\rho }}_{\mathrm{i}}{L}_{\mathrm{f}}{\displaystyle \frac{\partial {\mathit{\theta }}_{\mathrm{i}}}{\partial t}}\\ & ={\displaystyle \frac{\partial }{\partial z}}\left[k_{\mathrm{T}}{\displaystyle \frac{\partial T}{\partial z}}\right]-{\mathit{\rho }}_{\mathrm{l}}{c}_{\mathrm{l}}{\displaystyle \frac{\partial q_{\mathrm{l}}T}{\partial z}}-L_{\mathrm{v}}\left[{\displaystyle \frac{\partial q_{\mathrm{v}}}{\partial z}}+{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{v}}}{\partial t}}\right],\end{array}\end{array}$$

where T is temperature [Θ], θ_i is volumetric ice content [–], z is soil depth [L], t is time [T], ρ_i, ρ_l, and ρ_v are densities of ice, liquid water, and vapor [M L⁻³], respectively, L_f and L_v are latent heat of fusion and vaporization [L² T⁻²], respectively, q_l is liquid water flux [L T⁻¹], q_v is water vapor flux [M L⁻² T⁻¹], c_l is specific heat capacity of water [L² T⁻² Θ⁻¹], C_s is volumetric heat capacity of soil [M L⁻² T⁻² Θ⁻¹], and k_T is thermal conductivity of soil [M L² T⁻³ Θ⁻¹]. C_s of soil is determined by the sum of the volumetric heat capacities of the soil constituents:

$$\begin{array}{}\text{(2)}& C_{\mathrm{s}}=\sum {\mathit{\rho }}_j{c}_j{\mathit{\theta }}_j,\end{array}$$

and k_T of soil is calculated by the following equation proposed by de Vries (1963):

$$\begin{array}{}\text{(3)}& k_{\mathrm{T}}={\displaystyle \frac{\sum m_j{k}_j{\mathit{\theta }}_j}{\sum m_j{\mathit{\theta }}_j}},\end{array}$$

where ρ_j, c_j, θ_j, k_j, and m_j are the density, specific heat capacity, volumetric fraction, thermal conductivity, and weighting factor of the jth soil constituent (sand, silt, clay, water, ice, and air).

By considering ice content change and vapor flux, the water balance equation is written as

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial {\mathit{\theta }}_{\mathrm{l}}}{\partial t}}+{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{l}}}}{\displaystyle \frac{\partial {\mathit{\theta }}_{\mathrm{i}}}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left[K\left({\displaystyle \frac{\partial \mathit{\psi }}{\partial z}}+\mathrm{1}\right)\right]+{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{l}}}}{\displaystyle \frac{\partial q_{\mathrm{v}}}{\partial z}}+U,\end{array}$$

where θ_l is liquid water content [–], K is hydraulic conductivity [L T⁻¹], ψ is soil matric potential [L], and U is a source/sink term for water flux [T⁻¹]. The lateral inflow rate per unit volume, which is one form of source, is calculated using Darcy's law based on a constant horizontal hydraulic gradient at selected nodes within the saturated zone. Soil water characteristic curves, which describe the relationship between soil matric potential and liquid water content, are defined for both unfrozen and frozen soils. Here, the van Genuchten equation (van Genuchten, 1980) is selected, which is written as

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{{\mathit{\theta }}_{\mathrm{l}}-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}}={\left(\mathrm{1}+{\left|\mathit{\alpha }\mathit{\psi }\right|}^n\right)}^{-m},\end{array}$$

where θ_s and θ_r are saturated and residual water content [–], α [L⁻¹], n [–] and m [–] are empirical parameters. α equals the inverse of the air-entry value, n is the pore-size distribution index, and $m=\mathrm{1}-\mathrm{1}/n$. In the frozen zone with co-existence of liquid water and ice, the matric potential, which is strongly dependent on temperature, can be obtained by the following equation (Fuchs et al., 1978):

$$\begin{array}{}\text{(6)}& \mathit{\psi }={\displaystyle \frac{L_{\mathrm{f}}}{g}}{\displaystyle \frac{T}{T+\mathrm{273.15}}},\end{array}$$

where g is the acceleration of gravity [L T⁻²]. Equation (6) indicates that as the negative temperature increases, the soil suction also increases (matric potential becomes more negative).

In the unfrozen zone, the hydraulic conductivity is computed by

$$\begin{array}{}\text{(7)}& K=K_{\mathrm{s}}{S}_{\mathrm{e}}^l{\left[\mathrm{1}-{\left(\mathrm{1}-S_{\mathrm{e}}^{\mathrm{1}/m}\right)}^m\right]}^{\mathrm{2}},\end{array}$$

where K_s is the saturated hydraulic conductivity [L T⁻¹], S_e [–] is the effective saturation calculated by (${\mathit{\theta }}_{\mathrm{l}}-{\mathit{\theta }}_{\mathrm{r}})/({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}$), and l [–] is the pore-connectivity parameter which is assumed to be 0.5 (Mualem, 1976). In the frozen zone, occupation of ice in the pores would decrease the available porosity and thus hydraulic conductivity. When the available porosity is lower than 0.13, K is assumed to be 0; when the available porosity is above 0.13, K computed from Eq. (7) is reduced linearly with ice content assuming zero conductivity at an available porosity of 0.13 (Flerchinger, 2017).

The freezing–thawing cycle is usually accompanied by occasional or frequent snowfall events. In the SHAW model, precipitation is assumed to be snowfall if the air temperature is below 0 ^∘C. When snow falls on bare soil with a surface temperature of below 0 ^∘C, a snow cover layer is formed. The energy balance for the snow cover is written as follows (Flerchinger, 2017):

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& {\mathit{\rho }}_{\text{sp}}{c}_{\mathrm{i}}{\displaystyle \frac{\partial T}{\partial t}}+{\mathit{\rho }}_{\mathrm{l}}{L}_{\mathrm{f}}{\displaystyle \frac{\partial w_{\text{sp}}}{\partial t}}\\ & ={\displaystyle \frac{\partial }{\partial z}}\left[k_{\text{sp}}{\displaystyle \frac{\partial T}{\partial z}}\right]+{\displaystyle \frac{\partial R_{\mathrm{n}}}{\partial z}}-L_{\mathrm{s}}\left[{\displaystyle \frac{\partial q_{\mathrm{v}}}{\partial z}}+{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{v}}}{\partial t}}\right],\end{array}\end{array}$$

where ρ_sp is the density of snow [M L⁻³], w_sp is the liquid water content in the snow [–], c_i is the heat capacity of ice [L² T⁻² Θ⁻¹], L_s is the latent heat of sublimation [L² T⁻²], k_sp is the thermal conductivity of snow [M L² T⁻³ Θ⁻¹], and R_n is the net downward radiation flux within the snow [M L² T⁻⁴]. If the air temperature is increased to above 0 ^∘C, snow melted. Snowmelt in excess of the calculated interception could infiltrate into the subsurface medium when the maximum depth of ponding is set to be a positive value.

2.4 Model inputs

Based on the conceptual model shown in Fig. 1, we use SHAW to couple observed soil water and groundwater dynamics and wintertime atmospheric conditions during the freezing–thawing cycle from 28 November 2015 to 1 April 2016. The time step of numerical simulation is 1 h. Because the deepest water table depth is 143 cm and soil temperature at 150 cm below surface is available, the length of the 1D soil column is set to be 150 cm. The model domain is uniformly divided into 30 layers, which results in 31 nodes. The stable variables of the model include soil temperature and soil water content. The top node has an atmospheric boundary condition, while the bottom node has a zero-flux and specified-temperature boundary condition. The occurrence of lateral groundwater flow is realized by specifying a hydraulic gradient at a node within the saturated zone (the 30th node in the current study). In this way, the fluctuating groundwater level is an outcome of wintertime atmospheric conditions and lateral groundwater inflow, while the fluctuating soil water content is controlled by both wintertime atmospheric conditions and groundwater dynamics.

Figure 3 The maximum daily air temperature, the minimum daily air temperature, dew-point temperature, and precipitation during the freezing–thawing cycle from 28 November 2015 to 1 April 2016.

The inputs for the atmospheric boundary conditions include maximum air temperature, minimum air temperature, wind speed, precipitation, solar radiation, and dew-point temperature. The first four parameters are available from the observations at the Otak meteorological stations, the parameter of solar radiation is calculated by the Angström–Prescott equation (Yorukoglu and Celik, 2006), based on the sunshine duration, and the dew-point temperature is calculated by the Hyland–Wexter equation (Wexler et al., 1983), based on relative humidity, maximum air temperature, and minimum air temperature. Figure 3 shows the maximum and minimum air temperatures, dew-point temperature, and precipitation. The inputs for the lower boundary condition at the bottom include a zero hydraulic conductivity to obtain a no-flux boundary condition and the daily averaged soil temperature measured at 150 cm below surface for the specified-temperature boundary condition. The initial conditions of soil water content and soil temperature for spin-up are determined based on measured values on 29 October 2015, and the model spin-up is run for 30 d before the start of the freezing stage (28 November 2015).

Soil texture determines the thermal parameters and the initial hydraulic parameters. As shown in Table 1, the percentage of clay ranges between 5 % and 8 % in the layer between 70 and 100 cm and between 1 % and 4 % in other depths. Therefore, the model domain is divided into three layers. Based on average contents of clay, silt, and sand in each layer, the thermal parameters are calculated using Eqs. (2) and (3). The hydraulic parameters (θ_r, θ_s, α, n), which are initially estimated by the Rosetta pedotransfer function (Schaap and Leij, 1998; Zhang and Schaap et al., 2017), are further calibrated to fit the measured liquid soil water content (Table 1). Based on the possible range of lateral inflow rate varying between 0.96 and 1.16 mm d⁻¹ (Jiang et al., 2017), the lateral inflow rate is calibrated to be 1.03 mm d⁻¹, which leads to the lowest RMSE of water table depth.

Table 2 The different scenarios of sensitivity analysis with different initial water table depths and different conditions of groundwater supply.

* LFR represents the lateral flow rate.

To examine the role of lateral groundwater inflow and initial water table depth in freezing-induced groundwater migration, we build different scenarios with three different initial water table depths (120, 170, and 220 cm) and three different lateral flow rates (0, 0.51, and 1.03 mm d⁻¹) for sensitivity analysis. We also build one scenario with an initial water table depth of 250 cm but without lateral groundwater inflow and three scenarios with water table depths fixed at 120, 170, and 220 cm. The specific settings of the 13 scenarios are shown in Table 2. To ensure that the water table can be explicitly represented in the models, the length of model domain is extended to 250 cm. To exclude the influence of heterogeneity, the whole soil column is assumed to be filled with sand of the same soil texture as that of 0 to 70 cm in the field site. Because there is no temperature measurement at 250 cm below surface, the temperature at the lower boundary is estimated by the force-restore approach embedded in the SHAW model (Hirota, 2002), which can be written as

$$\begin{array}{}\text{(9)}& \left(\mathrm{1}+{\displaystyle \frac{\mathrm{2}z}{d_{\mathrm{d}}}}\right){\displaystyle \frac{\partial T}{\partial t}}={\displaystyle \frac{\mathrm{2}}{C_{\mathrm{s}}{d}_{\mathrm{d}}}}G-\mathit{\omega }\left(T-T_{\text{AVG}}\right),\end{array}$$

where z is the depth [L] below the surface, ω is the frequency [T⁻¹] of fluctuation period, d_d is damping depth [L] corresponding to ω, which is expressed as $d_{\mathrm{d}}={\left(\frac{\mathrm{2}{k}_{\mathrm{T}}}{C_{\mathrm{s}}\mathit{\omega }}\right)}^{\mathrm{1}/\mathrm{2}}$, and T_AVG is the average annual air temperature [Θ].

Figure 4 The observed hourly liquid water content and soil temperatures during the freezing–thawing cycle from 28 November 2015 to 1 April 2016.

3 Results and discussion

3.1 Field data showing freezing-induced water migration and water-level decline

The 0 ^∘C isothermal contour of soil temperature was shown in Fig. 2, which characterizes the gradual deepening of frost depth during the freezing stage from 28 November 2015 to 29 February 2016 well. Nevertheless, soil temperature alone is not able to demonstrate freezing-induced water migration. Figure 4 shows the evolutions of liquid water content and soil temperature at six depths (10, 20, 30, 50, 70, and 90 cm). The dates of the start of freezing, i.e., when soil temperature drops to 0 ^∘C, at 10, 20, 30, 50, and 70 cm, are 2 December, 15 December, 16 December, 9 January, and 23 January, respectively (Fig. 4). At 10, 20, and 30 cm below the surface, there is a sharp decrease in liquid water content when the soil temperature drops to 0 ^∘C, while at 50 and 70 cm below the surface, due to freezing-induced water migration toward the freezing front, there is a smooth decrease in liquid water content before reaching 0 ^∘C. At 90 cm, probably because the freezing front is close enough to the sensor, although the temperature at the sensor does not reach 0 ^∘C, there is a significant decrease in liquid water content on 29 February due to freezing-induced water redistribution. Therefore, the relationship between temperature and liquid water content recorded by 5TM sensors, which are based on the frequency domain reflectometer (FDR) method, reflects the effect of freezing on liquid water content well.

Due to the rising air temperature in spring, soil temperature at 10 cm increased from below 0 to above 0 ^∘C on 1 March, which is considered to be the start of the thawing stage in the current study. Because 1 April has the highest water level as a result of the downward infiltration of thawed water, this date is considered to be the end of the thawing period. For the sensors at 20 and 70 cm, the start date of the thawing stage is 2 March, while for the sensors at 30 and 50 cm, the start dates of the thawing stage are 3 and 5 March, implying the occurrence of bidirectional thawing. Different from the quick response of soil temperature to freezing, the response of soil temperature to thawing is a slow process. For the three sensors at 10, 20, and 30 cm, the duration from start to end of thawing is 13 d, while for the two sensors at 50 and 70 cm, the durations from start to end of thawing are 18 and 26 d, respectively. The liquid water content increases sharply at the end of thawing and is much higher than that before freezing as a result of the contribution of freezing-induced water gain.

As shown in Fig. 1, water table depths on 28 November and 1 April are 115 and 90 cm, respectively. Because snowfall in winter is only 11.7 mm, the much higher water level at the end of the thawing stage is mainly due to the occurrence of lateral groundwater inflow as reported in Jiang et al. (2017). The lowest water level of 143 cm below the surface occurs on 29 January. The trend of rising water level from 29 January to 29 February, when frozen soil has not been thawed yet, is attributed to lateral groundwater inflow. The occurrence of lateral groundwater inflow also alleviated the magnitude of freezing-induced groundwater-level decline from 28 November to 29 January; i.e., freezing-induced groundwater-level decline should be greater than 28 cm if there is no lateral inflow. To quantify the water budget during the freezing–thawing cycle, we applied the SHAW model to determine the fluxes of freezing-induced migration, lateral groundwater inflow, infiltration, and soil evaporation.

Figure 5 The observed and simulated liquid water content and simulated total water content (θ_t) (a₁-a₅) and soil temperature (b₁-b₅) at the five different depths within the frozen zone during the freezing–thawing cycle. The RMSEs of liquid water content and soil temperature are also listed.

Figure 6 The evolutions of simulated soil temperature (a), ice content (b), liquid water content (LWC) (c), and total water content (TWC) (d) at the field site. Also shown in the figures are evolutions of frost depth and groundwater level.

3.2 Numerical results of water redistribution during the freezing–thawing cycle

Based on the parameters listed in Table 1, and with a lateral groundwater inflow rate of 1.03 mm d⁻¹, there are good agreements between simulated and measured liquid water content in the frozen zone (Fig. 5a), between simulated and measured soil temperature (Fig. 5b), and between simulated and measured groundwater level (Fig. 6). Figure 5a also shows the simulated total water content, which coincides with liquid water content before freezing but increases rapidly after freezing and finally staying at a stable value. The rapid increase in total water content is due to freezing-induced water migration from the underlying soil.

Figure 6a shows the distribution of simulated soil temperature and frost depth. There is a good match between simulated and measured frost depth, the latter of which is determined by the dates of start of freezing and end of thawing as shown in Fig. 4. In the frozen zone, as the frost depth increases, ice content increases from the shallow to the deep, but liquid water content is generally maintained around the residual water content (Fig. 6b and c). In the node above the frost depth, the ice content is lower than the overlying node, and the liquid water content is higher than the residual water content. Note that when the frost depth reaches the middle layer with finer soils, there is a thicker zone with lower ice content and higher liquid water content, which is related to the higher residual water content of the middle layer. In the whole frozen zone, the trend of depth-increasing total water content toward the water table (Fig. 6d) indicates a trend of depth-increasing water gain controlled by the distance between the freezing front and the water table.

Figure 7 The simulated daily soil evaporation rate at the field site. Also shown are the mean evaporation rates for the three stages with different trends of daily soil evaporation rates.

3.3 Water budgets during the freezing–thawing cycle: implications for appropriate boundary conditions

During the freezing–thawing cycle, infiltration of snowmelt and lateral groundwater inflow are two processes of water gain in the soil column, while soil evaporation is the only form of water loss leaving the top of the soil column. At our field site, the cumulative infiltration of snowmelt is simulated to be 3.2 mm, while the cumulative evaporation is simulated to be 32.0 mm, both of which are restricted by the frozen soil. Based on the temporal evolution of evaporation rate (Fig. 7), we identify three stages of soil water evaporation. Stage I has a trend of decreasing evaporation rate, stage II has a limited evaporation rate, and stage III has a trend of increasing evaporation rate. Note that in stage II the atmospheric conditions associated with the snowfall event on 13 February result in a slightly higher evaporation rate. The mean evaporation rates in the three stages are 0.36, 0.18, and 0.59 mm d⁻¹, respectively. Similar three stages of temporal evolution of soil evaporation during the freezing–thawing cycle have been reported in several field studies (Kaneko et al., 2006; Wu et al., 2016; Xue et al., 2020). The high evaporation in stage III, which is a direct result of accumulation of liquid water during the thawing stage, explains the occurrence of soil salinization in spring in many regions of the world well (Liu et al., 2009; Lopez et al., 2010; Bechtold et al., 2011).

Table 3 Water budgets during the freezing–thawing cycle.

If the freezing–thawing cycle is divided into the freezing stage and the thawing stage, the balance of water in each stage can be written as

$$\begin{array}{}\text{(10a)}& {\displaystyle }{S}_{\mathrm{s}\mathrm{0}}+S_{\mathrm{u}\mathrm{0}}+q\cdot \mathrm{\Delta }{t}_{\mathrm{F}}+I_{\mathrm{F}}=S_{\mathrm{s}\mathrm{1}}+S_{\mathrm{u}\mathrm{1}}+{\mathrm{ET}}_{\mathrm{F}}\cdot \mathrm{\Delta }{t}_{\mathrm{F}}\text{(10b)}& {\displaystyle }{S}_{\mathrm{s}\mathrm{1}}+S_{\mathrm{u}\mathrm{1}}+q\cdot \mathrm{\Delta }{t}_{\mathrm{T}}+I_{\mathrm{T}}=S_{\mathrm{s}\mathrm{2}}+S_{\mathrm{u}\mathrm{2}}+{\mathrm{ET}}_{\mathrm{T}}\cdot \mathrm{\Delta }{t}_{\mathrm{T}},\end{array}$$

where S_s0 [L], S_s1 [L], and S_s2 [L] are cumulative water storages in the saturated zone at the beginning of the freezing stage, at the end of the freezing stage (the beginning of the thawing stage), and at the end of the thawing stage, respectively; S_u0 [L], S_u1 [L], and S_u2 [L] are cumulative water storages in the unsaturated zone at the beginning of the freezing stage, at the end of the freezing stage (the beginning of the thawing stage), and at the end of the thawing stage, respectively; q [L T⁻¹] is the mean rate of lateral groundwater inflow during the freezing–thawing cycle; Δt_F [T] and Δt_T [T] are the durations of the freezing and thawing stages, respectively; I_F [L] and I_T [L] are the cumulative snowmelt infiltration during the freezing and thawing stages, respectively; ET_F [L T⁻¹] and ET_T [L T⁻¹] are the mean soil evaporation rates during the freezing and the thawing stages, respectively. The values of each term in Eq. (10) are listed in Table 3.

Apparently, the change in groundwater storage, ΔS_s, during the freezing stage can be calculated by S_s1−S_s0, which equals −40.3 mm and corresponds to a water table decline of 10 cm. In fact, such a water table decline has been alleviated by the lateral groundwater inflow. The change in groundwater storage, ΔS_s, should be calculated by $S_{s\mathrm{1}}-S_{s\mathrm{0}}-q\cdot \mathrm{\Delta }{t}_{\mathrm{F}}$, which equals −136.1 mm; i.e., freezing-induced groundwater loss should be 136.1 mm. In a groundwater flow model without considering the unsaturated zone, this process can be characterized by a mean upward groundwater flux of 1.46 mm d⁻¹ with a duration of 93 d at the upper boundary. Similarly, the change in groundwater storage, ΔS_s, in the thawing stage should be calculated by $S_{\mathrm{s}\mathrm{2}}-S_{\mathrm{s}\mathrm{1}}-q\cdot \mathrm{\Delta }{t}_{\mathrm{T}}$, which equals 126.0 mm. In a groundwater flow model, this process can be characterized by a mean downward groundwater flux of 3.94 mm d⁻¹ at the upper boundary with a duration of 32 d.

As demonstrated in the pioneering studies of regional groundwater flow (Hubbert, 1940; Toth, 1962), shallow groundwater in the discharge area is dominated by upward flow due to strong evapotranspiration. The above calculation confirms that during $\mathrm{3}/\mathrm{4}$ of the freezing–thawing period, shallow groundwater is also dominated by upward flow. As demonstrated in Jiang et al. (2017), at our study site in the unfrozen period, the stage from late September 2014 to early October 2014 (with a duration of 18 d) is the only stage with net downward groundwater flux equaling 2.99 mm d⁻¹. The current study indicates that the thawing stage has an even larger downward groundwater flux and an even longer duration. Therefore, it will be interesting to examine the pattern of two-dimensional or three-dimensional groundwater flow induced by the freezing–thawing cycle using a fully coupled soil-groundwater model in the future.

Figure 8 The evolutions of simulated total water content, frost depth, and water table depth under three different initial water table depths (120, 170, and 220 cm) and four different groundwater supply conditions.

Table 4 A comparison of water table decline (ΔWTD), frost depth (FD), and water gain in the frozen zone under different scenarios.

Mean TWC is calculated by total water storage in the frozen zone (S_FZ) divided by frost depth, while ΔTWC is calculated by increased water storage in the frozen zone (ΔS_FZ) divided by frost depth.

3.4 The effect of water table depth and lateral groundwater inflow on freezing-induced water redistribution

In Sect. 3.1, we inferred that due to the occurrence of lateral groundwater inflow, freezing-induced groundwater-level decline is alleviated, and there is a trend of groundwater-level rise before the end of the freezing stage. Figure 8a shows the evolutions of total soil water content, frost depth, and groundwater level under four conditions of groundwater inflow when the initial water table depth (WTD₀) equals 120 cm. When there is no lateral groundwater inflow (Fig. 8a1), the lowest water level occurs at the end of the freezing period, corresponding to a water table decline of 45 cm. A comparison of water table fluctuations under three different rates of lateral groundwater inflow (Fig. 8a1–a3) confirmed that the water-level rise before the start of the thawing stage is caused by lateral groundwater inflow. Therefore, in the field, whether the timing of lowest water table corresponds to the end of the freezing stage can be used to infer the existence of lateral groundwater inflow–outflow. The higher water level caused by lateral groundwater inflow would be beneficial for freezing-induced water migration towards the freezing front. As lateral groundwater inflow rate changes from 0 to 0.51 to 1.03 mm d⁻¹, the increase in water storage in the frozen zone changes from 108.3 to 112.4 to 125.6 mm, while the increase in total water content in the frozen zone changes from 0.114 to 0.118 to 0.132 cm³ cm⁻³ (Table 4). As shown in Fig. 8a3, the total water content above the freezing front could be as high as 0.286 cm³ cm⁻³. We infer that the co-occurrence of shallow water table and lateral groundwater inflow would enhance the possibility of frost heave.

It is also interesting to examine the control of initial water table depth on freezing-induced groundwater migration. Figure 8b and c show the evolutions of total water content, frost depth, and groundwater level when WTD₀ equals 170 and 220 cm. If there is no lateral groundwater inflow, freezing-induced groundwater-level decline is 20 cm when WTD₀ equals 170 cm and is 0 when WTD₀ equals 220 and 250 cm. When WTD₀ equals 220 and 250 cm, the zero water table decline implies that groundwater is not directly involved in freezing-induced water redistribution, but there is still freezing-induced water migration from the unsaturated zone. When WTD₀ equals 250 cm, we find the lower soil water content below the freezing front led to an even smaller freezing-induced water migration, with the increased total water content as low as 0.010 cm³ cm⁻³ (Table 4). We also note that as WTD₀ increases, the smaller soil water content in the shallow part of the soil column also results in larger frost depths. When the lateral groundwater inflow is fixed at other rates, the increased water storage in the frozen zone also decreases with initial water table depth (Table 4).

When the water table depth is fixed at 120, 170, and 220 cm, the increase in total water content in the frozen zone is found to be 0.160, 0.083, and 0.022 cm³ cm⁻³, respectively. Compared with the cases without lateral groundwater inflow but with a dynamic water table (Scenarios A1, B1 and C1), a fixed water table leads to higher water gain in the frozen zone. This indicates that a fixed-head lower boundary condition would overestimate freezing-induced water migration as well as the risk of frost heave because such a lower boundary condition implies instantaneously replenishment of groundwater. This comparison also shows the necessity of properly characterizing water table fluctuations induced by freezing and groundwater inflow–outflow.

4 Conclusions

Based on field observations of soil temperature, liquid soil water content, and groundwater level at a site in the discharge area with shallow water table, it was found that the temporal variations of wintertime liquid water content and groundwater level are not only controlled by freezing and thawing, but also by lateral groundwater inflow. Therefore, we constructed a model to examine the responses of soil water and groundwater to wintertime climatic conditions and lateral groundwater inflow. The observed fluctuating soil water contents and groundwater level are well reproduced by the calibrated model, which increases our understanding of water balance as well as the influencing factors of freezing-induced groundwater migration.

Based on the groundwater budget during the freeing-thawing cycle, the mean flux of freezing-induced groundwater loss in the freezing stage with a duration of 3 months is 1.46 mm d⁻¹, and the mean flux of thawing-induced groundwater gain in the thawing stage with a duration of 1 month is 3.94 mm d⁻¹. As found out by Jiang et al. (2017), due to the semi-arid climate, water flux during the majority of the unfrozen period is controlled by evaporation, and the only stage with net recharge in the unfrozen period is from late September to early October. Therefore, combined with the fluxes obtained by Jiang et al. (2017) for the unfrozen period, the fluxes during the freezing–thawing cycle obtained in the current study can be useful for future studies on transient groundwater flow models to characterize the control of climatic conditions on deeper groundwater flow in semi-arid and seasonally frozen regions.

Based on a series of modeling with different initial water table depths, we found that when water table depth is maintained at 220 cm below surface, although part of water in the unsaturated zone can be migrated to the freezing front, groundwater in the saturated zone is not directly involved in freezing-induced water migration. However, when initial water table depths are shallower than 220 cm, freezing-induced groundwater migration could be significant, and the rise in groundwater level due to lateral groundwater inflow could increase the amount of groundwater migrated to the freezing front. Fan et al. (2013) reported that regions with water table depth shallower than 2 m account for around 31 % of the global land area, while Zhang et al. (2003) reported that seasonally frozen soils underlie approximately 24 % of the Northern Hemisphere exposed land surface. Moreover, wintertime water table recession has been reported in Canada (van der Camp et al., 2003; Ireson et al., 2013), China (Wu et al., 2016; Zhang et al., 2019), Russia (Vinnikov et al., 1996), Sweden (Stahli et al., 1999), and the USA (Drescher, 1955; Schneider, 1961; Daniel and Staricka, 2000), implying that the involvement of groundwater in freezing-induced water redistribution is a common physical process. Although the threshold water table depth in other study areas with different climate conditions remains unknown, our study demonstrated well the necessity of properly characterizing freezing-induced water table fluctuations to quantify freezing-induced groundwater migration and its effect on engineering problems and ecological processes in cold regions.