diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/pEqn.H
index b9029e4190..e352bbe914 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/pEqn.H
@@ -39,7 +39,7 @@ tmp<surfaceScalarField> phiF2;
     volScalarField D(fluid.D());
 
     // Phase-1 turbulent dispersion and particle-pressure flux
-    surfaceScalarField Df1
+    tmp<surfaceScalarField> DbyA1
     (
         fvc::interpolate
         (
@@ -48,7 +48,7 @@ tmp<surfaceScalarField> phiF2;
     );
 
     // Phase-2 turbulent dispersion and particle-pressure flux
-    surfaceScalarField Df2
+    tmp<surfaceScalarField> DbyA2
     (
         fvc::interpolate
         (
@@ -56,13 +56,6 @@ tmp<surfaceScalarField> phiF2;
         )
     );
 
-    // Cache the net diffusivity for implicit diffusion treatment in the
-    // phase-fraction equation
-    if (implicitPhasePressure)
-    {
-        fluid.pPrimeByA() = Df1 + Df2;
-    }
-
     // Lift and wall-lubrication forces
     volVectorField F(fluid.F());
 
@@ -72,16 +65,24 @@ tmp<surfaceScalarField> phiF2;
     // Phase-1 dispersion, lift and wall-lubrication flux
     phiF1 =
     (
-        Df1*snGradAlpha1
+        DbyA1()*snGradAlpha1
       + (fvc::interpolate(rAU1*F) & mesh.Sf())
     );
 
     // Phase-2 dispersion, lift and wall-lubrication flux
     phiF2 =
     (
-      - Df2*snGradAlpha1
+      - DbyA2()*snGradAlpha1
       - (fvc::interpolate(rAU2*F) & mesh.Sf())
     );
+
+    // Cache the phase diffusivities for implicit treatment in the
+    // phase-fraction equation
+    if (implicitPhasePressure)
+    {
+        phase1.DbyA(DbyA1);
+        phase2.DbyA(DbyA2);
+    }
 }
 
 
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H
index 8bd69f606c..161a08b06d 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H
@@ -69,11 +69,12 @@ tmp<surfaceScalarField> Ff2;
         fvc::interpolate(D + phase2.turbulence().pPrime())
     );
 
-    // Cache the net diffusivity for implicit diffusion treatment in the
+    // Cache the phase diffusivities for implicit treatment in the
     // phase-fraction equation
     if (implicitPhasePressure)
     {
-        fluid.pPrimeByA() = rAUf1*Df1 + rAUf2*Df2;
+        phase1.DbyA(rAUf1*Df1);
+        phase2.DbyA(rAUf2*Df2);
     }
 
     // Lift and wall-lubrication forces
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
index 01ddeb772f..cbefea2852 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
@@ -139,7 +139,7 @@ bool Foam::AnisothermalPhaseModel<BasePhaseModel>::compressible() const
 
 
 template<class BasePhaseModel>
-Foam::tmp<Foam::volScalarField>
+const Foam::volScalarField&
 Foam::AnisothermalPhaseModel<BasePhaseModel>::divU() const
 {
     return divU_;
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H
index 1ff8b79d8b..48b0cbd499 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H
@@ -90,7 +90,7 @@ public:
             virtual bool compressible() const;
 
             //- Phase dilatation rate (d(alpha)/dt + div(alpha*phi))
-            virtual tmp<volScalarField> divU() const;
+            virtual const volScalarField& divU() const;
 
             //- Set phase dilatation rate (d(alpha)/dt + div(alpha*phi))
             virtual void divU(const volScalarField& divU);
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
index 3d2e275080..aeeef98149 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
@@ -257,6 +257,31 @@ Foam::MovingPhaseModel<BasePhaseModel>::UEqn()
 }
 
 
+template<class BasePhaseModel>
+const Foam::surfaceScalarField&
+Foam::MovingPhaseModel<BasePhaseModel>::DbyA() const
+{
+    if (DbyA_.valid())
+    {
+        return DbyA_;
+    }
+    else
+    {
+        return surfaceScalarField::null();
+    }
+}
+
+
+template<class BasePhaseModel>
+void Foam::MovingPhaseModel<BasePhaseModel>::DbyA
+(
+    const tmp<surfaceScalarField>& DbyA
+)
+{
+    DbyA_ = DbyA;
+}
+
+
 template<class BasePhaseModel>
 Foam::tmp<Foam::volVectorField>
 Foam::MovingPhaseModel<BasePhaseModel>::U() const
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H
index f8d32d46b3..96e8942f6f 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H
@@ -88,6 +88,10 @@ class MovingPhaseModel
         //- Continuity error
         volScalarField continuityError_;
 
+        //- Phase diffusivity divided by the momentum coefficient.
+        //  Used for implicit treatment of the phase pressure and dispersion
+        tmp<surfaceScalarField> DbyA_;
+
 
     // Private static member functions
 
@@ -124,6 +128,16 @@ public:
         //- Return the momentum equation
         virtual tmp<fvVectorMatrix> UEqn();
 
+
+        // Implicit phase pressure and dispersion support
+
+            //- Return the phase diffusivity divided by the momentum coefficient
+            virtual const surfaceScalarField& DbyA() const;
+
+            //- Set the phase diffusivity divided by the momentum coefficient
+            virtual void DbyA(const tmp<surfaceScalarField>& DbyA);
+
+
         // Momentum
 
             //- Constant access the velocity
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C
index 5f290080d4..2697587b6a 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C
@@ -140,9 +140,10 @@ bool Foam::phaseModel::compressible() const
 }
 
 
-Foam::tmp<Foam::volScalarField> Foam::phaseModel::divU() const
+const Foam::volScalarField& Foam::phaseModel::divU() const
 {
-    return tmp<volScalarField>();
+    notImplemented("Foam::phaseModel::divU()");
+    return volScalarField::null();
 }
 
 
@@ -154,4 +155,18 @@ void Foam::phaseModel::divU(const volScalarField& divU)
 }
 
 
+const Foam::surfaceScalarField& Foam::phaseModel::DbyA() const
+{
+    return surfaceScalarField::null();
+}
+
+
+void Foam::phaseModel::DbyA(const tmp<surfaceScalarField>& DbyA)
+{
+    WarningIn("phaseModel::DbyA(const surfaceScalarField& DbyA)")
+        << "Attempt to set the dilatation rate of an incompressible phase"
+        << endl;
+}
+
+
 // ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
index 4fe4b9d1aa..de85779d55 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
@@ -169,13 +169,22 @@ public:
             //- Return true if the phase is compressible otherwise false
             virtual bool compressible() const;
 
-            //- Phase dilatation rate (d(alpha)/dt + div(alpha*phi))
-            virtual tmp<volScalarField> divU() const;
+            //- 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);
 
 
+        // Implicit phase pressure and dispersion support
+
+            //- Return the phase diffusivity divided by the momentum coefficient
+            virtual const surfaceScalarField& DbyA() const;
+
+            //- Set the phase diffusivity divided by the momentum coefficient
+            virtual void DbyA(const tmp<surfaceScalarField>& DbyA);
+
+
         // Thermo
 
             //- Return const access to the thermophysical model
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystem.C
index 825089539f..48b9703ebe 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystem.C
@@ -190,9 +190,9 @@ void Foam::twoPhaseSystem::solve()
         tdgdt =
         (
             alpha2.dimensionedInternalField()
-           *phase1_.divU()().dimensionedInternalField()
+           *phase1_.divU().dimensionedInternalField()
           - alpha1.dimensionedInternalField()
-           *phase2_.divU()().dimensionedInternalField()
+           *phase2_.divU().dimensionedInternalField()
         );
     }
     else if (phase1_.compressible())
@@ -200,7 +200,7 @@ void Foam::twoPhaseSystem::solve()
         tdgdt =
         (
             alpha2.dimensionedInternalField()
-           *phase1_.divU()().dimensionedInternalField()
+           *phase1_.divU().dimensionedInternalField()
         );
     }
     else if (phase2_.compressible())
@@ -208,7 +208,7 @@ void Foam::twoPhaseSystem::solve()
         tdgdt =
         (
           - alpha1.dimensionedInternalField()
-           *phase2_.divU()().dimensionedInternalField()
+           *phase2_.divU().dimensionedInternalField()
         );
     }
 
@@ -218,20 +218,18 @@ void Foam::twoPhaseSystem::solve()
     surfaceScalarField phic("phic", phi);
     surfaceScalarField phir("phir", phi1 - phi2);
 
-    tmp<surfaceScalarField> alpha1alpha2f;
+    tmp<surfaceScalarField> alphaDbyA;
 
-    if (pPrimeByA_.valid())
+    if (notNull(phase1_.DbyA()) && notNull(phase2_.DbyA()))
     {
-        alpha1alpha2f =
+        surfaceScalarField DbyA(phase1_.DbyA() + phase2_.DbyA());
+
+        alphaDbyA =
             fvc::interpolate(max(alpha1, scalar(0)))
-           *fvc::interpolate(max(alpha2, scalar(0)));
+           *fvc::interpolate(max(alpha2, scalar(0)))
+           *DbyA;
 
-        surfaceScalarField phiP
-        (
-            pPrimeByA_()*fvc::snGrad(alpha1, "bounded")*mesh_.magSf()
-        );
-
-        phir += phiP;
+        phir += DbyA*fvc::snGrad(alpha1, "bounded")*mesh_.magSf();
     }
 
     for (int acorr=0; acorr<nAlphaCorr; acorr++)
@@ -367,12 +365,12 @@ void Foam::twoPhaseSystem::solve()
             phase1_.alphaPhi() = alphaPhic1;
         }
 
-        if (pPrimeByA_.valid())
+        if (alphaDbyA.valid())
         {
             fvScalarMatrix alpha1Eqn
             (
                 fvm::ddt(alpha1) - fvc::ddt(alpha1)
-              - fvm::laplacian(alpha1alpha2f*pPrimeByA_(), alpha1, "bounded")
+              - fvm::laplacian(alphaDbyA, alpha1, "bounded")
             );
 
             alpha1Eqn.relax();
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystem.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystem.H
index 55b7378f4b..0ac971b790 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystem.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystem.H
@@ -63,9 +63,6 @@ protected:
         //- Phase model 2
         phaseModel& phase2_;
 
-        //- Optional dispersion diffusivity
-        tmp<surfaceScalarField> pPrimeByA_;
-
 
 public:
 
@@ -173,9 +170,6 @@ public:
 
             //- Return the virtual mass model for the given phase
             const virtualMassModel& virtualMass(const phaseModel& phase) const;
-
-            //- Return non-const access to the dispersion diffusivity
-            inline tmp<surfaceScalarField>& pPrimeByA();
 };
 
 
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystemI.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystemI.H
index 6a30bc73a3..4cf198a8cd 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystemI.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystemI.H
@@ -65,10 +65,4 @@ inline const Foam::phaseModel& Foam::twoPhaseSystem::otherPhase
 }
 
 
-inline Foam::tmp<Foam::surfaceScalarField>& Foam::twoPhaseSystem::pPrimeByA()
-{
-    return pPrimeByA_;
-}
-
-
 // ************************************************************************* //