diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C
index 9eddbe4efd..3db07ef804 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C
@@ -151,6 +151,16 @@ Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+template<class BasePhaseSystem>
+bool Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::transfersMass
+(
+    const phaseModel& phase
+) const
+{
+    return true;
+}
+
+
 template<class BasePhaseSystem>
 Foam::tmp<Foam::volScalarField>
 Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::dmdt
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.H
index d8aa29676a..3a0cbd4ba2 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.H
@@ -120,6 +120,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(const phasePairKey& key) const;
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C
index f4a1284db8..95636c3749 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C
@@ -92,20 +92,7 @@ Foam::HeatTransferPhaseSystem<BasePhaseSystem>::dmdt
     const Foam::phaseModel& phase
 ) const
 {
-    return tmp<volScalarField>
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                IOobject::groupName("dmdt", phase.name()),
-                this->mesh().time().timeName(),
-                this->mesh().time()
-            ),
-            this->mesh(),
-            dimensionedScalar("zero", dimDensity/dimTime, 0)
-        )
-    );
+    return tmp<volScalarField>(NULL);
 }
 
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
index f0453940c7..175a993e30 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
@@ -25,6 +25,7 @@ License
 
 #include "ThermalPhaseChangePhaseSystem.H"
 #include "fvCFD.H"
+#include "alphatPhaseChangeWallFunctionFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -39,7 +40,61 @@ ThermalPhaseChangePhaseSystem
     volatile_(this->lookup("volatile")),
     saturationModel_(saturationModel::New(this->subDict("saturationModel"))),
     massTransfer_(this->lookup("massTransfer"))
-{}
+{
+
+    forAllConstIter
+    (
+        phaseSystem::phasePairTable,
+        this->phasePairs_,
+        phasePairIter
+    )
+    {
+        const phasePair& pair(phasePairIter());
+
+        if (pair.ordered())
+        {
+            continue;
+        }
+
+        // Initially assume no mass transfer
+        iDmdt_.insert
+        (
+            pair,
+            new volScalarField
+            (
+                IOobject
+                (
+                    IOobject::groupName("iDmdt", pair.name()),
+                    this->mesh().time().timeName(),
+                    this->mesh(),
+                    IOobject::NO_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                this->mesh(),
+                dimensionedScalar("zero", dimDensity/dimTime, 0)
+            )
+        );
+
+        // Initially assume no mass transfer
+        wDmdt_.insert
+        (
+            pair,
+            new volScalarField
+            (
+                IOobject
+                (
+                    IOobject::groupName("wDmdt", pair.name()),
+                    this->mesh().time().timeName(),
+                    this->mesh(),
+                    IOobject::NO_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                this->mesh(),
+                dimensionedScalar("zero", dimDensity/dimTime, 0)
+            )
+        );
+    }
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -52,6 +107,14 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+template<class BasePhaseSystem>
+const Foam::saturationModel&
+Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::saturation() const
+{
+    return saturationModel_();
+}
+
+
 template<class BasePhaseSystem>
 Foam::autoPtr<Foam::phaseSystem::massTransferTable>
 Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::massTransfer() const
@@ -121,6 +184,9 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::massTransfer() const
 template<class BasePhaseSystem>
 void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
 {
+    typedef compressible::alphatPhaseChangeWallFunctionFvPatchScalarField
+        alphatPhaseChangeWallFunction;
+
     BasePhaseSystem::correctThermo();
 
     forAllConstIter
@@ -147,6 +213,8 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
         const volScalarField& he2(phase2.thermo().he());
 
         volScalarField& dmdt(*this->dmdt_[pair]);
+        volScalarField& iDmdt(*this->iDmdt_[pair]);
+        volScalarField& wDmdt(*this->wDmdt_[pair]);
 
         volScalarField& Tf = *this->Tf_[pair];
 
@@ -167,25 +235,28 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
 
             Tf = saturationModel_->Tsat(phase1.thermo().p());
 
-            dmdt =
-                (H1*(Tf - T1) + H2*(Tf - T2))
+            scalar iDmdtRelax(this->mesh().fieldRelaxationFactor("iDmdt"));
+
+            iDmdt =
+                (1 - iDmdtRelax)*iDmdt
+              + iDmdtRelax*(H1*(Tf - T1) + H2*(Tf - T2))
                /min
                 (
-                    (pos(dmdt)*he2 + neg(dmdt)*hef2)
-                  - (neg(dmdt)*he1 + pos(dmdt)*hef1),
+                    (pos(iDmdt)*he2 + neg(iDmdt)*hef2)
+                  - (neg(iDmdt)*he1 + pos(iDmdt)*hef1),
                     0.3*mag(hef2 - hef1)
                 );
 
-            Info<< "dmdt." << pair.name()
-                << ": min = " << min(dmdt.internalField())
-                << ", mean = " << average(dmdt.internalField())
-                << ", max = " << max(dmdt.internalField())
-                << ", integral = " << fvc::domainIntegrate(dmdt).value()
+            Info<< "iDmdt." << pair.name()
+                << ": min = " << min(iDmdt.internalField())
+                << ", mean = " << average(iDmdt.internalField())
+                << ", max = " << max(iDmdt.internalField())
+                << ", integral = " << fvc::domainIntegrate(iDmdt).value()
                 << endl;
         }
         else
         {
-            dmdt == dimensionedScalar("0", dmdt.dimensions(), 0);
+            iDmdt == dimensionedScalar("0", dmdt.dimensions(), 0);
         }
 
         volScalarField H1(this->heatTransferModels_[pair][pair.first()]->K());
@@ -200,12 +271,13 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
 
         volScalarField mDotL
         (
-            dmdt*
+            iDmdt*
             (
-                (pos(dmdt)*he2 + neg(dmdt)*hef2)
-              - (neg(dmdt)*he1 + pos(dmdt)*hef1)
+                (pos(iDmdt)*he2 + neg(iDmdt)*hef2)
+              - (neg(iDmdt)*he1 + pos(iDmdt)*hef1)
             )
         );
+
         Tf = (H1*T1 + H2*T2 + mDotL)/(H1 + H2);
 
         Info<< "Tf." << pair.name()
@@ -213,6 +285,68 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
             << ", mean = " << average(Tf.internalField())
             << ", max = " << max(Tf.internalField())
             << endl;
+
+        // Accumulate dmdt contributions from boundaries
+        if
+        (
+            phase2.mesh().foundObject<volScalarField>
+            (
+                "alphat." +  phase2.name()
+            )
+        )
+        {
+            scalar wDmdtRelax(this->mesh().fieldRelaxationFactor("wDmdt"));
+            wDmdt *= (1 - wDmdtRelax);
+
+            const volScalarField& alphat =
+                phase2.mesh().lookupObject<volScalarField>
+                (
+                    "alphat." +  phase2.name()
+                );
+
+            const fvPatchList& patches = this->mesh().boundary();
+            forAll(patches, patchi)
+            {
+                const fvPatch& currPatch = patches[patchi];
+
+                if
+                (
+                    isA<alphatPhaseChangeWallFunction>
+                    (
+                        alphat.boundaryField()[patchi]
+                    )
+                )
+                {
+                    const scalarField& patchDmdt =
+                        refCast<const alphatPhaseChangeWallFunction>
+                        (
+                            alphat.boundaryField()[patchi]
+                        ).dmdt();
+
+                    forAll(patchDmdt,facei)
+                    {
+                        label faceCelli = currPatch.faceCells()[facei];
+                        wDmdt[faceCelli] += wDmdtRelax*patchDmdt[facei];
+                    }
+                }
+            }
+
+            Info<< "wDmdt." << pair.name()
+                << ": min = " << min(wDmdt.internalField())
+                << ", mean = " << average(wDmdt.internalField())
+                << ", max = " << max(wDmdt.internalField())
+                << ", integral = " << fvc::domainIntegrate(wDmdt).value()
+                << endl;
+        }
+
+        dmdt = wDmdt + iDmdt;
+
+        Info<< "dmdt." << pair.name()
+            << ": min = " << min(dmdt.internalField())
+            << ", mean = " << average(dmdt.internalField())
+            << ", max = " << max(dmdt.internalField())
+            << ", integral = " << fvc::domainIntegrate(dmdt).value()
+            << endl;
     }
 }
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H
index edae5a6ff6..e67b8209fe 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H
@@ -73,6 +73,14 @@ protected:
         // Mass transfer enabled
         Switch massTransfer_;
 
+        //- Interfacial Mass transfer rate
+        HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>
+            iDmdt_;
+
+        //- Wall Mass transfer rate
+        HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>
+            wDmdt_;
+
 
 public:
 
@@ -88,6 +96,9 @@ public:
 
     // Member Functions
 
+        //- Return the saturationModel
+        const saturationModel& saturation() const;
+
         //- Return the mass transfer matrices
         virtual autoPtr<phaseSystem::massTransferTable> massTransfer() const;
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
index 78e015ff3a..8d705049f7 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
@@ -37,17 +37,7 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::AnisothermalPhaseModel
 )
 :
     BasePhaseModel(fluid, phaseName, index),
-    divU_
-    (
-        IOobject
-        (
-            IOobject::groupName("divU", this->name()),
-            fluid.mesh().time().timeName(),
-            fluid.mesh()
-        ),
-        fluid.mesh(),
-        dimensionedScalar("divU", dimless/dimTime, 0)
-    ),
+    divU_(NULL),
     K_
     (
         IOobject
@@ -105,7 +95,7 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
 
     volScalarField& he = this->thermo_->he();
 
-    return
+    tmp<fvScalarMatrix> tEEqn
     (
         fvm::ddt(alpha, this->rho(), he) + fvm::div(alphaRhoPhi, he)
       - fvm::Sp(contErr, he)
@@ -119,16 +109,23 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
            *fvc::interpolate(alphaEff),
             he
         )
-
-      + (
-            he.name() == this->thermo_->phasePropertyName("e")
-          ? fvc::ddt(alpha)*this->thermo().p()
-          + fvc::div(alphaPhi, this->thermo().p())
-          : -alpha*this->fluid().dpdt()
-        )
      ==
         this->Sh()
     );
+
+    // Add the appropriate pressure-work term
+    if (he.name() == this->thermo_->phasePropertyName("e"))
+    {
+        tEEqn() +=
+            fvc::ddt(alpha)*this->thermo().p()
+          + fvc::div(alphaPhi, this->thermo().p());
+    }
+    else if (this->thermo_->dpdt())
+    {
+        tEEqn() -= alpha*this->fluid().dpdt();
+    }
+
+    return tEEqn;
 }
 
 
@@ -140,7 +137,7 @@ bool Foam::AnisothermalPhaseModel<BasePhaseModel>::compressible() const
 
 
 template<class BasePhaseModel>
-const Foam::volScalarField&
+const Foam::tmp<Foam::volScalarField>&
 Foam::AnisothermalPhaseModel<BasePhaseModel>::divU() const
 {
     return divU_;
@@ -149,7 +146,10 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::divU() const
 
 template<class BasePhaseModel>
 void
-Foam::AnisothermalPhaseModel<BasePhaseModel>::divU(const volScalarField& divU)
+Foam::AnisothermalPhaseModel<BasePhaseModel>::divU
+(
+    const tmp<volScalarField>& divU
+)
 {
     divU_ = divU;
 }
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H
index 9d7a92b499..a8c4be694a 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H
@@ -55,7 +55,7 @@ class AnisothermalPhaseModel
     // Private data
 
         //- Dilatation rate
-        volScalarField divU_;
+        tmp<volScalarField> divU_;
 
         //- Kinetic energy
         volScalarField K_;
@@ -95,10 +95,10 @@ public:
             virtual bool compressible() const;
 
             //- Return the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
-            virtual const volScalarField& divU() const;
+            virtual const tmp<volScalarField>& divU() const;
 
             //- Set the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
-            virtual void divU(const volScalarField& divU);
+            virtual void divU(const tmp<volScalarField>& divU);
 
             //- Return the phase kinetic energy
             virtual const volScalarField& K() const;
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C
index 0849867944..eb2e871530 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C
@@ -165,16 +165,17 @@ bool Foam::phaseModel::compressible() const
 }
 
 
-const Foam::volScalarField& Foam::phaseModel::divU() const
+const Foam::tmp<Foam::volScalarField>& Foam::phaseModel::divU() const
 {
     notImplemented("Foam::phaseModel::divU()");
-    return volScalarField::null();
+    static tmp<Foam::volScalarField> divU_(NULL);
+    return divU_;
 }
 
 
-void Foam::phaseModel::divU(const volScalarField& divU)
+void Foam::phaseModel::divU(const tmp<volScalarField>& divU)
 {
-    WarningIn("phaseModel::divU(const volScalarField& divU)")
+    WarningIn("phaseModel::divU(const tmp<volScalarField>& divU)")
         << "Attempt to set the dilatation rate of an incompressible phase"
         << endl;
 }
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
index ae94353aeb..27d78390ca 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
@@ -216,10 +216,10 @@ public:
             virtual bool compressible() const;
 
             //- Return the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
-            virtual const volScalarField& divU() const;
+            virtual const tmp<volScalarField>& divU() const;
 
             //- Set the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
-            virtual void divU(const volScalarField& divU);
+            virtual void divU(const tmp<volScalarField>& divU);
 
             //- Return the phase kinetic energy
             virtual const volScalarField& K() const;
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/reactionThermo/hRefConstThermos.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/reactionThermo/hRefConstThermos.C
index 47ae733b14..76c2351109 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/reactionThermo/hRefConstThermos.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/reactionThermo/hRefConstThermos.C
@@ -120,6 +120,19 @@ constTransport
     >
 > constRefRhoConstEThermoPhysics;
 
+typedef
+constTransport
+<
+    species::thermo
+    <
+        hRefConstThermo
+        <
+            rhoConst<specie>
+        >,
+        sensibleEnthalpy
+    >
+> constRefRhoConstHThermoPhysics;
+
 
 // pureMixture, sensibleEnthalpy:
 
@@ -229,6 +242,18 @@ makeReactionMixtureThermo
 );
 
 
+// multiComponentMixture, sensibleEnthalpy:
+
+makeReactionMixtureThermo
+(
+    rhoThermo,
+    rhoReactionThermo,
+    heRhoThermo,
+    multiComponentMixture,
+    constRefRhoConstHThermoPhysics
+);
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index cf765d4c07..a79e9393b5 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
@@ -204,7 +204,7 @@ void Foam::multiphaseSystem::solveAlphas()
                 mesh_
             ),
             mesh_,
-            dimensionedScalar("Sp", phase.divU().dimensions(), 0.0)
+            dimensionedScalar("Sp", divU.dimensions(), 0.0)
         );
 
         volScalarField::DimensionedInternalField Su
@@ -220,9 +220,9 @@ void Foam::multiphaseSystem::solveAlphas()
             divU*min(alpha, scalar(1))
         );
 
-        if (phase.compressible())
+        if (phase.divU().valid())
         {
-            const scalarField& dgdt = phase.divU();
+            const scalarField& dgdt = phase.divU()();
 
             forAll(dgdt, celli)
             {
@@ -247,9 +247,9 @@ void Foam::multiphaseSystem::solveAlphas()
 
             if (&phase2 == &phase) continue;
 
-            if (phase2.compressible())
+            if (phase2.divU().valid())
             {
-                const scalarField& dgdt2 = phase2.divU();
+                const scalarField& dgdt2 = phase2.divU()();
 
                 forAll(dgdt2, celli)
                 {
@@ -515,6 +515,12 @@ 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 f4c185d9b7..a9390f2589 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H
@@ -175,6 +175,9 @@ public:
             const PtrList<volScalarField>& rAUs
         ) const = 0;
 
+        //- Return true if there is mass transfer for phase
+        virtual bool transfersMass(const phaseModel& phase) const;
+
         //- 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 e23746a3ad..230b650106 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
@@ -311,7 +311,7 @@ while (pimple.correct())
                     phasei,
                     (
                         (
-                            phase.continuityError() - fluid.dmdt(phase)
+                            phase.continuityError()
                           - fvc::Sp
                             (
                                 fvc::ddt(alpha) + fvc::div(phase.alphaPhi()),
@@ -340,7 +340,7 @@ while (pimple.correct())
                     phasei,
                     (
                         (
-                            phase.continuityError() - fluid.dmdt(phase)
+                            phase.continuityError()
                           - fvc::Sp
                             (
                                 (fvc::ddt(alpha) + fvc::div(phase.alphaPhi())),
@@ -352,6 +352,18 @@ while (pimple.correct())
                     ).ptr()
                 );
             }
+
+            if (fluid.transfersMass(phase))
+            {
+                if (pEqnComps.set(phasei))
+                {
+                    pEqnComps[phasei] -= fluid.dmdt(phase);
+                }
+                else
+                {
+                    pEqnComps.set(phasei, fvm::Su(-fluid.dmdt(phase), rho));
+                }
+            }
         }
     }
 
@@ -373,9 +385,7 @@ while (pimple.correct())
 
             forAll(phases, phasei)
             {
-                phaseModel& phase = phases[phasei];
-
-                if (phase.compressible())
+                if (pEqnComps.set(phasei))
                 {
                     pEqn += pEqnComps[phasei];
                 }
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
index 5ffc70ec84..c8a9e28288 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
@@ -321,12 +321,12 @@ while (pimple.correct())
         {
             fvScalarMatrix pEqn(pEqnIncomp);
 
-            if (phase1.compressible())
+            if (pEqnComp1.valid())
             {
                 pEqn += pEqnComp1();
             }
 
-            if (phase2.compressible())
+            if (pEqnComp2.valid())
             {
                 pEqn += pEqnComp2();
             }
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
index 161a08b06d..c3f7d472fa 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
@@ -227,7 +227,9 @@ while (pimple.correct())
             fvc::interpolate(psi2)*phi2
         );
 
-        pEqnComp1 =
+        if (phase1.compressible())
+        {
+            pEqnComp1 =
             (
                 phase1.continuityError() - fluid.dmdt()
               - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
@@ -237,10 +239,13 @@ while (pimple.correct())
                 psi1*fvm::ddt(p_rgh)
               + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
             );
-        deleteDemandDrivenData(pEqnComp1().faceFluxCorrectionPtr());
-        pEqnComp1().relax();
+            deleteDemandDrivenData(pEqnComp1().faceFluxCorrectionPtr());
+            pEqnComp1().relax();
+        }
 
-        pEqnComp2 =
+        if (phase2.compressible())
+        {
+            pEqnComp2 =
             (
                 phase2.continuityError() + fluid.dmdt()
               - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
@@ -250,24 +255,31 @@ while (pimple.correct())
                 psi2*fvm::ddt(p_rgh)
               + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
             );
-        deleteDemandDrivenData(pEqnComp2().faceFluxCorrectionPtr());
-        pEqnComp2().relax();
+            deleteDemandDrivenData(pEqnComp2().faceFluxCorrectionPtr());
+            pEqnComp2().relax();
+        }
     }
     else
     {
-        pEqnComp1 =
+        if (phase1.compressible())
+        {
+            pEqnComp1 =
             (
                 phase1.continuityError() - fluid.dmdt()
               - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
             )/rho1
           + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
+        }
 
-        pEqnComp2 =
+        if (phase2.compressible())
+        {
+            pEqnComp2 =
             (
                 phase2.continuityError() + fluid.dmdt()
               - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
             )/rho2
           + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh));
+        }
     }
 
     // Cache p prior to solve for density update
@@ -281,11 +293,25 @@ while (pimple.correct())
           - fvm::laplacian(rAUf, p_rgh)
         );
 
-        solve
-        (
-            pEqnComp1() + pEqnComp2() + pEqnIncomp,
-            mesh.solver(p_rgh.select(pimple.finalInnerIter()))
-        );
+        {
+            fvScalarMatrix pEqn(pEqnIncomp);
+
+            if (pEqnComp1.valid())
+            {
+                pEqn += pEqnComp1();
+            }
+
+            if (pEqnComp2.valid())
+            {
+                pEqn += pEqnComp2();
+            }
+
+            solve
+            (
+                pEqn,
+                mesh.solver(p_rgh.select(pimple.finalInnerIter()))
+            );
+        }
 
         if (pimple.finalNonOrthogonalIter())
         {
@@ -325,17 +351,23 @@ while (pimple.correct())
             fvOptions.correct(U2);
 
             // Set the phase dilatation rates
-            if (phase1.compressible())
+            if (pEqnComp1.valid())
             {
                 phase1.divU(-pEqnComp1 & p_rgh);
             }
-            if (phase2.compressible())
+            if (pEqnComp2.valid())
             {
                 phase2.divU(-pEqnComp2 & p_rgh);
             }
         }
     }
 
+    Info<< "min(p) = " << min(p_rgh + rho*gh).value() << endl;
+    if (min(p_rgh + rho*gh) < pMin)
+    {
+        Info<< "Clipping p" << endl;
+    }
+
     // Update and limit the static pressure
     p = max(p_rgh + rho*gh, pMin);
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/Make/files b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/Make/files
index d1cbb6eaeb..10a0e8679a 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/Make/files
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/Make/files
@@ -35,4 +35,6 @@ kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
 kineticTheoryModels/derivedFvPatchFields/JohnsonJacksonParticleTheta/JohnsonJacksonParticleThetaFvPatchScalarField.C
 kineticTheoryModels/derivedFvPatchFields/JohnsonJacksonParticleSlip/JohnsonJacksonParticleSlipFvPatchVectorField.C
 
+derivedFvPatchFields/alphatFixedDmdtWallBoilingWallFunction/alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField.C
+
 LIB = $(FOAM_LIBBIN)/libtwoPhaseReactingTurbulenceModels
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/Make/options b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/Make/options
index cd3f851fc8..f003ce2281 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/Make/options
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/Make/options
@@ -2,6 +2,7 @@ EXE_INC = \
     -I../twoPhaseSystem/lnInclude \
     -I../../phaseSystems/lnInclude \
     -I../../interfacialModels/lnInclude\
+    -I../../interfacialCompositionModels/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/transportModels/incompressible/transportModel \
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
index 401c4f622c..187ca4c238 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
@@ -203,30 +203,30 @@ void Foam::twoPhaseSystem::solve()
     // Construct the dilatation rate source term
     tmp<volScalarField::DimensionedInternalField> tdgdt;
 
-    if (phase1_.compressible() && phase2_.compressible())
+    if (phase1_.divU().valid() && phase2_.divU().valid())
     {
         tdgdt =
         (
             alpha2.dimensionedInternalField()
-           *phase1_.divU().dimensionedInternalField()
+           *phase1_.divU()().dimensionedInternalField()
           - alpha1.dimensionedInternalField()
-           *phase2_.divU().dimensionedInternalField()
+           *phase2_.divU()().dimensionedInternalField()
         );
     }
-    else if (phase1_.compressible())
+    else if (phase1_.divU().valid())
     {
         tdgdt =
         (
             alpha2.dimensionedInternalField()
-           *phase1_.divU().dimensionedInternalField()
+           *phase1_.divU()().dimensionedInternalField()
         );
     }
-    else if (phase2_.compressible())
+    else if (phase2_.divU().valid())
     {
         tdgdt =
         (
           - alpha1.dimensionedInternalField()
-           *phase2_.divU().dimensionedInternalField()
+           *phase2_.divU()().dimensionedInternalField()
         );
     }