The Qinghai Lake Basin, situated between 36°15′-38°20′N latitude and 97°50′-101°20′E longitude (see Figure 1), is located on the northeastern Qinghai-Tibet Plateau (QTP). This basin spans an area of 29,661 km², with elevations ranging from 3,265 to 5,280 m. Qinghai Lake, the largest inland saline lake in China, experiences low temperatures throughout the year. The freezing period begins from September to October and lasts until April to May of the following year. Meteorological data from stations in Gangcha, Haiyan, and Tianjun counties indicate that the mean annual temperature in the basin ranges between −1°C and 1°C (Zhang and Duan, 2021; Zhang, 2021). Precipitation is unevenly distributed, with 90–95% occurring from May to October, while winter snowfall is minimal, contributing less than 4% of the annual total. Observational data show that precipitation has fluctuated over the decades, increasing by 9 mm per decade from the 1950s to the early 1980s, then decreasing in the late 1980s, and subsequently rising by 25 mm per decade since the 1990s. The average temperature has risen significantly since the 1950s, at a rate of 0.244°C per decade (Zhang and Duan, 2021). Mountains cover over 68% of the basin, with the Buha and Shaliu Rivers, the two largest rivers, contributing more than 60% of the basin’s total runoff (Hu et al., 2022).

The meteorological data used in our study were derived from the China Meteorological Forcing Dataset (Yang et al., 2010; Kun et al., 2019; He et al., 2020). Initially, we created a one-year time series of air temperature by averaging daily temperatures from 1979 to 2000, which was used in model initialization to form a partly frozen soil as the initial condition. Subsequently, a typical-year time series for air temperature, relative humidity, wind speed, incoming shortwave radiation, and rain/snow data was generated by averaging values over 19 years (2001–2019; Supplementary Figure S1). Precipitation was classified as snow if the average air temperature was below 0°C and as rain if above 0°C. The typical-year dataset, which effectively represents seasonal variations in meteorological conditions, was employed during the spin-up process to establish a steady state for surface-subsurface hydro-thermal dynamics. Post spin-up, daily data on air temperature, relative humidity, wind speed, incoming shortwave radiation, and rain/snow from 2017 to 2019 (see Supplementary Figure S1) were used for simulations.

Simulated groundwater levels and soil temperatures at selected depths were compared with observation data (Figure 2) to evaluate model’s accuracy in representing the real system. Overall, the simulated water table reasonably matched the observations, with a root-mean-squared-error (RMSE) of 0.51 m. The simulated soil temperatures captured the major trends of the observed temperatures at different depths, with an average RMSE of 2.25°C and a Nash-Sutcliffe Efficiency (NSE) of 0.88.

In our study, we define the frozen state of SFG using “ice saturation > 0” instead of the zero-degree isotherm, which better characterizes the physical phase transition of water during the freezing and thawing process. Soil ice saturation at six selected dates were shown in Figure 3, illustrating the seasonal freezing–thawing process throughout a year. The surface soil begins to freeze downward in early November, leading to a notable increase in ice saturation in the shallow soil by mid-November (Figure 3F). In late March, the ground is frozen to a maximum depth of approximately 2 meters (Figure 3A). During the warm season, the frozen soil begins to thaw from the surface (Figure 3B) with rapid melting of ice in the upper soil and gradual reduction of ice saturation in the lower frozen soil (Figure 3C). By early June, only minimal ice content with lower saturation remains in the downhill lower soil, and shortly after, the SFG is completely thawed (Figures 3D,E).

Meteorological factors, particularly air temperature, significantly impact the freeze–thaw process of frozen ground by influencing the exchange of energy and mass between the atmosphere and the surface ground. The temporal variation of the ice saturation observed near the surface in the base case negatively correlates with the changes in air temperature (Figure 4C). As the air temperature drops below zero in November, soil water begins to freeze into ice, increasing ice saturation. This freezing process continues throughout winter, with ice saturation spreading downward. When the temperature rises above zero in late March, the soil ice starts to melt and eventually disappears completely in June.

Changes in annual temperature also influences the freezing depth and duration of SFG, as illustrated in Figure 4. The colder case demonstrates a deeper freezing depth (approximately 2.35 meters) and a longer duration (from late October to mid-July). In contrast, the warmer case shows a shallower freezing depth (around 1.7 meters) and a shorter duration (from early November to mid-May). The impact of these variations in freezing depth and duration on hydrological processes will be discussed in the following section.

The influence of SFG on rainfall infiltration and groundwater flow is evaluated by comparing the simulated liquid saturation at selected depths and groundwater discharge among four scenarios: the base case, colder case, warmer case and the No-SFG case. In the No-SFG case, fluctuations in liquid saturation correspond to rainfall in both shallow and deep soil (dashed lines in Figure 5). In contrast, in the base case, variations in liquid saturation at different depths are closely related to the freeze–thaw process of SFG. The sudden increase in liquid saturation observed in April is due to the phase transition from ice to liquid as the shallow soil thawed. Consequently, rainwater infiltrates the shallow soil, causing liquid saturation to vary according to the rainfall trend from May (blue lines in Figures 5A,B). However, remaining SFG in the deep soil hinders the infiltration process, preventing rainwater from infiltrating the deeper soil until the deep SFG thaws in June, allowing rainfall to start recharging groundwater (blue lines in Figures 5C,D). These effects vary for colder and warmer cases due to different freezing depths and durations of SFG. In the warmer case, deep soil SFG thaws in May, allowing rainwater to infiltrate and recharge groundwater earlier. In the colder case, SFG has a stronger and more prolonged influence on rainfall infiltration due to longer freezing duration in deep soil, delaying rainwater infiltration and groundwater recharge until July (Figures 5C,D).

The groundwater discharges measured at the downgradient outlet for the four cases are shown in Figure 6. In the No-SFG case, groundwater discharge variations closely follow rainfall patterns, with peaks and valleys corresponding to rainfall events. In the base case, residual SFG in deep soil initially obstructs rainwater infiltration, delaying groundwater recharge. Consequently, the increase in groundwater discharge occurs later than the early rainfall events and gradually reaches the first discharge peak as deep soil thaws in early June, approximately 2 months later than in the No-SFG case. The presence of SFG also causes rainwater to accumulate above the frozen layer, leading to a higher peak discharge once the SFG completely thaws. Additionally, as soil water freezes into ice from early November, reducing infiltration and subsequent groundwater recharge, groundwater discharge declines more rapidly. This decline is more pronounced compared to the No-SFG case. The arrival time of the first groundwater discharge peak and its magnitude also vary for the colder and warmer cases due to differences in SFG duration. In the colder case, the peak is delayed until late June and has the largest discharge among all cases, while in the warmer case, the peak occurs earlier in early May with a smaller discharge.

We simulated a steep case and a flat case to examine the effects on SFG distribution and groundwater flow. During warm periods, gravitational forces drive water flow from the uphill to the downhill, discharging into rivers via outlets on the downhill side. As a result, liquid saturation on the uphill decreases with increasing slope (Figure 7A). This also leads to decreased ice saturation on the uphill in the frozen period due to phase changes (Figure 7B). Consequently, the uphill side of the hillslope exhibits a decrease in liquid saturation in summer and ice saturation in winter with increasing slope, along with an increase in gas saturation (Figure 7C). Conversely, there is almost no difference in three-phase saturation on the downhill side as shown by the overlapping curves among three cases due to more significant flow exchange controlled by outlet boundary. The difference in ice saturation caused by gravity flow decreases along the hillslope. The greatest difference observed on the uphill side, gradually diminishing downgradient, and becoming negligible at the downhill side (Figure 8).

The steep hillslope results in less liquid saturation in the summer and ice saturation in the winter due to gravitational water flow. This leads to lower heat capacity in the steep case compared to the base case and the flat case. This effect is most obvious on the uphill slope and disappears gradually on the downhill slope. Lower heat capacity allows the freezing process to penetrate deeper in the steep case, resulting in a 5 cm increase in ice saturation distribution depth on average (Figures 8A–D), whereas the flat case shows a 5 cm decrease in ice saturation distribution (Figures 8G–J).

The simulated groundwater discharge and evapotranspiration for three different slope cases are shown in Figure 9. Variations in topographic slopes significantly influence evapotranspiration, which in turn affects groundwater discharge. The least evapotranspiration and the highest groundwater discharge are obtained for steep case, whereas the flat slope case shows the highest evapotranspiration and the lowest groundwater discharge. The variation in slope significantly influences both lateral and vertical soil water flow. As slope increases, increased hydraulic gradient accelerate groundwater movement downhill, reducing groundwater storage and leading to a decline in groundwater levels in uphill (Supplementary Figure S3A). During the later thawing period in steep case, the moisture gradient above and below the remaining ice-containing soil intensifies, enhancing vertical rainwater infiltration (Supplementary Figure S3B). Despite the deeper ice distribution in the steep case, this enhanced infiltration reduces its overall impact, resulting in only a slight delay in groundwater flow. Peak groundwater discharge occurs 3 days later in steep case compared to flat case.

SFG significantly influences hydrological processes by affecting rainfall infiltration and groundwater discharge in QLB. In winter, as the soil freezes downward, infiltration gradually decreases, leading to a reduction in groundwater discharge (Figure 10A). When rainfall begins, the shallow soil thaws, but rainwater infiltration is hindered by the remaining SFG in the deeper soil layers, obstructing groundwater recharge during the early rainy season (Figure 10B). Once the soil completely thaws, rainwater recharges the groundwater, resulting in increased groundwater discharge (Figure 10C). Our results align with previous studies, showing that the dynamic changes in groundwater level and discharge due to freeze–thaw processes are widespread in SFG regions (Zheng et al., 2023; Lyu et al., 2023). These changes are primarily driven by variations in the recharge from the unsaturated zone to the groundwater aquifer influenced by the freeze–thaw cycle.

The hydrological processes show unique characteristics compared to previous studies in Arctic and subarctic regions, where snowfall dominates annual precipitation (Evans et al., 2018; Appels et al., 2018; Wu et al., 2016). These differences arise from variations in the distribution of rain and snow, as well as snowpack conditions, driven by different precipitation patterns. Snow, with its low thermal conductivity, acts as thermal insulation, significantly affecting the ground energy balance. Typically, snow cover thickness inversely correlates with the depth and duration of frozen ground; thus, a thick snowpack can reduce or even prevent the development of frozen ground, even when air temperatures are below freezing (Iwata et al., 2018; Wang and Chen, 2024; Iwata et al., 2011). Consequently, regions with substantial snowpack often report minor or negligible hydrological effects from SFG (Rush et al., 2021; Rey et al., 2021). In contrast, in the QLB, summer rainfall dominates precipitation, resulting in a negligible snowpack. As a result, the freeze–thaw process of SFG is primarily controlled by air temperature, without snow influence.

Moreover, in regions with abundant snowfall, the main influence of SFG is during the spring snowmelt season, where a large amount of meltwater contributes to surface runoff, with less affecting groundwater discharge (Chang et al., 2018; Evans et al., 2018). Our study emphasizes the impact of SFG on groundwater recharge and discharge in the summer rainfall-dominant regions like QLB. In the early rainfall period, while the upper soil thaws, SFG in the deeper soil obstructs rainwater infiltration, delaying groundwater recharge and causing a two-month delay in the peak of groundwater discharge.

SFG exhibits a strong response to climate warming, including a reduction in maximum freezing depth and an extension of the thawing duration (Song et al., 2022; Peng et al., 2017). A primary influence of air temperature on the SFG freeze–thaw process was reported while snow cover partially affects this process (Wang and Chen, 2024; Peng et al., 2017). The freeze–thaw process of SFG in QLB is dominated by air temperature, as the influence of snow cover is negligible. In warmer case, the maximum SFG depth decreases by 0.65 meters, and the freezing period shortens by about 2 months compared to the colder case. Consequently, the impact of SFG on groundwater recharge and discharge diminishes significantly with a warming climate. In the colder case, the groundwater discharge peak is delayed by approximately 3 months, but this delay is reduced to 1 month in warmer case due to less overlap between the deep freeze period and the rainfall season. In the warmer case, rainfall begins in mid-April, and SFG thaws completely by mid-May, allowing minimal time for SFG to hinder rainfall infiltration. Thus, SFG has a limited effect on the groundwater recharge and discharge processes. Under predicted future climate warming, the freezing depth and period of SFG will further decrease (Zou et al., 2017; Song et al., 2022; Ji et al., 2020). Once the deeper frozen soil melts completely before the rainy season, the differences in rainfall infiltration and groundwater recharge and discharge processes can be essentially ignored.

Our results suggest that SFG can delay and reduce groundwater discharge, with significant variations observed under different climate scenarios and influenced by topographic slope. These results have critical implications for water resource management in regions such as QLB. The delay in peak groundwater discharge and the reduction in winter discharge significantly impact the spatial and temporal distribution of groundwater outflow. As a vital source of recharge for surface water bodies such as rivers and lakes, seasonal variations in groundwater discharge can affect surface water flow and volume (Rosenberry et al., 2015; Sebestyen et al., 2008; Abdelhalim et al., 2020). Therefore, water resource management must consider the effects of freeze–thaw processes and reassess the availability and allocation strategies. This may influence the operational planning of reservoirs, hydropower plants, and the scheduling of water use for agricultural irrigation.

Simulations with changed soil permeability and thermal conductivity were set to access the sensitivity of model output to the input data. For soil permeability, as soil permeability increases, the timing of groundwater discharge growth during the thawing period shifts earlier and becomes more gradual (Supplementary Figure S4A). For instance, in the k1 case, groundwater discharge starts to increase sharply in June, with the first peak occurring in mid-June. As permeability reaches k7, noticeable discharge growth begins in May, with the first peak advancing to early June, approximately half a month earlier. Starting from k8, this trend is more pronounced, with discharge growth initiating in late March, and the initial peak is disappeared on the groundwater discharge curve. This phenomenon arises because increased soil permeability enhances water flow, gradually reducing the obstructive effect of diminishing ground ice during the thawing period. In high-permeability soils, water can move more readily through ice-containing layers, recharging the groundwater aquifer more quickly and at an earlier stage, thereby advancing the increase of groundwater discharge.

Soil thermal conductivity significantly impacts the freeze–thaw processes. Higher thermal conductivity enables soil to freeze more deeply in winter and accelerates the thawing process in warmer months. As thermal conductivity increases, the timing of soil thawing advances, leading to an earlier increase of groundwater discharge (Supplementary Figure S4B). For instance, the timing of the first groundwater discharge peak is advanced about 1 week from the low thermal conductivity case (th1) the high thermal conductivity case (th10).

The simulation results demonstrate that air temperature and slope variations significantly influence freeze–thaw dynamics and groundwater discharge. Due to the limited range of temperature and slope values derived from available datasets, the input data may reflect averaged local conditions, which can introduce uncertainty into the model. This averaging effect may mask finer variations and limit the model’s accuracy in simulating local-scale processes. Additionally, soil permeability and thermal conductivity were identified as important parameters impacting outcomes. Given the assumption of two homogeneous soil layers in this study, accounting for soil heterogeneity in future work could further refine the model and improve its predictive accuracy.

Additionally, a 2D hillslope-scale model was employed. While the model captures freeze–thaw cycles and groundwater flow dynamics on a two-dimensional slope, the results cannot be directly extrapolated to a catchment scale due to the simplified geometry and limited hydrological processes. A 3D model would more effectively account for complex topography and better depict the intricate water exchanges between groundwater and surface water influenced by SFG at large scale. Currently, catchment-scale applications of physically-based thermal-hydrology models in SFG regions are scarce, largely due to data limitations and the challenges of modeling freeze–thaw processes at larger scales. Addressing these challenges is a key focus of our future research plans.