diff --git a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationDiffusionLimitedRate/COxidationDiffusionLimitedRate.C b/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationDiffusionLimitedRate/COxidationDiffusionLimitedRate.C
index f60e67099a94fece6951b1381b0a977350d3d390..1ac6efe24fcc153fbad4567fcb68d048af36279b 100644
--- a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationDiffusionLimitedRate/COxidationDiffusionLimitedRate.C
+++ b/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationDiffusionLimitedRate/COxidationDiffusionLimitedRate.C
@@ -157,12 +157,15 @@ Foam::scalar Foam::COxidationDiffusionLimitedRate<CloudType>::calculate
     dMassSRCarrier[O2GlobalId_] -= dmO2;
     dMassSRCarrier[CO2GlobalId_] += dmCO2;
 
-    const scalar HC = thermo.solids().properties()[CsLocalId_].H(T);
+    const scalar HsC = thermo.solids().properties()[CsLocalId_].Hs(T);
     const scalar HCO2 = thermo.carrier().H(CO2GlobalId_, T);
-    const scalar HO2 = thermo.carrier().H(O2GlobalId_, T);
+
+    // carrier enthalpy transfer handled by change in composition
+    // const scalar HsO2 = thermo.carrier().Hs(O2GlobalId_, T);
+    // dhsTrans -= dmO2*HsO2;
 
     // Heat of reaction [J]
-    return dmC*HC + dmO2*HO2 - dmCO2*HCO2;
+    return dmC*HsC - dmCO2*HCO2;
 }
 
 
diff --git a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationKineticDiffusionLimitedRate/COxidationKineticDiffusionLimitedRate.C b/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationKineticDiffusionLimitedRate/COxidationKineticDiffusionLimitedRate.C
index ef32b558f3525d5da7f15e3b0a3ffaaa160faa31..5f717530abb838483264d09e25fcc61c38839726 100644
--- a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationKineticDiffusionLimitedRate/COxidationKineticDiffusionLimitedRate.C
+++ b/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationKineticDiffusionLimitedRate/COxidationKineticDiffusionLimitedRate.C
@@ -160,12 +160,15 @@ Foam::scalar Foam::COxidationKineticDiffusionLimitedRate<CloudType>::calculate
     dMassSRCarrier[O2GlobalId_] -= dmO2;
     dMassSRCarrier[CO2GlobalId_] += dmCO2;
 
-    const scalar HC = thermo.solids().properties()[CsLocalId_].H(T);
+    const scalar HsC = thermo.solids().properties()[CsLocalId_].Hs(T);
     const scalar HCO2 = thermo.carrier().H(CO2GlobalId_, T);
-    const scalar HO2 = thermo.carrier().H(O2GlobalId_, T);
+
+    // carrier enthalpy transfer handled by change in composition
+    // const scalar HsO2 = thermo.carrier().Hs(O2GlobalId_, T);
+    // dhsTrans -= dmO2*HsO2;
 
     // Heat of reaction [J]
-    return dmC*HC + dmO2*HO2 - dmCO2*HCO2;
+    return dmC*HsC - dmCO2*HCO2;
 }
 
 
diff --git a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationMurphyShaddix/COxidationMurphyShaddix.C b/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationMurphyShaddix/COxidationMurphyShaddix.C
index a26b8242de19a396449664eeeffc2c504075105f..c2ba2899a7b809e56264284525f770fdb372c717 100644
--- a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationMurphyShaddix/COxidationMurphyShaddix.C
+++ b/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationMurphyShaddix/COxidationMurphyShaddix.C
@@ -220,12 +220,15 @@ Foam::scalar Foam::COxidationMurphyShaddix<CloudType>::calculate
     // Add to particle mass transfer
     dMassSolid[CsLocalId_] += dOmega*WC_;
 
-    const scalar HC = thermo.solids().properties()[CsLocalId_].H(T);
+    const scalar HsC = thermo.solids().properties()[CsLocalId_].Hs(T);
     const scalar HCO2 = thermo.carrier().H(CO2GlobalId_, T);
-    const scalar HO2 = thermo.carrier().H(O2GlobalId_, T);
+
+    // carrier enthalpy transfer handled by change in composition
+    // const scalar HsO2 = thermo.carrier().Hs(O2GlobalId_, T);
+    // dhsTrans -= dmO2*HsO2;
 
     // Heat of reaction
-    return dOmega*(WC_*HC + WO2_*HO2 - (WC_ + WO2_)*HCO2);
+    return dOmega*(WC_*HsC - (WC_ + WO2_)*HCO2);
 }
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
index fcad1a4ed3be5833a05df040c100dec90b4ff089..336b3ef705b897c1da725b6aaeb07f2109d6ec77 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
@@ -165,6 +165,11 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     const label cellI
 )
 {
+    typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
+    const CompositionModel<reactingCloudType>& composition =
+        td.cloud().composition();
+
+
     // Define local properties at beginning of timestep
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     const scalar np0 = this->nParticle_;
@@ -178,9 +183,9 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     const scalar pc = this->pc_;
 
     const scalarField& YMix = this->Y_;
-    const label idG = td.cloud().composition().idGas();
-    const label idL = td.cloud().composition().idLiquid();
-    const label idS = td.cloud().composition().idSolid();
+    const label idG = composition.idGas();
+    const label idL = composition.idLiquid();
+    const label idS = composition.idSolid();
 
 
     // Calc surface values
@@ -221,7 +226,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     scalar NCpW = 0.0;
 
     // Surface concentrations of emitted species
-    scalarField Cs(td.cloud().composition().carrier().species().size(), 0.0);
+    scalarField Cs(composition.carrier().species().size(), 0.0);
 
     // Calc mass and enthalpy transfer due to phase change
     this->calcPhaseChange
@@ -271,10 +276,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         Cs
     );
 
-    // Correct surface values due to emitted species
-    this->correctSurfaceValues(td, cellI, Ts, Cs, rhos, mus, Prs, kappas);
-    Res = this->Re(U0, d0, rhos, mus);
-
 
     // Surface reactions
     // ~~~~~~~~~~~~~~~~~
@@ -283,13 +284,9 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     scalarField dMassSRGas(YGas_.size(), 0.0);
     scalarField dMassSRLiquid(YLiquid_.size(), 0.0);
     scalarField dMassSRSolid(YSolid_.size(), 0.0);
-    scalarField dMassSRCarrier
-    (
-        td.cloud().composition().carrier().species().size(),
-        0.0
-    );
+    scalarField dMassSRCarrier(composition.carrier().species().size(), 0.0);
 
-    // Clac mass and enthalpy transfer due to surface reactions
+    // Calc mass and enthalpy transfer due to surface reactions
     calcSurfaceReactions
     (
         td,
@@ -313,6 +310,11 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     );
 
 
+    // Correct surface values due to emitted species
+    this->correctSurfaceValues(td, cellI, Ts, Cs, rhos, mus, Prs, kappas);
+    Res = this->Re(U0, d0, rhos, mus);
+
+
     // Update component mass fractions
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -324,6 +326,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         updateMassFractions(mass0, dMassGas, dMassLiquid, dMassSolid);
 
 
+
+
     // Heat transfer
     // ~~~~~~~~~~~~~
 
@@ -380,25 +384,33 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         // Transfer mass lost from particle to carrier mass source
         forAll(YGas_, i)
         {
-            label gid = td.cloud().composition().localToGlobalCarrierId(GAS, i);
+            label gid = composition.localToGlobalCarrierId(GAS, i);
             td.cloud().rhoTrans(gid)[cellI] += np0*dMassGas[i];
+//            td.cloud().hsTrans()[cellI] +=
+//                np0*dMassGas[i]*composition.carrier().Hs(gid, T0);
         }
         forAll(YLiquid_, i)
         {
-            label gid = td.cloud().composition().localToGlobalCarrierId(LIQ, i);
+            label gid = composition.localToGlobalCarrierId(LIQ, i);
             td.cloud().rhoTrans(gid)[cellI] += np0*dMassLiquid[i];
+//            td.cloud().hsTrans()[cellI] +=
+//                np0*dMassLiquid[i]*composition.carrier().Hs(gid, T0);
         }
 /*
         // No mapping between solid components and carrier phase
         forAll(YSolid_, i)
         {
-            label gid = td.cloud().composition().localToGlobalCarrierId(SLD, i);
+            label gid = composition.localToGlobalCarrierId(SLD, i);
             td.cloud().rhoTrans(gid)[cellI] += np0*dMassSolid[i];
+//            td.cloud().hsTrans()[cellI] +=
+//                np0*dMassSolid[i]*composition.carrier().Hs(gid, T0);
         }
 */
         forAll(dMassSRCarrier, i)
         {
             td.cloud().rhoTrans(i)[cellI] += np0*dMassSRCarrier[i];
+//            td.cloud().hsTrans()[cellI] +=
+//                np0*dMassSRCarrier[i]*composition.carrier().Hs(i, T0);
         }
 
         // Update momentum transfer
@@ -427,14 +439,12 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
             // Absorb parcel into carrier phase
             forAll(YGas_, i)
             {
-                label gid =
-                    td.cloud().composition().localToGlobalCarrierId(GAS, i);
+                label gid = composition.localToGlobalCarrierId(GAS, i);
                 td.cloud().rhoTrans(gid)[cellI] += np0*mass1*YMix[GAS]*YGas_[i];
             }
             forAll(YLiquid_, i)
             {
-                label gid =
-                    td.cloud().composition().localToGlobalCarrierId(LIQ, i);
+                label gid = composition.localToGlobalCarrierId(LIQ, i);
                 td.cloud().rhoTrans(gid)[cellI] +=
                     np0*mass1*YMix[LIQ]*YLiquid_[i];
             }
@@ -442,8 +452,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
             // No mapping between solid components and carrier phase
             forAll(YSolid_, i)
             {
-                label gid =
-                    td.cloud().composition().localToGlobalCarrierId(SLD, i);
+                label gid = composition.localToGlobalCarrierId(SLD, i);
                 td.cloud().rhoTrans(gid)[cellI] +=
                     np0*mass1*YMix[SLD]*YSolid_[i];
             }
@@ -508,6 +517,11 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
         return;
     }
 
+    typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
+    const CompositionModel<reactingCloudType>& composition =
+        td.cloud().composition();
+
+
     // Total mass of volatiles evolved
     td.cloud().devolatilisation().calculate
     (
@@ -535,10 +549,9 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
         // Note: hardcoded gaseous diffusivities for now
         // TODO: add to carrier thermo
         const scalar beta = sqr(cbrt(15.0) + cbrt(15.0));
-        const label id =
-            td.cloud().composition().localToGlobalCarrierId(GAS, i);
-        const scalar Cp = td.cloud().composition().carrier().Cp(id, Ts);
-        const scalar W = td.cloud().composition().carrier().W(id);
+        const label id = composition.localToGlobalCarrierId(GAS, i);
+        const scalar Cp = composition.carrier().Cp(id, Ts);
+        const scalar W = composition.carrier().W(id);
         const scalar Ni = dMassDV[i]/(this->areaS(d)*dt*W);
 
         // Dab calc'd using API vapour mass diffusivity function
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
index 855439b8939c627987e8aeae7ddf1104a94aac3e..11e61d8581d5ec1c187245d505f087e9142ffabb 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
@@ -252,6 +252,11 @@ void Foam::ReactingParcel<ParcelType>::calc
     const label cellI
 )
 {
+    typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
+    const CompositionModel<reactingCloudType>& composition =
+        td.cloud().composition();
+
+
     // Define local properties at beginning of time step
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     const scalar np0 = this->nParticle_;
@@ -301,7 +306,7 @@ void Foam::ReactingParcel<ParcelType>::calc
     scalar NCpW = 0.0;
 
     // Surface concentrations of emitted species
-    scalarField Cs(td.cloud().composition().carrier().species().size(), 0.0);
+    scalarField Cs(composition.carrier().species().size(), 0.0);
 
     // Calc mass and enthalpy transfer due to phase change
     calcPhaseChange
@@ -389,8 +394,10 @@ void Foam::ReactingParcel<ParcelType>::calc
         // Transfer mass lost from particle to carrier mass source
         forAll(dMassPC, i)
         {
-            label gid = td.cloud().composition().localToGlobalCarrierId(0, i);
+            label gid = composition.localToGlobalCarrierId(0, i);
             td.cloud().rhoTrans(gid)[cellI] += np0*dMassPC[i];
+//            td.cloud().hsTrans()[cellI] +=
+//                np0*dMassPC[i]*composition.carrier().Hs(gid, T0);
         }
 
         // Update momentum transfer
@@ -418,13 +425,12 @@ void Foam::ReactingParcel<ParcelType>::calc
             // Absorb parcel into carrier phase
             forAll(Y_, i)
             {
-                label gid =
-                    td.cloud().composition().localToGlobalCarrierId(0, i);
+                label gid = composition.localToGlobalCarrierId(0, i);
                 td.cloud().rhoTrans(gid)[cellI] += np0*mass1*Y_[i];
             }
             td.cloud().UTrans()[cellI] += np0*mass1*U1;
             td.cloud().hsTrans()[cellI] +=
-                np0*mass1*td.cloud().composition().H(0, Y_, pc_, T1);
+                np0*mass1*composition.H(0, Y_, pc_, T1);
         }
     }
 
@@ -434,7 +440,7 @@ void Foam::ReactingParcel<ParcelType>::calc
 
     else
     {
-        this->Cp_ = td.cloud().composition().Cp(0, Y_, pc_, T1);
+        this->Cp_ = composition.Cp(0, Y_, pc_, T1);
         this->T_ = T1;
         this->U_ = U1;
 
@@ -484,6 +490,11 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
         return;
     }
 
+    typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
+    const CompositionModel<reactingCloudType>& composition =
+        td.cloud().composition();
+
+
     // Calculate mass transfer due to phase change
     td.cloud().phaseChange().calculate
     (
@@ -511,19 +522,18 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
 
     forAll(YComponents, i)
     {
-        const label idc =
-            td.cloud().composition().localToGlobalCarrierId(idPhase, i);
-        const label idl = td.cloud().composition().globalIds(idPhase)[i];
+        const label idc = composition.localToGlobalCarrierId(idPhase, i);
+        const label idl = composition.globalIds(idPhase)[i];
 
         const scalar dh = td.cloud().phaseChange().dh(idc, idl, pc_, T);
         Sh -= dMassPC[i]*dh/dt;
 
         // Update particle surface thermo properties
         const scalar Dab =
-            td.cloud().composition().liquids().properties()[idl].D(pc_, Ts, Wc);
+            composition.liquids().properties()[idl].D(pc_, Ts, Wc);
 
-        const scalar Cp = td.cloud().composition().carrier().Cp(idc, Ts);
-        const scalar W = td.cloud().composition().carrier().W(idc);
+        const scalar Cp = composition.carrier().Cp(idc, Ts);
+        const scalar W = composition.carrier().W(idc);
         const scalar Ni = dMassPC[i]/(this->areaS(d)*dt*W);
 
         // Molar flux of species coming from the particle (kmol/m^2/s)
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
index 53a8e5545e703179369c18283446bbc500760a8b..6397e29437797466a8a629d8ac9800bd8b982318 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
@@ -208,8 +208,8 @@ void Foam::ThermoParcel<ParcelType>::calc
     // Sum Ni*Cpi*Wi of emission species
     scalar NCpW = 0.0;
 
-    // Calculate new particle velocity
-    scalar Cuh = 0.0;
+    // Calculate new particle temperature
+    scalar Sph = 0.0;
     scalar T1 =
         this->calcHeatTransfer
         (
@@ -226,7 +226,7 @@ void Foam::ThermoParcel<ParcelType>::calc
             NCpW,
             Sh,
             dhsTrans,
-            Cuh
+            Sph
         );
 
 
@@ -267,7 +267,7 @@ void Foam::ThermoParcel<ParcelType>::calc
         td.cloud().hsTrans()[cellI] += np0*dhsTrans;
 
         // Update sensible enthalpy coefficient
-        td.cloud().hsCoeff()[cellI] += np0*Cuh*this->areaS();
+        td.cloud().hsCoeff()[cellI] += np0*Sph;
     }
 
     // Set new particle properties
@@ -294,7 +294,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
     const scalar NCpW,
     const scalar Sh,
     scalar& dhsTrans,
-    scalar& Cuh
+    scalar& Sph
 )
 {
     if (!td.cloud().heatTransfer().active())
@@ -317,6 +317,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
 
     htc = max(htc, ROOTVSMALL);
     const scalar As = this->areaS(d);
+
     scalar ap = Tc_ + Sh/As/htc;
     scalar bp = 6.0*(Sh/As + htc*(Tc_ - T));
     if (td.cloud().radiation())
@@ -337,9 +338,9 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
 
     scalar Tnew = max(Tres.value(), td.cloud().constProps().TMin());
 
-    dhsTrans += dt*htc*As*(0.5*(T + Tnew) - Tc_);
+    Sph = dt*htc*As;
 
-    Cuh = dt*bp;
+    dhsTrans += Sph*(0.5*(T + Tnew) - Tc_);
 
     return Tnew;
 }
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
index e47e1da48a60ba7ca4153b9fa7629484fea7a5b2..c8fe14ad1a6eaa337f7c64f62d8c030bd3a7faf2 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
@@ -243,7 +243,7 @@ protected:
             const scalar NCpW,         // Sum of N*Cp*W of emission species
             const scalar Sh,           // explicit particle enthalpy source
             scalar& dhsTrans,          // sensible enthalpy transfer to carrier
-            scalar& Cuh                // linearised heat transfer coefficient
+            scalar& Sph                // linearised heat transfer coefficient
         );