diff --git a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C
index 8380d06ef7..3dea2ae910 100644
--- a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C
+++ b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C
@@ -77,6 +77,7 @@ Foam::surfaceFilmModels::standardPhaseChange::standardPhaseChange
     Tb_(readScalar(coeffs_.lookup("Tb"))),
     deltaMin_(readScalar(coeffs_.lookup("deltaMin"))),
     L_(readScalar(coeffs_.lookup("L"))),
+    TbFactor_(coeffs_.lookupOrDefault<scalar>("TbFactor", 1.1)),
     totalMass_(0.0),
     vapourRate_(0.0)
 {}
@@ -133,7 +134,7 @@ void Foam::surfaceFilmModels::standardPhaseChange::correct
             const scalar pc = pInf[cellI];
 
             // local temperature - impose lower limit of 200 K for stability
-            const scalar Tloc = min(liq.Tc(), max(200.0, T[cellI]));
+            const scalar Tloc = min(TbFactor_*Tb_, max(200.0, T[cellI]));
 
             // saturation pressure [Pa]
             const scalar pSat = liq.pv(pc, Tloc);
@@ -149,7 +150,9 @@ void Foam::surfaceFilmModels::standardPhaseChange::correct
                 const scalar qDotFilm = hFilm[cellI]*(T[cellI] - Tw[cellI]);
 
                 const scalar Cp = liq.Cp(pc, Tloc);
-                const scalar qCorr = availableMass[cellI]*Cp*(T[cellI] - Tb_);
+                const scalar Tcorr = max(0.0, T[cellI] - Tb_);
+                const scalar qCorr = availableMass[cellI]*Cp*(Tcorr);
+
                 dMass[cellI] =
                     dt*magSf[cellI]/hVap*(qDotInf + qDotFilm)
                   + qCorr/hVap;
@@ -185,7 +188,7 @@ void Foam::surfaceFilmModels::standardPhaseChange::correct
                 const scalar Sh = this->Sh(Re, Sc);
 
                 // mass transfer coefficient [m/s]
-                const scalar hm = Sh*Dab/L_;
+                const scalar hm = Sh*Dab/(L_ + ROOTVSMALL);
 
                 // add mass contribution to source
                 dMass[cellI] =
diff --git a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.H b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.H
index bb9e69e4c6..17138974c4 100644
--- a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.H
+++ b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.H
@@ -76,6 +76,10 @@ protected:
         //- Length scale / [m]
         const scalar L_;
 
+        //- Boiling temperature factor / []
+        //  Used to set max limit on temperature to Tb*TbFactor
+        const scalar TbFactor_;
+
         //- Total mass evolved / [kg]
         scalar totalMass_;
 
diff --git a/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayer.C b/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayer.C
index ec1391ec18..1977e339ad 100644
--- a/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayer.C
+++ b/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayer.C
@@ -209,12 +209,14 @@ void Foam::surfaceFilmModels::thermoSingleLayer::solveEnergy()
         massForPrimary_/magSf_/time_.deltaT()
     );
 
+//    dimensionedScalar hs0("SMALL", hs_.dimensions(), SMALL);
+
     solve
     (
         fvm::ddt(deltaRho_, hs_)
       + fvm::div(phi_, hs_)
      ==
-//        fvm::Sp(hsSp_/hs_, hs_)
+//        fvm::Sp(hsSp_/(hs_ + hs0), hs_)
         hsSp_
       + q(hs_)
       - fvm::Sp(mLossCoeff, hs_)