diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
index 864b9c7afb..21a38ea5a6 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
@@ -57,16 +57,12 @@
       + fvc::ddtPhiCorr(rAU2, alpha2, U2, phi2)
     );
 
-    phi = alpha1f*phiHbyA1 + alpha2f*phiHbyA2;
-    mrfZones.makeRelative(phi);
-
     phiHbyA1 +=
     (
         fvc::interpolate((1.0/rho1)*rAU1*dragCoeff)*phi2
       - phiP1
       + rAlphaAU1f*(g & mesh.Sf())
     );
-    mrfZones.makeRelative(phiHbyA1);
 
     phiHbyA2 +=
     (
@@ -74,8 +70,9 @@
       - phiP2
       + rAlphaAU2f*(g & mesh.Sf())
     );
-    mrfZones.makeRelative(phiHbyA2);
 
+    mrfZones.makeRelative(phiHbyA1);
+    mrfZones.makeRelative(phiHbyA2);
     mrfZones.makeRelative(phi1.oldTime());
     mrfZones.makeRelative(phi1);
     mrfZones.makeRelative(phi2.oldTime());
@@ -83,6 +80,16 @@
 
     surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2);
 
+    // Update phi BCs before pEqn
+    phi.boundaryField() =
+        mrfZones.relative
+        (
+            alpha1f.boundaryField()
+           *(mesh.Sf().boundaryField() & U1.boundaryField())
+          + alpha2f.boundaryField()
+           *(mesh.Sf().boundaryField() & U2.boundaryField())
+        );
+
     HbyA1 += (1.0/rho1)*rAU1*dragCoeff*U2;
     HbyA2 += (1.0/rho2)*rAU2*dragCoeff*U1;
 
@@ -171,11 +178,17 @@
             surfaceScalarField mSfGradp(pEqnIncomp.flux()/Dp);
 
             phi1.boundaryField() ==
-                (mesh.Sf().boundaryField() & U1.boundaryField());
+                mrfZones.relative
+                (
+                    mesh.Sf().boundaryField() & U1.boundaryField()
+                );
             phi1 = phiHbyA1 + rAlphaAU1f*mSfGradp/fvc::interpolate(rho1);
 
             phi2.boundaryField() ==
-                (mesh.Sf().boundaryField() & U2.boundaryField());
+                mrfZones.relative
+                (
+                    mesh.Sf().boundaryField() & U2.boundaryField()
+                );
             phi2 = phiHbyA2 + rAlphaAU2f*mSfGradp/fvc::interpolate(rho2);
 
             phi = alpha1f*phi1 + alpha2f*phi2;