diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index 2203d6e144..767600e394 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
@@ -888,7 +888,7 @@ void Foam::multiphaseSystem::solve()
         dimensionedScalar totalDeltaT = runTime.deltaT();
 
         PtrList<volScalarField> alpha0s(phases_.size());
-        PtrList<surfaceScalarField> phiSums(phases_.size());
+        PtrList<surfaceScalarField> alphaPhiSums(phases_.size());
 
         int phasei = 0;
         forAllIter(PtrDictionary<phaseModel>, phases_, iter)
@@ -902,7 +902,7 @@ void Foam::multiphaseSystem::solve()
                 new volScalarField(alpha.oldTime())
             );
 
-            phiSums.set
+            alphaPhiSums.set
             (
                 phasei,
                 new surfaceScalarField
@@ -936,7 +936,7 @@ void Foam::multiphaseSystem::solve()
             int phasei = 0;
             forAllIter(PtrDictionary<phaseModel>, phases_, iter)
             {
-                phiSums[phasei] += (runTime.deltaT()/totalDeltaT)*iter().phi();
+                alphaPhiSums[phasei] += iter().alphaPhi()/nAlphaSubCycles;
                 phasei++;
             }
         }
@@ -947,7 +947,7 @@ void Foam::multiphaseSystem::solve()
             phaseModel& phase = iter();
             volScalarField& alpha = phase;
 
-            phase.phi() = phiSums[phasei];
+            phase.alphaPhi() = alphaPhiSums[phasei];
 
             // Correct the time index of the field
             // to correspond to the global time
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index 5512fd6559..d92368a886 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
@@ -157,7 +157,7 @@ void Foam::multiphaseSystem::solveAlphas()
 
         MULES::limit
         (
-            1.0/mesh_.time().deltaT().value(),
+            1.0/mesh_.time().deltaT().value(), // ***HGW add support for LTS
             geometricOneField(),
             phase,
             phi_,
@@ -620,7 +620,7 @@ void Foam::multiphaseSystem::solve()
         dimensionedScalar totalDeltaT = runTime.deltaT();
 
         PtrList<volScalarField> alpha0s(phases().size());
-        PtrList<surfaceScalarField> phiSums(phases().size());
+        PtrList<surfaceScalarField> alphaPhiSums(phases().size());
 
         forAll(phases(), phasei)
         {
@@ -633,7 +633,7 @@ void Foam::multiphaseSystem::solve()
                 new volScalarField(alpha.oldTime())
             );
 
-            phiSums.set
+            alphaPhiSums.set
             (
                 phasei,
                 new surfaceScalarField
@@ -664,7 +664,7 @@ void Foam::multiphaseSystem::solve()
 
             forAll(phases(), phasei)
             {
-                phiSums[phasei] += phases()[phasei].phi();
+                alphaPhiSums[phasei] += phases()[phasei].alphaPhi();
             }
         }
 
@@ -673,7 +673,7 @@ void Foam::multiphaseSystem::solve()
             phaseModel& phase = phases()[phasei];
             volScalarField& alpha = phase;
 
-            phase.phi() = phiSums[phasei]/nAlphaSubCycles;
+            phase.alphaPhi() = alphaPhiSums[phasei]/nAlphaSubCycles;
 
             // Correct the time index of the field
             // to correspond to the global time