diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index 5dd28f50cc..ad564e6d89 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -39,7 +39,7 @@ void Foam::KinematicParcel<ParcelType>::updateCellQuantities
 )
 {
     rhoc_ = td.rhoInterp().interpolate(this->position(), cellI);
-    if (rhoc_ < SMALL)
+    if (rhoc_ < td.constProps().rhoMin())
     {
         WarningIn
         (
@@ -49,10 +49,17 @@ void Foam::KinematicParcel<ParcelType>::updateCellQuantities
                 "const scalar, "
                 "const label"
             ")"
-        )   << "Density < " << SMALL << " in cell " << cellI << nl << endl;
+        )   << "Limiting density in cell " << cellI << " to "
+            << td.constProps().rhoMin() <<  nl << endl;
+
+        rhoc_ = td.constProps().rhoMin();
     }
 
     Uc_ = td.UInterp().interpolate(this->position(), cellI);
+
+    // Apply correction to cell velocity to account for momentum transfer
+    Uc_ += td.cloud().UTrans()[cellI]/(massCell(cellI));
+
     muc_ = td.muInterp().interpolate(this->position(), cellI);
 
     // Apply dispersion components to carrier phase velocity
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
index f90c86b918..d5ccfd84b2 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
@@ -89,6 +89,9 @@ public:
             //- Constant properties dictionary
             const dictionary dict_;
 
+            //- Minimum density [kg/m3]
+            const scalar rhoMin_;
+
             //- Particle density [kg/m3] (constant)
             const scalar rho0_;
 
@@ -106,6 +109,9 @@ public:
             //- Return const access to the constant properties dictionary
             inline const dictionary& dict() const;
 
+            //- Return const access to the minimum density
+            inline scalar rhoMin() const;
+
             //- Return const access to the particle density
             inline scalar rho0() const;
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
index bf4545032a..2adc0e944d 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
@@ -33,6 +33,7 @@ inline Foam::KinematicParcel<ParcelType>::constantProperties::constantProperties
 )
 :
     dict_(parentDict.subDict("constantProperties")),
+    rhoMin_(dimensionedScalar(dict_.lookup("rhoMin")).value()),
     rho0_(dimensionedScalar(dict_.lookup("rho0")).value()),
     minParticleMass_
     (
@@ -99,6 +100,14 @@ Foam::KinematicParcel<ParcelType>::constantProperties::dict() const
 }
 
 
+template <class ParcelType>
+inline Foam::scalar
+Foam::KinematicParcel<ParcelType>::constantProperties::rhoMin() const
+{
+    return rhoMin_;
+}
+
+
 template <class ParcelType>
 inline Foam::scalar
 Foam::KinematicParcel<ParcelType>::constantProperties::rho0() const
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
index 7527690087..4529d4e077 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
@@ -40,7 +40,7 @@ void Foam::ReactingParcel<ParcelType>::updateCellQuantities
     ThermoParcel<ParcelType>::updateCellQuantities(td, dt, cellI);
 
     pc_ = td.pInterp().interpolate(this->position(), cellI);
-    if (pc_ < SMALL)
+    if (pc_ < td.constProps().pMin())
     {
         WarningIn
         (
@@ -50,8 +50,19 @@ void Foam::ReactingParcel<ParcelType>::updateCellQuantities
                 "const scalar, "
                 "const label"
             ")"
-        )   << "Pressure < " << SMALL << " in cell " << cellI << nl << endl;
+        )   << "Limiting pressure in cell " << cellI << " to "
+            << td.constProps().pMin() <<  nl << endl;
+
+        pc_ = td.constProps().pMin();
     }
+
+    // Apply correction to cell density to account for mass transfer
+    scalar addedMass = 0.0;
+    forAll(td.cloud().rhoTrans(), i)
+    {
+        addedMass += td.cloud().rhoTrans(i)[cellI];
+    }
+    this->rhoc_ += addedMass/td.cloud().pMesh().cellVolumes()[cellI];
 }
 
 
@@ -259,15 +270,16 @@ Foam::scalar Foam::ReactingParcel<ParcelType>::calcPhaseChange
         return 0.0;
     }
 
+    // Calculate mass transfer due to phase change
     td.cloud().phaseChange().calculate
     (
         dt,
         cellI,
         d,
-        min(T, td.constProps().Tbp()), // Limiting to boiling temperature
+        min(T, td.constProps().Tbp()), // Limit to boiling temperature
         pc_,
         this->Tc_,
-        this->muc_/this->rhoc_,
+        this->muc_/(this->rhoc_ + ROOTVSMALL),
         U - this->Uc_,
         dMassPC
     );
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
index f876fa6401..429ceffea3 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
@@ -80,6 +80,9 @@ public:
     {
         // Private data
 
+            //- Minimum pressure [Pa]
+            const scalar pMin_;
+
             //- Constant volume flag - e.g. during mass transfer
             Switch constantVolume_;
 
@@ -97,6 +100,9 @@ public:
 
         // Access
 
+            //- Return const access to the minimum pressure
+            inline scalar pMin() const;
+
             //- Return const access to the constant volume flag
             inline Switch constantVolume() const;
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
index 0cc40887d4..344eb09398 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
@@ -33,6 +33,7 @@ inline Foam::ReactingParcel<ParcelType>::constantProperties::constantProperties
 )
 :
     ThermoParcel<ParcelType>::constantProperties(parentDict),
+    pMin_(dimensionedScalar(this->dict().lookup("pMin")).value()),
     constantVolume_(this->dict().lookup("constantVolume")),
     Tbp_(dimensionedScalar(this->dict().lookup("Tbp")).value()),
     Tvap_(dimensionedScalar(this->dict().lookup("Tvap")).value())
@@ -106,6 +107,14 @@ inline Foam::ReactingParcel<ParcelType>::ReactingParcel
 
 // * * * * * * * * * constantProperties Member Functions * * * * * * * * * * //
 
+template<class ParcelType>
+inline Foam::scalar
+Foam::ReactingParcel<ParcelType>::constantProperties::pMin() const
+{
+    return pMin_;
+}
+
+
 template<class ParcelType>
 inline Foam::Switch
 Foam::ReactingParcel<ParcelType>::constantProperties::constantVolume() const
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
index e13420e56d..8b5a722d30 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
@@ -39,8 +39,15 @@ void Foam::ThermoParcel<ParcelType>::updateCellQuantities
 {
     KinematicParcel<ParcelType>::updateCellQuantities(td, dt, cellI);
 
+    cpc_ = td.cpInterp().interpolate(this->position(), cellI);
+
     Tc_ = td.TInterp().interpolate(this->position(), cellI);
-    if (Tc_ < SMALL)
+
+    // Apply correction to cell temperature to account for enthalpy transfer
+    scalar cpMean = td.cpInterp().psi()[cellI];
+    Tc_ += td.cloud().hTrans()[cellI]/(cpMean*this->massCell(cellI));
+
+    if (Tc_ < td.constProps().TMin())
     {
         WarningIn
         (
@@ -50,10 +57,11 @@ void Foam::ThermoParcel<ParcelType>::updateCellQuantities
                 "const scalar, "
                 "const label"
             ")"
-        )   << "Temperature < " << SMALL << " in cell " << cellI << nl << endl;
-    }
+        )   << "Limiting temperature in cell " << cellI << " to "
+            << td.constProps().TMin() <<  nl << endl;
 
-    cpc_ = td.cpInterp().interpolate(this->position(), cellI);
+        Tc_ = td.constProps().TMin();
+    }
 }
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
index c7ba92ab84..18a36dc024 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
@@ -85,6 +85,9 @@ public:
             //- Particle initial temperature [K]
             const scalar T0_;
 
+            //- Minimum temperature [K]
+            const scalar TMin_;
+
             //- Particle specific heat capacity [J/(kg.K)]
             const scalar cp0_;
 
@@ -104,18 +107,22 @@ public:
 
             // Access
 
-                //- Return const access to particle initial temperature [K]
+                //- Return const access to the particle initial temperature [K]
                 inline scalar T0() const;
 
-                //- Return const access to particle specific heat capacity
+                //- Return const access to minimum temperature [K]
+                inline scalar TMin() const;
+
+                //- Return const access to the particle specific heat capacity
                 //  [J/(kg.K)]
                 inline scalar cp0() const;
 
-                //- Return const access to Particle emissivity [] (radiation)
+                //- Return const access to the particle emissivity []
+                //  Active for radiation only
                 inline scalar epsilon0() const;
 
-                //- Return const access to Particle scattering factor []
-                //  (radiation)
+                //- Return const access to the particle scattering factor []
+                //  Active for radiation only
                 inline scalar f0() const;
     };
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H
index 1a2c0326d3..6e929b6e69 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H
@@ -34,6 +34,7 @@ inline Foam::ThermoParcel<ParcelType>::constantProperties::constantProperties
 :
     KinematicParcel<ParcelType>::constantProperties(parentDict),
     T0_(dimensionedScalar(this->dict().lookup("T0")).value()),
+    TMin_(dimensionedScalar(this->dict().lookup("TMin")).value()),
     cp0_(dimensionedScalar(this->dict().lookup("cp0")).value()),
     epsilon0_(dimensionedScalar(this->dict().lookup("epsilon0")).value()),
     f0_(dimensionedScalar(this->dict().lookup("f0")).value())
@@ -110,6 +111,14 @@ Foam::ThermoParcel<ParcelType>::constantProperties::T0() const
 }
 
 
+template <class ParcelType>
+inline Foam::scalar
+Foam::ThermoParcel<ParcelType>::constantProperties::TMin() const
+{
+    return TMin_;
+}
+
+
 template <class ParcelType>
 inline Foam::scalar
 Foam::ThermoParcel<ParcelType>::constantProperties::cp0() const