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

Freezing-induced water migration and groundwater level decline are widely observed in regions with shallow water table, but many existing studies trying to quantify freezing-induced groundwater migration do not account for water level fluctuations induced by freezing and thawing. Here, detailed field observations of liquid soil water content and groundwater level fluctuations at a site in the Ordos Plateau, China are combined with numerical modeling to show groundwater and soil water dynamics controlled by wintertime atmospheric conditions and topographically-driven lateral groundwater inflow. By comparing simulation results with and without lateral groundwater inflow, we find lateral groundwater inflow leads to an alleviated freezing-induced water level decline and enhanced freezing-induced water migration. At the field site with a lateral groundwater 20 inflow rate of 1.03 mm/d, compared with the case without lateral groundwater inflow, the water level decline decreases from 40 cm to 15 cm, and the increased total water content in the frozen zone enhances from 0.071 to 0.106. By calculating the budget of groundwater, the mean upward flux of freezing-induced groundwater loss is 1.46 mm/d for 93 days, and the mean flux of thawing-induced groundwater gain is as high as 3.94 mm/d for 32 days. The study enhances our understanding of the mechanisms controlling water redistribution between 25 saturated and unsaturated zones and the water budget in the freezing-thawing cycle. The fluxes of groundwater loss and gain in the freezing and thawing stages obtained in the current study can be useful for future studies on two- or three-dimensional transient groundwater flow in semi-arid regions with seasonally frozen soils. output. To obtain a clear understanding of wintertime water migration 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 thawing-induced water movement. In the current study, the field site reported in Jiang et al. (2017), which has shallow water table and lateral groundwater inflow, is used to numerically examine the interactions of soil water and groundwater and calculate the water budget during the freezing-thawing cycle. are 1 MAR, 2 MAR, 3 MAR, and 5 MAR, while for the sensor at 70 the start date of thawing stage is 2 MAR, implying the occurrence of bi-directional thawing. Different from the quick response of soil temperature to freezing, the response of soil temperature to thawing is a slow process. For the 3 sensors at 10 cm, 20 cm and 30 cm, the duration from start to end of thawing is 13 days, while for the 2 at 50 cm and 70 the durations from start to end of thawing 18 days and 26 days, respectively. The thawing of frozen soil is accompanied by a trend of increasing liquid water content. The liquid water at the end of thawing is much higher than that before freezing, indicating the contribution of freezing-induced water gain.

However, there is no numerical study on how groundwater and soil water dynamics respond to wintertime atmospheric conditions. The Wudu Lake catchment is one of the small catchments located to the southeast of the first-order water divide of the Ordos Plateau. The distribution of water table elevations in the Wudu lake catchment can be found 105 in Jiang et al. (2018) while the zone of flowing wells, i.e., with higher hydraulic head than the topography, can be found in Wang et al. (2015). There is a monitoring site of soil water content and groundwater level in the topographic lows of this catchment. The occurrence of several flowing wells near the monitoring site indicates that the site belong to the discharge area and could receive lateral groundwater flow from topographic highs. The water table depths are recorded by using water level loggers in a shallow well with a depth of 10 m while the 110 temperatures and liquid water content in the Quaternary deposits are recorded by using 5TM sensors produced by Decagon Devices at 8 depths (at around 10 cm, 20 cm, 30 cm, 50 cm, 70 cm, 90 cm, 110 cm, 150 cm). As shown in Fig. 1, water table depth increases from 115 cm on 28 NOV 2015 to 143 cm on 29 JAN 2016, which is a direct result of freezing-induced water migration from groundwater to the freezing front.
At the monitoring site, lateral groundwater flow had been found to be an indispensable component of 115 groundwater budget (Jiang et al., 2017). Fig. 2 shows the control of main hydrological and hydrogeological https://doi.org/10.5194/hess-2020-657 Preprint. Discussion started: 6 January 2021 c Author(s) 2021. CC BY 4.0 License. processes on groundwater level during the freezing-thawing cycle. The processes leading to water level decline include water migration induced by freezing and evaporation due to the dry climate, and the processes leading to water level rise is lateral groundwater flow throughout the cycle and infiltration of thawed water during the thawing stage. 120 Figure 1. The contour map showing the change of soil temperature at the study site during the freezing-thawing cycle. Also shown is the fluctuating groundwater level monitored at the study site. Figure 2 The conceptual models of hydrological processes and groundwater regime in semi-arid regions during 125 the seasonally freezing and thawing stages.

The SHAW model
The Simultaneous Heat and Water (SHAW) model (Flerchinger and Saxton, 1989) is one of the most robust models for freezing-induced water migration (Hayhoe, 1994;DeGaetano et al., 2001;Kahimba et al., 2009;Zhang et al., 2017). The contributions of liquid water, ice, and vapor contents have been considered in the water balance 130 equation, and the contributions of heat conduction, phase change, liquid flow, and vapor gas diffusion have been considered in the energy balance equation (Li et al., 2010). Because water flow in the unsaturated zone is predominantly vertical (Stephens, 1996;Romano et al., 1998;Van Dam and Feddes, 2000), the SHAW model only https://doi.org/10.5194/hess-2020-657 Preprint. Discussion started: 6 January 2021 c Author(s) 2021. CC BY 4.0 License. accounts for the effects of soil freezing and thawing on water redistribution within a one-dimensional (1D) soil column. Moreover, by adding a source term in the water balance equation, it has the ability to account for the 135 contribution of lateral groundwater flow to water level 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 By considering ice content change and vapor flux, the water balance equation is written as where θl is liquid water content [-], K is hydraulic conductivity [L T -1 ], ψ is soil matric potential [L], and U is a source/sink term for water flux [T -1 ]. Soil water characteristic curves, which describe the relationship between soil 150 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 where θs and θr are saturated and residual water content [-], α [L -1 ], n [-] and m [-] are empirical parameters. α equals the inverse of the air-entry value, n is a pore-size distribution index and m=1-1/n. In the frozen zone, due to 155 the co-existence of liquid water and ice, matric potential is strongly dependent on temperature, and the matric potential is obtained by the following equation (Fuchs et al., 1978): where g is the acceleration of gravity [L T -2 ]. Eq. (4) indicates that as the negative temperature increases, the soil suction also increases (matric potential becomes more negative).
where Ks is the saturated hydraulic conductivity [L T -1 ], Se [-] is the effective saturation calculated by (θl-θr)/(θs-θr), and l [-] is a pore-connectivity parameter which is assumed to be 0.5 in the original work of (Mualem, 1976). Note that in the saturated unfrozen zone, K equals Ks, but in the unsaturated unfrozen zone,

165
K equals unsaturated hydraulic conductivity (Ku). In the frozen zone, due to occupation of ice in the pores, the unsaturated hydraulic conductivity computed from Eq. (5) is reduced linearly with ice content and is assumed to be zero when the porosity is decreased to 0.13 (Flerchinger, 2000).

Model inputs
In the following, based on the conceptual model shown in Figure 2, we use SHAW to couple soil water and 170 groundwater dynamics in a 1D soil column during the freezing-thawing cycle from 28 NOV 2015 to 1 APR 2016, as such to understand how the coupled unsaturated-saturated processes respond to wintertime atmospheric conditions. To ensure that the lateral groundwater inflow occurs within the saturated zone, the length of the 1D soil column is set to be 155 cm. The model domain is uniformly divided into 31 layers, which results in 32 nodes.
The top node has an atmospheric boundary condition, while the bottom node has a zero-flux and 175 specified-temperature. Following Flerchinger (2000), the occurrence of lateral groundwater flow is realized by setting a specified hydraulic gradient at the 31 st node. The fluctuation of groundwater level is an outcome of wintertime atmospheric conditions and lateral groundwater flow, while the fluctuations of soil water content are controlled by both wintertime atmospheric conditions and groundwater dynamics.
The inputs for the atmospheric boundary condition include maximum air temperature, minimum air 180 temperature, wind speed, precipitation, solar radiation and dew-point temperature. The first four parameters 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 (data available from China Meteorological Data Service Center), and the dew-point temperature is calculated by the Hyland-Wexter equation (Wexler et al., 1983) based on relative humidity, maximum air temperature and 185 minimum air temperature.

195
Although soil samples from all depths belong to sand or loamy sand, soil texture differs in the profile. 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 (Zhao et al., 2020). Therefore, the model domain is divided into three layers. The parameters for the van Genuchten (1980) and Mualem (1976)

Field data showing freezing-induced water migration and water level decline
During the freezing stage from 28 NOV 2015 to 29 FEB 2016, there is a general trend of deepening frost depth (Fig. 1). The dates of the start of freezing, i.e., when soil temperature drops to 0 ℃, at 10 cm, 20 cm, 30 cm, 50 cm, and 70 cm are 02 DEC, 15 DEC, 16 DEC, 9 JAN and 23 JAN, respectively ( Fig. 4). At 10 cm, 20 cm and 220 30 cm below surface, there is a sharp decrease in liquid water content when the soil temperature drops to 0 ℃, while at 50 cm and 70 cm below surface, due to freezing-induced water migration toward the freezing front, there is a gradual decrease in liquid water content before reaching a temperature of 0 ℃. At 90 cm below surface, probably because the freezing front is close enough to the sensor, although the temperature at the sensor does not reach 0 ℃, there is a significant decrease in liquid water content on 29 FEB due to freezing-induced water 225 migration. Therefore, the relationship between temperature and liquid water content recorded by 5TM sensors, which are based on the frequency domain reflectometer (FDR) method, well reflects the effect of freezing on liquid water content.
Due to the rising air temperature in spring, soil temperature at 10 cm increased from below 0 ℃ to above 0 ℃ on 1 MAR, which is considered to be the start of the thawing stage in the current study. Because the water 230 level is highest on 1 APR as a result of the downward infiltration of thawed water, 1 APR is considered to be the end of the thawing period. For the four sensors at 10 cm, 20 cm, 30 cm and 50 cm, the start dates of thawing https://doi.org/10.5194/hess-2020-657 Preprint. Discussion started: 6 January 2021 c Author(s) 2021. CC BY 4.0 License.
stage are 1 MAR, 2 MAR, 3 MAR, and 5 MAR, while for the sensor at 70 cm, the start date of thawing stage is 2 MAR, implying the occurrence of bi-directional thawing. Different from the quick response of soil temperature to freezing, the response of soil temperature to thawing is a slow process. For the 3 sensors at 10 cm, 20 cm and 30 235 cm, the duration from start to end of thawing is 13 days, while for the 2 sensors at 50 cm and 70 cm, the durations from start to end of thawing are 18 days and 26 days, respectively. The thawing of frozen soil is accompanied by a trend of increasing liquid water content. The liquid water content at the end of thawing is much higher than that before freezing, indicating the contribution of freezing-induced water gain.
240 Figure 4 The observed hourly liquid water content and soil temperatures during the freezing-thawing cycle from 28 NOV 2015 to 1 APR 2016.
As shown in Fig. 1, water table depths on 28 NOV and 1 APR are 115 cm and 90 cm, respectively. Because snowfall in winter is limited, the much higher water level at the end of thawing stage is mainly due to the 245 occurrence of lateral groundwater inflow as reported in Jiang et al. (2017). The lowest water level equaling 143 cm below surface occurs on 29 JAN. From 29 JAN, there is a trend of rising water level, which is one month earlier than the occurrence of thawed water infiltration during the thawing stage. This phenomenon can also be attributed to lateral groundwater inflow. The occurrence of lateral groundwater inflow also alleviated the https://doi.org/10.5194/hess-2020-657 Preprint. Discussion started: 6 January 2021 c Author(s) 2021. CC BY 4.0 License.
magnitude of freezing-induced groundwater level decline from 28 NOV to 29 JAN, i.e., freezing-induced 250 groundwater level decline should be greater than 28 cm. To quantify the water budget during the freezing/thawing cycle, the abovementioned model is constructed to determine the fluxes of freezing-induced migration, lateral groundwater flow, infiltration and soil evaporation.

Freezing-induced water redistribution: the effect of lateral groundwater flow and soil heterogeneity
Based on the parameters of van Genuchten (1980) and Mualem (1976) Table 1, and the 255 lateral flow rate of 1.03 mm/d, there are good fits between simulated and measured liquid water content in the frozen zone, and also for soil temperature (Fig. 5) and water level (Fig. 6a). Fig. 5 also shows the simulated total water content, which coincides with liquid water content before freezing, but increases rapidly due to freezing-induced water migration from the underlying soil and finally maintains at a stable value corresponding to the completely frozen state. In the frozen zone, the magnitude of increase in total water content increases from 260 the shallow to the deep ( Fig.5 and 6a), indicating that the distance to the water table not only controls the initial liquid water content, but also the amount of freezing-induced water gain from underlying soils.

270
Soil texture and the corresponding hydraulic conductivity are also factors controlling the flux of freezing-induced water migration through the unfrozen zone. Because the soil texture in the layer from 70 cm to 100 cm is slightly different from overlying and underlying soils, we compare the distribution of total water content under the homogeneous case (Scenario B, Fig. 6b) with that of the base case (Scenario A, Fig. 6a). We 275 find the low Ks of the middle layer restricts water migration into the part of frozen zone overlying the middle layer. As shown in Table 3, the increase of total water content in the frozen zone is 0.126 for Scenario B (corresponding to an increased water storage of 119.7 mm), and is 0.106 for Scenario A (corresponding to an increased water storage of 96.1 mm). Scenario A also has a smaller magnitude of freezing-induced water table decline (Table 3). Moreover, due to the ability of clay to retain more water, before late January when the frost 280 depth is still away from the middle layer, the middle layer shown in Fig. 6a has more liquid water than that shown in Fig. 6b; from late January to late February when the freezing front in Fig. 6a is shallower than the bottom of the middle layer, the unfrozen part of the middle layer also has higher water content than the underlying soil.  ∆SFZ represents increased water storage in the frozen zone, ∆TWC represents increased total water content in the frozen zone, and Mean TWC represents the mean total water content in the frozen zone at the state of completely frozen.
Lateral groundwater flow is a component of water budget influencing water level. When there is no lateral 290 groundwater flow, in either the heterogeneous case (Scenario C, Fig. 6c) or the homogeneous case (Scenario D, Fig. 6d), the lowest water level occurs at the end of the freezing stage. By comparing the simulated groundwater level under a series of lateral groundwater inflow rates, we find the occurrence of lateral groundwater inflow alleviates the magnitude of freezing-induced water level decline, and leads to a trend of earlier water level rise before the start of the thawing stage (Fig. 7). Therefore, in the field, the existence of lateral groundwater inflow 295 can be determined by the earlier occurrence of groundwater level rise before the start of the thawing stage. Due to the occurrence of lateral groundwater inflow, the rise in water level from late January to late February also restricts the frost depth to some extent (Table 3). Figure 7 The temporal variations of simulated water table depth under different lateral groundwater inflow rates.
The higher water level caused by lateral groundwater flow would be beneficial for freezing-induced water migration towards the freezing front. In the heterogeneous cases, the increase in total water content in the frozen zone is found to be 0.106 in Scenario A, but is only 0.071 in Scenario C without lateral groundwater flow; in the homogeneous cases, the increase in total water content in the frozen zone is found to be 0.126 in Scenario B, but is only 0.094 in Scenario D without lateral groundwater flow (Table 3). Therefore, compared with scenarios 305 without lateral groundwater inflow, the occurrence of lateral groundwater inflow can inevitably lead to an increased freezing-induced water accumulation.
If a fixed head equaling 1.2 m below surface is assigned at the lower boundary, the increase in total water content in the frozen zone is 0.127 for the heterogeneous case, and is 0.150 for the homogeneous case (Table 3).
This confirms that a fixed head lower boundary condition would inevitably overestimate freezing-induced water 310 migration. As we pointed out in the Introduction part, a fixed head lower boundary condition implies instantaneously replenishment of groundwater. The results shown in Fig. 7 indicate that an even larger groundwater inflow is needed to maintain a stable water level.

Soil evaporation during the freezing-thawing cycle
By comparing model results with and without consideration of the three snowfall events, we find snowfall 315 did not infiltrate into the soil column due to the low permeability of frozen soil and got evaporated directly, implying that soil water evaporation is the only form of water loss leaving the top of the soil column during the freezing-thawing cycle. Therefore, we use the evaporation rate without considering snowfall events to represent actual soil evaporation. Note that although snowfall is not considered, the atmospheric conditions associated with the snowfall event on 13 FEB still result in higher evaporation rate (Fig. 8).

320
Based on the trend of evaporation rate versus time, we identify three stages of soil water evaporation. Stage Ⅰ has a trend of decreasing evaporation rate, stage Ⅱ has a limited evaporation rate, and stage Ⅲ has a trend of increasing evaporation rate. In both scenarios, the start of stage Ⅱ corresponds to a frost depth of around 20 cm.
In Scenario A, the mean evaporation rates in the three stages are 0.36, 0.18 and 0.59 mm/d, respectively (Fig. 8a), leading to a cumulative evaporation of 32.1 mm; in Scenario B, the mean evaporation rates in the three stages are 325 0.53, 0.21 and 1.55 mm/d (Fig. 8b), respectively, corresponding to a cumulative evaporation of 55.3 mm. The high evaporation rate in stage Ⅲ , which is a direct result of accumulation of liquid water during the thawing stage, well explains the occurrence of soil salinization in spring (Liu et al., 2009;Lopez et al., 2010;Bechtold et al., 2011). Similar three stages of temporal evolution of soil evaporation during the freezing-thawing cycle have https://doi.org/10.5194/hess-2020-657 Preprint. Discussion started: 6 January 2021 c Author(s) 2021. CC BY 4.0 License. been reported in several field studies (Kaneko et al., 2006;Wu et al., 2016;Xue et al., 2020).

Water budgets during the freezing/thawing cycle
Because almost all snowfall are evaporated without entering the subsurface, soil evaporation and lateral groundwater inflow are the only two processes influencing the balance of water in the soil column. 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  Table 4. 355 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 days at the upper boundary. This is consistent with the dominance of upward flux at the upper boundary in the unfrozen period at the study site (Jiang et al., 2017). Therefore, the current study confirms that shallow groundwater in the discharge area of a semi-arid catchment is dominated by upward flow, which is consistent with the direction of regional 365 groundwater flow in the discharge area (Hubbert, 1940;Toth, 1962).
Similarly, in the thawing stage, groundwater gain from downward infiltration should be calculated by

4 Conclusions
Based on field observations of soil temperature, liquid soil water content and groundwater level in the discharge area of a catchment, it was found that the observed 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, a model is built to examine the responses of soil water and groundwater to wintertime climatic 380 conditions and occurrence of lateral groundwater inflow in the saturated part of the soil column. The calibrated model increases our understanding of water balance as well as freezing-induced water migration during the freezing-thawing cycle.
Our model results showed that almost all snowfall is evaporated without infiltrating into the soil column due to the low permeability of frozen soil. Therefore, soil evaporation is the only form of water loss, while lateral 385 groundwater inflow is the only source of water gain. By comparing models with and without lateral groundwater inflow, we found the occurrence of lateral groundwater inflow led to much smaller magnitude of freezing-induced water level decline, and a trend of rising water level before the start of the thawing stage.
Consequently, compared with the case without lateral groundwater inflow, the higher water level due to lateral groundwater inflow resulted in much larger freezing-induced water migration into the frozen zone. We also found 390 the middle layer with lower saturated hydraulic conductivity restricts both evaporation rates and freezing-induced water migration.
Based on the budget of groundwater during the freeing-thawing cycle, the mean flux of freezing-induced groundwater loss in the freezing stage with a duration of three months is 1.46 mm/d, and the mean flux of thawing-induced groundwater gain in the thawing stage with a duration of one month is found to be 3.94 mm/d.

395
As found by Jiang et al. (2017), due to the semi-arid climate, 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 freeing-thawing cycle obtained in the current study can be useful for future transient groundwater flow models to characterize the control of climatic conditions on groundwater flow in semi-arid regions with seasonally freezing and thawing processes.