A Review of Hydrological Models Applied in the Permafrost-Dominated Arctic Region

Abstract

The Arctic region is the most sensitive region to climate change. Hydrological models are fundamental tools for climate change impact assessment. However, due to the extreme weather conditions, specific hydrological process, and data acquisition challenges in the Arctic, it is crucial to select suitable hydrological model(s) for this region. In this paper, a comprehensive review and comparison of different models is conducted based on recently available studies. The functionality, limitations, and suitability of the potential hydrological models for the Arctic hydrological process are analyzed, including: (1) The surface hydrological models Topoflow, DMHS (deterministic modeling hydrological system), HBV (Hydrologiska Byråns Vattenbalansavdelning), SWAT (soil and water assessment tool), WaSiM (water balance simulation model), ECOMAG (ecological model for applied geophysics), and CRHM (cold regions hydrological model); and (2) the cryo-hydrogeological models ATS (arctic terrestrial simulator), CryoGrid 3, GEOtop, SUTRA-ICE (ice variant of the existing saturated/unsaturated transport model), and PFLOTRAN-ICE (ice variant of the existing massively parallel subsurface flow and reactive transport model). The review finds that Topoflow, HBV, SWAT, ECOMAG, and CRHM are suitable for studying surface hydrology rather than other processes in permafrost environments, whereas DMHS, WaSiM, and the cryo-hydrogeological models have higher capacities for subsurface hydrology, since they take into account the three phase changes of water in the near-surface soil. Of the cryo-hydrogeological models reviewed here, GEOtop, SUTRA-ICE, and PFLOTRAN-ICE are found to be suitable for small-scale catchments, whereas ATS and CryoGrid 3 are potentially suitable for large-scale catchments. Especially, ATS and GEOtop are the first tools that couple surface/subsurface permafrost thermal hydrology. If the accuracy of simulating the active layer dynamics is targeted, DMHS, ATS, GEOtop, and PFLOTRAN-ICE are potential tools compared to the other models. Further, data acquisition is a challenging task for cryo-hydrogeological models due to the complex boundary conditions when compared to the surface hydrological models HBV, SWAT, and CRHM, and the cryo-hydrogeological models are more difficult for non-expert users and more expensive to run compared to other models.

1. Introduction

1.1. Extreme Global Climate Change in the Arctic Region

Global climate change (GCC) is more intensive in the Arctic region than in other parts of the world [1,2]. The annual average temperature in the Arctic has increased at twice the rate of that in the rest of the world since 1980 [1]. Since 2005, the surface air temperature in the Arctic has been higher than any recorded five-year period from 1880 [1]. During the 1971–2017 period, the annual mean air temperature in the Arctic increased by around 2.7 °C [2]. It is predicted that the air temperature in autumn and winter in the Arctic region will continue to increase by around 4 °C over the next two decades [2]. Intensive climate change in the Arctic has a significant impact on the hydrological processes in the region.

1.2. The Presence of Permafrost and Its Relation to Hydrological Processes in the Arctic Region

Permafrost accounts for approximately 24% of the exposed land area in the Northern Hemisphere [3]. Permafrost regions have different hydrology and hydrology-related conditions compared to non-permafrost regions [4,5]. Permafrost can affect many hydrological processes in Arctic and sub-Arctic environments [6], for example, surface and subsurface water fluxes [4,7,8,9]. Especially, permafrost accelerates the initiation of runoff [10] and shortens the time of response to rainfall [4]. Additionally, the evapotranspiration process from vegetation is limited in permafrost environments since permafrost prevents the downward growth of roots and therefore limits the uptake of water for evapotranspiration [4,7]. Other processes, including microclimatology and thermal regimes (related to the evapotranspiration process) [11,12], water storage processes [4,13,14], and energy and water balances [13,14] are also affected by permafrost. The thin soil layer overlying permafrost is the active layer that seasonally freezes and thaws [4,5]. This active layer in the Arctic varies from several centimeters to one or two meters in depth [6] and most of the hydrological and biogeochemical processes occur in this layer [7,15]. The active layer determines the conditions for plant growth, gas fluxes, groundwater flow regimes, and soil formation [13]. While the active layer supports many hydrological and biogeochemical processes, permafrost beneath the active layer limits the amount of soil water percolation and subsurface storage of water [16]. Unlike non-permafrost soils, where a groundwater system is available and deep, the subsurface movement of water in permafrost-affected soil is mostly confined to the shallow active layer [6]. Because of this limitation of the vertical movement of subsurface flow in permafrost soil, subsurface flow in the horizontal direction is therefore important [6,17]. Moreover, the presence of permafrost has influenced the magnitude of the specific base flow in some Arctic basins. According to the study by McNamara et al., the specific base flow in a permafrost basin is lower than that in a non-permafrost basin [10].

1.3. The Impacts of Permafrost Thawing on the Arctic Hydrological Processes

Climate change is expected to alter the hydrological processes in the Arctic [6,18], e.g., through thawing of the permafrost layer, which has been observed from the field measurement data obtained during the last decades [19,20,21,22,23,24]. Global warming is expected to lead to permafrost degradation through changes in the three-dimensional distribution of permafrost (Figure 1), e.g., via changes in the active layer thickness (ALT) [8,25,26,27,28], spatial extent [26], open vertical taliks [29], and lateral taliks [30]. The thawing of permafrost is expected to change both surface and subsurface hydrology in the Arctic [31]. Permafrost thawing leads to alterations in: (1) Water fluxes and flow paths, including increasing soil drainage, increasing suprapermafrost flow (SUPRA-PF) [32], runoff variation [33], increasing evapotranspiration (ET), increasing the exchange of water flow between lakes and groundwater [34], increasing subpermafrost flow (SUB-PF), and increasing baseflow [35,36,37,38]; (2) secondly, water storage and ecosystem responses including water distribution, e.g., variation in soil moisture [39], lakes and wetlands [40,41], groundwater (GW) storage [42,43]), aufeis (icing) volume [44], decreasing winter river ice thickness [45], and ecosystem variables, e.g., increasing vegetation [46,47], variation in surface water connectivity [48], increasing subsurface connectivity, decreasing streamflow, and temperature seasonality [49]; and (3) also alterations in greenhouse gasses (GHGs) and the surface energy flux [50]. It is projected that, because of permafrost degradation under global climate change, surface water systems will be transferred to more groundwater-based systems in large-scale assessments of the Arctic region [31].

1.4. Importance of Choosing the Suitable Modeling Tools for the Arctic Region

It has been demonstrated that the presence of permafrost has a high impact on hydrological processes in the Arctic. The changes of Arctic hydrology (including surface and subsurface hydrology) are expected to be more complicated under the context of global climate change, which causes permafrost degradation and changes in other related processes. Hydrological models are state-of-the-art tools for the investigation of global climate change impacts, and numerous models have been developed in recent decades. Both surface and subsurface hydrological models have demonstrated their capacities to simulate permafrost hydrology. For example, many surface hydrological models have approached analytical solutions (by using a simple heat transfer equation, e.g., Stefan’s equation) and numerical solutions (e.g., finite difference, finite element, and finite volume methods) to simulate the seasonal freezing–thawing process. Many subsurface hydrological models have coupled a three-dimensional (3D) equation for water flow (e.g., the 3D Richards equation) and a three-dimensional equation for heat transfer, especially considering the three phase changes of water in near-surface soil. However, both surface and subsurface hydrological models still have their limitations when dealing with permafrost hydrology. For example, one-dimensional (vertical direction) heat transfer, used in surface hydrological models, cannot be used to simulate multidecadal and multidimensional changes. Additionally, surface hydrological models still lack important processes in permafrost environments, such as consideration of heat capacity, thermodynamic equilibrium, and the three phase changes (ice, liquid, and gas) of water in near-surface soils. Subsurface hydrological models do not feature land surface schemes in their structures and have complex boundary conditions, and they are also difficult for non-expert users to use, etc. Moreover, sparse data of the Arctic make it more difficult to collect the necessary input data for the models. Therefore, finding suitable models for the Arctic is a challenge for hydrological modelers. The current paper presents a review and analysis of the functions, advantages, and disadvantages of different hydrological models and evaluates their suitability for simulating hydrological processes in the Arctic region. The selection of suitable models is carried out via answering the following pertinent questions:

• 1.. Do the models consider the important processes in permafrost environments, including the following factors:

  • Surface energy balance;

  • Snow processes, snow insulation, and snow melt;

  • Infiltration processes;

  • The dynamics of soil thermal and soil moisture fluxes;

  • Soil heterogeneities;

  • The dynamics (seasonal thawing) of the active layer;

  • Subsidence;

  • A three-phase change of water (ice, liquid, and gas) during the freezing and thawing of near-surface soil.

• 2.. Can the models be widely applied for Arctic permafrost, particularly considering the following requirements:

  • Requirement for input data, i.e., large or small requirement;

  • Requirement for computation processes, i.e., strong or low requirement;

  • Ability to be applied with different sizes of watersheds, i.e., small-scale and/or large-scale.

2. Some Well-Known Hydrological Models Applied in the Arctic

2.1. Topoflow Model

Topoflow is a spatially distributed and process-based hydrological model that was firstly designed for Arctic and sub-Arctic basins [6]. Topoflow has the capacity to simulate hydrological processes in permafrost environments, including the surface energy balance, snowmelt (using simple degree day or full energy balance methods), infiltration (using Green–Ampt, Smith–Parlange, or 1D Richards’ equation with three layers), and volumetric soil moisture content (using Darcian theory with multiple uniform layers). Moreover, Topoflow is able to reasonably simulate the ALT with a relatively simple method. However, such a method should be further improved in order to analyze the dynamics of the active layer more accurately, since the active layer has a high impact on the hydrological processes in permafrost environments [51]. The simulation of subsidence and a three-phase change of water (ice, liquid, and gas) during the freezing–thawing process of the near-surface soil is not mentioned in the Topoflow model.

The Topoflow model was applied in a small-scale watershed in the Arctic, Imnavait Creek, (around 2 km²) in Alaska, United States, which was underlain by continuous permafrost [6]. The results showed that the model has good performance for simulation of the hydrological cycle, including evapotranspiration, snowmelt, infiltration, runoff, and energy balances in the Arctic [6]. The hydrology change of the Imnavait Creek watershed under different climate change scenarios was also simulated [6], including evaluation of the performance of the model. The model has some limitations, such as the spatial variability of the active layer’s depth (not presented in the model) and the simplification of the complex soil moisture heterogeneity, etc.

2.2. DMHS Model

The DMHS (deterministic modeling hydrological system) model was developed for both mountainous and flat topographies, including basins with different climate zones, regardless of their scales [52]. This model can be applied to any geographical area in the world [53]. The DMHS model is a kind of physically-based, semi-distributed model for runoff estimation. The model can simulate the important processes in permafrost-affected regions, such as the surface energy balance, snow accumulation, snowmelt, sublimation, infiltration, soil heterogeneity, heat and water dynamics, and phase changes in soil layers [52,53,54,55,56]. Additionally, DMHS can also simulate the dynamics of the ALT with reliable results [57].

Vinogradov et al. [52] applied the DMHS model to six mountainous watersheds of different sizes, with areas varying from 40 to 2.4 million km² across eastern Siberia and inside the Lena River basin in the Arctic region, showing promising results. However, the DMHS model has limitations regarding the routing scheme, which is especially not applicable to rivers which have a backwater phenomenon [58]. Another limitation of the DMHS model is the difficulty in the acquisition of soil profile properties in a required format for the model [52]. This could be a limiting factor for its wide application in the Arctic, where the available data are usually limited.

2.3. HBV Model

The HBV (Hydrologiska Byråns Vattenbalansavdelning) model is a rainfall–runoff model that was first successfully applied in 1972 [59,60,61,62,63]. The HBV model has developed significantly since then, and nowadays it is considered to be a standard tool for an increasing number of applications, such as flood forecasting, forecasts of inflow to reservoirs of hydropower dams, assessment of the impacts of climate change on water resources, or for the simulation of hazards in designing high hydropower dams, etc. [64,65,66,67]. The HBV model can be used as a semi-distributed, conceptual model. The HBV model considers the important processes of permafrost hydrology, such as the surface energy balance, snow routine, snow melt, soil moisture, and infiltration. Especially, the model uses an accumulated degree day coefficient, which is set up based on field measurement, to simulate the ALT [68]. In some case, the HBV model is coupled with other thermal hydrological models, e.g., the CryoGrid model, for simulation of the active layer dynamics [67].

The HBV model has been applied in more than 40 countries with different climate conditions [69]. In Nordic countries, the HBV model is currently used for flood forecasting and several other purposes, such as the simulation of design flood for supporting the design for spillway structures [70], water resource evaluation [71,72], and nutrient load estimates [73]. In Norway, the HBV model has been applied since 1983 [74] and has recently been considered as an important tool for water management by NVE (the Norwegian Water Resources and Energy Directorate). The model has been used for flood forecasting under the impacts of climate change [75].

The advantage of the HBV model is that it requires less input data and computer facilities, which is suitable for the sparse data in the Arctic region. Especially, the model could be a potential tool for hydrological simulation in ungauged river basins. However, obtaining the optimum model parameters for the HBV is not an easy task [64]. In addition, HBV model simplifies the equation for water storage in the catchment. Moreover, the HBV model requires the most basic input data, i.e., only daily precipitation and average daily air temperature, and these data are spatially averaged over the watersheds by Thiessen polygons. However, the Thiessen polygon method does not consider the elevation and temperature gradients in the basin. Hence, the interpolation method used in the HBV model could result in unsatisfactory modeling results.

2.4. SWAT Model

The SWAT (soil and water assessment tool) model [76] is a robust watershed modeling tool which was designed by the USDA (United States Department of Agriculture) Agricultural Research Service, the USDA Natural Resources Conservation Service, and the Texas A&M University. The model is currently used in more than 100 countries around the world. The SWAT is a physically-based, semi-distributed, watershed-scale model which operates on a daily time step [77,78]. The model was developed to simulate the impacts of land use management activities on water resources, sediment, and agricultural chemical yields [79] in large, complex, and ungauged watersheds with variation in soils, land uses, and management conditions over long periods of time [80,81]. Moreover, the model supports the understanding of complex ecosystems, climate change, and agricultural production issues around the world [82,83,84,85,86,87,88,89,90]. The SWAT model is capable of modeling river basins for thousands of square miles. Especially, it is able to simulate watersheds without requiring monitoring data. Regarding permafrost hydrology, the SWAT is able to simulate the surface energy balance, snow cover, snow melt, infiltration, soil moisture, soil heterogeneity, and lateral subsurface flow from the soil profile. Additionally, the SWAT takes into account the presence and development of the ALT in the model’s structure but only via using the average values [91]. Therefore, further development in the model’s structure to obtain a higher resolution for the spatiotemporal variation of the ALT is recommended, or otherwise coupling SWAT with other models to return better modeling results for the active layer dynamics. Moreover, simulation of soil heat transfer, snow insulation effects, and subsidence has not been developed yet for the SWAT model.

The SWAT model has been applied in many river basins in Europe and some regions in the Arctic. For instance, the SWAT model was applied to assess the variation of the hydrological regime and water quality under the impacts of human activities and climate change for the whole European territory [92]. In another study, the SWAT model was coupled with carbon modules to assess organic carbon exportation in an Arctic watershed via examination of the Yenisei River [91]. The advantage of the SWAT model is its capacity to model the temporal and spatial variations of hydrological process in large-scale Arctic watersheds [91]. The NRCS (Natural Resources Conservation Service) runoff curve number method in the SWAT provides a relatively easy way for the model to be adapted to a wide range of hydrological conditions [79].

2.5. WaSiM

The WaSiM (water balance simulation model) is a deterministic, spatially distributed hydrological model that was firstly developed during 1994-1996 by Jörg Schulla [93]. The model was designed to simulate water cycles in both surface and subsurface areas in river basins with variations in the spatial and temporal scales [93]. The WaSiM is able to simulate the hydrological processes for the catchments with sizes varying from <1 up to >100,000 km² [93]. The WaSiM is able to solve several issues in river basins, including assessing the impacts of climate change and land use change on water resources, runoff forecasting, groundwater recharge, soil water, substance transport, etc. [93]. The model runs with time steps from minutes to several days [93]. The WaSiM can be used for both short-term (e.g., floods) and long-term (e.g., water balance) simulations [93]. The WaSiM is able to simulate the important hydrological processes in Arctic permafrost, including snow accumulation, snow melt, infiltration, soil moisture, and soil heterogeneities. Additionally, the model can simulate the ALT dynamics in permafrost regions. However, such simulation is relatively simple and is only based on empirical approaches. Particularly, the WaSiM calculates the thaw depth based on the simple formula given as follows [94]:

(1)${\mathrm{d}}_{\mathrm{thaw}}=\mathsf{\alpha } \sqrt{{\mathrm{n}}_{\mathrm{sf}}},$

where:

• ${\mathrm{d}}_{\mathrm{thaw}}$ is the thaw depth (m);

• $\mathsf{\alpha }$ is an empirical coefficient (~0.02, …, 0.05);

• ${\mathrm{n}}_{\mathrm{sf}}$ is the number of snow-free days.

Since the calculation approach in Equation (1) is relatively simple, in order to simulate the freezing–thawing process of the active layer more accurately, a physically-based heat transfer model is recommended [94].

The model is able to deal with soil heat flux. It provides the variations of soil temperatures in three dimensions. In addition, the WaSiM considers phase changes, which allows the model to simulate the freezing and thawing of an active layer in permafrost environments [95]. However, the drawbacks of the WaSiM is its high sensitivity to temporal and spatial resolutions. Therefore, the model is not suitable for transfer to other scales without recalibration. Moreover, the WaSiM requires intensive input data and model parameters for each grid cell of the model domain since it is a fully distributed model. This is considered as a big challenge in light of the scarce data in the Arctic region.

The WaSiM was first applied in the Arctic watershed by Liljedahl [93]. The long-term water balance during 1999–2009 was evaluated and its changes under global climate change in the permafrost-dominated watersheds were projected, such as the two vegetated drained thaw lake basins near Utqiaġvik, formerly known as Barrow, in the north of Alaska. In another study, WaSiM was applied to assess the influence of permafrost degradation on the Arctic tundra hydrology [96].

2.6. ECOMAG Model

The ECOMAG (ecological model for applied geophysics) model is a physically-based, distributed hydrological model for the simulation of hydrological cycles and water quality transformation in catchments in cold climate regions [97,98]. The model includes two separate submodels, e.g., a hydrological submodel and a water quality submodel, which operate at a daily time step. The hydrological submodel describes several processes occurring in the catchments, including surface runoff, evapotranspiration, infiltration into soil layers, soil moisture, and subsurface flow. It is able to simulate the hydrothermal processes, which are important in permafrost regions, including the formation of snow cover, snowmelt rate (using the degree day method), ALT dynamics, infiltration of snowmelt into the unfrozen and frozen soils by integrating the governing equations of basic hydrodynamic and thermodynamic of water, heat vertical transfer, horizontal water flow, etc. [99]. The water quality submodel simulates the pollution transformation process from point sources and non-point sources in the catchment, including geochemical processes and the biochemical degradation process of dissolved organic pollutants.

The ECOMAG model was tested for hydrological simulation in numerous river basins in cold climate regions such as Canada, Russia, Norway, and Sweden. The model has been applied in the large Arctic basins, including the Mackenzie River basin in Canada, with satisfactory performance [100], and the Lena River basin in Russia, also with satisfactory performance [101]. Additionally, the ECOMAG was approached to investigate climate change impacts on the water regimes of those river basins [102,103].

2.7. CRHM Model

The CRHM (cold regions hydrological model) is a physically-based, semi-distributed model developed for hydrological cycle simulation over small to medium river basins in cold climate regions [104]. The CRHM is a flexible, object-oriented modeling system with the capacity to simulate a wide range of important permafrost-related processes in the cold regions, including snow processes (e.g., snow redistribution by wind, snow interception, sublimation, snowmelt, and infiltration of snowmelt into unfrozen and frozen soils), glacier melts, actual evaporation and evapotranspiration, radiation exchange, and soil moisture balance, etc. The CRHM is able to simulate the ALT dynamics by approaching Stefan’s heat flow equation [105] as follows:

(2)$\mathsf{\xi }=\sqrt{\frac{2\mathrm{kF}}{\mathrm{Lw}\mathsf{\rho }}},$

where:

• $\mathsf{\xi }$ is the frost/thaw front depth (m);

• k is the thermal conductivity of the soil (W m⁻¹ k⁻¹);

• F is the surface freeze/thaw index (°C degree days);

• L is the latent heat of fusion (J kg⁻¹);

• w is volumetric water content (m³ m⁻³);

• $\mathsf{\rho }$ is the bulk density of the soil (kg m⁻³).

The CRHM is flexible in terms of its spatial solutions (from lumped to distributed) and model structures (from conceptual to physically-based), depending on the objectives of the studies and available input data for catchments. This is considered as a benefit for use in the Arctic region, which features sparse data. However, there is no required calibration for the model, and the model parameters are normally obtained based on the expert knowledge of the modeled catchment. Therefore, the modeling results have high uncertainty.

2.8. ATS Model

The ATS (arctic terrestrial simulator) model, which was developed from Amanzi code [106], is an integrated tool of the permafrost-related process model and the physically-based model [107]. The model couples the surface energy balance model [108] and the three-dimensional (3D) subsurface thermal hydrology model [109,110] for multidimensional simulations in permafrost-affected regions. The model employs the diffusion wave equation to simulate the surface hydrology, energy transport, and phase change. However, simulation of the dynamic topography due to the influence of thaw-induced subsidence has not been developed yet.

The ATS model is an open source model and one of the first tools for addressing fully-coupled surface/subsurface permafrost thermal hydrology in multidimensional simulations, which is considered as the existing challenge for hydrological models in permafrost environments. By approaching a novel multiphysics management system [111] to increase the complexity of the model, the ATS model can meet the challenge of operating many numerical models in permafrost regions, i.e., performing nonlinear constitutive modeling, phase change modeling, the coupling of several processes in different spatial domains, the ability for simulating transitions among different states on the land surface, and the possibility for using an unstructured grid according to the given topographical features. Moreover, the ATS model has a fine-scale structure and is suitable for the specific topography of the polygonal Arctic tundra, which requires a high spatial resolution, large spatial model domain, and large grid size because of long-range surface flow. Therefore, the incorporation of the extension ATS code with good parallel scaling of the Amanzi code becomes a potential tool for investigating the responses of small-scale topographic features in permafrost environments in the context of global climate change. Moreover, the ATS is able to simulate comprehensive snow processes, such as an increasing snow density with snow age, snow insulation, and snow thermal conduction [107].

However, the ATS model has limitations regarding (1) solving the convergence of nonlinear systems around the transition between freezing and thawing [112] and (2) simulating topography change by subsurface ice melting requiring the movement of grid/mesh cells but having to maintain the water and energy balance inside the moving grid cells.

2.9. CryoGrid 3 Model

The CryoGrid 3 model is a new, simple, and one-dimensional land surface model developed to simulate ground surface temperatures in permafrost-affected regions. The model was built based on the thermal permafrost model CryoGrid 2 [113] by developing further calculation of the surface energy balance and modifying the snow scheme, which are important factors in permafrost hydrology processes. The CryoGrid 3 model is considered as the upper boundary condition of the CryoGrid 2 model. The surface energy balance describes the processes of energy transfer between the atmosphere and the ground. Such processes include the radiation balance, the exchange of sensible heat, evaporation, etc. [114], as described in the following equation:

(3)$\frac{\partial \mathrm{E}}{\partial \mathrm{t}}={\mathrm{S}}_{\mathrm{in}}+{\mathrm{S}}_{\mathrm{out}}+{\mathrm{L}}_{\mathrm{in}}+{\mathrm{L}}_{\mathrm{out}}+{\mathrm{Q}}_{\mathrm{h}}+{\mathrm{Q}}_{\mathrm{e}}+{\mathrm{Q}}_{\mathrm{g}},$

where:

• ${\mathrm{S}}_{\mathrm{in}},{\mathrm{S}}_{\mathrm{out}}$ are the short-wave radiation input and output, respectively (W m⁻²);

• ${\mathrm{L}}_{\mathrm{in}},{\mathrm{L}}_{\mathrm{out}}$ are the long-wave radiation input and output, respectively (W m⁻²);

• ${\mathrm{Q}}_{\mathrm{h}},{\mathrm{Q}}_{\mathrm{e}},{\mathrm{Q}}_{\mathrm{g}}$ are the sensible, latent, and ground heat fluxes, respectively (W m⁻²).

The CryoGrid 3 model is able to simulate other important processes in permafrost-affected regions, including subsurface heat transfer (considering the phase change of soil water), energy transfer, and mass balance of the snowpack. The subsurface thermal scheme is described by the concept of conductive heat transfer via Fourier’s law as the following equation:

(4)${\mathrm{c}}_{\mathrm{eff}}\left(\mathrm{z},\mathrm{T}\right)\frac{\partial \mathrm{T}}{\partial \mathrm{t}}-\frac{\partial }{\partial \mathrm{z}}\left(\mathrm{k}\left(\mathrm{z},\mathrm{T}\right)\frac{\partial \mathrm{T}}{\partial \mathrm{z}}\right)=0,$

where:

• ${\mathrm{c}}_{\mathrm{eff}}\left(\mathrm{z},\mathrm{T}\right)$ is the effective volume capacity (J m⁻³ K⁻¹);

• $\mathrm{k}\left(\mathrm{z},\mathrm{T}\right)$ is the thermal conductivity (W m⁻¹ K⁻¹).

The energy transfer within the snowpack is conductive heat transfer, which is similar to the soil, as in the following equation:

(5)${\mathrm{c}}_{\mathrm{snow}}\left(\mathrm{z},\mathrm{T}\right)\frac{\partial {\mathrm{T}}^{*}}{\partial \mathrm{t}}-\frac{\partial }{\partial \mathrm{z}}\left({\mathrm{k}}_{\mathrm{snow}}\left(\mathrm{z},\mathrm{T}\right)\frac{\partial {\mathrm{T}}^{*}}{\partial \mathrm{t}}\right)=0,$

where:

• ${\mathrm{c}}_{\mathrm{snow}}\left(\mathrm{z},\mathrm{T}\right)$ is the snow heat capacity (J m⁻³ K⁻¹);

• ${\mathrm{k}}_{\mathrm{snow}}\left(\mathrm{z},\mathrm{T}\right)$ is the thermal conductivity of the snow (W m⁻¹ K⁻¹);

• ${\mathrm{T}}^{*}$ is the snow temperature (°C).

The CryoGrid 3 model is run in the MATLAB programming environment. The model has simple code and is open for modification. Therefore, the CryoGrid 3 model can be considered as a platform to integrate further processes in permafrost environments.

The CryoGrid 3 model has been tested in Arctic conditions with the large Lena River basin and showed satisfactory results to simulate the surface temperature, surface energy balance, ground temperature, ALT, and ground subsidence [115]. However, it is not guaranteed that the model could perform well in other permafrost basins. Instead, further considerations should be made before applying the model to a wider range of permafrost-affected regions. The CryoGrid 3 model has three major challenges [115]: (1) The model considers relatively simple snow processes (e.g., the assumption of a constant density of snowfall); (2) secondly, it is unclear that the model could perform well with the simulation of energy transfer and ground thermal regimes in regions with high vegetation cover; and, finally, (3) the simulation of the water balance needs further improvements, particularly for the simulation of seasonal changes.

2.10. GEOtop Model

The GEOtop model is a physically-based, distributed hydrological model that couples the water and energy balance [116,117]. The model was developed specifically for small catchments and complex mountain terrains. The GEOtop model combines the strengths of both land surface models and flood forecasting models [117]. In the GEOtop model, the interaction of topography with radiation is treated in detail, which is normally not considered in many hydrological models. The model simulates not only the energy balance (e.g., evapotranspiration and heat transfer) but also the water cycle (e.g., cycles of water, snow, and glaciers). The energy and mass balance is calculated based on a 3D Richards equation [118] and a 1D energy equation [112]. Vegetation, which contributes to the turbulent fluxes, is calculated based on a double layer scheme [119]. The snow processes (e.g., snow accumulation and snowmelt) are simulated by a multilayer discretization of the snowpack [120]. Additionally, a blowing snow module [121,122] is also included in the model to calculate the accumulation and blowing of snow because of wind.

The GEOtop model has been applied in a wide range of studies, including studies of soil water content [123], evaporation from soil [124,125], transpiration from vegetation [119], snow in basins [120,126,127], surface temperature [125], the temperature of soil and rock under freezing conditions [128], the mass balance of glaciers [129], interactions between the ground water table and thaw depth [118], and discharge at basin outlets [116].

The GEOtop model is a fully distributed model and requires intensive inputs, model parameters, field measurements, and experiments, as well as intensive computation. Additionally, modification of the code in the GEOtop model is not allowed, unlike other cryo-hydrogeological models (e.g., CryoGrid 3 or SUTRA-ICE).

2.11. SUTRA-ICE Model

The SUTRA-ICE [130] model was developed by modification of the existing SUTRA (saturated/unsaturated transport) numerical model [131,132] in order to simulate the processes of subsurface ice formation and melting, which is important for heat transport, groundwater, and biological activities in permafrost environments [130]. Basically, the SUTRA model is a finite element numerical model for the simulation of saturated/unsaturated groundwater flow and solute energy transport. Regarding the working mechanisms, the SUTRA model approaches the problem via a two-dimensional hybrid method (e.g., finite element and finite difference methods) to describe two interdependent processes, including: (1) Fluid density-dependent saturated or unsaturated groundwater flow; and (2) the transport of a solute or thermal energy in groundwater flow, as well as in the solid matrix of the aquifer. Especially, the SUTRA model has the capacity to model the thermal regimes, thermal energy storage, and thermal pollution in aquifers, and additionally subsurface heat conductivity, geothermal reservoirs, and natural hydrogeological convection. In addition, the SUTRA model has been applied to investigate the impacts of climate change on permafrost thawing, ALT dynamics, and groundwater flow in cold climate regions [133,134,135,136,137,138]. Besides the functions existing in the SUTRA version, the modified SUTRA-ICE (with 2-dimensional and 3-dimensional versions) model can deal with the subsurface, saturated/unsaturated freezing processes (including phase changes), latent heat, permeability, heat capacity, thermal conductivity, and liquid porosity [139]. The SUTRA-ICE model has demonstrated its high capacity to simulate the formation and melting of near-surface ice and subsurface temperature distribution in cold climate region in previous studies [140]. The SUTRA-ICE model is a free model (applying for version 3.0) and it allows users to modify the code for further development. This is considered as an advantage of the model. However, the model has some challenges [139]: (1) The boundary conditions are complex; (2) it is still a problem to simulate freezing in unsaturated zones; (3) the model is not accurate for the simulation of active layer dynamics when the pore space is filled with both liquid and ice; (4) intensive computation is required; (5) the model is not able to utilize massive parallel computing hardware [109]; and (6) field verification is needed.

2.12. PFLOTRAN-ICE Model

The PFLOTRAN-ICE is a non-isothermal, single component (water), three-phase (ice, liquid, and gas) numerical model that has been recently developed to simulate subsurface hydrology in permafrost-affected regions [141]. The model calculates the balance of mass and energy for water components in three phases (ice, liquid, and gas) via the following equations [142]:

(6)$\frac{\partial }{\partial \mathrm{t}}\left[\varnothing \left({\mathrm{s}}_{\mathrm{g}}{\mathsf{\eta }}_{\mathrm{g}}{\mathrm{X}}_{\mathrm{w}}^{\mathrm{g}}+{\mathrm{s}}_{\mathrm{l}}{\mathsf{\eta }}_{\mathrm{l}}+{\mathrm{s}}_{\mathrm{i}}{\mathsf{\eta }}_{\mathrm{i}}\right)\right]+\nabla .\left[{\mathrm{v}}_{\mathrm{l}}{\mathsf{\eta }}_{\mathrm{l}}\right]-\nabla .\left[\varnothing {\mathrm{s}}_{\mathrm{g}}{\mathsf{\tau }}_{\mathrm{g}}{\mathsf{\eta }}_{\mathrm{g}}{\mathrm{D}}_{\mathrm{g}}\nabla {\mathrm{X}}_{\mathrm{w}}^{\mathrm{g}}\right]={\mathrm{Q}}_{\mathrm{W}},$

(7)$\frac{\partial }{\partial \mathrm{t}}\left[\varnothing \left({\mathrm{s}}_{\mathrm{l}}{\mathsf{\eta }}_{\mathrm{l}}{\mathrm{U}}_{\mathrm{l}}+{\mathrm{s}}_{\mathrm{g}}{\mathsf{\eta }}_{\mathrm{g}}{\mathrm{U}}_{\mathrm{g}}+{\mathrm{s}}_{\mathrm{i}}{\mathsf{\eta }}_{\mathrm{i}}{\mathrm{U}}_{\mathrm{i}}\right)+\left(1-\varnothing \right){\mathsf{\rho }}_{\mathrm{r}}{\mathrm{c}}_{\mathrm{r}}\mathrm{T}\right]+\nabla .\left[{\mathrm{v}}_{\mathrm{l}}{\mathsf{\eta }}_{\mathrm{l}}{\mathrm{H}}_{\mathrm{l}}\right]-\nabla .\left[\mathrm{k}\nabla \mathrm{T}\right]={\mathrm{Q}}_{\mathrm{e}},$

(8)${\mathrm{v}}_{\mathrm{l}}=\frac{{\mathrm{k}}_{\mathrm{rl}}\mathrm{k}}{{\mathsf{\mu }}_{\mathrm{l}}}\nabla \left[{\mathrm{p}}_{\mathrm{l}}+{\mathsf{\rho }}_{\mathrm{l}}\mathrm{gz}\right],$

where:

• Subscripts l, g, and i, are the liquid, gas, and ice phases, respectively;

• ∅ is the porosity (-);

• ${\mathrm{s}}_{\mathrm{l}},{\mathrm{s}}_{\mathrm{g}},{\mathrm{s}}_{\mathrm{i}}$ (constraint: ${\mathrm{s}}_{\mathrm{l}}+{\mathrm{s}}_{\mathrm{g}}+{\mathrm{s}}_{\mathrm{i}}$ = 1) are saturation indices of the liquid, gas, and ice phases, respectively (m³ m⁻³);

• ${\mathsf{\eta }}_{\mathrm{l}},{\mathsf{\eta }}_{\mathrm{g}},{\mathsf{\eta }}_{\mathrm{i}}$ are the molar densities of the liquid, gas, and ice phases, respectively (kmol m⁻³);

• ${\mathrm{X}}_{\mathrm{w}}^{\mathrm{g}}$ is the mole fraction of H₂O in the gas phase (-);

• ${\mathsf{\tau }}_{\mathrm{g}}$ is tortuosity of the gas phase (-);

• ${\mathrm{D}}_{\mathrm{g}}$ is the diffusion coefficient in the gas phase (-);

• T is the temperature (it is assumed that all the phases and soil are in thermal equilibrium) (K);

• ${\mathrm{c}}_{\mathrm{r}}$ is the heat capacity of the soil (J K⁻¹);

• ${\mathsf{\rho }}_{\mathrm{r}}$ is the density of the soil (kg m⁻³);

• ${\mathrm{U}}_{\mathrm{l}},{\mathrm{U}}_{\mathrm{g}},{\mathrm{U}}_{\mathrm{i}}$ are the molar internal energies of the liquid, gas, and ice phases, respectively (kJ mol⁻¹);

• ${\mathrm{H}}_{\mathrm{l}}$ is the molar enthalpy of the liquid phase (kJ mol⁻¹);

• ${\mathrm{Q}}_{\mathrm{W}}$ is the mass source of H₂O (kmol m⁻³ s⁻¹);

• ${\mathrm{Q}}_{\mathrm{e}}$ is the heat source (kmol m⁻³ s⁻¹);

• $\nabla ()$ is the gradient operator (-);

• $\nabla . ()$ is the divergence operator (-);

• ${\mathrm{v}}_{\mathrm{l}}$ is Darcy velocity of the liquid phase (m s⁻¹);

• ${\mathrm{k}}_{\mathrm{rl}}$ is the relative permeability of the liquid phase (-);

• $\mathrm{k}$ is the absolute permeability (m²);

• ${\mathsf{\rho }}_{\mathrm{l}}$ is mass density of the liquid phase (kg m⁻³);

• ${\mathrm{p}}_{\mathrm{l}}$ is the partial pressure of the liquid phase (Pa);

• $\mathrm{g}$ is acceleration because of gravity (m s⁻²);

• $\mathrm{z}$ is the vertical distance from a reference datum (m).

Together with the above balance equations, further requirements for the constitutive relationships (e.g., the mole fraction of water vapor, saturations of the phases, thermal conductivity, relative permeability, and water vapor diffusion coefficient) for the simulation of non-isothermal elements and the multiple phases of water are also described in the PFLOTRAN-ICE model.

The PFLOTRAN-ICE model makes use of the capacity in the PFLOTRAN code [143] for finding highly-scalable, parallel subsurface multiphysics solutions. Therefore, the model can simulate the degradation of ice wedge polygon bogs, which requires elucidation of three-phase change and relatively large model domains [144]. Although the PFLOTRAN-ICE model considers one component (water), the model is able to produce similar results to the more complicated two-component model [142] for the same application. This could be an advantage for the simulation of permafrost hydrology by using fewer demanding components (e.g., a single component of water substance). Moreover, the PFLOTRAN-ICE model is suitable for the large-scale range of model domains (i.e., kilometer scales).

However, the PFLOTRAN-ICE model does not consider surface flows, the surface energy balance, and topography dynamics because of permafrost thawing and the melting of ground ice. Such processes are expected to be coupled with subsurface hydrology for the comprehensive modeling of hydrology in permafrost-dominated regions [144].

2.13. The Present Capacities and Challenges of Hydrological Models to Deal with Permafrost Hydrology in the Arctic

2.13.1. Surface Hydrological Models

Surface hydrological models have demonstrated their capacity to simulate the seasonal freeze–thaw process in soil by approaching analytical and numerical solutions [25]. Regarding the analytical solutions, many physically-based models have incorporated a simple heat transfer equation (e.g., Stefan’s equation [145]) into their structure to deal with the soil freeze–thaw process. Because of the simplification of soil freeze–thaw algorithms, simulation for large-scale basins is an advantage [146,147]. Additionally, Stefan’s equation may be easily modified to include other processes, such as the temporal variation of soil moisture [148], spatial variation of the moisture content, thermal properties [149], the freezing–thawing process in two directions [150], heat advection [151], and the soil heat capacity [152]. Beside the analytical solutions, some surface hydrological models, especially distributed models [153,154,155,156], have approached numerical solutions (e.g., via finite difference, finite element, and finite volume methods) to simulate the ground freezing–thawing process. Such numerical approaches are able to deal with complex conditions (e.g., varying soil heterogeneities and complex temperature boundaries) or complex processes (e.g., coupling heat and water transfer, discontinuous freezing–thawing process, and the temporal variation of thermal factors). Compared to analytical solutions, numerical solutions perform better in the context of the freezing–thawing process in the ground [147]. For example, numerical solutions simulate the freezing and thawing process in soil over a wide range of temperatures, unlike the assumption of sharp change by analytical solutions. Therefore, such a freezing range permits some subsurface water flow at sub-zero temperatures.

However, there are still numerous challenges for surface hydrological models in permafrost environments [25]. First of all, many surface hydrological models coupled with numerical soil freezing models only consider a vertical direction of heat transfer but ignore the lateral direction, which is crucial to the isolated and discontinuous permafrost bodies located in the lowland regions [157,158], as well as in steep and alpine regions [159,160]. Therefore, one-dimensional heat transfer is only suitable for simulating the seasonal freeze–thaw process of permafrost, but not for multidecadal and multidimensional changes. Secondly, many surface hydrological models do not consider the important factors influencing the rate of the freezing–thawing process in soil, such as the heat capacity and soil layering. A lack of those factors could result in an inaccurate calculation of frost, depths of thawing, subsurface water storage, and groundwater routing. Thirdly, many surface hydrological models or numerical models of soil freezing–thawing processes do not represent the thermo-dynamic equilibrium well, particularly the disequilibrium phase change processes occurring during freezing and thawing, as well as disequilibrium pressure during the infiltration of snowmelt into partially frozen soils. Fourthly, surface hydrological models, even complex physically-based models, cannot accurately simulate the ground surface temperature, which controls the seasonal freezing–thawing process of soil under snowpack and snow-free conditions. This is because ground thermal regimes can be highly sensitive to model parameters [140,161,162]. Lastly, in permafrost environments, surface hydrological models are still unable to clearly represent the relationship between the hydraulic conductivity (K) of soil and the ice, liquid, and gas existing in partially frozen and unsaturated soils. For example, it is unclear whether the K factor of partially frozen soils is decreased or not by resistance factors.

2.13.2. Subsurface Hydrological Models/Groundwater Models/Cryo-Hydrogeological Models

Several one-dimensional groundwater models, coupled with energy transport models, have been developed to simulate the freezing and thawing process in soil since the 1970s. However, the simulation is only limited to the vertical direction. In the last decade, many multi-dimensional models have been developed to support groundwater simulation in permafrost environments [109,130,141,163,164]. Such models work by coupling a three-dimensional equation for water flow (the Richards equation) into a three-dimensional equation for heat transfer (e.g., heat conduction, heat advection, and thermal dispersion). Such models are also known as cryo-hydrogeological models, and they take into account the pore water phase change [165] and the impacts of the latent heat of it on the effectiveness of subsurface heat capacity, as well the decreasing hydraulic conductivity of soil because of the formed pore ice. Moreover, cryo-hydrogeological models have been approached in numerous studies of climate change impacts (e.g., the increase of the ALT) on groundwater (e.g., increase of baseflow and the groundwater exchange between supra-permafrost and sub-permafrost aquifers) in permafrost basins [134,135,136,137,138,157,166,167,168,169,170].

Beside the existing achievements, the cryo-hydrogeological models have several limitations. First of all, most of the models have been applied in ideal environmental conditions, except for a few studies conducting investigations under field conditions [108,157,171]. This is due to the lack of hydrogeological data at high altitude or latitude regions. Secondly, most of the cryo-hydrogeological models do not include the land surface scheme and thermal processes in their structure. Therefore, the boundary conditions for such models must be subsurface conditions, such as the groundwater table, groundwater recharge rate, and soil temperature. However, such input data are limited and they are normally simplified and assumed. Thirdly, the models are limited into two dimensions with simplified structures, only considering small-scale areas. Additionally, the models require a fine grid structure and a small time step for application. Therefore, intensive computational processes are normally expected for the high spatial and temporal resolutions required. Finally, most of cryo-hydrogeological models are difficult for non-expert users to use because of their complexities.

2.14. Model Comparison Regarding the Capacities of the Models to Deal with Permafrost Hydrology in the Arctic

This section aims to summarize and identify the capabilities of the twelve well-known hydrological models for simulating Arctic permafrost hydrology. Permafrost is an important feature affecting hydrological processes in the Arctic. The main factor dominating water storage and transmission in permafrost basins is the ground thermal regime, which has been identified as the typical challenge of hydrological models in permafrost environments [31]. In addition to the ground thermal regime, other important processes related to permafrost hydrology and the Arctic conditions should be considered, such as the surface energy balance, snow process, snow insulation (influencing the air–soil temperature relationships), snowmelt, infiltration, soil heterogeneities, soil moisture regime, the three-phase change of water (ice, liquid, and gas) during the freezing–thawing process of near-surface soil, as well as the dynamics of the ALT. Therefore, the most suitable models for application in permafrost conditions should include as many of the important processes as possible. Table 1 summarizes the important processes for permafrost hydrology that are considered in each model’s structure.

According to Table 1, the important processes in permafrost environments considered in the model include: (1) Surface energy balance; (2) snow process; (3) snow insulation; (4) snowmelt; (5) infiltration; (6) soil thermal; (7) soil moisture; (8) soil heterogeneities; (9) active layer thickness (ALT) dynamics; (10) subsidence; and (11) three-phase change of water during the freezing–thawing process in near-surface soils. Additionally, the symbol of √ denotes available and the abbreviation of n/a defines not available, unclear, or no information.

The twelve hydrological models may be classified into two major groups. The first group are surface hydrological models, including Topoflow, DMHS, HBV, SWAT, WaSiM, ECOMAG, and CRHM. The second group are cryo-hydrogeological models, including ATS, CryoGrid 3, GEOtop, SUTRA-ICE, and PFLOTRAN-ICE. Generally, each model group has its own competences and limitations in dealing with permafrost hydrology, as addressed in prior sections. The first model group performs well for the simulation of surface hydrology, but the capacities of these models to simulate sophisticated ground thermal processes, groundwater flow, or other processes related to permafrost degradation are low. For instance, most of models in group 1 are not able to take into account the influences of snow insulation on air–soil temperature relationships, as well as a required three-phase change of water in permafrost soil (except for DMHS and WaSiM). Some models exclude the important soil thermal process in their structure, such as Topoflow, HBV, and SWAT. None of the surface hydrological models are able to simulate the subsidence as a result of permafrost thawing. On the contrary, many cryo-hydrogeological models are able to simulate subsurface hydrological processes well.

Most of the cryo-hydrogeological models consider a required three-phase change in permafrost environments. Exceptionally, the SUTRA-ICE and PFLOTRAN-ICE models do not include land surface schemes and only simulate the subsurface hydrology. These two models are also not able to simulate the subsidence. Interestingly, almost all of the twelve models in the two groups have capacities to simulate the dynamics of the ALT. However, the accuracy or reliability of the simulation results are not the same for each model. Topoflow, HBV, SWAT, WaSiM, and SUTRA-ICE have low model performance, since they approach a relatively simple method for the simulation of the active layer dynamics, and they lack consideration of spatiotemporal variation.

2.15. Model(s) Selection for the Arctic

Concerning the specific hydrological processes and the sparse data of the Arctic region, the selection of a suitable model should fulfil two main criteria: (1) The capability of the model to deal with unique permafrost hydrological processes; and (2) the possibility to be applied with moderate amounts of data. Table 2, Table 3, Table 4 and Table 5 summarize the main criteria as a basis to choose suitable models for application in the Arctic region. Since the freezing and thawing process of the ALT strongly impacts hydrological processes in the Arctic, the ALT was prioritized among other factors listed in Table 1 for selection of the suitable models for the Arctic.

First of all, if a study aims to only evaluate the surface hydrology process, Topoflow, HBV, SWAT, ECOMAG, and CRHM could be solutions, since they have been verified in many different case studies. Particularly, the conceptual model HBV and the physically-based, semi-distributed SWAT and CRHM models fit best for catchments with moderate data requirements, while the physically-based, distributed models Topoflow and ECOMAG require intensive inputs for every cell of the catchment.

Secondly, if the accuracy of subsurface hydrology simulation is prioritized, DMHS, WaSiM, and the cryo-hydrogeology models are good options, since such models take into account the three phase changes of water, which are important in permafrost environments for analyzing water during the freezing and thawing process of near-surface soils. Other models also have capacities to simulate the subsurface hydrology, but their performances are lower. Further, some cryo-hydrogeology models, e.g., GEOtop, SUTRA-ICE, and PFLOTRAN-ICE, are only suitable for applications in small-scale study areas, while ATS and CryoGrid 3 can be applied to larger basins. Since SUTRA-ICE and PFLOTRAN-ICE do not include a land surface scheme in their structure, these models should be coupled with other land surface models for applications in the Arctic region.

Further, among the five cryo-hydrogeology models, the ATS model has demonstrated its capacity to comprehensively simulate snow processes. It takes into account the snow thermal conductivity, increasing snow density by snow age, and snow insulation. While other models, e.g., CryoGrid 3, have relatively simple snow simulation processes, the ATS and GEOtop models are the first tools that have coupled surface/subsurface permafrost thermal hydrology. Although, such models still have very few applications in the Arctic, so they should be widely tested.

Although the cryo-hydrogeology models appear to be potential tools for researching permafrost hydrology, they could be difficult for non-expert users because of their complicated model structures. Workshops or training should be helpful for new users. Moreover, it is normally difficult to collect enough input data for subsurface hydrological models in Arctic regions. Thus, coupling surface hydrological models as upper boundary conditions for subsurface hydrological models is recommended as one of the solutions, where the outputs from surface hydrological models would become inputs for the subsurface hydrological models. For example, outputs from the SWAT model, e.g., groundwater recharge, soil deep percolation, and stream stage, etc., could be used as inputs for groundwater models.

Thirdly, if the accuracy of ALT dynamics is targeted, DMHS, ATS, GEOtop, and PFLOTRAN-ICE are desirable choices. The other remaining models are not recommended for the study of ALT dynamics since their approaches are relatively simple, resulting in highly uncertainty for the simulation results.

Fourthly, choosing a proper model also depends on the available input data of the catchments and the study conditions (e.g., funding, etc.). Although the cryo-hydrogeology models have high capacities to simulate permafrost hydrology, they still present challenges in terms of gathering the required input data, especially for information about the subsurface boundary conditions in the Arctic. Therefore, it is very expensive to run such simulations.

Lastly, successful application of a model also depends on the knowledge and experience of the hydrologist users to the given model. For example, although HBV is a kind of conceptual, semi-distributed model, it has been widely used in Nordic regions, i.e., cold climate region. In addition, the ECOMAG and CRHM models are well-known in Russia and Canada, respectively (Figure 2). However, the application of CRHM requires accumulated knowledge and data of the study catchments, since there is no calibration procedure or testing cases for the model. This could present a challenge for non-expert or less experienced users.

3. Conclusions

The hydrological processes in the Arctic are more complicated than in other regions. Finding a suitable hydrological model to simulate and predict changes in the hydrological processes is therefore a challenging task for modelers. This paper has reviewed the functions, structures, operational mechanisms, and performances of twelve hydrological models that have previously been applied or have potential for application in the Arctic. Of them, the DMHS, ATS, GEOtop, and PFLOTRAN-ICE models are desirable options for the simulation of ALT dynamics, which has strong impacts on surface/subsurface hydrology. The three cryo-hydrogeology models GEOtop, SUTRA-ICE, and PFLOTRAN-ICE are only suitable for small-scale case studies, while the ATS and CryoGrid 3 models could be applied to large-scale catchments. The SUTRA-ICE and PFLOTRAN-ICE models appear to be suitable subsurface hydrological models. However, since they lack land surface schemes in their structures, they should be applied in tandem with other land surface models for the comprehensive simulation of permafrost hydrology. A significant revolution of hydrological modeling is the recently developed ATS and GEOtop models, which are the first tools that couple surface/subsurface permafrost thermal hydrology. It would be valuable to test the new models and to verify their performances in the Arctic condition. In the situation of only studying the surface hydrology, Topoflow, HBV, SWAT, ECOMAG, and CRHM could be good choices. Particularly, the semi-distributed models HBV, SWAT, and CRHM fit best for catchments with moderate data requirements, whereas the fully distributed models Topoflow and ECOMAG require intensive input data and model parameters. It is worth mentioning that since there is no required calibration procedure for the CRHM and Topoflow models, application of these models requires expert knowledge of the catchments, otherwise the modeling results could feature high uncertainty. Further, if the accuracy of simulating the subsurface hydrology in permafrost environments is targeted, DMHS, WaSiM, and the cryo-hydrogeological models are good options, since they consider the three phase changes of water during the freezing and thawing of near-surface soils. Finally, the selection of the suitable models for Arctic permafrost regions should be based on several factors, besides simulating the permafrost hydrology. Other factors, such as data acquisition, the required research period, and funding could influence the model selection and reliability of the results.

Author Contributions

M.T.B. initiated the idea, planned the scope, collected the information and data, structured and produced a draft manuscript, and performed revisions. J.L. supervised the progress, supported the idea, and conducted the review and editing. L.N. contributed to the reformulation of the paper structure, HBV hydrological model, language revision, and co-supervision of the review and revision processes. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

The authors would like to thank the Department of Technology and Safety, University of Tromsø—The Arctic University of Norway, for their financial support and administration for this study under project number 301184. Additionally, the publication charges for this article have been funded by a grant from the publication fund of UiT The Arctic University of Norway under project 551011. The collaboration with CSDI WaterTech is also highly appreciated. We especially thank the three anonymous reviewers for their constructive comments that helped to reasonably improve our manuscript. Finally, we also thank MDPI for English language editing, especially the English editor Max Kooij.

Conflicts of Interest

The authors declare no conflict of interest.

Figures and Tables

Figure 1. A sketch of impacts and responses of permafrost thawing on water fluxes and distribution. The symbols of ∆, ↗, ↘, ↕, and ? denote change, increase, decrease, variation, and unknown changes, respectively. Modified from [25].

Figure 2. Number of recent publications on applications of hydrological/cryo-hydrogeological models in the Arctic (found via https://www.scopus.com/).

Table 1

Model comparison regarding the model’s capacities to simulate the important processes in permafrost environments.

Model	Important Processes in Permafrost Environments Considered in the Model
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)
Topoflow	√	n/a	n/a	√	√	n/a	√	√	√	n/a	n/a
DMHS	√	√	n/a	√	√	√	√	√	√	n/a	√
HBV	√	√	√	√	√	n/a	√	n/a	√	n/a	n/a
SWAT	√	√	n/a	√	√	n/a	√	√	√	n/a	n/a
WaSiM	√	√	n/a	√	√	√	√	√	√	n/a	√
ECOMAG	√	√	n/a	√	√	√	√	√	√	n/a	n/a
CRHM	√	√	n/a	√	√	√	√	√	√	n/a	n/a
ATS	√	√	√	√	√	√	√	√	√	√	√
CryoGrid 3	√	√	√	√	√	√	√	√	√	√	√
GEOtop	√	√	√	√	√	√	√	√	√	√	√
SUTRA-ICE	n/a	n/a	n/a	n/a	√	√	√	√	√	n/a	√
PFLOTRAN-ICE	n/a	n/a	n/a	n/a	√	√	√	√	√	n/a	√

Table 2

Summarization of model comparisons.

Model	Data Requirements	Time Step	Simulating the ALT Dynamics	Study Area (km2)	Ease-of-Use	Model Availability
Topoflow	Spatial data: Digital elevation model (DEM) and soil.Meteorological data: Precipitation, air temperature, air pressure, wind speed, wind direction, relative humidity, solar radiation, and soil temperature.	Seconds to minutes	Using a relatively simple method. Spatial variability of the ALT is not presented.	<250	Simply learned command syntax; no calibration procedure; requires expert knowledge of the given catchment.	Open code
WaSiM	Spatial data: DEM, land use, and soil.Meteorological data: Precipitation, air temperature, wind speed, vapor pressure, and solar radiation.	Minutes to days	Using a relatively simple method only based on empirical parameters.	<1 to >100,000	The model allows various model configurations depending on the targets of studies, available input data, and quality of input data.The model can be operated with various spatial and temporal discretization solutions.	Open software
ECOMAG	Spatial data: DEM, land use, and soil;Meteorological data (for hydrological submodel): Precipitation, air temperature, and humidity.	Daily	Solving the thermodynamic equations and heat vertical transfer.	<250 to >2500	The model structure can be flexibly adjusted according to the available input data.	Requires a licensed ArcView platform

Table 3

Summarization of model comparisons (continued).

Model	Data Requirements	Time Step	Simulating the ALT Dynamics	Study Area (km2)	Ease-of-Use	Model Availability
DMHS	Spatial data: DEM, land use, and soil.Meteorological data: Precipitation, air temperature, and relative humidity.	Daily/sub-daily	Using a heat transfer analytical solution via considering phase change in the soil profile.	<250 to >2500	Less effort for weather data collection, but difficult for acquisition of the soil profile properties in a suitable format for the model.	n/a
HBV	Spatial data: DEM, land use, and soil;Meteorological data: Precipitation, air temperature, and estimates of potential evapotranspiration.	Daily	Using an accumulated degree day coefficient from field measurements.	<250 to >2500	Requires little time to learn and run the model.	Open software
SWAT	Spatial data: DEM, land use, and soil.Meteorological data: Precipitation, max. and min. air temperature, wind speed, relative humidity, and solar radiation.	Daily	Only using the average values of the ALT.	<250 to >2500	Time consuming for data collection and processing, calibration, and validation.	Requires a licensed ArcGIS platform (applying for ArcSWAT)
CRHM	Spatial data: DEM, land use, and soil.Meteorological data: Precipitation, air temperature, wind speed, relative humidity, and short- and long-wave radiation.	Daily	Solving Stefan’s heat flow equation.	<250 to 2500	No calibration procedure, but requires expert knowledge of the catchment.	Open software

Table 4

Summarization of model comparisons (continued).

Model	Data Requirements	Time Step	Simulating the ALT Dynamics	Study Area (km2)	Ease-of-Use	Model Availability
ATS	Meteorological data: Air temperature, snow precipitation, rain precipitation, wind speed, relative humidity, incoming short-wave and long-wave radiation, and water table elevations.	Daily	Solving equations for the coupled surface (using the diffusion wave equation) and subsurface (using a three phase (ice, liquid, and gas) Richards-like equation), including the energy and flow (using an advection-diffusion equation) and surface energy balance, including snow.	<250 to >2500	Preparing the XML input code is a challenge.	Open code
CryoGrid 3	Meteorological data: Air temperature, relative or absolute humidity, wind speed, incoming short-wave and long-wave radiation, air pressure, and rates of snowfall and rainfall.	Daily	Using a simple 1D parameterization.	<250 to >2500	The source code is simple and modifiable.	Open code
GEOtop	Spatial data: Elevation (DTM), land use, and soil.Meteorological data: Precipitation intensity, wind velocity, wind direction, relative humidity, air temperature, dew temperature, air pressure, short-wave solar global radiation, short-wave solar direct radiation, short-wave solar diffuse radiation, short-wave solar net radiation, and long-wave incoming radiation.	Hourly	Solving the energy and mass balance equations dealing with phase change based on the globally convergent Newton scheme.	<250	The source code is modifiable. High effort is required for data acquisition and processing of the hourly forcing data.	Open code

Table 5

Summarization of model comparisons (continued).

Model	Data Requirements	Time Step	Simulating the ALT Dynamics	Study Area (km2)	Ease-of-Use	Model Availability
SUTRA-ICE	Flow data (specified pressures, specified flows and fluid sources) and energy or solute data (specified temperatures or concentrations, diffusive fluxes of energy, or solute mass at boundaries).	<1 to <12 h	Approaching a two-zone (frozen and thawed) analytical solution to simulate ice forming and melting in porous media, also ignoring a mushy zone containing both ice and water (as considered in the three-zone analytical solution of Lunardini).	Applicable for small-scale study areas (approx. a few hundred square meters).	The source code can be easily modified to add new processes by users.	Open code
PFLOTRAN-ICE	River stage, river chemistry, groundwater recharge, specified infiltration rate, temperature, gas pressure, and infiltration chemistry.	Hours to days	Solving the energy and mass balance equation for soil water in three-phase (ice, liquid, and gas) change.	The size of the study area can be up to a few kilometers.	The source code can be easily modified and further developed by users.	Open code