diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index 2b06f1a7f452f17f9bf6801e1d6deabb66f384ef..2ec0b1590739110aa2424b6ef17e3f70b8fb8482 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
@@ -71,6 +71,8 @@ void Foam::multiphaseSystem::calcAlphas()
 
 void Foam::multiphaseSystem::solveAlphas()
 {
+    bool LTS = fv::localEulerDdt::enabled(mesh_);
+
     PtrList<surfaceScalarField> alphaPhiCorrs(phases().size());
     forAll(phases(), phasei)
     {
@@ -155,14 +157,11 @@ void Foam::multiphaseSystem::solveAlphas()
             }
         }
 
-        if (fv::localEulerDdt::enabled(mesh_))
+        if (LTS)
         {
-            const volScalarField& rDeltaT =
-                fv::localEulerDdt::localRDeltaT(mesh_);
-
             MULES::limit
             (
-                rDeltaT,
+                fv::localEulerDdt::localRDeltaT(mesh_),
                 geometricOneField(),
                 phase,
                 phi_,
@@ -636,8 +635,6 @@ void Foam::multiphaseSystem::solve()
                 fv::localEulerDdt::localRSubDeltaT(mesh_, nAlphaSubCycles);
         }
 
-        dimensionedScalar totalDeltaT = runTime.deltaT();
-
         PtrList<volScalarField> alpha0s(phases().size());
         PtrList<surfaceScalarField> alphaPhiSums(phases().size());
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
index 7b7ddee6f49d622ff4068791a0bc7b3d3ec6b2fa..36e73662cdc05a3043055e09f18c3b96d213a6bd 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
@@ -186,18 +186,17 @@ Foam::twoPhaseSystem::dmdt() const
 
 void Foam::twoPhaseSystem::solve()
 {
-    const fvMesh& mesh = this->mesh();
-    const Time& runTime = mesh.time();
+    const Time& runTime = mesh_.time();
 
     volScalarField& alpha1 = phase1_;
     volScalarField& alpha2 = phase2_;
 
-    const dictionary& alphaControls = mesh.solverDict(alpha1.name());
+    const dictionary& alphaControls = mesh_.solverDict(alpha1.name());
 
     label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
     label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr")));
 
-    bool LTS = fv::localEulerDdt::enabled(mesh);
+    bool LTS = fv::localEulerDdt::enabled(mesh_);
 
     word alphaScheme("div(phi," + alpha1.name() + ')');
     word alpharScheme("div(phir," + alpha1.name() + ')');
@@ -264,9 +263,9 @@ void Foam::twoPhaseSystem::solve()
             (
                 "Sp",
                 runTime.timeName(),
-                mesh
+                mesh_
             ),
-            mesh,
+            mesh_,
             dimensionedScalar("Sp", dimless/dimTime, 0.0)
         );
 
@@ -276,7 +275,7 @@ void Foam::twoPhaseSystem::solve()
             (
                 "Su",
                 runTime.timeName(),
-                mesh
+                mesh_
             ),
             // Divergence term is handled explicitly to be
             // consistent with the explicit transport solution
@@ -345,7 +344,7 @@ void Foam::twoPhaseSystem::solve()
             if (LTS)
             {
                 trSubDeltaT =
-                    fv::localEulerDdt::localRSubDeltaT(mesh, nAlphaSubCycles);
+                    fv::localEulerDdt::localRSubDeltaT(mesh_, nAlphaSubCycles);
             }
 
             for
@@ -420,7 +419,7 @@ void Foam::twoPhaseSystem::solve()
             fvc::interpolate(phase2_.rho())*phase2_.alphaPhi();
 
         Info<< alpha1.name() << " volume fraction = "
-            << alpha1.weightedAverage(mesh.V()).value()
+            << alpha1.weightedAverage(mesh_.V()).value()
             << "  Min(alpha1) = " << min(alpha1).value()
             << "  Max(alpha1) = " << max(alpha1).value()
             << endl;
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
index f3bb5432a52912232cf8ed3de76ad504c99773d2..61e0d885504e8f55a3b4330a289aa05853a82911 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
@@ -89,8 +89,17 @@ void Foam::MULES::explicitSolve
 )
 {
     const fvMesh& mesh = psi.mesh();
-    const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
-    explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su);
+
+    if (fv::localEulerDdt::enabled(mesh))
+    {
+        const volScalarField& rDeltaT = fv::localEulerDdt::localRDeltaT(mesh);
+        explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su);
+    }
+    else
+    {
+        const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
+        explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su);
+    }
 }