  In this paper we consider a nonlinear partial differential equation
describing heat flow with ice-water phase transition in permafrost soils. Such
models and their numerical approximations have been well explored in the
applications literature. In this paper we describe a new direction in which the
allow relaxation and hysteresis of the phase transition which introduce
additional nonlinear terms and complications for the analysis. We present
numerical algorithms as well as analysis of the well-posedness and convergence
of the fully implicit iterative schemes. The analysis we propose handles the
equilibrium, non-equilibrium, and hysteresis cases in a unified way.
  We also illustrate with numerical examples for a model ODE and PDE.
