Ground heat flux (G0) is a key component of the land-surface energy balance of high-latitude regions. Despite its crucial role in controlling permafrost degradation due to global warming, G0 is sparsely measured and not well represented in the outputs of global scale model simulation. In this study, an analytical heat transfer model is tested to reconstruct G0 across seasons using soil temperature series from field measurements, Global Climate Model, and climate reanalysis outputs. The probability density functions of ground heat flux and of model parameters are inferred using available G0 data (measured or modeled) for snow-free period as a reference. When observed G0 is not available, a numerical model is applied using estimates of surface heat flux (dependent on parameters) as the top boundary condition. These estimates (and thus the corresponding parameters) are verified by comparing the distributions of simulated and measured soil temperature at several depths. Aided by state-of-the-art uncertainty quantification methods, the developed G0 reconstruction approach provides novel means for assessing the probabilistic structure of the ground heat flux for regional permafrost change studies.