diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/Make/options b/applications/solvers/multiphase/twoLiquidMixingFoam/Make/options
index c8ce69c074b6f87f924a391ae8778b3f19e946ee..6e20519ec98d90ae8f38ba26587a05c72e39374c 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/Make/options
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/Make/options
@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I../interFoam \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H
index d152baba0bb44a1f2cd042596f94030157e0f713..c68a0e56c5c794542ad1925514cef31ce555d4ec 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H
@@ -1,22 +1,22 @@
 {
-    fvScalarMatrix alpha1Eqn
+    word alphaScheme("div(phi,alpha1)");
+
+    surfaceScalarField phiAlpha
     (
-        fvm::ddt(alpha1)
-      + fvm::div(phi, alpha1)
-      - fvm::laplacian
+        fvc::flux
         (
-            Dab + alphatab*turbulence->nut(), alpha1,
-            "laplacian(Dab,alpha1)"
+            phi,
+            alpha1,
+            alphaScheme
         )
     );
 
-    alpha1Eqn.solve();
+    MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
 
-    rhoPhi = alpha1Eqn.flux()*(rho1 - rho2) + phi*rho2;
-    rho = alpha1*rho1 + (scalar(1) - alpha1)*rho2;
+    rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
 
-    Info<< "Phase 1 volume fraction = "
-        << alpha1.weightedAverage(mesh.V()).value()
+    Info<< "Phase-1 volume fraction = "
+        << alpha1.weightedAverage(mesh.Vsc()).value()
         << "  Min(alpha1) = " << min(alpha1).value()
         << "  Max(alpha1) = " << max(alpha1).value()
         << endl;
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C b/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
index 8a2b20e153be50b07fd86ae745f5b4e7776503f4..22c709c960afbe9c910bac8523cc81a800f3f833 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
@@ -32,6 +32,8 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "MULES.H"
+#include "subCycle.H"
 #include "twoPhaseMixture.H"
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
@@ -60,19 +62,21 @@ int main(int argc, char *argv[])
     {
         #include "readTimeControls.H"
         #include "CourantNo.H"
+        #include "alphaCourantNo.H"
         #include "setDeltaT.H"
 
         runTime++;
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
+        twoPhaseProperties.correct();
+
+        #include "alphaEqnSubCycle.H"
+        #include "alphaDiffusionEqn.H"
+
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
-            twoPhaseProperties.correct();
-
-            #include "alphaEqn.H"
-
             #include "UEqn.H"
 
             // --- Pressure corrector loop
diff --git a/tutorials/multiphase/twoLiquidMixingFoam/lockExchange/system/controlDict b/tutorials/multiphase/twoLiquidMixingFoam/lockExchange/system/controlDict
index 13e3d95facdd911b10dca02fd2236d3355211e63..d8e11662d8eab7c6d1a2a7254f51d78ad931b714 100644
--- a/tutorials/multiphase/twoLiquidMixingFoam/lockExchange/system/controlDict
+++ b/tutorials/multiphase/twoLiquidMixingFoam/lockExchange/system/controlDict
@@ -48,6 +48,7 @@ runTimeModifiable yes;
 adjustTimeStep  on;
 
 maxCo           0.5;
+maxAlphaCo      0.5;
 maxDeltaT       1;
 
 // ************************************************************************* //
diff --git a/tutorials/multiphase/twoLiquidMixingFoam/lockExchange/system/fvSolution b/tutorials/multiphase/twoLiquidMixingFoam/lockExchange/system/fvSolution
index a4e4e913ec6feecc616f1d46f18592c27713a6e5..19d27d99becf5739284643fa38934c9dc77f7cc8 100644
--- a/tutorials/multiphase/twoLiquidMixingFoam/lockExchange/system/fvSolution
+++ b/tutorials/multiphase/twoLiquidMixingFoam/lockExchange/system/fvSolution
@@ -21,7 +21,7 @@ solvers
     {
         solver          smoothSolver;
         smoother        GaussSeidel;
-        tolerance       1e-7;
+        tolerance       1e-9;
         relTol          0;
         nSweeps         1;
     }
@@ -66,6 +66,7 @@ PIMPLE
     nOuterCorrectors    1;
     nCorrectors         2;
     nNonOrthogonalCorrectors 0;
+    nAlphaSubCycles     2;
     pRefValue           0;
     pRefPoint           (0.1 0.1 1);
 }