diff --git a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/noPhaseChange/noPhaseChange.C b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/noPhaseChange/noPhaseChange.C
index 67094e226d..8e023485b5 100644
--- a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/noPhaseChange/noPhaseChange.C
+++ b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/noPhaseChange/noPhaseChange.C
@@ -61,6 +61,7 @@ Foam::surfaceFilmModels::noPhaseChange::~noPhaseChange()
 void Foam::surfaceFilmModels::noPhaseChange::correct
 (
     const scalar,
+    scalarField&,
     scalarField&
 )
 {
diff --git a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/noPhaseChange/noPhaseChange.H b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/noPhaseChange/noPhaseChange.H
index 22a2a17e47..e75a9f78c4 100644
--- a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/noPhaseChange/noPhaseChange.H
+++ b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/noPhaseChange/noPhaseChange.H
@@ -84,7 +84,12 @@ public:
         // Evolution
 
             //- Correct
-            virtual void correct(const scalar dt, scalarField& dMass);
+            virtual void correct
+            (
+                const scalar dt,
+                scalarField& dMass,
+                scalarField& dEnergy
+            );
 };
 
 
diff --git a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/phaseChangeModel/phaseChangeModel.H b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/phaseChangeModel/phaseChangeModel.H
index 8421520d98..63f0c227c1 100644
--- a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/phaseChangeModel/phaseChangeModel.H
+++ b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/phaseChangeModel/phaseChangeModel.H
@@ -137,7 +137,12 @@ public:
         // Evolution
 
             //- Correct
-            virtual void correct(const scalar dt, scalarField& dMass) = 0;
+            virtual void correct
+            (
+                const scalar dt,
+                scalarField& dMass,
+                scalarField& dEnergy
+            ) = 0;
 };
 
 
diff --git a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C
index 9011d02fca..4619a4a76a 100644
--- a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C
+++ b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C
@@ -90,7 +90,8 @@ Foam::surfaceFilmModels::standardPhaseChange::~standardPhaseChange()
 void Foam::surfaceFilmModels::standardPhaseChange::correct
 (
     const scalar dt,
-    scalarField& dMass
+    scalarField& dMass,
+    scalarField& dEnergy
 )
 {
     const thermoSingleLayer& film = refCast<const thermoSingleLayer>(owner_);
@@ -111,11 +112,11 @@ void Foam::surfaceFilmModels::standardPhaseChange::correct
     const scalarField& TInf = film.TPrimary();
     const scalarField& rho = film.rho();
     const scalarField& mu = film.mu();
-    const scalarField& hs = film.hs();
     const scalarField& magSf = film.magSf();
     const scalarField hInf = film.htcs().h();
     const scalarField hFilm = film.htcw().h();
     const vectorField dU = film.UPrimary() - film.Us();
+    const scalarField availableMass = delta*rho*magSf;
 
     // Reynolds number
     const scalarField Re = rho*mag(dU)*L_/mu;
@@ -136,18 +137,16 @@ void Foam::surfaceFilmModels::standardPhaseChange::correct
             // saturation pressure
             const scalar pSat = liq.pv(pc, Ts[cellI]);
 
+            // latent heat
+            const scalar hVap = liq.hl(pc, Ts[cellI]);
+
             // calculate mass transfer
             if (pSat > pc)
             {
                 // boiling
                 const scalar qDotInf = hInf[cellI]*(TInf[cellI] - T[cellI]);
                 const scalar qDotFilm = hFilm[cellI]*(T[cellI] - Tw[cellI]);
-                dMass[cellI] +=
-                    max
-                    (
-                        0.0,
-                        dt*magSf[cellI]/hs[cellI]*(qDotInf - qDotFilm)
-                    );
+                dMass[cellI] = dt*magSf[cellI]/hVap*(qDotInf - qDotFilm);
             }
             else
             {
@@ -171,12 +170,11 @@ void Foam::surfaceFilmModels::standardPhaseChange::correct
 
                 // add mass contribution to source
                 dMass[cellI] =
-                    max
-                    (
-                        0.0,
-                        dt*magSf[cellI]*rhoAve*hm*(Ys - YInf[cellI])/(1.0 - Ys)
-                    );
+                    dt*magSf[cellI]*rhoAve*hm*(Ys - YInf[cellI])/(1.0 - Ys);
             }
+
+            dMass[cellI] = min(availableMass[cellI], max(0.0, dMass[cellI]));
+            dEnergy[cellI] = dMass[cellI]*hVap;
         }
     }
 }
diff --git a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.H b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.H
index dd16488e03..726f36e2cc 100644
--- a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.H
+++ b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.H
@@ -105,7 +105,12 @@ public:
         // Evolution
 
             //- Correct
-            virtual void correct(const scalar dt, scalarField& dMass);
+            virtual void correct
+            (
+                const scalar dt,
+                scalarField& dMass,
+                scalarField& dEnergy
+            );
 };
 
 
diff --git a/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayer.C b/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayer.C
index 2d0b4a2fff..480759783c 100644
--- a/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayer.C
+++ b/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayer.C
@@ -131,7 +131,7 @@ void Foam::surfaceFilmModels::thermoSingleLayer::updateSurfaceTemperatures()
     Ts_ =
         (
             // qRad
-          - massPhaseChangeForPrimary_*hs_/(time_.deltaT()*magSf_)
+          - energyPhaseChangeForPrimary_/(time_.deltaT()*magSf_)
           + TPrimary_*htcs_->h()
           + kappaDeltaBy2*T_
         )
@@ -176,7 +176,14 @@ void Foam::surfaceFilmModels::thermoSingleLayer::updateSubmodels()
 
     // Update phase change
     massPhaseChangeForPrimary_ == dimensionedScalar("zero", dimMass, 0.0);
-    phaseChange_->correct(time_.deltaTValue(), massPhaseChangeForPrimary_);
+    energyPhaseChangeForPrimary_ == dimensionedScalar("zero", dimEnergy, 0.0);
+
+    phaseChange_->correct
+    (
+        time_.deltaTValue(),
+        massPhaseChangeForPrimary_,
+        energyPhaseChangeForPrimary_
+    );
     massPhaseChangeForPrimary_.correctBoundaryConditions();
     totalMassPhaseChange_ += sum(massPhaseChangeForPrimary_).value();
 
@@ -187,7 +194,8 @@ void Foam::surfaceFilmModels::thermoSingleLayer::updateSubmodels()
     const dimensionedScalar deltaT = time_.deltaT();
     rhoSp_ -= massPhaseChangeForPrimary_/magSf_/deltaT;
     USp_ -= massPhaseChangeForPrimary_*U_/magSf_/deltaT;
-    hsSp_ -= (massForPrimary_ + massPhaseChangeForPrimary_)*hs_/magSf_/deltaT;
+    hsSp_ -=
+        (massForPrimary_*hs_ + energyPhaseChangeForPrimary_)/magSf_/deltaT;
 }
 
 
@@ -380,7 +388,21 @@ Foam::surfaceFilmModels::thermoSingleLayer::thermoSingleLayer
         heatTransferModel::New(*this, coeffs_.subDict("lowerSurfaceModels"))
     ),
     phaseChange_(phaseChangeModel::New(*this, coeffs_)),
-    totalMassPhaseChange_(0.0)
+    totalMassPhaseChange_(0.0),
+    energyPhaseChangeForPrimary_
+    (
+        IOobject
+        (
+            "energyPhaseChangeForPrimary",
+            time_.timeName(),
+            filmRegion_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        filmRegion_,
+        dimensionedScalar("zero", dimEnergy, 0),
+        zeroGradientFvPatchScalarField::typeName
+    )
 {
     if (thermo_.hasMultiComponentCarrier())
     {
diff --git a/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayer.H b/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayer.H
index a7255a9994..a5c5b892e7 100644
--- a/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayer.H
+++ b/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayer.H
@@ -150,6 +150,9 @@ protected:
             //- Total mass transferred to primary region [kg]
             scalar totalMassPhaseChange_;
 
+            //- Film energy evolved via phase change
+            volScalarField energyPhaseChangeForPrimary_;
+
 
     // Protected member functions