diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/EEqns.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/EEqns.H
index 9d91ada887..d9f1e85e6f 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/EEqns.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/EEqns.H
@@ -29,37 +29,25 @@ for (int Ecorr=0; Ecorr<nEnergyCorrectors; Ecorr++)
     }
 
     {
-        tmp<fvScalarMatrix> E2eqn(phase2.heEqn());
+        tmp<fvScalarMatrix> E2Eqn(phase2.heEqn());
 
-        if (E2eqn.valid())
+        if (E2Eqn.valid())
         {
-            E2eqn =
+            E2Eqn =
             (
-                E2eqn
+                E2Eqn
              ==
                *heatTransfer[phase2.name()]
               + alpha2*rho2*(U2&g)
               + fvOptions(alpha2, rho2, phase2.thermo().he())
             );
 
-            E2eqn->relax();
-            fvOptions.constrain(E2eqn.ref());
-            E2eqn->solve();
+            E2Eqn->relax();
+            fvOptions.constrain(E2Eqn.ref());
+            E2Eqn->solve();
         }
     }
 
     fluid.correctThermo();
     fluid.correct();
 }
-
-Info<< phase1.name() << " min/max T "
-    << min(phase1.thermo().T()).value()
-    << " - "
-    << max(phase1.thermo().T()).value()
-    << endl;
-
-Info<< phase2.name() << " min/max T "
-    << min(phase2.thermo().T()).value()
-    << " - "
-    << max(phase2.thermo().T()).value()
-    << endl;
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/createRDeltaTf.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/createRDeltaTf.H
new file mode 100644
index 0000000000..5a4c77eacf
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/createRDeltaTf.H
@@ -0,0 +1,21 @@
+tmp<surfaceScalarField> trDeltaTf;
+
+if (LTS && faceMomentum)
+{
+    trDeltaTf = tmp<surfaceScalarField>
+    (
+        new surfaceScalarField
+        (
+            IOobject
+            (
+                fv::localEulerDdt::rDeltaTfName,
+                runTime.timeName(),
+                mesh,
+                IOobject::READ_IF_PRESENT,
+                IOobject::AUTO_WRITE
+            ),
+            mesh,
+            dimensionedScalar("one", dimless/dimTime, 1)
+        )
+    );
+}
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/DDtU.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/DDtU.H
index 63b61cd61a..c7bccf66ea 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/DDtU.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/DDtU.H
@@ -1,9 +1,2 @@
-    ddtPhi1 =
-    (
-        (phi1 - phi1.oldTime())/runTime.deltaT()
-    );
-
-    ddtPhi2 =
-    (
-        (phi2 - phi2.oldTime())/runTime.deltaT()
-    );
+ddtPhi1 = fvc::ddt(phi1);
+ddtPhi2 = fvc::ddt(phi2);
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/createDDtU.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/createDDtU.H
index 7445d4af23..14e66e47a8 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/createDDtU.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/createDDtU.H
@@ -1,9 +1,2 @@
-    surfaceScalarField ddtPhi1
-    (
-        (phi1 - phi1.oldTime())/runTime.deltaT()
-    );
-
-    surfaceScalarField ddtPhi2
-    (
-        (phi2 - phi2.oldTime())/runTime.deltaT()
-    );
+surfaceScalarField ddtPhi1(fvc::ddt(phi1));
+surfaceScalarField ddtPhi2(fvc::ddt(phi2));
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
index 67374f26ee..b76932216f 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
@@ -27,15 +27,17 @@ surfaceScalarField Kdf("Kdf", fluid.Kdf());
 // Virtual-mass coefficient
 surfaceScalarField Vmf("Vmf", fluid.Vmf());
 
+const surfaceScalarField& rDeltaTf = fv::localEulerDdt::localRDeltaTf(mesh);
+
 surfaceScalarField rAUf1
 (
     IOobject::groupName("rAUf", phase1.name()),
     1.0
    /(
-        (alphaRhof10 + Vmf)/runTime.deltaT()
+        rDeltaTf*(alphaRhof10 + Vmf)
       + fvc::interpolate(U1Eqn.A())
       + Kdf
-   )
+    )
 );
 
 surfaceScalarField rAUf2
@@ -43,10 +45,10 @@ surfaceScalarField rAUf2
     IOobject::groupName("rAUf", phase2.name()),
     1.0
    /(
-        (alphaRhof20 + Vmf)/runTime.deltaT()
+        rDeltaTf*(alphaRhof20 + Vmf)
       + fvc::interpolate(U2Eqn.A())
       + Kdf
-   )
+    )
 );
 
 
@@ -157,7 +159,7 @@ while (pimple.correct())
         rAUf1
        *(
             (alphaRhof10 + Vmf)
-           *MRF.absolute(phi1.oldTime())/runTime.deltaT()
+           *rDeltaTf*MRF.absolute(phi1.oldTime())
           + fvc::flux(U1Eqn.H())
           + Vmf*ddtPhi2
           + Kdf*MRF.absolute(phi2)
@@ -175,7 +177,7 @@ while (pimple.correct())
         rAUf2
        *(
             (alphaRhof20 + Vmf)
-           *MRF.absolute(phi2.oldTime())/runTime.deltaT()
+           *rDeltaTf*MRF.absolute(phi2.oldTime())
           + fvc::flux(U2Eqn.H())
           + Vmf*ddtPhi1
           + Kdf*MRF.absolute(phi1)
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C
index a2ca33d51f..6e9d4e4be4 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C
@@ -73,6 +73,7 @@ int main(int argc, char *argv[])
         )
     );
 
+    #include "createRDeltaTf.H"
     #include "pUf/createDDtU.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -91,6 +92,10 @@ int main(int argc, char *argv[])
         if (LTS)
         {
             #include "setRDeltaT.H"
+            if (faceMomentum)
+            {
+                #include "setRDeltaTf.H"
+            }
         }
         else
         {
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/setRDeltaTf.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/setRDeltaTf.H
new file mode 100644
index 0000000000..1c0366dc3b
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/setRDeltaTf.H
@@ -0,0 +1 @@
+trDeltaTf.ref() = fvc::interpolate(fv::localEulerDdt::localRDeltaT(mesh));
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
index 3ac77ae34c..2fd8acaea7 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
@@ -317,6 +317,33 @@ EulerDdtScheme<Type>::fvcDdt
 }
 
 
+template<class Type>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+EulerDdtScheme<Type>::fvcDdt
+(
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& sf
+)
+{
+    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+
+    IOobject ddtIOobject
+    (
+        "ddt("+sf.name()+')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    return tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+    (
+        new GeometricField<Type, fvsPatchField, surfaceMesh>
+        (
+            ddtIOobject,
+            rDeltaT*(sf - sf.oldTime())
+        )
+    );
+}
+
+
 template<class Type>
 tmp<fvMatrix<Type>>
 EulerDdtScheme<Type>::fvmDdt
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.H
index 4cc92a40e2..b8ce3488dd 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.H
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.H
@@ -124,6 +124,11 @@ public:
             const GeometricField<Type, fvPatchField, volMesh>& psi
         );
 
+        tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> fvcDdt
+        (
+            const GeometricField<Type, fvsPatchField, surfaceMesh>&
+        );
+
         tmp<fvMatrix<Type>> fvmDdt
         (
             const GeometricField<Type, fvPatchField, volMesh>&
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C
index 41024f831c..98d0ec1057 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C
@@ -132,6 +132,22 @@ tmp<fvMatrix<Type>> ddtScheme<Type>::fvmDdt
 }
 
 
+template<class Type>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> ddtScheme<Type>::fvcDdt
+(
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& sf
+)
+{
+    NotImplemented;
+
+    return tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+    (
+        GeometricField<Type, fvsPatchField, surfaceMesh>::null()
+    );
+}
+
+
+
 template<class Type>
 tmp<surfaceScalarField> ddtScheme<Type>::fvcDdtPhiCoeff
 (
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.H
index f475ed24f4..9a5861426f 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.H
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.H
@@ -167,6 +167,11 @@ public:
             const GeometricField<Type, fvPatchField, volMesh>&
         ) = 0;
 
+        virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> fvcDdt
+        (
+            const GeometricField<Type, fvsPatchField, surfaceMesh>&
+        );
+
         virtual tmp<fvMatrix<Type>> fvmDdt
         (
             const GeometricField<Type, fvPatchField, volMesh>&
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdt.C b/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdt.C
index 0d607c86b4..3d09e85917 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdt.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdt.C
@@ -29,6 +29,7 @@ License
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 Foam::word Foam::fv::localEulerDdt::rDeltaTName("rDeltaT");
+Foam::word Foam::fv::localEulerDdt::rDeltaTfName("rDeltaTf");
 Foam::word Foam::fv::localEulerDdt::rSubDeltaTName("rSubDeltaT");
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -53,6 +54,18 @@ const Foam::volScalarField& Foam::fv::localEulerDdt::localRDeltaT
 }
 
 
+const Foam::surfaceScalarField& Foam::fv::localEulerDdt::localRDeltaTf
+(
+    const fvMesh& mesh
+)
+{
+    return mesh.objectRegistry::lookupObject<surfaceScalarField>
+    (
+        rDeltaTfName
+    );
+}
+
+
 Foam::tmp<Foam::volScalarField> Foam::fv::localEulerDdt::localRSubDeltaT
 (
     const fvMesh& mesh,
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C
index 2dcdb83bb6..cad1b1aa90 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C
@@ -46,6 +46,13 @@ const volScalarField& localEulerDdtScheme<Type>::localRDeltaT() const
 }
 
 
+template<class Type>
+const surfaceScalarField& localEulerDdtScheme<Type>::localRDeltaTf() const
+{
+    return localEulerDdt::localRDeltaTf(mesh());
+}
+
+
 template<class Type>
 tmp<GeometricField<Type, fvPatchField, volMesh>>
 localEulerDdtScheme<Type>::fvcDdt
@@ -333,6 +340,33 @@ localEulerDdtScheme<Type>::fvcDdt
 }
 
 
+template<class Type>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+localEulerDdtScheme<Type>::fvcDdt
+(
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& sf
+)
+{
+    const surfaceScalarField& rDeltaT = localRDeltaTf();
+
+    IOobject ddtIOobject
+    (
+        "ddt("+sf.name()+')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    return tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+    (
+        new GeometricField<Type, fvsPatchField, surfaceMesh>
+        (
+            ddtIOobject,
+            rDeltaT*(sf - sf.oldTime())
+        )
+    );
+}
+
+
 template<class Type>
 tmp<fvMatrix<Type>>
 localEulerDdtScheme<Type>::fvmDdt
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.H
index 73ba7777a8..3e07cda8a2 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.H
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.H
@@ -69,6 +69,9 @@ public:
     //- Name of the reciprocal local time-step field
     static word rDeltaTName;
 
+    //- Name of the reciprocal local face time-step field
+    static word rDeltaTfName;
+
     //- Name of the reciprocal local sub-cycling time-step field
     static word rSubDeltaTName;
 
@@ -88,6 +91,10 @@ public:
         //  looked-up from the objectRegistry
         static const volScalarField& localRDeltaT(const fvMesh& mesh);
 
+        //- Return the reciprocal of the local face time-step
+        //  looked-up from the objectRegistry
+        static const surfaceScalarField& localRDeltaTf(const fvMesh& mesh);
+
         //- Calculate and return the reciprocal of the local sub-cycling
         //  time-step
         static tmp<volScalarField> localRSubDeltaT
@@ -119,6 +126,9 @@ class localEulerDdtScheme
         //- Return the reciprocal of the local time-step
         const volScalarField& localRDeltaT() const;
 
+        //- Return the reciprocal of the local face time-step
+        const surfaceScalarField& localRDeltaTf() const;
+
 
 public:
 
@@ -178,6 +188,11 @@ public:
             const GeometricField<Type, fvPatchField, volMesh>& psi
         );
 
+        tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> fvcDdt
+        (
+            const GeometricField<Type, fvsPatchField, surfaceMesh>&
+        );
+
         tmp<fvMatrix<Type>> fvmDdt
         (
             const GeometricField<Type, fvPatchField, volMesh>&
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcDdt.C b/src/finiteVolume/finiteVolume/fvc/fvcDdt.C
index 391da38be3..6087724834 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcDdt.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcDdt.C
@@ -125,6 +125,21 @@ ddt
 }
 
 
+template<class Type>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+ddt
+(
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& sf
+)
+{
+    return fv::ddtScheme<Type>::New
+    (
+        sf.mesh(),
+        sf.mesh().ddtScheme("ddt(" + sf.name() + ')')
+    ).ref().fvcDdt(sf);
+}
+
+
 template<class Type>
 tmp<GeometricField<Type, fvPatchField, volMesh>>
 ddt
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcDdt.H b/src/finiteVolume/finiteVolume/fvc/fvcDdt.H
index c164445ebf..0c6ee4b35a 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcDdt.H
+++ b/src/finiteVolume/finiteVolume/fvc/fvcDdt.H
@@ -88,6 +88,12 @@ namespace fvc
         const GeometricField<Type, fvPatchField, volMesh>&
     );
 
+    template<class Type>
+    tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> ddt
+    (
+        const GeometricField<Type, fvsPatchField, surfaceMesh>&
+    );
+
     template<class Type>
     tmp<GeometricField<Type, fvPatchField, volMesh>> ddt
     (