The present paper proposes an iterative implicit numerical method for simulating the thaw depth of permafrost soil. For this purpose, the enthalpy-porosity model was used for the phase change process, and the finite difference scheme FTCS (Forward Time Centered Space) was used for discretization. An artificial mushy zone was maintained with the same thickness by keeping the regularization parameter proportional to the temperature gradient. In doing so, we made the scheme more stable and convergence occurred faster. The model accuracy was validated by comparing the numerical results with the analytical Stefan solution and with the results of a derived numerical model, based on an explicit scheme. The model performance was also tested against observation data collected on four different landscapes with different soil profiles and located on a basin underlain by continuous permafrost. It was found that the proposed model matched noticeably well the analytical solution for a volumetric liquid fraction (phi) equal to 0.5 regardless of the grid resolution. Furthermore, compared with the observation data, the model reproduced the annual maximum thaw depth with an absolute error lying between 0.7 and 7.7%. In addition, the designed algorithm allowed the model to converge after a maximum of eight iterations, reducing the computational time by around 75% compared to the explicit model. The results were so encouraging that the model can be included in a hydrological modeling of permafrost watersheds or cold regions in general.

In regions characterized by continuous permafrost, hydrological modeling remains a complex activity, primarily due to constraints related to the prevailing climatic conditions and the specific behavior of the active layer. High-latitude regions receive less solar radiation; thus, most creeks are active only during summertime and stay frozen in the winter. To realistically simulate watersheds underlain by continuous permafrost, the heat transfer through the soil needs to be accounted for in the modeling process. In this study, a watershed located in a continuous permafrost zone in Russia is investigated. A model is proposed to integrate this heat transfer into an existing conceptual rain-flow transformation model, Hydrologiska Byråns Vattenbalansavdelning (HBV), to calculate the seasonal thaw depth and determine the components of water balance. The proposed integration is a novelty compared to the standard model, as it enables the physical and thermal properties of the soil to be taken into account. It was found that the proposed model, HBV-Heat, performs better than the stand-alone HBV model. Specifically, the average Nash–Sutcliffe efficiency (NSE) increases by 30% for the whole calibration period. In terms of the water balance components, the results are consistent with previous studies, showing that surface runoff represents 64% of the observed precipitation.