diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
index 21e52385e9..77a86260d5 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
@@ -68,7 +68,8 @@ ThermalPhaseChangePhaseSystem
 )
 :
     BasePhaseSystem(mesh),
-    volatile_(this->template lookupOrDefault<word>("volatile", "none"))
+    volatile_(this->template lookupOrDefault<word>("volatile", "none")),
+    dmdt0s_(this->phases().size())
 {
     this->generatePairsAndSubModels
     (
@@ -325,9 +326,15 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::heatTransfer() const
         const volScalarField dmdtf21(posPart(dmdtf));
         const volScalarField dmdtf12(negPart(dmdtf));
 
-        *eqns[phase1.name()] += - fvm::Sp(dmdtf21, he1) + dmdtf21*(K2 - K1);
+        *eqns[phase1.name()] +=
+            fvm::Sp(dmdt0s_[phase1.index()],he1)
+          - fvm::Sp(dmdtf21, he1)
+          + dmdtf21*K2 + dmdtf12*K1;
 
-        *eqns[phase2.name()] -= - fvm::Sp(dmdtf12, he2) + dmdtf12*(K1 - K2);
+        *eqns[phase2.name()] -=
+          - fvm::Sp(dmdt0s_[phase2.index()],he2)
+          - fvm::Sp(dmdtf12, he2)
+          + dmdtf12*K1 + dmdtf21*K2;
 
         if (this->saturationModels_.found(phasePairIter.key()))
         {
@@ -376,14 +383,12 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::heatTransfer() const
         else
         {
             *eqns[phase1.name()] += dmdtf21*he2;
-
             *eqns[phase2.name()] -= dmdtf12*he1;
         }
 
         if (this->nDmdtLfs_.found(phasePairIter.key()))
         {
             *eqns[phase1.name()] += negPart(*this->nDmdtLfs_[pair]);
-
             *eqns[phase2.name()] -= posPart(*this->nDmdtLfs_[pair]);
         }
 
@@ -398,7 +403,6 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::heatTransfer() const
             *eqns[phase2.name()] -=
                 phase2.thermo().p()*dmdtf/phase2.thermo().rho();
         }
-
     }
 
     return eqnsPtr;
@@ -457,6 +461,34 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::specieTransfer() const
 }
 
 
+template<class BasePhaseSystem>
+void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::
+correctContinuityError()
+{
+    dmdt0s_ = PtrList<volScalarField>(this->phases().size());
+
+    forAllConstIter(phaseSystem::dmdtfTable, dmdtfs_, dmdtfIter)
+    {
+        const phasePair& pair = this->phasePairs_[dmdtfIter.key()];
+        const volScalarField& dmdtf = *dmdtfIter();
+
+        addField(pair.phase1(), "dmdt", dmdtf, dmdt0s_);
+        addField(pair.phase2(), "dmdt", - dmdtf, dmdt0s_);
+    }
+
+    forAllConstIter(phaseSystem::dmdtfTable, nDmdtfs_, nDmdtfIter)
+    {
+        const phasePair& pair = this->phasePairs_[nDmdtfIter.key()];
+        const volScalarField& nDmdtf = *nDmdtfIter();
+
+        addField(pair.phase1(), "dmdt", nDmdtf, dmdt0s_);
+        addField(pair.phase2(), "dmdt", - nDmdtf, dmdt0s_);
+    }
+
+    BasePhaseSystem::correctContinuityError();
+}
+
+
 template<class BasePhaseSystem>
 void
 Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctInterfaceThermo()
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H
index e40a63d85e..756bffe50b 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H
@@ -91,16 +91,9 @@ class ThermalPhaseChangePhaseSystem
         //- Name of the volatile specie
         word volatile_;
 
-        //- The saturation model used to evaluate Tsat = Tf
-        //autoPtr<saturationModel> saturationModel_;
-
         //- Saturation models used to evaluate Tsat = Tf
         saturationModelTable saturationModels_;
 
-        //- PhasePairs between which phaseChange occurs, e.g.,
-        //  "((gasI and liquid) (gasII and liquid))"
-        //List<phasePairKey> phaseChangePairKeys_;
-
         //- Mass transfer rates
         phaseSystem::dmdtfTable dmdtfs_;
 
@@ -110,6 +103,10 @@ class ThermalPhaseChangePhaseSystem
         //- Nucleate thermal energy transfer rates
         phaseSystem::dmdtfTable nDmdtLfs_;
 
+        //- Previous continuity error update phase dmdts for the heat transfer
+        //  function
+        PtrList<volScalarField> dmdt0s_;
+
 
     // Private member functions
 
@@ -153,6 +150,9 @@ public:
         virtual autoPtr<phaseSystem::specieTransferTable>
             specieTransfer() const;
 
+        //- Store phase dmdts at the during the continuity error update
+        virtual void correctContinuityError();
+
         //- Correct the interface thermodynamics
         virtual void correctInterfaceThermo();