diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaDiffusionEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaDiffusionEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..139646e43c6023c1486e3ee53dfdb0a14c63e158
--- /dev/null
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaDiffusionEqn.H
@@ -0,0 +1,18 @@
+{
+    fvScalarMatrix alpha1Eqn
+    (
+        fvm::ddt(alpha1)
+      - fvc::ddt(alpha1)
+      - fvm::laplacian
+        (
+            volScalarField("Dab", Dab + alphatab*turbulence->nut()),
+            alpha1
+        )
+    );
+
+    alpha1Eqn.solve();
+
+    rhoPhi += alpha1Eqn.flux()*(rho1 - rho2);
+}
+
+rho = alpha1*rho1 + (scalar(1) - alpha1)*rho2;
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H
new file mode 100644
index 0000000000000000000000000000000000000000..6c82f940499ea39567db1e400274832dfa354b0f
--- /dev/null
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H
@@ -0,0 +1,26 @@
+label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
+label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
+
+if (nAlphaSubCycles > 1)
+{
+    dimensionedScalar totalDeltaT = runTime.deltaT();
+    surfaceScalarField rhoPhiSum(0.0*rhoPhi);
+
+    for
+    (
+        subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+        !(++alphaSubCycle).end();
+    )
+    {
+        #include "alphaEqn.H"
+        rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
+    }
+
+    rhoPhi = rhoPhiSum;
+}
+else
+{
+    #include "alphaEqn.H"
+}
+
+rho == alpha1*rho1 + (scalar(1) - alpha1)*rho2;