diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C
index ba153f166a87addf99eefeaaac431c37ea7784f3..ab62b151362510d4ea34c248948ef4f41119be45 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C
@@ -288,6 +288,38 @@ Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::heatTransfer() const
         }
     }
 
+    // Source term due to mass trasfer
+    forAllConstIter
+    (
+        phaseSystem::phasePairTable,
+        this->phasePairs_,
+        phasePairIter
+    )
+    {
+        const phasePair& pair(phasePairIter());
+
+        if (pair.ordered())
+        {
+            continue;
+        }
+
+        const volScalarField& he1(pair.phase1().thermo().he());
+        const volScalarField& he2(pair.phase2().thermo().he());
+
+        const volScalarField& K1(pair.phase1().K());
+        const volScalarField& K2(pair.phase2().K());
+
+        const volScalarField dmdt(this->dmdt(pair));
+        const volScalarField dmdt12(dmdt*pos(dmdt));
+        const volScalarField dmdt21(dmdt*neg(dmdt));
+
+        *eqns[pair.phase1().name()] +=
+            fvm::Sp(dmdt21, he1) + dmdt21*K1 - dmdt21*(he2 + K2);
+
+        *eqns[pair.phase2().name()] +=
+            dmdt12*(he1 + K1) - fvm::Sp(dmdt12, he2) - dmdt12*K2;
+    }
+
     return eqnsPtr;
 }
 
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
index cbefea2852e42763878b35f7f2e31279a2b91da8..3a2786d75c1851d2d1d7d5c0de0fc67addcaadba 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
@@ -154,4 +154,12 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::divU(const volScalarField& divU)
 }
 
 
+template<class BasePhaseModel>
+const Foam::volScalarField&
+Foam::AnisothermalPhaseModel<BasePhaseModel>::K() const
+{
+    return K_;
+}
+
+
 // ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H
index 48b0cbd499118bcd674eeab7556c97e69832fbe2..e7a76bf496e24375cb3956492abc254ee800ccd4 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H
@@ -89,11 +89,14 @@ public:
             //- Return true if the phase is compressible otherwise false
             virtual bool compressible() const;
 
-            //- Phase dilatation rate (d(alpha)/dt + div(alpha*phi))
+            //- Return the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
             virtual const volScalarField& divU() const;
 
-            //- Set phase dilatation rate (d(alpha)/dt + div(alpha*phi))
+            //- Set the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
             virtual void divU(const volScalarField& divU);
+
+            //- Return the phase kinetic energy
+            virtual const volScalarField& K() const;
 };
 
 
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
index aeeef98149cddb0e0d06c51a92ec9615796fb0f2..b1869b9fbfe2ac7f0e85b3bbdbd3523e999f7951 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
@@ -314,14 +314,6 @@ Foam::MovingPhaseModel<BasePhaseModel>::continuityError() const
 }
 
 
-template<class BasePhaseModel>
-Foam::volScalarField&
-Foam::MovingPhaseModel<BasePhaseModel>::continuityError()
-{
-    return continuityError_;
-}
-
-
 template<class BasePhaseModel>
 Foam::tmp<Foam::surfaceScalarField>
 Foam::MovingPhaseModel<BasePhaseModel>::phi() const
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H
index 96e8942f6fcc380c231f8268115ab7031d7cc343..14e9c057b7eaaaf2f92db29eea3c0b73abab3bf4 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H
@@ -152,9 +152,6 @@ public:
             //- Constant access the continuity error
             virtual tmp<volScalarField> continuityError() const;
 
-            //- Access the continuity error
-            virtual volScalarField& continuityError();
-
             //- Constant access the volumetric flux
             virtual tmp<surfaceScalarField> phi() const;
 
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C
index 2697587b6a56fa303e4cb41ad4abbf39ed47d94d..61ea75e586c4c1c27e64bb3e88866d87e12aad6f 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C
@@ -155,6 +155,13 @@ void Foam::phaseModel::divU(const volScalarField& divU)
 }
 
 
+const Foam::volScalarField& Foam::phaseModel::K() const
+{
+    notImplemented("Foam::phaseModel::K()");
+    return volScalarField::null();
+}
+
+
 const Foam::surfaceScalarField& Foam::phaseModel::DbyA() const
 {
     return surfaceScalarField::null();
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
index de85779d550c474bc6427f1551d689a711ced063..082bb6c019169e69038b8fe7870ef5519601d609 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
@@ -175,6 +175,9 @@ public:
             //- Set the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
             virtual void divU(const volScalarField& divU);
 
+            //- Return the phase kinetic energy
+            virtual const volScalarField& K() const;
+
 
         // Implicit phase pressure and dispersion support
 
@@ -218,9 +221,6 @@ public:
             //- Constant access the continuity error
             virtual tmp<volScalarField> continuityError() const = 0;
 
-            //- Access the continuity error
-            virtual volScalarField& continuityError() = 0;
-
             //- Constant access the volumetric flux
             virtual tmp<surfaceScalarField> phi() const = 0;