diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
index 03a8c01aa6eb2436e0acf233ea360bf17074bf6f..aacd51f04ae380f413dfa941bc407abdff95ce46 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
@@ -162,7 +162,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     scalarField dMassSRGas(YGas_.size(), 0.0);
     scalarField dMassSRLiquid(YLiquid_.size(), 0.0);
     scalarField dMassSRSolid(YSolid_.size(), 0.0);
-    scalarField dMassSRc(td.cloud().gases().size(), 0.0);
+    scalarField dMassSRCarrier(td.cloud().gases().size(), 0.0);
 
 
     // Phase change in liquid phase
@@ -177,8 +177,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     // ~~~~~~~~~~~~~~~~
 
     // Return enthalpy source and calc mass transfer due to devolatilisation
-    scalar ShDV =
-        calcDevolatilisation(td, dt, T0, mass0, idG, YMix, dMassDV);
+    scalar ShDV = calcDevolatilisation(td, dt, T0, mass0, idG, YMix, dMassDV);
 
 
     // Surface reactions
@@ -196,7 +195,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
             dMassSRGas,
             dMassSRLiquid,
             dMassSRSolid,
-            dMassSRc,
+            dMassSRCarrier,
             dhTrans
         );
 
@@ -257,9 +256,9 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
 //            label id = td.composition().localToGlobalGasId(SOLID, i);
 //            td.cloud().rhoTrans(id)[cellI] += np0*dMassSolid[i];
 //        }
-        forAll(dMassSRc, i)
+        forAll(dMassSRCarrier, i)
         {
-            td.cloud().rhoTrans(i)[cellI] += np0*dMassSRc[i];
+            td.cloud().rhoTrans(i)[cellI] += np0*dMassSRCarrier[i];
         }
 
         // Update momentum transfer
@@ -270,7 +269,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
 
         // Update enthalpy transfer
         // - enthalpy of lost solids already accounted for
-        td.cloud().hTrans()[cellI] += np0*dhTrans;
+        td.cloud().hTrans()[cellI] += np0*(dhTrans + Sh);
 
         // Coefficient to be applied in carrier phase enthalpy coupling
         td.cloud().hCoeff()[cellI] += np0*htc*this->areaS();
@@ -286,25 +285,36 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
 
         if (td.cloud().coupled())
         {
-            // Absorb particle(s) into carrier phase
+            // Absorb parcel into carrier phase
             forAll(YGas_, i)
             {
                 label id = td.composition().localToGlobalGasId(GAS, i);
-                td.cloud().rhoTrans(id)[cellI] += np0*dMassGas[i];
+                td.cloud().rhoTrans(id)[cellI] += np0*mass1*YMix[GAS]*YGas_[i];
             }
             forAll(YLiquid_, i)
             {
                 label id = td.composition().localToGlobalGasId(LIQUID, i);
-                td.cloud().rhoTrans(id)[cellI] += np0*dMassLiquid[i];
+                td.cloud().rhoTrans(id)[cellI] +=
+                    np0
+                   *mass1
+                   *YMix[LIQUID]
+                   *YLiquid_[i];
             }
             // No mapping between solid components and carrier phase
 //            forAll(YSolid_, i)
 //            {
 //                label id = td.composition().localToGlobalGasId(SOLID, i);
-//                td.cloud().rhoTrans(id)[cellI] += np0*dMassSolid[i];
+//                td.cloud().rhoTrans(id)[cellI] +=
+//                    np0
+//                   *mass1
+//                   *YMix[SOLID]
+//                   *YSolid_[i];
 //            }
 
-            td.cloud().hTrans()[cellI] += np0*mass1*HEff(td, pc, T1, idG, idL, idS);
+            td.cloud().hTrans()[cellI] +=
+                np0
+               *mass1
+               *HEff(td, pc, T1, idG, idL, idS);
             td.cloud().UTrans()[cellI] += np0*mass1*U1;
         }
     }
@@ -397,7 +407,7 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
     scalarField& dMassSRGas,
     scalarField& dMassSRLiquid,
     scalarField& dMassSRSolid,
-    scalarList& dMassSRc,
+    scalarList& dMassSRCarrier,
     scalar& dhTrans
 )
 {
@@ -426,7 +436,7 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
         dMassSRGas,
         dMassSRLiquid,
         dMassSRSolid,
-        dMassSRc
+        dMassSRCarrier
     );
 
     // Heat of reaction divided between particle and carrier phase by the
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
index e213ddace8f2fc92dd1586898c618e52a40a593d..58e2ef60dc8664a96b764f1f97b4232b5b459c1c 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
@@ -66,9 +66,15 @@ class ReactingMultiphaseParcel
 :
     public ReactingParcel<ParcelType>
 {
-
 public:
 
+    // IDs of phases in ReacingParcel phase list (Y)
+
+        static const label GAS;
+        static const label LIQUID;
+        static const label SOLID;
+
+
     //- Class to hold reacting particle constant properties
     class constantProperties
     :
@@ -185,13 +191,6 @@ protected:
 
     // Protected data
 
-        // IDs of phases in parent phase list
-
-            static const label GAS;
-            static const label LIQUID;
-            static const label SOLID;
-
-
         // Parcel properties
 
             //- Mass fractions of gases []
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
index 8b59c5017a81b966394fd1f4e1f7d134883d106c..676e0f10b2cde314df6442e0532a3ed586d30565 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
@@ -166,7 +166,7 @@ void Foam::ReactingParcel<ParcelType>::calc
 
         if (td.cloud().coupled())
         {
-            // Absorb particle(s) into carrier phase
+            // Absorb parcel into carrier phase
             forAll(Y_, i)
             {
                 label id = td.composition().localToGlobalGasId(0, i);
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
index b22a8e57e67362beffa7b9d7908f8f404a7f2060..276fcd0b59940da15569b0cad31bbaba25e84c05 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
@@ -135,13 +135,12 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
     const label cellI,
     const scalar Sh,
     scalar& htc,
-    scalar& ShHT
+    scalar& dhTrans
 )
 {
     if (!td.cloud().heatTransfer().active())
     {
         htc = 0.0;
-        ShHT = 0.0;
         return T_;
     }
 
@@ -190,7 +189,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
 
     // Enthalpy transfer
     // - Using average particle temperature
-    ShHT = dt*Ap*htc*(Tres.average() - Tc_);
+    dhTrans = dt*Ap*htc*(Tres.average() - Tc_);
 
     return Tres.value();
 }
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
index 584588d695857a379ef255fb8cd55e13bfe9fbea..830c4d7f1d2f703ac4ac90345c1b8d08951fbbd0 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
@@ -29,44 +29,6 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-template <class CloudType>
-Foam::label Foam::LiquidEvaporation<CloudType>::carrierSpecieId
-(
-    const word& specieName
-) const
-{
-    forAll (this->owner().carrierThermo().composition().Y(), i)
-    {
-        if
-        (
-            this->owner().carrierThermo().composition().Y()[i].name()
-         == specieName
-        )
-        {
-            return i;
-        }
-    }
-
-    wordList species(this->owner().carrierThermo().composition().Y().size());
-    forAll (this->owner().carrierThermo().composition().Y(), i)
-    {
-        species[i] = this->owner().carrierThermo().composition().Y()[i].name();
-    }
-
-    FatalErrorIn
-    (
-        "Foam::label Foam::LiquidEvaporation<CloudType>::carrierSpecieId"
-        "("
-            "const word&"
-        ") const"
-    )   << "Could not find " << specieName << " in species list" << nl
-        <<  "Avialable species:" << nl << species
-        << nl << exit(FatalError);
-
-    return -1;
-}
-
-
 template<class CloudType>
 Foam::scalarField Foam::LiquidEvaporation<CloudType>::calcXc
 (
@@ -136,35 +98,18 @@ Foam::LiquidEvaporation<CloudType>::LiquidEvaporation
     }
 
     // Determine mapping between liquid and carrier phase species
+    label idLiquid = owner.composition().idLiquid();
     forAll(activeLiquids_, i)
     {
-        liqToGasMap_[i] = carrierSpecieId(activeLiquids_[i]);
+        liqToGasMap_[i] =
+            owner.composition().globalId(idLiquid, activeLiquids_[i]);
     }
 
     // Determine mapping between local and global liquids
     forAll(activeLiquids_, i)
     {
-        forAll(liquids_->components(), j)
-        {
-            if (liquids_->components()[j] == activeLiquids_[i])
-            {
-                liqToLiqMap_[i] = j;
-            }
-        }
-
-        if (liqToLiqMap_[i] < 0)
-        {
-            FatalErrorIn
-            (
-                "Foam::LiquidEvaporation<CloudType>::LiquidEvaporation"
-                "("
-                    "const dictionary& dict, "
-                    "CloudType& owner"
-                ")"
-            )   << "Unable to find liquid species " << activeLiquids_[i]
-                << " in global liquids list. Avaliable liquids:" << nl
-                << liquids_->components() << nl << exit(FatalError);
-        }
+        liqToLiqMap_[i] =
+            owner.composition().localId(idLiquid, activeLiquids_[i]);
     }
 }
 
@@ -239,7 +184,7 @@ void Foam::LiquidEvaporation<CloudType>::calculate
         scalar Ni = max(kc*(Cs - Cinf), 0.0);
 
         // mass transfer [kg]
-        dMass[gid] += Ni*A*liquids_->properties()[lid].W()*dt;
+        dMass[lid] += Ni*A*liquids_->properties()[lid].W()*dt;
     }
 }
 
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H
index 0aa6953521697aca6edd2082b74f7820d23cdcae..c70c081e90dc424e92aef6d3ba940e0fab5b3b79 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H
@@ -68,9 +68,6 @@ protected:
 
     // Protected member functions
 
-        //- Return the id of species specieName in the carrier phase
-        label carrierSpecieId(const word& specieName) const;
-
         //- Sherwood number as a function of Reynolds and Schmidt numbers
         scalar Sh(const scalar Re, const scalar Sc) const;
 
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H b/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H
index 505902f8c8687ef5f220c5545e772c452889546c..6566f60456c9f7b06b3bc150170919840cee083c 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H
@@ -92,7 +92,7 @@ public:
             scalarField& dMassGas,
             scalarField& dMassLiquid,
             scalarField& dMassSolid,
-            scalarList& dMassSRc
+            scalarList& dMassSRCarrier
         ) const;
 };
 
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H b/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H
index 1a3b664d4a25b33f381f3ed9365f805a827a291d..269c37de37e2c317c71fc270920fb97ec41ad293 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H
@@ -151,7 +151,7 @@ public:
             scalarField& dMassGas,
             scalarField& dMassLiquid,
             scalarField& dMassSolid,
-            scalarList& dMassSRc
+            scalarList& dMassSRCarrier
         ) const = 0;
 };
 
diff --git a/src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/RanzMarshall/RanzMarshall.C b/src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/RanzMarshall/RanzMarshall.C
index 3cfe1f712c78f92c238bf23aee85ffa3fe6440fe..b728282001a0dc7a3f9a8b454fe052fae6d9336f 100644
--- a/src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/RanzMarshall/RanzMarshall.C
+++ b/src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/RanzMarshall/RanzMarshall.C
@@ -65,7 +65,7 @@ Foam::scalar Foam::RanzMarshall<CloudType>::Nu
     const scalar Pr
 ) const
 {
-    return 2.0 + 0.6*pow(Re, 1.0/2.0)*pow(Pr, 1.0/3.0);
+    return 2.0 + 0.6*pow(Re, 0.5)*pow(Pr, 0.333);
 }