diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C
index 95636c3749..07b41eb771 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C
@@ -57,6 +57,16 @@ Foam::HeatTransferPhaseSystem<BasePhaseSystem>::~HeatTransferPhaseSystem()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+template<class BasePhaseSystem>
+bool Foam::HeatTransferPhaseSystem<BasePhaseSystem>::transfersMass
+(
+    const phaseModel& phase
+) const
+{
+    return false;
+}
+
+
 template<class BasePhaseSystem>
 Foam::tmp<Foam::volScalarField>
 Foam::HeatTransferPhaseSystem<BasePhaseSystem>::dmdt
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.H
index cef3b696bd..0c908239fd 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.H
@@ -89,6 +89,9 @@ public:
 
     // Member Functions
 
+        //- Return true if there is mass transfer for phase
+        virtual bool transfersMass(const phaseModel& phase) const;
+
         //- Return the interfacial mass flow rate
         virtual tmp<volScalarField> dmdt
         (
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/reactionThermo/hRefConstThermos.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/reactionThermo/hRefConstThermos.C
index 76c2351109..f43303a2b7 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/reactionThermo/hRefConstThermos.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/reactionThermo/hRefConstThermos.C
@@ -254,6 +254,16 @@ makeReactionMixtureThermo
 );
 
 
+makeReactionMixtureThermo
+(
+    rhoThermo,
+    rhoReactionThermo,
+    heRhoThermo,
+    multiComponentMixture,
+    constRefGasHThermoPhysics
+);
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index a79e9393b5..5512fd6559 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
@@ -515,12 +515,6 @@ Foam::multiphaseSystem::~multiphaseSystem()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-bool Foam::multiphaseSystem::transfersMass(const phaseModel& phase) const
-{
-    return false;
-}
-
-
 Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseSystem::surfaceTension
 (
     const phaseModel& phase1
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H
index a9390f2589..e33bf8c487 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H
@@ -176,7 +176,7 @@ public:
         ) const = 0;
 
         //- Return true if there is mass transfer for phase
-        virtual bool transfersMass(const phaseModel& phase) const;
+        virtual bool transfersMass(const phaseModel& phase) const = 0;
 
         //- Return the total interfacial mass transfer rate for phase
         virtual tmp<volScalarField> dmdt(const phaseModel& phase) const = 0;
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
index 230b650106..2703d75ef1 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
@@ -357,11 +357,15 @@ while (pimple.correct())
             {
                 if (pEqnComps.set(phasei))
                 {
-                    pEqnComps[phasei] -= fluid.dmdt(phase);
+                    pEqnComps[phasei] -= fluid.dmdt(phase)/rho;
                 }
                 else
                 {
-                    pEqnComps.set(phasei, fvm::Su(-fluid.dmdt(phase), rho));
+                    pEqnComps.set
+                    (
+                        phasei,
+                        fvm::Su(-fluid.dmdt(phase)/rho, p_rgh)
+                    );
                 }
             }
         }
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
index c8a9e28288..1af167cc40 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
@@ -247,7 +247,7 @@ while (pimple.correct())
 
             pEqnComp1 =
                 (
-                    phase1.continuityError() - fluid.dmdt()
+                    phase1.continuityError()
                   - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
                 )/rho1
               + (alpha1/rho1)*correction
@@ -270,7 +270,7 @@ while (pimple.correct())
 
             pEqnComp2 =
                 (
-                    phase2.continuityError() + fluid.dmdt()
+                    phase2.continuityError()
                   - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
                 )/rho2
                + (alpha2/rho2)*correction
@@ -288,7 +288,7 @@ while (pimple.correct())
         {
             pEqnComp1 =
                 (
-                    phase1.continuityError() - fluid.dmdt()
+                    phase1.continuityError()
                   - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
                 )/rho1
               + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
@@ -298,13 +298,34 @@ while (pimple.correct())
         {
             pEqnComp2 =
                 (
-                    phase2.continuityError() + fluid.dmdt()
+                    phase2.continuityError()
                   - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
                 )/rho2
               + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh));
         }
     }
 
+    if (fluid.transfersMass())
+    {
+        if (pEqnComp1.valid())
+        {
+            pEqnComp1() -= fluid.dmdt()/rho1;
+        }
+        else
+        {
+            pEqnComp1 = fvm::Su(-fluid.dmdt()/rho1, p_rgh);
+        }
+
+        if (pEqnComp2.valid())
+        {
+            pEqnComp2() += fluid.dmdt()/rho2;
+        }
+        else
+        {
+            pEqnComp2 = fvm::Su(fluid.dmdt()/rho2, p_rgh);
+        }
+    }
+
     // Cache p prior to solve for density update
     volScalarField p_rgh_0(p_rgh);
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
index c3f7d472fa..e2edf14bed 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
@@ -231,7 +231,7 @@ while (pimple.correct())
         {
             pEqnComp1 =
             (
-                phase1.continuityError() - fluid.dmdt()
+                phase1.continuityError()
               - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
             )/rho1
           + (alpha1/rho1)*correction
@@ -247,7 +247,7 @@ while (pimple.correct())
         {
             pEqnComp2 =
             (
-                phase2.continuityError() + fluid.dmdt()
+                phase2.continuityError()
               - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
             )/rho2
           + (alpha2/rho2)*correction
@@ -265,7 +265,7 @@ while (pimple.correct())
         {
             pEqnComp1 =
             (
-                phase1.continuityError() - fluid.dmdt()
+                phase1.continuityError()
               - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
             )/rho1
           + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
@@ -275,13 +275,34 @@ while (pimple.correct())
         {
             pEqnComp2 =
             (
-                phase2.continuityError() + fluid.dmdt()
+                phase2.continuityError()
               - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
             )/rho2
           + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh));
         }
     }
 
+    if (fluid.transfersMass())
+    {
+        if (pEqnComp1.valid())
+        {
+            pEqnComp1() -= fluid.dmdt()/rho1;
+        }
+        else
+        {
+            pEqnComp1 = fvm::Su(-fluid.dmdt()/rho1, p_rgh);
+        }
+
+        if (pEqnComp2.valid())
+        {
+            pEqnComp2() += fluid.dmdt()/rho2;
+        }
+        else
+        {
+            pEqnComp2 = fvm::Su(fluid.dmdt()/rho2, p_rgh);
+        }
+    }
+
     // Cache p prior to solve for density update
     volScalarField p_rgh_0("p_rgh_0", p_rgh);
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
index 187ca4c238..7b7ddee6f4 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
@@ -168,6 +168,12 @@ Foam::twoPhaseSystem::D() const
 }
 
 
+bool Foam::twoPhaseSystem::transfersMass() const
+{
+    return transfersMass(phase1());
+}
+
+
 Foam::tmp<Foam::volScalarField>
 Foam::twoPhaseSystem::dmdt() const
 {
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H
index deca55fa6e..3c4c98f589 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H
@@ -78,6 +78,9 @@ class twoPhaseSystem
         //  Multiplies the phase-fraction gradient
         virtual tmp<volScalarField> D(const phasePairKey& key) const = 0;
 
+        //- Return true if there is mass transfer for phase
+        virtual bool transfersMass(const phaseModel& phase) const = 0;
+
         //- Return the interfacial mass flow rate for phase pair
         virtual tmp<volScalarField> dmdt(const phasePairKey& key) const = 0;
 
@@ -192,6 +195,9 @@ public:
         //  Multiplies the phase-fraction gradient
         tmp<volScalarField> D() const;
 
+        //- Return true if there is mass transfer
+        bool transfersMass() const;
+
         //- Return the interfacial mass flow rate
         tmp<volScalarField> dmdt() const;