diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
index 3d67afa89b7ca401d9c5b3aaf9740ca1cfd1643e..385fe5ffbcdf5fcd2e92bcd3f3f2df5108b89986 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
@@ -108,18 +108,17 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::updateMassFractions
 {
     scalarField& YMix = this->Y_;
 
-    scalar dMassGasTot = sum(dMassGas);
-    scalar dMassLiquidTot = sum(dMassLiquid);
-    scalar dMassSolidTot = sum(dMassSolid);
+    scalar massGas =
+        this->updateMassFraction(mass0*YMix[GAS], dMassGas, YGas_);
+    scalar massLiquid =
+        this->updateMassFraction(mass0*YMix[LIQ], dMassLiquid, YLiquid_);
+    scalar massSolid =
+        this->updateMassFraction(mass0*YMix[SLD], dMassSolid, YSolid_);
 
-    this->updateMassFraction(mass0*YMix[GAS], dMassGas, YGas_);
-    this->updateMassFraction(mass0*YMix[LIQ], dMassLiquid, YLiquid_);
-    this->updateMassFraction(mass0*YMix[SLD], dMassSolid, YSolid_);
+    scalar massNew = max(massGas + massLiquid + massSolid, ROOTVSMALL);
 
-    scalar massNew = mass0 - (dMassGasTot + dMassLiquidTot + dMassSolidTot);
-
-    YMix[GAS] = (mass0*YMix[GAS] - dMassGasTot)/massNew;
-    YMix[LIQ] = (mass0*YMix[LIQ] - dMassLiquidTot)/massNew;
+    YMix[GAS] = massGas/massNew;
+    YMix[LIQ] = massLiquid/massNew;
     YMix[SLD] = 1.0 - YMix[GAS] - YMix[LIQ];
 
     return massNew;
@@ -189,6 +188,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
             d0,
             T0,
             U0,
+            mass0,
             idL,
             YMix[LIQ],
             YLiquid_,
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
index 1db979c6da654a75e2722d8a271a7abea384a2ca..cb7e83c267e66ec6a2b172bcefe5f7d3d9efffce 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
@@ -67,19 +67,25 @@ void Foam::ReactingParcel<ParcelType>::updateCellQuantities
 
 
 template<class ParcelType>
-void Foam::ReactingParcel<ParcelType>::updateMassFraction
+Foam::scalar Foam::ReactingParcel<ParcelType>::updateMassFraction
 (
     const scalar mass0,
     const scalarField& dMass,
     scalarField& Y
-)
+) const
 {
-    scalar mass1 = mass0 + sum(dMass);
+    scalar mass1 = mass0 - sum(dMass);
 
-    forAll(Y, i)
+    // only update the mass fractions if the new particle mass is finite
+    if (mass1 > ROOTVSMALL)
     {
-        Y[i] = (Y[i]*mass0 - dMass[i])/mass1;
+        forAll(Y, i)
+        {
+            Y[i] = (Y[i]*mass0 - dMass[i])/mass1;
+        }
     }
+
+    return mass1;
 }
 
 
@@ -116,13 +122,10 @@ void Foam::ReactingParcel<ParcelType>::calc
 
     // Return enthalpy source and calc mass transfer due to phase change
     scalar ShPC =
-        calcPhaseChange(td, dt, cellI, d0, T0, U0, 0, 1.0, Y_, dMassPC);
-
-    // Update particle component mass fractions
-    updateMassFraction(mass0, dMassPC, Y_);
+        calcPhaseChange(td, dt, cellI, d0, T0, U0, mass0, 0, 1.0, Y_, dMassPC);
 
-    // Update mass
-    scalar mass1 = mass0 - sum(dMassPC);
+    // Update particle component mass and mass fractions
+    scalar mass1 = updateMassFraction(mass0, dMassPC, Y_);
 
 
     // Heat transfer
@@ -218,6 +221,7 @@ Foam::scalar Foam::ReactingParcel<ParcelType>::calcPhaseChange
     const scalar d,
     const scalar T,
     const vector& U,
+    const scalar mass,
     const label idPhase,
     const scalar YPhase,
     const scalarField& YComponents,
@@ -248,6 +252,9 @@ Foam::scalar Foam::ReactingParcel<ParcelType>::calcPhaseChange
         dMassPC
     );
 
+    // Limit phase change mass by availability of each specie
+    dMassPC = min(mass*YPhase*YComponents, dMassPC);
+
     scalar dMassTot = sum(dMassPC);
 
     // Add to cumulative phase change mass
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
index 429ceffea33d833ff429d831cdfea603ac7f865b..f7b85d22c21bcada609f4453227b9484d13cd01d 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
@@ -196,6 +196,7 @@ protected:
             const scalar d,            // diameter
             const scalar T,            // temperature
             const vector& U,           // velocity
+            const scalar mass,         // mass
             const label idPhase,       // id of phase involved in phase change
             const scalar YPhase,       // total mass fraction
             const scalarField& YComponents, // component mass fractions
@@ -203,12 +204,12 @@ protected:
         );
 
         //- Update mass fraction
-        void updateMassFraction
+        scalar updateMassFraction
         (
             const scalar mass0,
             const scalarField& dMass,
             scalarField& Y
-        );
+        ) const;
 
 
 public: