We study heat conduction in permafrost soil at the Darcy scale. These involve nonlinear parabolic PDEs accounting for the phase transitions between frozen and unfrozen water phase transitions. We extend previous work of [Bigler, Vohra, Slugg, Mathangadeera and Peszynska] on the Stefan (ST) problem and adapted permafrost (P) models to treat the most general non-adapted case called (P*). Mathematically and computationally, (P*) is most challenging, while it is also most flexible from the applications point of view as it allows, in particular, immobile pockets of air and macro-pores. The model is only posed in the sense of distributions, and we develop a new analytical solution to the (P*)-model.

The main challenge of the (P*)-model is that it involves the enthalpy-temperature relationship, $w \in \alpha(\theta)$, which is a multivalued graph at the freezing point. We show that it can be reformulated using a single-valued bijective map from $(\theta,\chi)$ to $w$, where $\chi$ is the unfrozen water content. After the PDE is discretized, we require both a local and global nonlinear solver to solve for the primary unknowns. We propose a new accelerated method with Newton-Anderson iteration as the global nonlinear solver. We also suggest alternative options for local nonlinear solvers and study other methods to accelerate the convergence including the use of look-up tables and/or extrapolating initial guesses to aid in the convergence of the scheme. Finally, we test the solvers with uncertain data through non-intrusive Monte-Carlo simulations in order to provide better insight when modeling permafrost applications.