diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C
index a553a01b5e..cd60d2cbb0 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C
@@ -248,4 +248,30 @@ bool Foam::phaseModel::read(const dictionary& phaseProperties)
 }
 
 
+void Foam::phaseModel::correctInflowFlux(surfaceScalarField& alphaPhi) const
+{
+    surfaceScalarField::Boundary& alphaPhiBf = alphaPhi.boundaryFieldRef();
+
+    // Ensure that the flux at inflow BCs is preserved
+    forAll(alphaPhiBf, patchi)
+    {
+        fvsPatchScalarField& alphaPhip = alphaPhiBf[patchi];
+
+        if (!alphaPhip.coupled())
+        {
+            const scalarField& phip = phi().boundaryField()[patchi];
+            const scalarField& alphap = boundaryField()[patchi];
+
+            forAll(alphaPhip, facei)
+            {
+                if (phip[facei] < SMALL)
+                {
+                    alphaPhip[facei] = alphap[facei]*phip[facei];
+                }
+            }
+        }
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.H
index 5fb0245513..a5ae107c90 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.H
@@ -319,6 +319,9 @@ public:
             return alphaRhoPhi_;
         }
 
+        //- Ensure that the flux at inflow BCs is preserved
+        void correctInflowFlux(surfaceScalarField& alphaPhi) const;
+
         //- Correct the phase properties
         //  other than the thermodynamics and turbulence
         //  which have special treatment
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
index 7c2323a030..808fe255f9 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
@@ -444,28 +444,7 @@ void Foam::twoPhaseSystem::solve()
             )
         );
 
-        surfaceScalarField::Boundary& alphaPhic1Bf =
-            alphaPhic1.boundaryFieldRef();
-
-        // Ensure that the flux at inflow BCs is preserved
-        forAll(alphaPhic1Bf, patchi)
-        {
-            fvsPatchScalarField& alphaPhic1p = alphaPhic1Bf[patchi];
-
-            if (!alphaPhic1p.coupled())
-            {
-                const scalarField& phi1p = phi1.boundaryField()[patchi];
-                const scalarField& alpha1p = alpha1.boundaryField()[patchi];
-
-                forAll(alphaPhic1p, facei)
-                {
-                    if (phi1p[facei] < 0)
-                    {
-                        alphaPhic1p[facei] = alpha1p[facei]*phi1p[facei];
-                    }
-                }
-            }
-        }
+        phase1_.correctInflowFlux(alphaPhic1);
 
         if (nAlphaSubCycles > 1)
         {
@@ -537,6 +516,7 @@ void Foam::twoPhaseSystem::solve()
 
         phase2_.alphaPhi() = phi_ - phase1_.alphaPhi();
         alpha2 = scalar(1) - alpha1;
+        phase2_.correctInflowFlux(phase2_.alphaPhi());
         phase2_.alphaRhoPhi() =
             fvc::interpolate(phase2_.rho())*phase2_.alphaPhi();