diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H
index c27a3e752e..d152baba0b 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H
@@ -3,7 +3,11 @@
     (
         fvm::ddt(alpha1)
       + fvm::div(phi, alpha1)
-      - fvm::laplacian(turbulence->nut() + Dab, alpha1)
+      - fvm::laplacian
+        (
+            Dab + alphatab*turbulence->nut(), alpha1,
+            "laplacian(Dab,alpha1)"
+        )
     );
 
     alpha1Eqn.solve();
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
index a38135dc19..3b5e064372 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
@@ -50,6 +50,9 @@
 
     dimensionedScalar Dab(twoPhaseProperties.lookup("Dab"));
 
+    // Read the reciprocal of the turbulent Schmidt number
+    dimensionedScalar alphatab(twoPhaseProperties.lookup("alphatab"));
+
     // Need to store rho for ddt(rho, U)
     volScalarField rho("rho", alpha1*rho1 + (scalar(1) - alpha1)*rho2);
     rho.oldTime();