Integrated Hydrologic Modelling of Groundwater-Surface Water Interactions in Cold Regions

Yang, Xiaofan; Hu, Jinhua; Ma, Rui; Sun, Ziyong

doi:10.3389/feart.2021.721009

REVIEW article

Front. Earth Sci., 01 December 2021

Sec. Cryospheric Sciences

Volume 9 - 2021 | https://doi.org/10.3389/feart.2021.721009

This article is part of the Research Topic Ecological Impacts of Degrading Permafrost View all 19 articles

Integrated Hydrologic Modelling of Groundwater-Surface Water Interactions in Cold Regions

Xiaofan Yang¹*

Jinhua Hu¹

Rui Ma²

Ziyong Sun²

¹State Key Laboratory of Earth Surface Processes and Resource Ecology, Faculty of Geographical Science, Beijing Normal University, Beijing, China
²School of Environmental Studies and State Key Laboratory of Biogeology and Environmental Geology, China University of Geoscience, Wuhan, China

Groundwater-surface water (GW-SW) interaction, as a key component in the cold region hydrologic cycle, is extremely sensitive to seasonal and climate change. Specifically, the dynamic change of snow cover and frozen soil bring additional challenges in observing and simulating hydrologic processes under GW-SW interactions in cold regions. Integrated hydrologic models are promising tools to simulate such complex processes and study the system behaviours as well as its responses to perturbations. The cold region integrated hydrologic models should be physically representative and fully considering the thermal-hydrologic processes under snow cover variations, freeze-thaw cycles in frozen soils and GW-SW interactions. Benchmarking and integration with scarce field observations are also critical in developing cold region integrated hydrologic models. This review summarizes the current status of hydrologic models suitable for cold environment, including distributed hydrologic models, cryo-hydrogeologic models, and fully-coupled cold region GW-SW models, with a specific focus on their concepts, numerical methods, benchmarking, and applications across scales. The current research can provide implications for cold region hydrologic model development and advance our understanding of altered environments in cold regions disturbed by climate change, such as permafrost degradation, early snow melt and water shortage.

Introduction

Cold regions are headwaters of many prominent rivers around the world and considered as the water towers of inland river basins in arid/semi-arid regions (Yao et al., 2012; Qin et al., 2017; Immerzeel et al., 2020). As unique and critical elements of the terrestrial cryosphere, frozen soil and snow cover not only have significant influences on the hydrologic processes in cold regions (Ding et al., 2020; Wang et al., 2020; Gao et al., 2021a; McKenzie et al., 2021), but they are also extremely sensitive to seasonal (short-term, from day to month) and climate changes (long-term, from year to decade) (Kang et al., 2020; IPCC, 2019 SROCC). The pronounced global warming has led to permafrost degradation due to increased ground temperature, permafrost thawing, thickening of the active layer, shortening of the freezing period, and prolonged snowmelting period, which directly affect groundwater recharge, runoff, and promote groundwater-surface water (GW-SW) interactions (Nitze et al., 2018; Cheng et al., 2019; Teufel and Sushama, 2019; Zhao et al., 2019; Evans et al., 2020; Lemieux et al., 2020). Meanwhile, the snow-induced changes, such as shortened snow duration, reduction of the snow cover area, earlier snowmelt, and more frequent extreme snowfall events, have significantly altered the surface albedo and runoff, and indirectly affect the groundwater hydrothermal conditions, the GW-SW interactions, and the energy exchange by changing the surface temperature and soil water infiltration (Barnett et al., 2005; Barnhart et al., 2016; Wu et al., 2018; Che et al., 2019). Therefore, GW-SW interaction is one of the key components in the cold region hydrologic cycle (Ge et al., 2011; Liljedahl et al., 2016). Comprehensive research on GW-SW interactions in cold regions under seasonal and climate changes, their driving factors, and the associated mechanisms is essential to understanding the source and mechanisms of runoff and the response of the hydrologic processes in cold regions. Those insights would provide fundamental evidence for the scientific management of inland basin water resources and promote inter-disciplinary research (i.e., cryosphere science, atmospheric science, hydrology, hydrogeology, etc.).

Seasonal and climate changes significantly affect the hydrologic and energy cycles in cold regions. The dynamic flow and heat transfer processes in such a complex system are driven by the interplay between snow cover, soil freeze-thaw cycles (FTC), surface water and groundwater flow, which should be fully considered in an integrated framework (Flerchinger and Saxton, 1989; Winter et al., 1998; Atchley et al., 2016; Hubbard et al., 2018). However, due to the harsh environment, the measuring range, frequency, and precision of hydrometeorological and geophysical field observations are quite limited (Schilling et al., 2019a; Che et al., 2019). Yet it is difficult to predict the hydrologic processes through field/site control experiments (i.e., soil warming experiments and snow cover control experiments, Melillo et al., 2002; Fu et al., 2018) and indoor snowmelt and freeze-thaw experiments (Waldner et al., 2004; Costa and Pomeroy, 2019; Riviere et al., 2019; Deprez et al., 2020) due to constrained spatio-temporal scales (Tokunaga et al., 2019; Chiasson-Poirier et al., 2020). Therefore, under the scenario of seasonal and climate changes, establishing an integrated hydrologic model to quantitatively simulate and predict the coupled processes of flow and heat transfer, their influencing factors and mechanisms, and most importantly, investigating the dynamics of GW-SW interactions under freeze-thaw cycles and snow cover variations in cold regions are emerging and Frontier topics of cryosphere, hydrology, and hydrogeology research. This paper provides an advanced review on the current research status of integrated hydrologic models of cold regions with future perspectives. The rest of the paper is organized as follows. The research progress in numerical models, benchmarking, and simulations are summarized in Research Progress in Integrated Hydrologic Modeling of Groundwater-Surface Water Interactions in Cold Regions, followed by future research and perspectives on model development and potential applications in Challenges and Future Research Needs. A brief summary is given in Summary.

Research Progress in Integrated Hydrologic Modeling of Groundwater-Surface Water Interactions in Cold Regions

Surface water (SW) and groundwater (GW), which are essential elements in the global and regional hydrologic cycles, are strongly coupled and dynamically interactive. However, surface water and groundwater are often considered to be separate entities and are investigated individually because the physical, chemical, and biological properties of surface water and groundwater are largely varied (Yang et al., 2012). Under seasonal and climate changes, snow cover and frozen soil poses additional challenges in studying GW-SW interactions in cold regions, and thus, the fundamental processes become more unique, complex, and uncertain (Figure 1, Ireson et al., 2013). The temperature variations caused by seasonal and climate changes influence the snow accumulation and melting processes first. Decreasing temperature in autumn and winter provide feasible conditions for snowfall and promote the accumulation of snow cover. Deep snow cover has a heat conservation effect on the soil while shallow snow cover decreases the soil temperature (Rasouli et al., 2019; Mazzotti et al., 2020). Rising air temperature also leads to enlarged areas of snow melting, increasing the surface water flow velocity and water depth as the snowmelt flows into surface runoff, leading to more surface water infiltration and accelerated interaction with groundwater (Barnhart et al., 2016; Chen et al., 2017). Since the surface temperature and evapotranspiration in spring are low in some mountainous and plateau areas, the snow melt even becomes an important source of groundwater recharge (Carroll et al., 2019; Li et al., 2019). Snow melting, radiation, reflection, and blowing snow also strongly alter the local and regional surface energy balance, thus affecting the thermal conditions of the groundwater (Kinar and Pomeroy, 2015). Furthermore, permafrost can act as a barrier in the soil system, which affects the spatial-temporal heterogeneity of the aquifer, thus regulating the preferential flow path of groundwater and influencing groundwater recharge, runoff, and outflow (Xu et al., 2001; Bense et al., 2009). In addition, hydrometeorological factors, such as precipitation, temperature, and evapotranspiration, evidently have impacts on the hydrothermal conditions and moisture transport in the active layer, leading to periodic thawing (Woo, 2012; McKenzie and Voss, 2013). Soil freeze-thaw processes not only result in the phase change between ice and liquid water in groundwater systems, but they also lead to the frequent transport of groundwater in different layers of aquifers and between subsurface and surface systems, affecting the amount of water involved in the regional hydrologic cycle. Moreover, this also changes the structure of the aquifers, thus affecting the preferential flow path of groundwater and making the GW-SW interactions more complex in cold regions (Walvoord and Kurylyk, 2016).

FIGURE 1

FIGURE 1. Cold region GW-SW interactions before (A) and after (B) climate change.

With the rapid development of numerical algorithms, computers, and scientific computing techniques, numerical models have become important tools for investigating the hydrologic cycle in cold regions under seasonal and climate changes (Kurylyk et al., 2014a; Ding et al., 2017; Yang et al., 2019). As the connection between the subsurface and surface components of the hydrologic cycle in cold regions, the major problem of modelling GW-SW interactions in cold regions is to accurately describe the coupled thermal-hydrologic (TH) processes in such a complex system consisting of snow cover, surface water, groundwater, and frozen soil under disturbances. Under the influence of temperature, precipitation, and evapotranspiration, the flow and heat transport within the snow cover and the phase change between water and ice are dynamically coupled and complex, the key processes of which include snow accumulation, snowmelt, snow sublimation, blowing snow, and the radiation and reflection from the snow cover (Domine et al., 2019). Surface hydrothermal conditions also affect the flow and heat transport processes in frozen soils, including multiphase flow, heat conduction and convection, water vapour diffusion, and phase changes between ice, liquid water, and vapour. It is noted that different phases and moisture distributions in the soil pores are strongly correlated with temperature during the FTCs, which has a significant impact on the soil physical properties, such as the porosity and permeability. As a result, the constitutive equations describing the coupled TH processes are highly nonlinear in groundwater-permafrost systems (Kurylyk and Watanabe, 2013). Despite these difficulties, the dynamic interactions between the surface and groundwater make it more difficult for hydrologic modeling and simulations (Maxwell et al., 2015). Therefore, it is extremely challenging to establish an integrated GW-SW model in such a way that all the aforementioned processes are fully considered and precisely described based on fundamental phenomena, including, but not limited to, multi-phase flow in porous media, phase changes, and the coupled thermal-hydrologic processes and interactions between groundwater and surface water (Coon et al., 2016; Painter et al., 2016).

There are a growing number of theoretical and numerical models that are capable of describing the integrated surface and subsurface flow and heat processes in cold regions. These models have been validated through benchmark cases, laboratory experiments, and field observations, and have been applied to predict future trends and responses of the hydrologic system under different seasonal and climate change scenarios (Figure 2). Studies of integrated hydrologic models and numerical simulations of cold regions can be divided into three categories: 1) numerical model development, 2) model benchmarking and validation, and 3) model simulation and application.

FIGURE 2

FIGURE 2. The current research status of modelling GW-SW interactions in cold regions.

Numerical Model Development

Most existing hydrologic models are developed for general purposes, including: 1) distributed hydrologic models based on terrestrial hydrology (Yang et al., 2015), 2) groundwater/hydrogeological models based on hydrogeology (Yao et al., 2014; Tian et al., 2018), and 3) physically-based and fully-coupled GW-SW models (Kuffour et al., 2020). In order to fundamentally understand and precisely describe the hydrologic cycle in cold regions, it is natural to modify/improve selected previously developed hydrologic models in order to facilitate the characteristics of the study area and the research needs, leading to three groups of cold region integrated hydrologic models (Table 1): 1) distributed hydrologic models of cold regions, 2) cryo-hydrogeological models, and 3) cold region GW-SW models.

TABLE 1

TABLE 1. Summary of cold region integrated hydrologic models.

The concept of modular-based hydrologic models and distributed hydrologic models based on terrestrial hydrology in cold regions is to establish sub-models/modules for all the sub-processes/subunits in the hydrologic cycle then integrate them into a unified modelling framework. In distributed hydrologic models, such as CBHM (Chen et al., 2018), CRHM (Marsh et al., 2020), GBEHM (Li et al., 2019), VIC (Hamman et al., 2018), HydroSiB2-SF (Wang G. X. et al., 2017), WEB-DHM-pF (Song et al., 2020), different sub-models/modules are “stitched” and tuned before being used to simulate and predict certain processes in selected research areas. For example, the module-based cold regions hydrologic model (CRHM, Marsh et al., 2020), which was developed at the University of Saskatchewan in Canada, integrates the most complete snow model composed of a series of snow cover-related modules, including radiation, reflection, snow sublimation, snow melting, and blowing snow. Although such module-based distributed hydrologic models integrate all the possible components of the hydrologic cycle in cold regions, they mainly focus on the surface hydrologic process and the water and energy balances, which subsequently simplifies the groundwater by adopting a parameterized scheme or 1D bucket model, ignoring processes such as lateral flow and heat transfer. Moreover, the spatial resolution of these models is usually low (at the kilometer level), so they cannot be used to investigate smaller-scale processes. At present, distributed and module-based hydrologic models of cold regions are mostly used to simulate the effects of seasonal and climate changes on frozen soil, surface hydrologic processes, and vegetation in large-scale watersheds (Zhou et al., 2014; Chen et al., 2018; Gao et al., 2018; Qi et al., 2019).

In theory, the main idea of the cryo-hydrogeological model, which considers the groundwater, hydrogeology, and permafrost, is to establish the constitutive equations of multiphase groundwater flow and heat transfer (based on the theories of flow and heat transfer in porous media), which characterizes the effects of the FTCs on the pore water, soil porosity, and permeability using the soil freeze-thaw curve or a physical process-based constitutive relationship (such as the Clausius-Clapyeron equation). Furthermore, they simulate the soil moisture and temperature, groundwater runoff, and its recharge from surface runoff in the basin (Kurylyk and Watanabe, 2013). If considering surface hydrologic processes in cold regions, such models mostly use temperature, precipitation, evapotranspiration, and snowmelt as the boundary conditions and forcing data, or one-way coupling (weak coupling) with surface hydrologic process models (such as simultaneous heat and water, SHAW models) (Voss et al., 2018). In particular, the cryo-hydrogeological models use snow accumulation and ablation (key elements in snow hydrology) as the upper boundary of the simulation domain by altering the surface temperature, while ignoring the physical processes of moisture and energy transport during snowmelt and its impact on the runoff and GW-SW interactions. Therefore, hydrogeological models in cold regions are mostly used to study the hydrogeological response of frozen soil or the effects of seasonal and annual changes on surface processes (i.e., soil moisture, temperature, and dynamic changes in the groundwater) (Evans et al., 2018; Jafarov et al., 2018). The models that have been reported for hydrogeological research in cold regions include but not limited to SUTRA (or modified SUTRA, SUTRA-ICE, McKenzie et al., 2007; McKenzie and Voss, 2013), Cast3M (www-cast3m.cea.fr), GSFLOW (Markstrom et al., 2008), FEFLOW (DHI, 2016), and PFLOTRAN-ICE (Karra et al., 2014), DarcyTools (Svensson et al., 2010; Svensson and Ferry, 2014), GEOAN (Holmén, 2019), Ginette (Rivière et al., 2019), MELT (Frederick and Buffet, 2014), SMOKER (Molson and Frind, 2019), FlexPDE (Bense et al., 2009), PermaFOAM (Orgogozo et al., 2019), COMSOL (Scheidegger et al., 2017).

In recent years, to accurately and efficiently describe the hydrologic processes in cold regions, the cold region GW-SW model has gradually emerged on the basis of the original GW-SW models and cryo-hydrogeological models (Jan et al., 2020). Such models are considered as the extensions of the cryo-hydrogeological models, while the groundwater-frozen soil modules still follow or directly adopt the cryo-hydrogeological models. The difference is that the surface hydrologic processes (i.e., snowfall, snow accumulation and melt, runoff, and evapotranspiration) are no longer simply presented as boundary conditions, but rather, they are included in the form of governing equations (Schilling et al., 2019b). In this type of models, the constitutive equations of each hydrologic process are simultaneously solved and coupled in an interactive way (two-way coupling or strong coupling). Such models constitute the transient partial differential equations (PDEs) that are iteratively solved using numerical discretization and high-performance computing techniques. Moreover, soil hydrologic properties (e.g., the soil water content and infiltration rate) are directly used to help identify the water flow and energy transport processes under saturated/unsaturated conditions (Painter et al., 2016; Jan et al., 2018). Recently developed models include but not limited to the integrated cryohydrologic model GeoTop (Rigon et al., 2006; Endrizzi et al., 2014), the Advanced Terrestrial Simulator (ATS; Atchley et al., 2015; Painter et al., 2016; Coon et al., 2019)—a hydro-ecological model for cold regions developed by the U.S. Department of Energy (US DOE), and the HydroGeoSphere (HGS; Aquanty, 2017) - an integrated hydrologic model developed by the University of Waterloo in Canada. The cold region GW-SW models combine the advantages of the previous two types of models. Physical processes in frozen soil and snow are both been taken into consideration in the cold region GW-SW models, which use a two-way coupling numerical method to connect and simultaneously simulate the dynamic GW-SW processes. In addition, their spatial resolutions are quite high (∼cm), which is suitable for simulating small-scale processes. Although such models have a great potential for future development and applications, the model structure, coupling algorithms, and parallel computing efficiency still need to be improved. In addition, most of the current cold region GW-SW models are mostly in-house developed or commercialized software, which makes it difficult for community development that greatly limits the scope and promotion of such models. It is also noted that the aforementioned models are often assigned with prior assumptions on permafrost hydrological behavior as well as the subsequent impacts on catchment hydrology (Walvoord and Kurylyk, 2016; Gao et al., 2021a). Such models follow a “bottom-up” modeling approach with most of the process understanding obtained from in-situ observations, which have limited spatial and invariably limited temporal coverage (Lindstrom et al., 2002; Niu and Yang, 2006; Brutsaert and Hiyama, 2012). In recent years, “top-down” modeling methodology learned from data has emerged, such as soil temperature–based water saturation function (Wang L. et al., 2017), HBV (Osuch et al., 2019), FLEX-Topo (Gao et al., 2021b), etc. Undoubtedly, these “top-down” modeling approaches provide new insights on the impact of permafrost on catchment hydrology, which also have strong potential in future model development. In this paper, we are mainly focused on “bottom-up” modeling approach.

Model Benchmarking and Intercomparison

The representation of the physical processes, parameterization schemes, numerical discretization schemes, and grid types of the same type of numerical models are often diverse, and their accuracy and physical realism is quite uncertain. This motives research on model testing, intercomparison, and validation using benchmark cases in order to comprehensively evaluate the physical accuracy, completeness, calculation accuracy, and efficiency of these models (Kurylyk and Watanabe, 2013). A summary of the model benchmarking and intercomparison studies is given in Table 2.

TABLE 2

TABLE 2. Summary of the model benchmarking studies.

With the recent advancement of integrated GW-SW modelling frameworks, Maxwell et al. (2014) identified the need for a formalized Integrated Hydrologic Model Intercomparison Project (IH-MIP), and to this end, they established standard procedures and benchmark test cases for coupled surface-subsurface models and presented the results of an intercomparison study of seven coupled surface-subsurface models, which included the CATHY (Bixio et al., 2002; Camporese et al., 2010), HGS, OGS (Delfs et al., 2009; Delfs et al., 2012), PIHM (Kumar et al., 2009), PAWS (Shen and Phanikumar, 2010), ParFlow (Kollet and Maxwell, 2006) and tRIBS-VEGGIE (Ivanov et al., 2004) models. These seven integrated hydrologic models were intercompared using a standard set of benchmark test cases (i.e., excess infiltration, excess saturation, tilted V-catchments, slabs, and return flow). Kollet et al. (2017) introduced the results of the second phase of the Integrated Hydrologic Model Intercomparison Project (IH-MIP2). ATS, Cast3M, CATHY (Bixio et al., 2002; Camporese et al., 2010), GEOtop, HGS, MIKE-SHE (Abbott et al., 1986; Butts et al., 2004), and ParFlow (Kollet and Maxwell, 2006) were investigated in the IH-MIP2. Compared with the first phase, the complexity of the benchmarking tests was increased in IH-MIP2, which include a tilted V-catchment with 3D subsurface, a superslab case with additional horizontal subsurface heterogeneity, and the Borden site rainfall-runoff experiment. Parameters such as saturated, unsaturated, and ponded storages, discharge, as well as velocity profiles in the superslab and Borden site benchmarking tests were intercompared among seven integrated hydrologic models.

McKenzie et al. (2007) first considered an analytical solution of the 1D Stefan heat transfer equation in model benchmarking and proposed a 2D thawing case in an ideal frozen soil medium. Kurylyk et al. (2014b) revisited the 1D analytical solution and suggested that heat conduction and phase change should be considered for cryo-hydrogeological code benchmarking. In order to comprehensively verify and optimize the integrated hydrologic model of cold regions, the InterFrost project (wiki.lsce.ipsl.fr/interfrost) was launched at the end of 2014. The objective of the InterFrost project was to propose a series of simple to complex, practical, representative benchmarking cases for model comparison and benchmark verification and to further optimize and develop these models. Identical input and setups were configured in each of the models investigated, and the numerical simulation results were inter-compared with the 1D analytical solution, results from 2D benchmark cases, measured data from laboratory experiments, and observation data collected from cold region watersheds. By using the 1D analytical solution and the 2D benchmark cases, Rühaak et al. (2015) attempted to compare the evolution of temperature simulated by 4 of the codes assessed by the InterFrost project. Grenier et al. (2018) introduced 13 cryo-hydrogeology codes/softwares that participated in the first phase of the InterFrost project, comprehensively evaluated and analysed the performance of the models using the analytical solution of the 1D Stefan equation and results from 2D ideal benchmark cases (frozen inclusion thaw and talik opening/closure) using a series of indicators, and finally, they proposed a path to develop and verify cryo-hydrogeology models (Figure 3). It is worth noting that the InterFrost project is still actively advancing, and the next two phases will continue to discuss more practical benchmarking verification and evaluation based on laboratory freezing-thawing experiments and field observations in cold regions.

FIGURE 3

FIGURE 3. The flowchart and examples of cold region hydrologic model benchmarking.

Nevertheless, the current research on integrated hydrologic models of cold regions largely lacks benchmarking verification, also the sub-models (e.g., snow cover and frozen soil representation) and algorithms (e.g., GW-SW coupling and phase change) still need to be verified. In particular, the constitutive equations, parameter selections, and calibrations need to be assessed for groundwater flow and heat transport in frozen soil and the snow accumulation and ablation processes. The integration and coupling method of these sub-models also needs to be evaluated. In addition, the current benchmarking problems are too simple and mostly focus on the verification of freeze-thaw models. It is essential to evaluate models based on field data (e.g., soil moisture and temperature, snow depth and snow water equivalent, Schilling et al., 2019b; Chen et al., 2017; Gao et al., 2017; Gao et al., 2020) and control experiments.

Model Simulation and Prediction

Distributed hydrologic models of cold regions are mostly used to simulate and predict surface hydrologic processes in larger-scale watersheds, but they may not be feasible to simulate the dynamic processes of GW-SW interactions in cold regions due to their over-simplification of subsurface processes. Cryo-hydrogeological models are mainly used to simulate and predict the hydrogeological response of frozen soil affected by climate and seasonal changes, without fully considering the physical processes of snow cover and its influence on the interaction between groundwater and surface water in cold regions. The fully-coupled cold region GW-SW model can accurately simulate the dynamic processes of GW-SW interactions; however, it is still in development and has limited applications, especially in the field. In recent years, research on modelling GW-SW interactions in cold regions included, but was not limited to, describing the dynamic processes of GW-SW interactions at multiple scales (e.g., hillslope scale (Woo, 2012), catchment scale (Evans et al., 2015), watershed scale (Lemieux et al., 2020)) and predicting their temporal and spatial changes, quantitatively analyzing groundwater lateral flow and heat transfer, local-scale permafrost heterogeneity (e.g., taliks), micro-surface processes (e.g., micro-topography), the impact of snow accumulation and melting, frozen soil on the hydrologic processes, and GW-SW interactions in cold regions under seasonal and climate changes (Table 3).

TABLE 3

TABLE 3. Summary of cold region integrated hydrologic models and their applications.

Hillslope Scale

At the hillslope scale, SUTRA (including SUTRA-ICE and modified versions of SUTRA), developed by the US Geological Survey (USGS) has been extensively utilized in various studies. Ge et al. (2011) established a two-dimensional conceptual model of the slope using SUTRA-ICE, which motivated the study of GW-SW interactions in cold regions. The model was used to simulate the subsurface TH processes in the permafrost region of the Qinghai-Tibetan Plateau (QTP) and to calculate the subsurface runoff and its recharge of the surface runoff in the basin. Wellman et al. (2013) applied SUTRA-ICE to address the development of sub-lake and supra-permafrost talik under climate change. A 2D cylindrical cross section (1800 m in radius and 500 m in height) in Yukon Flats, Alaska (United States) was designed to study the evolution of sub-lake talik evolution and supra-permafrost near a lake under various environmental conditions over 1,000 years. Evans and Ge (2017) applied SUTRA to establish a 2D slope model in the Alaska permafrost area, the northern QTP, the seasonal permafrost in the Colorado Rockies, United States, and New Brunswick, Canada. Different hydrogeological effects of permafrost and seasonal frozen ground degradation on groundwater flow were comprehensively explored under climate change scenarios. Simulated results revealed that there is greater groundwater flow in hillslopes with permafrost than seasonal frozen ground under warming scenario. Lamontagne-Hallé et al. (2018) used SUTRA and an improved set of boundary conditions to simulate the distribution of the frozen soil and the dynamic groundwater discharge in permafrost regions under climate warming. Zipper et al. (2018) employed the modified version of SUTRA to study groundwater-permafrost interactions under the influence of wildfire in a 2D transect (200 m width, 20–25 m height) that simply represented the Anaktuvuk River field sites. To examine active-layer dynamics and perennial thaw zones in boreal hillslopes, Walvoord et al. (2019) applied the enhanced version of SUTRA to an idealized 2D hillslope model (100 m long) under fully-saturated condition. Various idealized boundary conditions subjected to seasonal air-temperature variation in a long-term warming trend were superimposed. The enhanced version of SUTRA could simulate coupled fluid-flow and heat transport incorporating the nonlinear phase change between ice and liquid water. Model results indicated that perennial thaw zones serve as a detectable precursor to accelerated permafrost degradation. Rushlow et al. (2020) investigated the effects of snow cover and air temperature on the active layer, the freeze-thaw states, and the groundwater flow near water tracks using SUTRA.

The Advanced Terrestrial Simulator (ATS), developed by Los Alamos and Oak Ridge National Laboratories (USA), is an ecosystem hydrology code focusing on solving problems of integrated surface/subsurface hydrology that involve with vegetation, surface energy balance, snow, and other environmental interactions. HEATFLOW/SMOKER, developed by Université Laval and University of Waterloo (Canada), is a finite element numerical model for solving complex density-dependent groundwater flow, contaminant transport, groundwater age, and thermal energy transport problems. Both codes are open-source and have emerged as powerful tools in solving flow and heat transport problems within partially frozen ground in cold regions. Sjöberg et al. (2016) employed the ATS to estimate the thermal effects of groundwater flow in permafrost environments. Numerical simulations were configured using 1D vertical columns and 2D vertical (10 m long × 5 m deep) cross-sections of the subsurface using data from a subarctic fen in northern Sweden. Shojae Ghias et al. (2017) developed a conceptual model in the Iqaluit Airport taxiway, Nunavut, Canada. A corresponding 100 m wide and 50 m high vertical cross-section was calibrated via HEATFLOW/SMOKER. Simulated results were compared with the observed ground temperature, and the impacts of future climate-warming on the thermal regime and flow system under climate scenarios proposed by the Intergovernmental Panel on Climate Change (IPCC) were analyzed. To better understand the thermal-hydrologic controls and consequences of taliks, Jafarov et al. (2018) used the ATS to simulate the evolution of taliks developed under preferentially distributed snow. The numerical simulations were conducted in a hillslope transect (45 m deep, 1 km wide, 2D domain) on the southern Seward Peninsula. Dagenais et al. (2020) adopted the cryo-hydrogeological numerical model HEATFLOW to study the flow and heat transfer processes in an ice-rich permafrost mound. First, the HEATFLOW model was verified using a vertical cross-section (600 m long × 200 m high) extracted from a 3D geological model (0.8 km long × 2.6 km wide × 200 m high) of the Tasiapik valley, near Umiujaq, Nunavik (Québec), Canada. Numerical simulations were then performed for future projection under a regional climate-change scenario. Based on the “Umiujaq model” provided in Dagenais et al. (2020), Albers et al. (2020) performed HEATFLOW-SMOKER to investigate the sensitivity of the model to changes in the thermal and hydraulic properties of the system. It was demonstrated that the model was quite sensitive to the near-surface thermal and hydraulic parameters, and to those parameters which define the characteristics of the active layer. Sjöberg et al. (2021) applied ATS to investigate if the impact of permafrost on flow path depth could cause the same pattern in temperatures of groundwater discharging from hillslopes to streams. The model domain was 1,000 m in the horizontal (x) direction, 8 m thick (z), and had a slope of 0.1, which represents the slopes near the streams in the Agashashok basins and the Cutler-Imelyak basins located in the Noatak National Preserve in northwestern Alaska. Hamm and Frampton (2021) used ATS to simulate two idealized hillslopes (50 m in x-direction, 1 m in y-direction, and 20 m in z-direction) with inclinations (22°, 11°) that can be found in Adventdalen, Svalbard, and compare them to a baseline case (0°). Results show that hillslope inclination causes differences in ground temperature uphill and downhill and highlight the relevance of considering lateral flow of water in the subsurface combined with heat flux for modelling arctic catchments with permafrost.

Other in-house developed or commercial software were also used in simulating groundwater flow and heat transport in cold regions. Bense et al. (2012) studied the climate change-induced hydraulic evolution of aquifers and the consequent increase in baseflow in idealized permafrost environments using the FlexPDE (developed using the finite element method). Wicky and Hauck (2017) used the commercially available GeoStudio software to model the air flow and heat transfer in a natural talus slope and studied the ground thermal regime. The entire domain represented an idealized 170 m long talus slope with a talus thickness of 17 m and a constant slope angle of 33°. Orgogozo et al. (2019) developed the cryo-hydrogeological numerical model (permaFOAM) in the open-source CFD environment, OpenFOAM (https://openfoam.com). Based on permafrost observation in central Siberia in Russia, seasonal changes in soil moisture and temperature were simulated in a 2D transect. The permaFOAM accurately described the phase change between the pore water and ice in frozen soils and quantitatively analysed the influence of surface evapotranspiration on the coupled TH processes in the groundwater in frozen soils. Scheidegger et al. (2019) applied the commercial software COMSOL to simulate the impact of long-term climatic variability on groundwater flow and permafrost dynamics in two contrasting geological settings (20 km long × 3 km high and 300 km long ×1.5 km high) in Great Britain. The sensitivity of simulated permafrost thickness to a variety of climatic and subsurface conditions was also evaluated. Using FEFLOW, Huang et al. (2020) established an ideal 2D slope model of the central part of the QTP and predicted the hydrogeological processed in high altitude areas in the context of future climate change. The impacts of the slope aspect, land surface temperature, FTCs, and permeability were all analysed, and finally, the land surface temperature was identified as the main factor controlling the seasonal dynamics of the shallow groundwater systems in permafrost regions.

Catchment Scale

At the catchment scale, Endrizzi et al. (2014) applied GEOtop 2.0 to a synthetic catchment (two 160 m long and 200 m wide hillslopes forming a convergent topography) and simulated the GW-SW interactions with frozen soil and snow cover under different meteorological forcings. Evans et al. (2015) used SUTRA to describe the groundwater flow in mountainous headwater catchments and evaluated the impact of permafrost on groundwater. The model was applied to a representative mountainous Hulugou catchment (8 km long, 3 km wide. elevation range from 2,970 m to 4,470 m) on the QTP, China, and then, they obtained predictions under a warming scenario with a mean annual surface temperature increase of 2°C. Kumar et al. (2016) applied a multiphase subsurface thermal hydrology model (PFLOTRAN; http://www.pflotran.org/) to study the thermal regimes of 3D ice-wedge polygons extracted from four sites near Barrow, Alaska. A high-resolution (25 cm) LiDAR (light detection and ranging) digital elevation model (DEM) was used to characterize the microtopographic features of the landscape, which were characterized and represented as a high-resolution grid. Painter et al. (2016) developed the ATS model by combing a surface energy balance model to simulate the surface water flows and snow distribution in the microtopography. The tilted open book catchment was used to conduct fully integrated thermal hydrology simulations, and a cross section across a low-centered, ice-wedge polygon in Barrow, Alaska, in a fine-scale 100-years projected future warming climate was also presented. Schilling et al. (2019b) improved the GW-SW flow model HydroGeoSphere (HGS) by adding a surface flow model with snowmelt and the degree-day model in an attempt to simultaneously balance the computational efficiency and accuracy. First, the analytical solution of the 1D Stefan heat transfer equation was simulated to verify the accuracy of the freeze-thaw parameterization scheme in the integrated GW-SW model, and then the classical 3D V-catchment benchmark case was calculated to determine the impact of the snow melt and FTCs on the GW-SW interactions. Furthermore, the model was directly used to simulate the GW-SW interactions of the Borden site in Canada, whose feasibility was evaluated based on the snow depth, snow water equivalent, and groundwater and surface water flow under different precipitation and temperature conditions. Evans et al. (2020) employed SUTRA in a 3D domain (550 m × 100 m × 60 m) with a horizontal spatial resolution of 3 m × 3 m and a vertical resolution of 20 cm. The model was representing a small catchment adopted from the Upper Kuparuk River watershed on the North Slope of Alaska, United States. The effect of surface temperature variation on groundwater discharge with water tracks was investigated. Jan et al. (2020) evaluated the GW-SW interaction in the ice-wedge polygon area of Alaska, United States, and analysed the effects of the microtopography on the TH processes using the ATS. The model was driven by on-site meteorological data with a 25 cm horizontal spatial resolution, and the simulation results exhibited a high degree of consistency with the snow depth, water table, and soil temperature observed in the Next Generation Ecosystem Experiments (NGEE) Arctic Program of the U.S. Department of Energy (US DOE), which demonstrated the accuracy of the fully-coupled GW-SW models in cold regions.

Watershed Scale

At the watershed scale, Wang et al. (2009) incorporated a frozen soil parameterization into a distributed hydrological model (WEB-DHM) and applied the model to the Binngou watershed (30.48 km²) of the upper reaches of the Heihe River Basin for evaluation with the in-situ observations. WEB-DHM was subsequently developed into WEB-DHM-SF by incorporating the enthalpy-based coupled snow and frozen ground physics (Qi et al., 2019; Song et al., 2020). Holmén et al. (2011) developed the GEOAN model to study the impact of climate and permafrost on groundwater in the entire Paris Basin (approx. 163,400 km²). This large-scale 3D hydrogeological model had the highest numerical resolution (333 m × 333 m) in the cross-cutting plane, and the vertical resolution ranged from a few meters up to more than 100 m. It included geological layers from the Triassic to the Tertiary (a total of 34 geological layers represented by 58 numerical layers). The model was calibrated using the actual TH conditions, and then, the future groundwater flow for the next 130,000 years was calculated. Zhang et al. (2013) used SHAWDHM in the Binggou watershed and the results showed that the model was able to predict soil freezing and thawing, unfrozen soil water content, and snow depth reasonably well. Yang et al. (2015) developed a distributed scheme (GBHM) for eco-hydrological simulation in the upper Heihe River, which has been subsequently developed into the distributed Geomorphology-Based Eco-hydrological Model (GBEHM) (Gao et al., 2016) and a physical process and machine learning combined hydrological model GBHM-ANN-CA-CV (Yang et al., 2020). The results identified that the soil freezing and thawing process would significantly influence the runoff generation. Yao et al. (2017) used regional groundwater model MODFLOW-NWT (Niswonger et al., 2011) to explore the partitioning between baseflow and mountain block recharge in the QTP. The model was first calibrated using mountain stream baseflows and then run with multiple scenarios in the Qilian Mountains. Carroll et al. (2019) comprehensively used high-quality field snow observations and the GSFLOW model to study the GW-SW interactions under seasonal and climate changes in the Copper Creek watershed (24 km²) of the upper Colorado River in the United States from 1987 to 2018. The driving and verification data for GSFLOW was obtained from two SNOTEL site observations (i.e., snow area, snow depth, snow water equivalent, during snow, and snowmelting periods), high-precision (3 m) LiDAR-based Airborne Snow Observatory data, and hydrometeorological data. Their results indicate that seasonal snowmelt and vegetation water consumption are the main factors controlling the hydrologic processes in the low-altitude areas of the basin, while the runoff is completely controlled by snowmelt in the high-altitude mountains. Cochand et al. (2019) used a modified the HGS model to assess the potential impact of climate change on the GW-SW interactions in the 553 km² Saint Charles River Basin in Quebec, Canada. First, the model was developed and calibrated to reproduce the observed streamflow and the hydraulic heads from 1970 to 2000, and then, it was used with three different temperature and precipitation scenarios to simulate the surface flow and groundwater dynamics from 2070 to 2100. The simulations predict that the warmer winters will cause more precipitation (rain and snow), which will increase the winter stream discharges and the water table. Conversely, the summer stream discharges will decrease due to an increase in the evapotranspiration, demonstrating that winter processes play a key role in the seasonal modifications anticipated for surface and subsurface flow dynamics. Langford et al. (2020) used the piFreeze plug-in to one-way couple the surface water heat balance SHAW model with the FEFLOW groundwater model based on the Scotty Creek Basin in Northwestern Territories of Canada. This study not only explored the importance of the initial conditions but also demonstrated the critical role that the variations in the microtopography of the land surface and permafrost table and talik development play in permafrost degradation. Based on the model CRHM developed by Pomeroy et al. (2007), Marsh et al. (2020) presented a new distributed model, the Canadian Hydrological Model (CHM). The model was further used in the Marmot Creek Research Basin of the Kananaskis River Valley of the Canadian Rockies and focused on the role of cold region processes in hydrology. Yao et al. (2021) used MODFLOW-NWT to evaluate the annual and seasonal contributions of groundwater discharge to total river flow in the Yarlung Zangbo basin on the northern Himalayan mountainous region. Results indicated that the groundwater recharge and discharge would increase with projected increasing ratio of recharge to precipitation, and thawing permafrost would increase groundwater recharge and discharge.

In summary, due to the complexity in simulating cold region hydrologic cycles and lack of high-quality data, although some cutting-edge research has emerged in recent years, many scientific issues are still unclear, and more advanced model developments are needed for quantitative analysis and accurate prediction. For now, due to computational limitations and complexities, 2D slope models have been more commonly employed than 3D catchment and watershed models. In addition, conceptual domains are often adopted instead of realistic terrains in cold regions, especially in the 2D slope scale. Also, the assumptions and setups of the models are quite simple (e.g., ignoring microtopography and heterogeneities, simplified subsurface aquifer layers, the assignment of hydraulic and thermal boundary conditions) that cannot accurately reflect the actual situations. Therefore, even the aforementioned models have improved our theoretical understanding of cold region hydrologic processes, but there are still limitations and uncertainties in solving practical problems (Lamontagne-Hallé et al., 2020). Future model development should focus on describing the dynamic processes on multiple scales (i.e., 3D groundwater flow, heat transfer, groundwater-frozen soil, snowmelt-groundwater, and GW-SW interactions), quantitatively analyzing the influencing mechanisms and factors, and predicting the temporal and spatial evolution of GW-SW interactions in cold regions under seasonal and climate change.

Challenges and Future Research Needs

GW-SW Interactions in Cold Regions

With projections of rapid and abrupt warming for the foreseeable future, surface water-groundwater interactions are considered as critical elements in cold-region hydrologic processes. As a key link in the hydrologic cycle in cold regions, the interactions between groundwater and surface water are unique, complex, and uncertain. The core issue is the dynamic processes of water and energy transport between the surface and the subsurface, especially in regions with snow cover and frozen soil. Therefore, the physics of surface water-groundwater interactions should be accurately represented in hydrologic models.

Future Model Development

However, the existing integrated hydrologic models of cold regions still have room to develop in several aspects:

a) The physical processes are not completely considered in the existing integrated hydrologic models of cold regions. In the distributed hydrologic models based on terrestrial hydrology, the surface hydrologic processes are well considered, but the groundwater processes are often simplified, and even lateral flow and heat transfer processes are directly neglected. In contrast, in the cryo-hydrogeological models, the interactions between soil and groundwater are accurately simulated, but the surface hydrologic processes are often described by simple parameterization or one-way coupling.

b) Compared with the previous models, the coupled GW-SW models of cold regions consider both the surface and subsurface hydrologic processes. Since incorporating more processes can potentially increase the uncertainty and computational time, the coupling algorithms and parallel computing efficiency of the coupled GW-SW models of cold regions still need to be improved. Moreover, most coupled GW-SW models of cold regions are not open source, which inhibits its further development and community applications.

c) In addition, various models have been developed and applied in recent years, while rigorous and systematic benchmark verifications are relatively rare. The physical processes and solution strategies behind various models vary among codes, and their accuracies and physical realism are unknown, which limits their practical application and prediction capabilities.

Perspectives

Therefore, it is vital to develop the next generation integrated GW-SW models. Not only should a complete description of GW-SW processes be ensured, but the integrated GW-SW models should achieve a better balance between accuracy and computational efficiency. Furthermore, to evaluate and optimize the present models of cold regions, it is urgent that more code validation and intercomparison efforts should be initiated across multi-disciplinaries. The next-generation GW-SW models are expected to fundamentally explain the hydrologic process in cold regions and address the challenges of climate change.

Summary

Integrated hydrologic models of cold regions have emerged as promising tools to simulate the impacts of pronounced climate warming and seasonal changes on the groundwater-surface water (GW-SW) interactions and groundwater resources. This review revealed the current state-of-the-art integrated hydrologic models of cold regions and subsequently prospected the future development. We first introduced different types of numerical models and explained the cons and pros of simulating GW-SW interactions if using these models. Selected intercomparison project and benchmarking cases for testing the accuracy of the numerical models were also reviewed. It should be highlighted that the research on rigorous and systematic benchmark validations should be largely promoted for benefiting the next generation models. Recent research on modelling GW-SW interactions in cold regions were summarized from hillslope scale to watershed scale, with aspects in temporal and spatial changes of the hydrologic cycle, and quantitatively analyzing groundwater lateral flow and heat transfer, local-scale permafrost heterogeneity, micro-surface roughness, the impact of snow accumulation and melting and frozen soil on the GW-SW interactions in cold regions, and the prediction under seasonal and climate changes. In conclusion, more efforts are absolutely needed to further develop integrated hydrologic models, specifically for applications in cold regions, as suggested in future research perspectives.

Author Contributions

XY designed the research; JH and XY wrote the initial draft of the manuscript; RM, ZS, XY, and JH revised the manuscript; all authors agreed to its submission.

Funding

The current research is supported by the National Natural Science Foundation of China (NSFC, Grant No. 42077172, 41730854), the Second Tibetan Plateau Scientific Expedition and Research Program (STEP, Grant No. 2019QZKK0306), and the Strategic Priority Research Program of Chinese Academy of Sciences (Grant No. XDA20100104).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Abbott, M. B., Bathurst, J. C., Cunge, J. A., O'Connell, P. E., and Rasmussen, J. (1986). An Introduction to the European Hydrological System – Systeme Hydrologique Europeen, “She”. 1: History and Philosophy of a Physically-Based, Distributed Modeling System. J. Hydrol. 87 (1-2), 45–59. doi:10.1016/0022-1694(86)90114-9