diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
index 9da399d914..eda03ec719 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
@@ -41,6 +41,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "interfaceCompression.H"
 #include "CMULES.H"
 #include "EulerDdtScheme.H"
 #include "localEulerDdtScheme.H"
diff --git a/applications/solvers/multiphase/interFoam/interFoam.C b/applications/solvers/multiphase/interFoam/interFoam.C
index 256e6213c8..3a6e011722 100644
--- a/applications/solvers/multiphase/interFoam/interFoam.C
+++ b/applications/solvers/multiphase/interFoam/interFoam.C
@@ -33,6 +33,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "interfaceCompression.H"
 #include "CMULES.H"
 #include "EulerDdtScheme.H"
 #include "localEulerDdtScheme.H"
diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C b/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
index 2b0a04851a..b6c1641d0b 100644
--- a/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
@@ -32,6 +32,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "interfaceCompression.H"
 #include "CMULES.H"
 #include "localEulerDdtScheme.H"
 #include "subCycle.H"
diff --git a/etc/caseDicts/postProcessing/solvers/scalarTransport b/etc/caseDicts/postProcessing/solvers/scalarTransport
index 7131d7885f..7a8ee73732 100644
--- a/etc/caseDicts/postProcessing/solvers/scalarTransport
+++ b/etc/caseDicts/postProcessing/solvers/scalarTransport
@@ -21,5 +21,9 @@ Description
 field           <fieldName>;    // Name of the transported scalar
 schemesField    $field;         // Name of the field from which to use schemes
                                 // and solvers settings
+diffusion       viscosity;
+
+alphal          1;
+alphat          1;
 
 // ************************************************************************* //
diff --git a/src/functionObjects/solvers/Make/options b/src/functionObjects/solvers/Make/options
index 4e63d7d26c..06e6f834ef 100644
--- a/src/functionObjects/solvers/Make/options
+++ b/src/functionObjects/solvers/Make/options
@@ -4,6 +4,7 @@ EXE_INC = \
     -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
 
@@ -13,6 +14,7 @@ LIB_LIBS = \
     -lmomentumTransportModels \
     -lincompressibleMomentumTransportModels \
     -lcompressibleMomentumTransportModels \
+    -ltwoPhaseMixture \
     -lspecie \
     -lfiniteVolume \
     -lmeshTools \
diff --git a/src/functionObjects/solvers/scalarTransport/scalarTransport.C b/src/functionObjects/solvers/scalarTransport/scalarTransport.C
index 9d19d9dd4f..7823091e93 100644
--- a/src/functionObjects/solvers/scalarTransport/scalarTransport.C
+++ b/src/functionObjects/solvers/scalarTransport/scalarTransport.C
@@ -26,13 +26,23 @@ License
 #include "scalarTransport.H"
 #include "surfaceFields.H"
 #include "fvmDdt.H"
+#include "fvcDdt.H"
 #include "fvmDiv.H"
 #include "fvmLaplacian.H"
 #include "fvmSup.H"
+#include "fvcFlux.H"
 #include "fvModels.H"
 #include "fvConstraints.H"
 #include "incompressibleMomentumTransportModel.H"
 #include "compressibleMomentumTransportModel.H"
+
+#include "CMULES.H"
+#include "EulerDdtScheme.H"
+#include "localEulerDdtScheme.H"
+#include "CrankNicolsonDdtScheme.H"
+#include "subCycle.H"
+#include "interfaceCompression.H"
+
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -53,16 +63,33 @@ namespace functionObjects
 }
 
 
+template<>
+const char* Foam::NamedEnum
+<
+    Foam::functionObjects::scalarTransport::diffusionType,
+    3
+>::names[] =
+{
+    "none",
+    "constant",
+    "viscosity"
+};
+
+const Foam::NamedEnum
+<
+    Foam::functionObjects::scalarTransport::diffusionType,
+    3
+> Foam::functionObjects::scalarTransport::diffusionTypeNames_;
+
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::tmp<Foam::volScalarField> Foam::functionObjects::scalarTransport::D
-(
-    const surfaceScalarField& phi
-) const
+Foam::tmp<Foam::volScalarField>
+Foam::functionObjects::scalarTransport::D() const
 {
     const word Dname("D" + s_.name());
 
-    if (constantD_)
+    if (diffusion_ == diffusionType::constant)
     {
         return volScalarField::New
         (
@@ -71,13 +98,7 @@ Foam::tmp<Foam::volScalarField> Foam::functionObjects::scalarTransport::D
             dimensionedScalar(Dname, dimViscosity, D_)
         );
     }
-    else if
-    (
-        mesh_.foundObject<momentumTransportModel>
-        (
-            momentumTransportModel::typeName
-        )
-    )
+    else
     {
         const momentumTransportModel& turbulence =
             mesh_.lookupObject<momentumTransportModel>
@@ -85,16 +106,7 @@ Foam::tmp<Foam::volScalarField> Foam::functionObjects::scalarTransport::D
                 momentumTransportModel::typeName
             );
 
-        return alphaD_*turbulence.nu() + alphaDt_*turbulence.nut();
-    }
-    else
-    {
-        return volScalarField::New
-        (
-            Dname,
-            mesh_,
-            dimensionedScalar(Dname, phi.dimensions()/dimLength, 0)
-        );
+        return alphal_*turbulence.nu() + alphat_*turbulence.nut();
     }
 }
 
@@ -110,6 +122,7 @@ Foam::functionObjects::scalarTransport::scalarTransport
 :
     fvMeshFunctionObject(name, runTime, dict),
     fieldName_(dict.lookupOrDefault<word>("field", "s")),
+    diffusion_(diffusionType::none),
     D_(0),
     nCorr_(0),
     s_
@@ -123,9 +136,49 @@ Foam::functionObjects::scalarTransport::scalarTransport
             IOobject::AUTO_WRITE
         ),
         mesh_
-    )
+    ),
+    MULES_(false),
+    deltaN_
+    (
+        "deltaN",
+        1e-8/pow(average(mesh_.V()), 1.0/3.0)
+    ),
+    sRestart_(false)
 {
     read(dict);
+
+    const dictionary& controls = mesh_.solverDict(s_.name());
+
+    if (controls.found("nSubCycles"))
+    {
+        MULES_ = true;
+
+        typeIOobject<surfaceScalarField> sPhiHeader
+        (
+            IOobject::groupName("sPhi", s_.group()),
+            runTime.timeName(),
+            mesh_,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        );
+
+        sRestart_ = sPhiHeader.headerOk();
+
+        if (sRestart_)
+        {
+            Info << "Restarting s" << endl;
+        }
+
+        const surfaceScalarField& phi =
+            mesh_.lookupObject<surfaceScalarField>(phiName_);
+
+        // Scalar volumetric flux
+        tsPhi_ = new surfaceScalarField
+        (
+            sPhiHeader,
+            phi*fvc::interpolate(s_)
+        );
+    }
 }
 
 
@@ -145,9 +198,22 @@ bool Foam::functionObjects::scalarTransport::read(const dictionary& dict)
     rhoName_ = dict.lookupOrDefault<word>("rho", "rho");
     schemesField_ = dict.lookupOrDefault<word>("schemesField", fieldName_);
 
-    constantD_ = dict.readIfPresent("D", D_);
-    alphaD_ = dict.lookupOrDefault("alphaD", 1.0);
-    alphaDt_ = dict.lookupOrDefault("alphaDt", 1.0);
+    diffusion_ = diffusionTypeNames_.read(dict.lookup("diffusion"));
+
+    switch(diffusion_)
+    {
+        case diffusionType::none:
+            break;
+
+        case diffusionType::constant:
+            dict.lookup("D") >> D_;
+            break;
+
+        case diffusionType::viscosity:
+            dict.lookup("alphal") >> alphal_;
+            dict.lookup("alphat") >> alphat_;
+            break;
+    }
 
     dict.readIfPresent("nCorr", nCorr_);
 
@@ -168,11 +234,7 @@ bool Foam::functionObjects::scalarTransport::execute()
     const surfaceScalarField& phi =
         mesh_.lookupObject<surfaceScalarField>(phiName_);
 
-    // Calculate the diffusivity
-    volScalarField D("D" + s_.name(), this->D(phi));
-
-    word divScheme("div(phi," + schemesField_ + ")");
-    word laplacianScheme("laplacian(" + D.name() + "," + schemesField_ + ")");
+    const word divScheme("div(phi," + schemesField_ + ")");
 
     // Set under-relaxation coeff
     scalar relaxCoeff = 0.0;
@@ -189,25 +251,38 @@ bool Foam::functionObjects::scalarTransport::execute()
 
     if (phi.dimensions() == dimVolume/dimTime)
     {
-        for (int i=0; i<=nCorr_; i++)
+        if (MULES_)
         {
-            fvScalarMatrix sEqn
-            (
-                fvm::ddt(s_)
-              + fvm::div(phi, s_, divScheme)
-              - fvm::laplacian(D, s_, laplacianScheme)
-             ==
-                fvModels.source(s_)
-            );
-
-            sEqn.relax(relaxCoeff);
-
-            fvConstraints.constrain(sEqn);
-
-            sEqn.solve(schemesField_);
+            subCycleMULES();
 
             fvConstraints.constrain(s_);
         }
+        else
+        {
+            for (int i=0; i<=nCorr_; i++)
+            {
+                fvScalarMatrix sEqn
+                (
+                    fvm::ddt(s_)
+                  + fvm::div(phi, s_, divScheme)
+                 ==
+                    fvModels.source(s_)
+                );
+
+                if (diffusion_ != diffusionType::none)
+                {
+                    sEqn -= fvm::laplacian(D(), s_);
+                }
+
+                sEqn.relax(relaxCoeff);
+
+                fvConstraints.constrain(sEqn);
+
+                sEqn.solve(schemesField_);
+
+                fvConstraints.constrain(s_);
+            }
+        }
     }
     else if (phi.dimensions() == dimMass/dimTime)
     {
@@ -220,11 +295,15 @@ bool Foam::functionObjects::scalarTransport::execute()
             (
                 fvm::ddt(rho, s_)
               + fvm::div(phi, s_, divScheme)
-              - fvm::laplacian(rho*D, s_, laplacianScheme)
              ==
                 fvModels.source(rho, s_)
             );
 
+            if (diffusion_ != diffusionType::none)
+            {
+                sEqn -= fvm::laplacian(rho*D(), s_);
+            }
+
             sEqn.relax(relaxCoeff);
 
             fvConstraints.constrain(sEqn);
@@ -248,6 +327,303 @@ bool Foam::functionObjects::scalarTransport::execute()
 }
 
 
+void Foam::functionObjects::scalarTransport::subCycleMULES()
+{
+    const dictionary& controls = mesh_.solverDict(s_.name());
+    const label nSubCycles(controls.lookup<label>("nSubCycles"));
+    const bool LTS = fv::localEulerDdt::enabled(mesh_);
+
+    if (nSubCycles > 1)
+    {
+        tmp<volScalarField> trSubDeltaT;
+
+        if (LTS)
+        {
+            trSubDeltaT =
+                fv::localEulerDdt::localRSubDeltaT(mesh_, nSubCycles);
+        }
+
+        for
+        (
+            subCycle<volScalarField> sSubCycle(s_, nSubCycles);
+            !(++sSubCycle).end();
+        )
+        {
+            solveMULES();
+        }
+    }
+    else
+    {
+        solveMULES();
+    }
+
+
+    // Apply the diffusion term separately to allow implicit solution
+    // and boundedness of the explicit advection
+    if (diffusion_ != diffusionType::none)
+    {
+        fvScalarMatrix sEqn
+        (
+            fvm::ddt(s_) - fvc::ddt(s_)
+          - fvm::laplacian(D(), s_)
+        );
+
+        sEqn.solve(controls.subDict("diffusion"));
+
+        Info<< s_.name() << " volume fraction = "
+            << s_.weightedAverage(mesh_.V()).value()
+            << "  Min(" << s_.name() << ") = " << min(s_).value()
+            << "  Max(" << s_.name() << ") = " << max(s_).value()
+            << endl;
+    }
+}
+
+
+void Foam::functionObjects::scalarTransport::solveMULES()
+{
+    const dictionary& controls = mesh_.solverDict(s_.name());
+    const label nCorr(controls.lookup<label>("nCorr"));
+    const label nSubCycles(controls.lookup<label>("nSubCycles"));
+    const bool MULESCorr(controls.lookupOrDefault<Switch>("MULESCorr", false));
+
+    // Apply the compression correction from the previous iteration
+    // Improves efficiency for steady-simulations but can only be applied
+    // once the s field is reasonably steady, i.e. fully developed
+    const bool applyPrevCorr
+    (
+        controls.lookupOrDefault<Switch>("applyPrevCorr", false)
+    );
+
+    const bool LTS = fv::localEulerDdt::enabled(mesh_);
+
+    const word divScheme("div(phi," + schemesField_ + ")");
+
+    const surfaceScalarField& phi =
+        mesh_.lookupObject<surfaceScalarField>(phiName_);
+
+    surfaceScalarField& sPhi_ = tsPhi_.ref();
+
+    const word sScheme(mesh_.divScheme(divScheme)[1].wordToken());
+
+    // If a compressive convection scheme is used
+    // the interface normal must be cached
+    tmp<surfaceScalarField> nHatf;
+
+    if (compressionSchemes.found(sScheme))
+    {
+        const surfaceVectorField gradsf(fvc::interpolate(fvc::grad(s_)));
+
+        nHatf = new surfaceScalarField
+        (
+            IOobject
+            (
+                "nHatf",
+                s_.time().timeName(),
+                mesh_
+            ),
+            gradsf/(mag(gradsf) + deltaN_) & mesh_.Sf()
+        );
+    }
+
+    // Set the off-centering coefficient according to ddt scheme
+    scalar ocCoeff = 0;
+    {
+        tmp<fv::ddtScheme<scalar>> tddtS
+        (
+            fv::ddtScheme<scalar>::New
+            (
+                mesh_,
+                mesh_.ddtScheme("ddt(s)")
+            )
+        );
+        const fv::ddtScheme<scalar>& ddtS = tddtS();
+
+        if
+        (
+            isType<fv::EulerDdtScheme<scalar>>(ddtS)
+         || isType<fv::localEulerDdtScheme<scalar>>(ddtS)
+        )
+        {
+            ocCoeff = 0;
+        }
+        else if (isType<fv::CrankNicolsonDdtScheme<scalar>>(ddtS))
+        {
+            if (nSubCycles > 1)
+            {
+                FatalErrorInFunction
+                    << "Sub-cycling is not supported "
+                       "with the CrankNicolson ddt scheme"
+                    << exit(FatalError);
+            }
+
+            if
+            (
+                sRestart_
+             || mesh_.time().timeIndex() > mesh_.time().startTimeIndex() + 1
+            )
+            {
+                ocCoeff =
+                    refCast<const fv::CrankNicolsonDdtScheme<scalar>>(ddtS)
+                   .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);
+
+    tmp<surfaceScalarField> phiCN(phi);
+
+    // Calculate the Crank-Nicolson off-centred volumetric flux
+    if (ocCoeff > 0)
+    {
+        phiCN = surfaceScalarField::New
+        (
+            "phiCN",
+            cnCoeff*phi + (1.0 - cnCoeff)*phi.oldTime()
+        );
+    }
+
+    if (MULESCorr)
+    {
+        fvScalarMatrix sEqn
+        (
+            (
+                LTS
+              ? fv::localEulerDdtScheme<scalar>(mesh_).fvmDdt(s_)
+              : fv::EulerDdtScheme<scalar>(mesh_).fvmDdt(s_)
+            )
+          + fv::gaussConvectionScheme<scalar>
+            (
+                mesh_,
+                phiCN,
+                upwind<scalar>(mesh_, phiCN)
+            ).fvmDiv(phiCN, s_)
+        );
+
+        sEqn.solve();
+
+        Info<< s_.name() << " volume fraction = "
+            << s_.weightedAverage(mesh_.Vsc()).value()
+            << "  Min(" << s_.name() << ") = " << min(s_).value()
+            << "  Max(" << s_.name() << ") = " << max(s_).value()
+            << endl;
+
+        tmp<surfaceScalarField> tsPhiUD(sEqn.flux());
+        sPhi_ = tsPhiUD();
+
+        if (applyPrevCorr && tsPhiCorr0_.valid())
+        {
+            Info<< "Applying the previous iteration compression flux" << endl;
+            MULES::correct
+            (
+                geometricOneField(),
+                s_,
+                sPhi_,
+                tsPhiCorr0_.ref(),
+                oneField(),
+                zeroField()
+            );
+
+            sPhi_ += tsPhiCorr0_();
+        }
+
+        // Cache the upwind-flux
+        tsPhiCorr0_ = tsPhiUD;
+    }
+
+    for (int sCorr=0; sCorr<nCorr; sCorr++)
+    {
+        // Split operator
+        tmp<surfaceScalarField> tsPhiUn
+        (
+            fvc::flux
+            (
+                phiCN(),
+                (cnCoeff*s_ + (1.0 - cnCoeff)*s_.oldTime())(),
+                mesh_.divScheme(divScheme)
+            )
+        );
+
+        if (MULESCorr)
+        {
+            tmp<surfaceScalarField> tsPhiCorr(tsPhiUn() - sPhi_);
+            volScalarField s0("s0", s_);
+
+            MULES::correct
+            (
+                geometricOneField(),
+                s_,
+                tsPhiUn(),
+                tsPhiCorr.ref(),
+                oneField(),
+                zeroField()
+            );
+
+            // Under-relax the correction for all but the 1st corrector
+            if (sCorr == 0)
+            {
+                sPhi_ += tsPhiCorr();
+            }
+            else
+            {
+                s_ = 0.5*s_ + 0.5*s0;
+                sPhi_ += 0.5*tsPhiCorr();
+            }
+        }
+        else
+        {
+            sPhi_ = tsPhiUn;
+
+            MULES::explicitSolve
+            (
+                geometricOneField(),
+                s_,
+                phiCN,
+                sPhi_,
+                oneField(),
+                zeroField()
+            );
+        }
+    }
+
+    if (applyPrevCorr && MULESCorr)
+    {
+        tsPhiCorr0_ = sPhi_ - tsPhiCorr0_;
+        tsPhiCorr0_.ref().rename("sPhiCorr0");
+    }
+    else
+    {
+        tsPhiCorr0_.clear();
+    }
+
+    if
+    (
+        word(mesh_.ddtScheme("ddt(s)"))
+     == fv::CrankNicolsonDdtScheme<scalar>::typeName
+    )
+    {
+        if (ocCoeff > 0)
+        {
+            // Calculate the end-of-time-step s flux
+            sPhi_ = (sPhi_ - (1.0 - cnCoeff)*sPhi_.oldTime())/cnCoeff;
+        }
+    }
+
+    Info<< s_.name() << "volume fraction = "
+        << s_.weightedAverage(mesh_.Vsc()).value()
+        << "  Min(" << s_.name() << ") = " << min(s_).value()
+        << "  Max(" << s_.name() << ") = " << max(s_).value()
+        << endl;
+}
+
+
 bool Foam::functionObjects::scalarTransport::write()
 {
     return true;
diff --git a/src/functionObjects/solvers/scalarTransport/scalarTransport.H b/src/functionObjects/solvers/scalarTransport/scalarTransport.H
index 6c892e76a5..5fa03666e9 100644
--- a/src/functionObjects/solvers/scalarTransport/scalarTransport.H
+++ b/src/functionObjects/solvers/scalarTransport/scalarTransport.H
@@ -30,16 +30,58 @@ Description
     - To specify the field name set the \c field entry
     - To employ the same numerical schemes as another field set
       the \c schemesField entry,
-    - A constant diffusivity may be specified with the \c D entry,
-
-    - Alternatively if a turbulence model is available a turbulent diffusivity
-      may be constructed from the laminar and turbulent viscosities using the
-      optional diffusivity coefficients \c alphaD and \c alphaDt (which default
-      to 1):
+    - The \c diffusivity entry can be set to \c none, \c constant, \c viscosity
+    - A constant diffusivity is specified with the \c D entry,
+    - If a momentum transport model is available and the \c viscosity
+      diffusivety option specified an effective diffusivity may be constructed
+      from the laminar and turbulent viscosities using the diffusivity
+      coefficients \c alphal and \c alphat:
       \verbatim
-          D = alphaD*nu + alphaDt*nut
+          D = alphal*nu + alphat*nut
       \endverbatim
 
+    Example:
+    \verbatim
+        #includeFunc scalarTransport(T, alphaD=1, alphaDt=1)
+    \endverbatim
+
+    For incompressible flow the passive scalar may optionally be solved with the
+    MULES limiter and sub-cycling or semi-implicit in order to maintain
+    boundedness, particularly if a compressive, PLIC or MPLIC convection
+    scheme is used.
+
+    Example:
+    \verbatim
+        #includeFunc scalarTransport(tracer, diffusion=none)
+
+    with scheme specification:
+        div(phi,tracer)     Gauss interfaceCompression vanLeer 1;
+
+    and solver specification:
+        tracer
+        {
+            nCorr      1;
+            nSubCycles 3;
+
+            MULESCorr       no;
+            nLimiterIter    5;
+            applyPrevCorr   yes;
+
+            solver          smoothSolver;
+            smoother        symGaussSeidel;
+            tolerance       1e-8;
+            relTol          0;
+
+            diffusion
+            {
+                solver          smoothSolver;
+                smoother        symGaussSeidel;
+                tolerance       1e-8;
+                relTol          0;
+            }
+        }
+    \endverbatim
+
 See also
     Foam::functionObjects::fvMeshFunctionObject
 
@@ -69,6 +111,19 @@ class scalarTransport
 :
     public fvMeshFunctionObject
 {
+public:
+
+        //- Enumeration defining the type of the diffusion
+        enum class diffusionType
+        {
+            none,
+            constant,
+            viscosity
+        };
+
+
+private:
+
     // Private Data
 
         //- Name of field to process
@@ -80,17 +135,20 @@ class scalarTransport
         //- Name of density field (optional)
         word rhoName_;
 
-        //- Diffusion coefficient (optional)
+        //- Diffusion type names
+        static const NamedEnum<diffusionType, 3> diffusionTypeNames_;
+
+        //- The type of diffusion
+        diffusionType diffusion_;
+
+        //- Constant diffusion coefficient (optional)
         scalar D_;
 
-        //- Flag to indicate whether a constant, uniform D_ is specified
-        bool constantD_;
-
         //- Laminar diffusion coefficient (optional)
-        scalar alphaD_;
+        scalar alphal_;
 
         //- Turbulent diffusion coefficient (optional)
-        scalar alphaDt_;
+        scalar alphat_;
 
         //- Number of corrector iterations (optional)
         int nCorr_;
@@ -101,11 +159,30 @@ class scalarTransport
         //- The scalar field
         volScalarField s_;
 
+        //- Switch for MULES limited solution
+        bool MULES_;
+
+        //- Stabilisation for normalisation of the interface normal
+        //  needed if a compressive convection scheme is used
+        const dimensionedScalar deltaN_;
+
+        //- Scalar volumetric flux
+        tmp<surfaceScalarField> tsPhi_;
+
+        //- MULES Correction
+        tmp<surfaceScalarField> tsPhiCorr0_;
+
+        //- Switch to indicate that s has been restarted for Crank-Nicolson
+        bool sRestart_;
+
 
     // Private Member Functions
 
         //- Return the diffusivity field
-        tmp<volScalarField> D(const surfaceScalarField& phi) const;
+        tmp<volScalarField> D() const;
+
+        void subCycleMULES();
+        void solveMULES();
 
 
 public:
diff --git a/src/twoPhaseModels/twoPhaseMixture/VoF/alphaScheme.H b/src/twoPhaseModels/twoPhaseMixture/VoF/alphaScheme.H
index 76b2b37def..075751866b 100644
--- a/src/twoPhaseModels/twoPhaseMixture/VoF/alphaScheme.H
+++ b/src/twoPhaseModels/twoPhaseMixture/VoF/alphaScheme.H
@@ -1,13 +1,3 @@
-static const wordHashSet compressionSchemes
-{
-    "interfaceCompression",
-    "noInterfaceCompression",
-    "PLIC",
-    "PLICU",
-    "MPLIC",
-    "MPLICU"
-};
-
 static const word divAlphaName("div(phi,alpha)");
 
 const word alphaScheme(mesh.divScheme(divAlphaName)[1].wordToken());
diff --git a/src/twoPhaseModels/twoPhaseMixture/interfaceCompression/interfaceCompression.C b/src/twoPhaseModels/twoPhaseMixture/interfaceCompression/interfaceCompression.C
index e1248426a0..7956f00994 100644
--- a/src/twoPhaseModels/twoPhaseMixture/interfaceCompression/interfaceCompression.C
+++ b/src/twoPhaseModels/twoPhaseMixture/interfaceCompression/interfaceCompression.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2020-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,6 +35,16 @@ namespace Foam
     surfaceInterpolationScheme<scalar>::
         addMeshFluxConstructorToTable<interfaceCompressionNew>
         addinterfaceCompressionScalarMeshFluxConstructorToTable_;
+
+    const wordHashSet compressionSchemes
+    {
+        "interfaceCompression",
+        "noInterfaceCompression",
+        "PLIC",
+        "PLICU",
+        "MPLIC",
+        "MPLICU"
+    };
 }
 
 
diff --git a/src/twoPhaseModels/twoPhaseMixture/interfaceCompression/interfaceCompression.H b/src/twoPhaseModels/twoPhaseMixture/interfaceCompression/interfaceCompression.H
index ab567c1bd4..671733878b 100644
--- a/src/twoPhaseModels/twoPhaseMixture/interfaceCompression/interfaceCompression.H
+++ b/src/twoPhaseModels/twoPhaseMixture/interfaceCompression/interfaceCompression.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2020-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -143,6 +143,9 @@ public:
 };
 
 
+extern const wordHashSet compressionSchemes;
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/tutorials/incompressible/pimpleFoam/RAS/flowWithOpenBoundary/system/controlDict b/tutorials/incompressible/pimpleFoam/RAS/flowWithOpenBoundary/system/controlDict
index 0d878bd7ae..3d55a89d0d 100644
--- a/tutorials/incompressible/pimpleFoam/RAS/flowWithOpenBoundary/system/controlDict
+++ b/tutorials/incompressible/pimpleFoam/RAS/flowWithOpenBoundary/system/controlDict
@@ -45,7 +45,7 @@ runTimeModifiable true;
 
 functions
 {
-    #includeFunc scalarTransport(T, alphaD=1, alphaDt=1)
+    #includeFunc scalarTransport(T, alphal=1, alphat=1)
     #includeFunc time
 }