diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.C
index d5d805c085..0bb29c1bea 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.C
@@ -26,6 +26,8 @@ License
 #include "turbulentDispersionModel.H"
 #include "phasePair.H"
 #include "fvcGrad.H"
+#include "surfaceInterpolate.H"
+#include "fvcSnGrad.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -66,4 +68,11 @@ Foam::turbulentDispersionModel::F() const
 }
 
 
+Foam::tmp<Foam::surfaceScalarField>
+Foam::turbulentDispersionModel::Ff() const
+{
+    return fvc::interpolate(D())*fvc::snGrad(pair_.dispersed());
+}
+
+
 // ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.H b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.H
index 932fdb3801..56e72ebcac 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.H
@@ -120,6 +120,9 @@ public:
 
         //- Turbulent dispersion force
         virtual tmp<volVectorField> F() const;
+
+        //- Turbulent dispersion force on faces
+        virtual tmp<surfaceScalarField> Ff() const;
 };
 
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
index dc9335f7ea..f43a3a0354 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
@@ -38,6 +38,7 @@ License
 #include "fvmDiv.H"
 #include "fvmSup.H"
 #include "fvcDiv.H"
+#include "fvcSnGrad.H"
 #include "fvMatrix.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -575,24 +576,75 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Fs() const
         Fs[pair.phase2().index()] -= F;
     }
 
+    return tFs;
+}
+
+
+template<class BasePhaseSystem>
+Foam::autoPtr<Foam::PtrList<Foam::surfaceScalarField> >
+Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::phiDs
+(
+    const PtrList<volScalarField>& rAUs
+) const
+{
+    autoPtr<PtrList<surfaceScalarField> > tphiDs
+    (
+        new PtrList<surfaceScalarField>(this->phases().size())
+    );
+
+    PtrList<surfaceScalarField>& phiDs = tphiDs();
+
+    forAll(phiDs, phasei)
+    {
+        phiDs.set
+        (
+            phasei,
+            new surfaceScalarField
+            (
+                IOobject
+                (
+                    liftModel::typeName + ":F",
+                    this->mesh_.time().timeName(),
+                    this->mesh_,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                this->mesh_,
+                dimensionedScalar
+                (
+                    "zero",
+                    dimTime*dimArea*turbulentDispersionModel::dimF/dimDensity,
+                    0
+                )
+            )
+        );
+    }
+
     // Add the turbulent dispersion force
-    // forAllConstIter
-    // (
-    //     turbulentDispersionModelTable,
-    //     turbulentDispersionModels_,
-    //     turbulentDispersionModelIter
-    // )
-    // {
-    //     const volVectorField F(turbulentDispersionModelIter()->F<vector>());
-
-    //     const phasePair&
-    //         pair(this->phasePairs_[turbulentDispersionModelIter.key()]);
-
-    //     *eqns[pair.phase1().name()] += F;
-    //     *eqns[pair.phase2().name()] -= F;
-    // }
+    forAllConstIter
+    (
+        turbulentDispersionModelTable,
+        turbulentDispersionModels_,
+        turbulentDispersionModelIter
+    )
+    {
+        const phasePair&
+            pair(this->phasePairs_[turbulentDispersionModelIter.key()]);
 
-    return tFs;
+        const volScalarField D(turbulentDispersionModelIter()->D());
+        const surfaceScalarField snGradAlpha1
+        (
+            fvc::snGrad(pair.phase1())*this->mesh_.magSf()
+        );
+
+        phiDs[pair.phase1().index()] +=
+            fvc::interpolate(rAUs[pair.phase1().index()]*D)*snGradAlpha1;
+        phiDs[pair.phase2().index()] -=
+            fvc::interpolate(rAUs[pair.phase2().index()]*D)*snGradAlpha1;
+    }
+
+    return tphiDs;
 }
 
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H
index 194602dfc9..980f854a9e 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H
@@ -175,6 +175,12 @@ public:
         //- Return the combined force (lift + wall-lubrication)
         virtual autoPtr<PtrList<volVectorField> > Fs() const;
 
+        //- Return the turbulent dispersion force on faces for phase pair
+        virtual autoPtr<PtrList<surfaceScalarField> > phiDs
+        (
+            const PtrList<volScalarField>& rAUs
+        ) const;
+
         //- Return the combined face-force (lift + wall-lubrication)
         virtual tmp<surfaceScalarField> Ff(const phasePairKey& key) const;
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H
index 9b7e3da3c1..f4c185d9b7 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H
@@ -169,6 +169,12 @@ public:
         //- Return the combined force (lift + wall-lubrication) for phase pair
         virtual autoPtr<PtrList<Foam::volVectorField> > Fs() const = 0;
 
+        //- Return the turbulent dispersion force on faces for phase pair
+        virtual autoPtr<PtrList<Foam::surfaceScalarField> > phiDs
+        (
+            const PtrList<volScalarField>& rAUs
+        ) 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 9fd728ec68..136d595418 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
@@ -34,10 +34,11 @@ forAll(phases, phasei)
     );
 }
 
-// Turbulent diffusion, particle-pressure, lift and wall-lubrication fluxes
+// Lift, wall-lubrication and turbulent diffusion fluxes
 PtrList<surfaceScalarField> phiFs(phases.size());
 {
     autoPtr<PtrList<volVectorField> > Fs = fluid.Fs();
+    autoPtr<PtrList<surfaceScalarField> > phiDs = fluid.phiDs(rAUs);
 
     forAll(phases, phasei)
     {
@@ -50,6 +51,7 @@ PtrList<surfaceScalarField> phiFs(phases.size());
             (
                 IOobject::groupName("phiF", phase.name()),
                 (fvc::interpolate(rAUs[phasei]*Fs()[phasei]) & mesh.Sf())
+              + phiDs()[phasei]
             )
         );
     }