diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index 620499a2e7b0997f0d1ff0fd09c60c60b7b0bcc0..b55ffed9f7f58c4d372c334358c87df3faf394cb 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
@@ -58,8 +58,6 @@ void Foam::multiphaseSystem::calcAlphas()
 
 void Foam::multiphaseSystem::solveAlphas()
 {
-    surfaceScalarField phic(mag(phi_/mesh_.magSf()));
-
     PtrList<surfaceScalarField> phiAlphaCorrs(phases_.size());
     int phasei = 0;
 
@@ -94,6 +92,11 @@ void Foam::multiphaseSystem::solveAlphas()
 
             if (&phase2 == &phase1) continue;
 
+            surfaceScalarField phic
+            (
+                (mag(phi_) + mag(phase1.phi() - phase2.phi()))/mesh_.magSf()
+            );
+
             surfaceScalarField phir
             (
                 (phase1.phi() - phase2.phi())
@@ -186,32 +189,6 @@ void Foam::multiphaseSystem::solveAlphas()
 }
 
 
-Foam::dimensionedScalar Foam::multiphaseSystem::sigma
-(
-    const phaseModel& phase1,
-    const phaseModel& phase2
-) const
-{
-    scalarCoeffTable::const_iterator sigma
-    (
-        sigmas_.find(interfacePair(phase1, phase2))
-    );
-
-    if (sigma == sigmas_.end())
-    {
-        FatalErrorIn
-        (
-            "multiphaseSystem::sigma(const phaseModel& phase1,"
-            "const phaseModel& phase2) const"
-        )   << "Cannot find interface " << interfacePair(phase1, phase2)
-            << " in list of sigma values"
-            << exit(FatalError);
-    }
-
-    return dimensionedScalar("sigma", dimSigma_, sigma());
-}
-
-
 Foam::scalar Foam::multiphaseSystem::cAlpha
 (
     const phaseModel& phase1,
@@ -263,10 +240,10 @@ Foam::dimensionedScalar Foam::multiphaseSystem::Cvm
 
     FatalErrorIn
     (
-        "multiphaseSystem::sigma"
+        "multiphaseSystem::Cvm"
         "(const phaseModel& phase1, const phaseModel& phase2) const"
     )   << "Cannot find interface " << interfacePair(phase1, phase2)
-        << " in list of sigma values"
+        << " in list of Cvm values"
         << exit(FatalError);
 
     return Cvm()*phase2.rho();
@@ -729,52 +706,56 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::dragCoeff
 }
 
 
-Foam::tmp<Foam::surfaceScalarField>
-Foam::multiphaseSystem::surfaceTensionForce() const
+Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseSystem::surfaceTension
+(
+    const phaseModel& phase1
+) const
 {
-    tmp<surfaceScalarField> tstf
+    tmp<surfaceScalarField> tSurfaceTension
     (
         new surfaceScalarField
         (
             IOobject
             (
-                "surfaceTensionForce",
+                "surfaceTension",
                 mesh_.time().timeName(),
                 mesh_
             ),
             mesh_,
             dimensionedScalar
             (
-                "surfaceTensionForce",
+                "surfaceTension",
                 dimensionSet(1, -2, -2, 0, 0),
-                0.0
+                0
             )
         )
     );
 
-    surfaceScalarField& stf = tstf();
-
-    forAllConstIter(PtrDictionary<phaseModel>, phases_, iter1)
+    forAllConstIter(PtrDictionary<phaseModel>, phases_, iter)
     {
-        const phaseModel& phase1 = iter1();
-
-        PtrDictionary<phaseModel>::const_iterator iter2 = iter1;
-        ++iter2;
+        const phaseModel& phase2 = iter();
 
-        for (; iter2 != phases_.end(); ++iter2)
+        if (&phase2 != &phase1)
         {
-            const phaseModel& phase2 = iter2();
+            scalarCoeffTable::const_iterator sigma
+            (
+                sigmas_.find(interfacePair(phase1, phase2))
+            );
 
-            stf += sigma(phase1, phase2)
-               *fvc::interpolate(K(phase1, phase2))*
-                (
-                    fvc::interpolate(phase2)*fvc::snGrad(phase1)
-                  - fvc::interpolate(phase1)*fvc::snGrad(phase2)
-                );
+            if (sigma != sigmas_.end())
+            {
+                tSurfaceTension() +=
+                    dimensionedScalar("sigma", dimSigma_, sigma())
+                   *fvc::interpolate(K(phase1, phase2))*
+                    (
+                        fvc::interpolate(phase2)*fvc::snGrad(phase1)
+                      - fvc::interpolate(phase1)*fvc::snGrad(phase2)
+                    );
+            }
         }
     }
 
-    return tstf;
+    return tSurfaceTension;
 }
 
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H
index 4c3e52e0ea548db8647d3a625e5a143133b7d638..4e2d0193850a28f66927f96239d0aee0756e06cf 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H
@@ -196,12 +196,6 @@ private:
 
         void solveAlphas();
 
-        dimensionedScalar sigma
-        (
-            const phaseModel& phase1,
-            const phaseModel& phase2
-        ) const;
-
         scalar cAlpha
         (
             const phaseModel& phase1,
@@ -299,7 +293,7 @@ public:
             const dragCoeffFields& dragCoeffs
         ) const;
 
-        tmp<surfaceScalarField> surfaceTensionForce() const;
+        tmp<surfaceScalarField> surfaceTension(const phaseModel& phase) const;
 
         //- Indicator of the proximity of the interface
         //  Field values are 1 near and 0 away for the interface.
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
index 6dc67ce353c5ca4b2080cac6cf28ef4663084492..9ed3fecce9c41af14c013fe65390d0b795703f9c 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
@@ -74,7 +74,11 @@
         (
             (fvc::interpolate(HbyAs[phasei]) & mesh.Sf())
           + fvc::ddtPhiCorr(rAUs[phasei], alpha, phase.U(), phase.phi())
-          + rAlphaAUfs[phasei]*(g & mesh.Sf())
+          + rAlphaAUfs[phasei]
+           *(
+               fluid.surfaceTension(phase)*mesh.magSf()/phase.rho()
+             + (g & mesh.Sf())
+            )
         );
 
         mrfZones.relativeFlux(phiHbyAs[phasei]);