diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/UEqns.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/UEqns.H
index 82b9295486..31bd20ff7a 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/UEqns.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/UEqns.H
@@ -18,8 +18,7 @@ volScalarField Kd(fluid.Kd());
         U1Eqn =
         (
             fvm::ddt(alpha1, rho1, U1) + fvm::div(alphaRhoPhi1, U1)
-          - fvm::Sp(phase1.continuityError() - dmdt21, U1)
-          - dmdt21*U2
+          - fvm::Sp(phase1.continuityError() + dmdt21, U1) + dmdt21*U2
           + MRF.DDt(alpha1*rho1, U1)
           + phase1.turbulence().divDevRhoReff(U1)
          ==
@@ -44,8 +43,7 @@ volScalarField Kd(fluid.Kd());
         U2Eqn =
         (
             fvm::ddt(alpha2, rho2, U2) + fvm::div(alphaRhoPhi2, U2)
-          - fvm::Sp(phase2.continuityError() + dmdt12, U2)
-          + dmdt12*U1
+          - fvm::Sp(phase2.continuityError() - dmdt12, U2) - dmdt12*U1
           + MRF.DDt(alpha2*rho2, U2)
           + phase2.turbulence().divDevRhoReff(U2)
          ==
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/UEqns.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/UEqns.H
index 0303eea106..eb1456bc15 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/UEqns.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/UEqns.H
@@ -21,10 +21,14 @@ fvVectorMatrix U2Eqn(U2, rho2.dimensions()*U2.dimensions()*dimVol/dimTime);
       + MRF.DDt(U2)
     );
 
+    const volScalarField dmdt12(posPart(fluid.dmdt()));
+    const volScalarField dmdt21(negPart(fluid.dmdt()));
+
     {
         U1Eqn =
         (
             fvm::div(alphaRhoPhi1, U1) - fvm::Sp(fvc::div(alphaRhoPhi1), U1)
+          + fvm::Sp(dmdt12, U1) + dmdt21*U2
           + MRF.DDt(alpha1*rho1, U1)
           + phase1.turbulence().divDevRhoReff(U1)
           + Vm*(UgradU1 - (UgradU2 & U2))
@@ -39,6 +43,7 @@ fvVectorMatrix U2Eqn(U2, rho2.dimensions()*U2.dimensions()*dimVol/dimTime);
         U2Eqn =
         (
             fvm::div(alphaRhoPhi2, U2) - fvm::Sp(fvc::div(alphaRhoPhi2), U2)
+          - fvm::Sp(dmdt21, U2) - dmdt12*U1
           + MRF.DDt(alpha2*rho2, U2)
           + phase2.turbulence().divDevRhoReff(U2)
           + Vm*(UgradU2 - (UgradU1 & U1))