An evaluation of three numerical models used in simulations of the active layer and permafrost temperature regimes

Abstract

Three numerical models (designated the Goodrich, Guymon/Hromadka, and Seregina models) used for calculations of the ground thermal regime which are based on different numerical methods and employ different treatments of freezing and thawing were compared with each other, with analytical solutions, and with measured temperature data. Comparisons of the models with the Neumann solution show differences generally less than 0.2°C between calculated temperatures using a wide range of time and depth steps. The Goodrich and Guymon/Hromadka models have been shown to predict temperature field dynamics reliably in the active layer and permafrost using small time and depth steps. However, comparisons of the models with each other using large time and depth steps and field data for the surface boundary condition showed significant differences between them (RMS deviations exceeding 1°C) and, in addition, the development of a non-physical feature (thaw bulb after freeze-up). Therefore, with large time and depth steps, the models cannot reproduce the temperature field dynamics in the active layer and permafrost. Consequently, agreement with the Neumann solution is necessary but not sufficient to qualify the models for calculations of real temperature fields. The Goodrich model requires a time step not longer than 1 h and depth step in the upper 1 m not larger than 0.02 m to reproduce the temperature regime with reasonable accuracy. However, the choice of optimum time and depth steps appears to be specific to the application. Using the Guymon/Hromadka model, similar accuracy can be obtained with a 1 h time step and 0.1 m space step within the upper 1 m depth or a 1 day time step and 0.01 m space step. However, the use of larger steps does not necessarily decrease the calculational time compared to the Goodrich model. For the case with unfrozen water present in the frozen soil, the results of calculations using the numerical models were compared with an analytical solution and were found to agree within 0.02°C.