diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PhaseTransferPhaseSystem/PhaseTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PhaseTransferPhaseSystem/PhaseTransferPhaseSystem.C
index 9d66b4d2b..06e04aff4 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PhaseTransferPhaseSystem/PhaseTransferPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PhaseTransferPhaseSystem/PhaseTransferPhaseSystem.C
@@ -162,8 +162,8 @@ Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::massTransfer() const
         const phaseModel& otherPhase = pair.phase2();
 
         const volScalarField dmdt(this->rDmdt(pair));
-        const volScalarField dmdt12(posPart(dmdt));
-        const volScalarField dmdt21(negPart(dmdt));
+        const volScalarField dmdt12(negPart(dmdt));
+        const volScalarField dmdt21(posPart(dmdt));
 
         const PtrList<volScalarField>& Yi = phase.Y();
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C
index f7f464d02..08e596def 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C
@@ -193,8 +193,8 @@ Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::massTransfer() const
         const phaseModel& otherPhase = pair.phase2();
 
         const volScalarField dmdt(this->pDmdt(pair));
-        const volScalarField dmdt12(posPart(dmdt));
-        const volScalarField dmdt21(negPart(dmdt));
+        const volScalarField dmdt12(negPart(dmdt));
+        const volScalarField dmdt21(posPart(dmdt));
 
         const PtrList<volScalarField>& Yi = phase.Y();
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
index dbe044c81..25b01d94f 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
@@ -283,11 +283,6 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::massTransfer() const
 
         forAll(Yi, i)
         {
-            if (Yi[i].member() != volatile_)
-            {
-                continue;
-            }
-
             const word name
             (
                 IOobject::groupName(volatile_, phase.name())
@@ -299,12 +294,15 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::massTransfer() const
             );
 
             const volScalarField dmdt(this->iDmdt(pair) + this->wDmdt(pair));
-            const volScalarField dmdt12(posPart(dmdt));
-            const volScalarField dmdt21(negPart(dmdt));
 
-            *eqns[name] += fvm::Sp(dmdt21, eqns[name]->psi()) - dmdt21;
-            *eqns[otherName] +=
-                dmdt12 - fvm::Sp(dmdt12, eqns[otherName]->psi());
+            *eqns[name] -= fvm::Sp(dmdt, eqns[name]->psi());
+            *eqns[otherName] += fvm::Sp(dmdt, eqns[otherName]->psi());
+
+            if (Yi[i].member() == volatile_)
+            {
+                *eqns[name] += dmdt;
+                *eqns[otherName] -= dmdt;
+            }
         }
     }