Advancing hydrological understanding in cold regions: development and application of the WEP model for lateral flow estimation in the Great Lakes Depression of Mongolia

Analyzing and measuring water flow is crucial for studying and understanding cycles, water migration and impacts, sediment movement, pollutant transport, water resource evolution, ecological processes, and related socioeconomic aspects (Dong et al. 2022; Dong et al. 2023). Hydrological cycle research in watersheds establishes the scientific and technological foundation for long-term development. Specifically, analyzing the hydrological cycle at the river basin scale is fundamental for addressing resource and environmental challenges. It offers a solid scientific foundation for research on flood forecasting, water resources, and aquatic ecology (Jia et al. 2009; Hao et al. 2024). The significance of modeling hydrological systems has been increasingly recognized. Over the past decade, hydrologic modeling has evolved substantially (Clark et al. 2015; Palla & Gnecco 2015; Strauch et al. 2017; Anand et al. 2018), facilitating in-depth investigations of hydrologic processes and their modifications at the watershed level. The sources of uncertainty in hydrological modeling are diverse, encompassing parameters, model structure, input data (including forcing), and assessment data (such as discharge) (Kauffeldt et al. 2016). These uncertainties extend beyond the realm of meteorological forcing. The accessibility of high-resolution meteorological data, water gauge stations, land cover, and soil data has notably enhanced the convenience and reliability of hydrologic modeling. This wealth of data has become invaluable in designing and implementing effective hydrological processing and modeling for water resource planning (Ignatius & Jones 2018).

Ongoing research and development in hydrologic modeling remain imperative. Scientists have developed various effective hydrological models to simulate the spatial variability of water and energy processes in watersheds with different ecosystems. The WEP model is one of the most influential models. This grid-based model replicates the spatial variability of energy and hydrological processes in watersheds characterized by intricate ecosystem patterns (Jia et al. 2001a, b, 2006, 2009). The WEP model has been successfully applied to numerous watersheds in China, Japan, Korea, and Sri Lanka, where a variety of climatic and geographical conditions have been successfully executed by the WEP model (Jia et al. 2001a, b, 2002), its updated version WEP-L (WEP in Large basins) (Wang et al. 2021, 2023), and distributed WEP modeling (Jia et al. 2006, 2009). The WEP model operates on a spatial calculation unit composed of square or rectangular grids (Avissar & Pielke 1989; Jia et al. 2002, 2009; Kim et al. 2005; Qin et al. 2008; Noh 2011; Dahanayake & Rajapakse 2016). Runoff routing on slopes and rivers uses a one-dimensional kinematic wave technique that routes the flow downstream along the river course. For both plains and mountainous regions, numerical simulation of multilayered aquifers is carried out independently, considering streamflow, soil moisture, and groundwater interaction with surface water (Jia et al. 2009). The model takes into account the water cycle and its thermal impact in the calculation of hydrological processes.

The simulated hydrological processes include evapotranspiration, precipitation, saturation, infiltration, subsurface runoff, surface runoff, groundwater outflow, river flow, overland flow, snow melting, frozen soil consideration, and water use (Jia et al. 2001a, 2009). The hydrological cycle commences as water vapor condenses and falls as precipitation due to gravity. This process incorporates all aspects of the water cycle and employs widely recognized computational equations (Jia & Tamai 1998; Jia et al. 2003; Gao et al. 2023). Many scientists have studied research related to groundwater processes and modeling (Yang et al. 2021b, 2024). What differentiates our study from previous studies is the research focus on improving models that are theoretically sound, capable of capturing detailed hydrological processes, and suitable for sensitive permafrost regions. Previous WEP model calculations have considered surface heat exchange in comprehensive detail, but the effects of permafrost on deep soils in cold regions have yet to be well considered.

Permafrost plays a unique role in the water cycle of cold regions. Permafrost occurs in some of the colder regions of the Central Asian continent with extreme climates. The presence of permafrost significantly impacts hydrological processes (Pomeroy et al. 2022). Permafrost often affects groundwater recharge and migration, limiting groundwater flow, water movement, and groundwater storage and drainage (Walvoord & Kurylyk 2016; Gao et al. 2021). Furthermore, permafrost frequently obstructs surface water–groundwater exchanges. As a result, groundwater dynamics and interactions with surface water are greatly influenced by permafrost (Gao et al. 2021). The distribution of permafrost influences groundwater flow patterns and distribution (Liao & Zhuang 2017; Ahmed et al. 2024). Permafrost deterioration promotes infiltration, expands the groundwater reservoir, and results in a sluggish river discharge recession (Ye et al. 2009; Xie et al. 2024). All of these imply that in river basins with considerable permafrost covering, hydrological processes are significantly impacted by permafrost degradation (Niu et al. 2011; Kuang et al. 2024; Lemieux et al. 2024). Moreover, surface and subsurface thermal effects are essential in permafrost dispersion.

Permafrost degradation, driven by both natural and human activities and often linked to global warming, is characterized by reduced permafrost extent, increased active layer thickness, and the expansion of thermokarst areas (Sharkhuu et al. 2007; Kopp et al. 2017; Zhang et al. 2017). This degradation improves hydraulic links between supra-, intra-, and sub-permafrost fluids, indicating a likely transition in the hydrological system toward a groundwater-dominated regime (Lawrence & Slater 2005).

This paper addresses several unanswered questions concerning the changing environmental conditions in GLDRM. Mainly, it is essential to consider how precipitation penetrates the soil in permafrost areas, how water saturation occurs in the soil, how impervious water becomes excess on the surface due to saturation, and how it affects floods and river discharge. The vast Mongolian plateau, characterized by distinct environmental gradients primarily influenced by temperature and precipitation, is highly sensitive to water-related processes, especially those involving permafrost (Dorjsuren et al. 2023). A better understanding of hydrological regimes in arid and semi-arid regions having four distinct seasons, extreme continental climates, and diverse natural zones across the belt is essential to anticipate potential problems. This study focuses on developing and applying the groundwater process, lateral flow, and permafrost computation within the framework of the WEP model. The approach used in this research has been instrumental in estimating permafrost in the semi-arid and cold zones of the region.

The input data for the WEP model are organized into five main groups: (1) meteorology, (2) hydrology, (3) remote sensing, (4) soil, and (5) permafrost. Table 1 lists the fundamental information gathered and serves as the basis for the WEP input data.

The National Centers for Environmental Information (NCEI), hosted by NOAA (https://ngdc.noaa.gov/), and the Mongolian Information and Research Institute of Meteorology, Hydrology, and Environment (IRIMHE) are the sources for meteorological and hydrological data. The worldwide digital elevation model (DEM) is derived from remote sensing data, while the USGS National Elevation Dataset (NED) provides DEM data. Primary geographic data, including elevation, topography, slope, and overland flow direction, are derived from the 30 m resolution DEM. The soil map at a 1:1,000,000 scale serves as the source for soil types and their characteristic parameters. The statistical profiles of soil types in Mongolia provide details on the thickness and composition of each soil type. Permafrost data, including borehole records and characteristic parameters, are obtained from the IRIMHE and the Global Terrestrial Network for Permafrost (GTN-P) of the Global Climate Observing System (GCOS) (http://gtnpdatabase.org).

The WEP model's primary parameters are soil, permafrost, groundwater aquifer hydraulic conductivity, specific yield, lateral flow, and permafrost layer infiltration coefficient. Most of the model's parameters are default and do not require calibration. However, certain modifications are anticipated by comparing the simulated river discharge with observed levels over the entire simulation period (1970–2021), which was used for calibration according to the methodology recommended by (Shen et al. 2022). A comprehensive list of the calibrated parameters is provided in Table 2. It includes descriptions of the parameters, their range, and the specific calibrated values for the three gauging stations (Khovd-Myangad, Khovd-Bayannuur, and Khovd-Ulgii). Table 2 also presents the default values for reference, as summarized below:

Table 2 ensures clarity and transparency regarding each gauging station's calibration process and parameter adjustments.

As a soil phenomenon, permafrost is challenging to analyze and monitor since it depends on numerous microclimatic elements such as snow, surface radiation budget, soil temperature, physical qualities, moisture, and vegetation (Williams & Smith 1993). Calculating the rates at which the thickness of the active layer increases in boreholes is a complicated process that significantly depends on the sites' temperature, moisture content, and soil characteristics (Sharkhuu et al. 2007). The permafrost parameters used in the model are selected to calculate infiltration, saturation, and lateral flow based on borehole data, as presented in Table 4.

Permafrost classifications, as outlined by the International Permafrost Association map, encompass various extents and typologies, such as isolated, sporadic, discontinuous, and continuous permafrost. These classifications are distinguished by varying ground-ice contents and land type (Table 5) (Brown et al. 1998; Lawrence & Slater 2005).

Table 5

Permafrost typology

View Large

The permafrost extent (percentage of area) and the imposed proportion and rate are calculated using a proportion of equal sharing method (coefficient). This represents the region's proportion of permafrost area, and the residual portion, or percentage difference, can reveal the water's permeability.

Calculating the water balance facilitates the computation of underground lateral flow on the permafrost. The lateral flow (Q_lf) is calculated using Darcy's Law, and Ying Fan's water table dynamics (Fan et al. 2007).

Accurately modeling hydrological processes in permafrost regions is critical due to the unique interactions between soil, ice, and water (Yang et al. 2021a). Traditional hydrological models often struggle to represent the complex dynamics of lateral flow over permafrost, leading to uncertainties in predicting water balance components in cold regions (Sedaghatkish et al. 2024). Our introduction of a new formula for calculating lateral flow over permafrost tackles these challenges by integrating permafrost-specific parameters and processes into the WEP model framework.

One of the primary advantages of our improved formula is its ability to capture the influence of permafrost on subsurface hydrological processes more realistically. Unlike conventional models that may oversimplify or neglect the impermeable nature of permafrost layers, our formula explicitly accounts for the reduced infiltration and enhanced lateral flow above the permafrost table. This leads to a more accurate representation of soil moisture dynamics, especially during thawing periods when active layers deepen and lateral flow becomes more significant.

The natural ecological processes in Central Asia's arid and semi-arid regions, characterized by extreme continental climates and diverse natural zones, are highly fragile and sensitive (Klinge et al. 2021). The hydrological system of the GLDRM in this area is unique, and the water source is the highlands with permafrost and glaciers in the high mountains (Dashtseren 2021). However, the tributary traverses various natural zones and belts, eventually flowing into several large lakes. The hydrological regime of this region is strongly influenced by permafrost. With four distinct seasons in the area, the hydrological cycle is influenced by the cold season and soil temperature fluctuations (Munkhjargal et al. 2020). Therefore, a good study of the impact of permafrost remains important in hydrological management considerations. For example, it is essential to calculate the lateral flow of water on the underground permafrost to calculate the saturation of soil infiltration during rainfall, to assess the flood risk and the water balance correctly (Liu et al. 2023). To accurately assess the total water flow at the basin level, it is crucial to calculate and model the entire water flow process.

This study estimates the lateral flow of water over the permafrost in the GLDRM basin using river flow, precipitation, and glacial period data spanning 52 years (1970–2021). Calibration and validation criteria included (1) reducing river flow modeling error, (2) improving Nash–Sutcliffe flow efficiency, and (3) enhancing the correlation coefficient between the modeled and observed flow (Jang et al. 2018; Althoff & Rodrigues 2021). Since groundwater flow is much slower than surface or river flow, several factors must be considered such as borehole data from nearby water gauging stations, ice age data, and soil infiltration (Cook & Böhlke 2000; Somers & McKenzie 2020).

Simulating the intricate dynamics of hydrological processes in cool zone basins is challenging (Gao et al. 2021). Hydrological processes in cold regions varied considerably throughout the year, with water flow decreasing in cold winter and snow cover beginning to stabilize. Significant differences were observed between measured and simulated flows at the three locations from 1970 to 1990, while the differences were relatively minor from 1990 to 2021. Notably, from 1970 to 1990, the measured values at the Khovd-Myangad (Figure 5a) station (NSE = 0.68) and Khovd-Bayannuur (Figure 5b) station (NSE = 0.64) showed significant discrepancies from the model values. For instance, there was a notable discrepancy between observed and model values from 1984 to 1986 at the Khovd-Bayannuur gauging station. A period of high river overflow was observed, despite relatively low precipitation. However, the best results were obtained at the Khovd-Ulgii (Figure 5c) hydrological station, with an NSE = 0.75.

Based on these calculations, although there were slight discrepancies between the measured and modeled values of river water flow at some water measurement stations at the beginning of the study period, the difference was minimal toward the end of the period, from 2000 to 2020. Assessing lateral water flow beneath the soil and over the permafrost in regions with an extreme continental climate, such as Mongolia, presents several challenges. This is primarily due to the differences in regions with sporadic and isolated patches of frost, especially during spring and fall when soil thawing and freezing occur slowly relative to precipitation inputs. In such cases, model estimates may deviate from actual measurements.

When assessing the accuracy of hydrological models, the Nash–Sutcliffe model efficiency (NSE) coefficient is utilized, which ranges from 0.64 to 0.75 for the three stations during the study period, indicating that the model results meet the requirements. This suggests that Formulas (7) and (8) used in the WEP model can accurately calculate the lateral flow of the ground with the soil, demonstrating its full potential for use in future research.

While the numerical improvements in model performance indicators such as the NSE were modest, the enhanced physical realism of the model is a significant step forward (Homan et al. 2024). The new formula provides a mechanistic understanding of lateral flow in permafrost terrains, which is crucial for predicting hydrological responses to climatic variations. For instance, the model successfully simulated the delayed response of soil temperature changes at different depths, capturing the thermal inertia characteristic of permafrost soils. This temporal lag between surface and subsurface temperature variations has important implications for seasonal hydrological processes, including runoff and groundwater recharge timing.

Accurately calculating changes in soil temperature both at the surface and at depth is essential for understanding aquifer hydrological processes and improving water modeling (Condon et al. 2021). The quantity and duration of soil freezing and thawing temperatures in the groundwater process must be calculated with precision (Cui et al. 2020; Xie et al. 2021). The topsoil temperatures and temperatures at 100 centimeters below the surface were calculated and modeled using soil temperature data from three permafrost borehole monitoring stations in the study area. The measured topsoil temperatures and the simulated temperature results at the three boreholes align closely, with NSE coefficients ranging from 0.95 to 0.98 (Figure 6).

It is observed that when the soil temperature reaches a depth of 100 cm, there is a discrepancy between the actual measurements and the model values. For example, at the Khongor Ulun (Figure 7a) well, the soil temperature remains positive from early June to mid-October, indicating that the temperature change at a depth of 100 cm occurs approximately 1 month after the surface soil temperature changes. A similar 1-month lag is seen in the temperature changes at both the upper and lower soil layers. In other words, the transfer of heat and cold in the soil does not follow the same dynamics as the warm and cold air and the precipitation input. The transfer of heat and water in the atmosphere occurs rapidly, while in the soil, the process is slower, leading to a delayed response and resulting in a difference between the two. This leads to consistency in the results of the model.

Precipitation infiltration into the soil, soil absorption and water seepage, along with lateral flow over snow are significantly lower than surface water migration and river flow (Wu et al. 2023; Li et al. 2024). This will result in a discrepancy compared to the actual observed value. However, in such areas, winter and summer hydrological processes are influenced by consistent environmental factors, leading to no significant difference between the two. The calculations in winter and summer are correct and effective. Therefore, new Equations (8) and (9), which account for the lateral flow of water over subsurface permafrost, are novel and important for water modeling in extreme continental climate regions, so this ‘new’ ‘equation’ can be used. It is recommended that the equations from the WEP model and those from similar water models be utilized in future studies.

Another key benefit of our approach is its applicability to diverse permafrost conditions without extensive recalibration. The formula incorporates parameters such as permafrost extent and ground-ice content, allowing it to adapt to various permafrost typologies outlined by the International Permafrost Association. This flexibility makes the model suitable for large-scale applications across different cold regions, enhancing its utility for hydrological forecasting and water resource management.

The improved formula also addresses soil heterogeneity and anisotropy limitations in permafrost regions. By considering the soil's hydraulic conductivity and layering above the permafrost, the model captures the preferential flow paths that often occur due to soil texture variations (De Vrese et al. 2023). This results in more accurate estimates of subsurface flow contributions to the total streamflow, which are particularly important during periods of active layer development.

Furthermore, our work highlights the importance of integrating permafrost data, such as borehole measurements and ground temperature profiles, into hydrological modeling. The alignment between modeled and observed soil temperatures at various depths reinforces the model's robustness in simulating thermal and hydrological processes concurrently (Mihalakakou et al. 2024). This integration is essential for assessing the impacts of permafrost thaw on hydrological regimes, especially under the influence of climate change.

In terms of practical implications, the enhanced model can improve flood risk assessments and water resource planning in permafrost regions. By accurately simulating lateral flow and soil moisture dynamics, stakeholders can make more informed decisions regarding infrastructure development, ecosystem conservation, and disaster mitigation (Kumar et al. 2021).

Despite the advancements, our study has limitations. The modest numerical improvements suggest that while the new formula enhances physical realism, other factors may influence model performance. These may include data quality issues, such as gaps or inaccuracies in meteorological and hydrological records or the need for finer-scale spatial resolution to capture local heterogeneities.

Future research should focus on integrating remote sensing data to refine spatial inputs like soil moisture and active layer thickness. Additionally, coupling the model with climate change projections could provide insights into how permafrost degradation might alter hydrological processes. Long-term monitoring and data collection efforts are also essential to validate and further enhance the model's predictive capabilities.

The WEP model perfectly simulated the lateral flow of water in permafrost regions. This study focused on the circulation of the hydrological system in the permafrost zone of the Khovd River-Khar-Us Lake basin, located within the GLDRM. Regarding the computational aspect of the WEP model, two new Formulas (7) and (8) were developed to accurately calculate the lateral flow of water in the groundwater aquifer. These formulas were applied to measured and modeled data spanning 52 years, from 1970 to 2021. Additionally, the study included estimating and modeling the soil temperature at the surface and a depth of 100 centimeters below ground.

Formulas (7) and (8) were applied to estimate the lateral flow of water below the ground surface and above the permafrost in cold permafrost regions over 52 years, which produced promising results. The calculations for the GLDRM, specifically in the Khovd River-Khar-Us Lake Sub Basin, were executed with considerable accuracy. The NSE coefficient was estimated to be between 0.64 and 0.75, indicating that the model is plausible and suitable for further research and application.

Moreover, modeling soil temperatures at the surface (0 cm) and a depth of 100 cm below the ground revealed a temporal discrepancy. The temperature difference between these two depths was about 1 month, with the temperature decrease occurring almost simultaneously at both depths. The NSE coefficients for these temperature simulations were notably high, ranging from 0.95 to 0.98 for the upper soil layer. However, the NSE coefficient values decreased to 0.35–0.80 at a depth of 100 cm below the ground. Nevertheless, there was observed a good agreement between measured and modeled values for the simulated soil temperature.

It would be beneficial to apply the newly formulated Formulas (7) and (8) to assess the lateral flow of water over permafrost under various conditions in other basins of diverse topography. Such extensive testing will provide deeper insights and enhance the robustness of the model, thereby contributing significantly to the understanding and management of hydrological processes in permafrost regions.

B.D., and D.Y. contributed to conceptualization; B.D. performed methodology; Y.Y. did software analysis; Y.Y., and D.N. did validation; D.S. did formal analysis; B.Z. did investigation; H.G. collected resources; Y.Y. did data curation; D.B. contributed to writing – original draft preparation; D.N. contributed to writing – review and editing; H.Z. and H.G. contributed to visualization; D.Y. did supervision; D.Y. did project administration; B.D. acquired funds. All authors have read and agreed to the published version of the manuscript.

This research was funded by the Ministerial Innovation Scholarship for Postdoctoral Research (grant: 19ХХ04DI208), the National Key Research and Development Project (Grant No. 2016YFA0601503), the Mongolian Science and Technology Foundation (grant CHN-2022/274), and the National Key Research and Development Project of China (grant 2022YFE0119400).

The authors declare there is no conflict.

All relevant data are included in the paper or its Supplementary Information.