diff --git a/applications/solvers/multiphase/VoF/alphaEqn.H b/applications/solvers/multiphase/VoF/alphaEqn.H
index 570442a30..abae00fb3 100644
--- a/applications/solvers/multiphase/VoF/alphaEqn.H
+++ b/applications/solvers/multiphase/VoF/alphaEqn.H
@@ -2,46 +2,57 @@
     word alphaScheme("div(phi,alpha)");
     word alpharScheme("div(phirb,alpha)");
 
-    tmp<fv::ddtScheme<scalar>> ddtAlpha
-    (
-        fv::ddtScheme<scalar>::New
-        (
-            mesh,
-            mesh.ddtScheme("ddt(alpha)")
-        )
-    );
-
     // Set the off-centering coefficient according to ddt scheme
     scalar ocCoeff = 0;
-    if
-    (
-        isType<fv::EulerDdtScheme<scalar>>(ddtAlpha())
-     || isType<fv::localEulerDdtScheme<scalar>>(ddtAlpha())
-    )
     {
-        ocCoeff = 0;
-    }
-    else if (isType<fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha()))
-    {
-        if (nAlphaSubCycles > 1)
+        tmp<fv::ddtScheme<scalar>> tddtAlpha
+        (
+            fv::ddtScheme<scalar>::New
+            (
+                mesh,
+                mesh.ddtScheme("ddt(alpha)")
+            )
+        );
+        const fv::ddtScheme<scalar>& ddtAlpha = tddtAlpha();
+
+        if
+        (
+            isType<fv::EulerDdtScheme<scalar>>(ddtAlpha)
+         || isType<fv::localEulerDdtScheme<scalar>>(ddtAlpha)
+        )
+        {
+            ocCoeff = 0;
+        }
+        else if (isType<fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha))
+        {
+            if (nAlphaSubCycles > 1)
+            {
+                FatalErrorInFunction
+                    << "Sub-cycling is not supported "
+                       "with the CrankNicolson ddt scheme"
+                    << exit(FatalError);
+            }
+
+            if
+            (
+                alphaRestart
+             || mesh.time().timeIndex() > mesh.time().startTimeIndex() + 1
+            )
+            {
+                ocCoeff =
+                    refCast<const fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha)
+                   .ocCoeff();
+            }
+        }
+        else
         {
             FatalErrorInFunction
-                << "Sub-cycling is not supported "
-                   "with the CrankNicolson ddt scheme"
+                << "Only Euler and CrankNicolson ddt schemes are supported"
                 << exit(FatalError);
         }
-
-        ocCoeff =
-            refCast<const fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha())
-           .ocCoeff();
-    }
-    else
-    {
-        FatalErrorInFunction
-            << "Only Euler and CrankNicolson ddt schemes are supported"
-            << exit(FatalError);
     }
 
+    // Set the time blending factor, 1 for Euler
     scalar cnCoeff = 1.0/(1.0 + ocCoeff);
 
     // Standard face-flux compression coefficient
@@ -136,8 +147,8 @@
         (
             fvc::flux
             (
-                phi,
-                alpha1,
+                phiCN(),
+                cnCoeff*alpha1 + (1.0 - cnCoeff)*alpha1.oldTime(),
                 alphaScheme
             )
           + fvc::flux
@@ -148,13 +159,6 @@
             )
         );
 
-        // Calculate the Crank-Nicolson off-centred alpha flux
-        if (ocCoeff > 0)
-        {
-            talphaPhiUn =
-                cnCoeff*talphaPhiUn + (1.0 - cnCoeff)*alphaPhi.oldTime();
-        }
-
         if (MULESCorr)
         {
             tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);
diff --git a/applications/solvers/multiphase/VoF/createAlphaFluxes.H b/applications/solvers/multiphase/VoF/createAlphaFluxes.H
new file mode 100644
index 000000000..e75c2fe0b
--- /dev/null
+++ b/applications/solvers/multiphase/VoF/createAlphaFluxes.H
@@ -0,0 +1,20 @@
+IOobject alphaPhiHeader
+(
+    "alphaPhi",
+    runTime.timeName(),
+    mesh,
+    IOobject::READ_IF_PRESENT,
+    IOobject::AUTO_WRITE
+);
+
+const bool alphaRestart = alphaPhiHeader.headerOk();
+
+// MULES flux from previous time-step
+surfaceScalarField alphaPhi
+(
+    alphaPhiHeader,
+    phi*fvc::interpolate(alpha1)
+);
+
+// MULES Correction
+tmp<surfaceScalarField> talphaPhiCorr0;
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
index f82665a49..896af439a 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
@@ -64,6 +64,7 @@ int main(int argc, char *argv[])
     #include "initContinuityErrs.H"
     #include "createControl.H"
     #include "createFields.H"
+    #include "createAlphaFluxes.H"
     #include "createFvOptions.H"
     #include "createUf.H"
     #include "createControls.H"
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
index f6ec5d20c..ce6446d0e 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
@@ -61,6 +61,7 @@ int main(int argc, char *argv[])
     #include "createControl.H"
     #include "createTimeControls.H"
     #include "createFields.H"
+    #include "createAlphaFluxes.H"
     #include "createFvOptions.H"
 
     volScalarField& p = mixture.p();
diff --git a/applications/solvers/multiphase/compressibleInterFoam/createFields.H b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
index 9bac4d825..46168e6b8 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
@@ -101,21 +101,4 @@ autoPtr<compressible::turbulenceModel> turbulence
 Info<< "Creating field kinetic energy K\n" << endl;
 volScalarField K("K", 0.5*magSqr(U));
 
-// MULES flux from previous time-step
-surfaceScalarField alphaPhi
-(
-    IOobject
-    (
-        "alphaPhi",
-        runTime.timeName(),
-        mesh,
-        IOobject::READ_IF_PRESENT,
-        IOobject::AUTO_WRITE
-    ),
-    phi*fvc::interpolate(alpha1)
-);
-
-// MULES Correction
-tmp<surfaceScalarField> talphaPhiCorr0;
-
 #include "createMRF.H"
diff --git a/applications/solvers/multiphase/interFoam/createFields.H b/applications/solvers/multiphase/interFoam/createFields.H
index 4a82afbd2..dca977548 100644
--- a/applications/solvers/multiphase/interFoam/createFields.H
+++ b/applications/solvers/multiphase/interFoam/createFields.H
@@ -121,21 +121,4 @@ if (p_rgh.needReference())
 mesh.setFluxRequired(p_rgh.name());
 mesh.setFluxRequired(alpha1.name());
 
-// MULES flux from previous time-step
-surfaceScalarField alphaPhi
-(
-    IOobject
-    (
-        "alphaPhi",
-        runTime.timeName(),
-        mesh,
-        IOobject::READ_IF_PRESENT,
-        IOobject::AUTO_WRITE
-    ),
-    phi*fvc::interpolate(alpha1)
-);
-
-// MULES Correction
-tmp<surfaceScalarField> talphaPhiCorr0;
-
 #include "createMRF.H"
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
index df826479f..6a6da9d2d 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
@@ -60,6 +60,7 @@ int main(int argc, char *argv[])
     #include "createTimeControls.H"
     #include "createDyMControls.H"
     #include "createFields.H"
+    #include "createAlphaFluxes.H"
     #include "createFvOptions.H"
 
     volScalarField rAU
diff --git a/applications/solvers/multiphase/interFoam/interFoam.C b/applications/solvers/multiphase/interFoam/interFoam.C
index 93d70dea5..7b755b6d1 100644
--- a/applications/solvers/multiphase/interFoam/interFoam.C
+++ b/applications/solvers/multiphase/interFoam/interFoam.C
@@ -63,6 +63,7 @@ int main(int argc, char *argv[])
     #include "createTimeControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
+    #include "createAlphaFluxes.H"
     #include "createFvOptions.H"
     #include "correctPhi.H"
 
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
index 956898e4b..0c7fe053b 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
@@ -198,7 +198,7 @@ scalar CrankNicolsonDdtScheme<Type>::coef_
     const DDt0Field<GeoField>& ddt0
 ) const
 {
-    if (mesh().time().timeIndex() - ddt0.startTimeIndex() > 0)
+    if (mesh().time().timeIndex() > ddt0.startTimeIndex())
     {
         return 1 + ocCoeff_;
     }
@@ -216,7 +216,7 @@ scalar CrankNicolsonDdtScheme<Type>::coef0_
     const DDt0Field<GeoField>& ddt0
 ) const
 {
-    if (mesh().time().timeIndex() - ddt0.startTimeIndex() > 1)
+    if (mesh().time().timeIndex() > ddt0.startTimeIndex() + 1)
     {
         return 1 + ocCoeff_;
     }